如何首先 group_by() 然后通过列迭代 lm()？答案

【问题标题】：How to first group_by() then iterate lm() through columns?如何首先 group_by() 然后通过列迭代 lm()？
【发布时间】：2017-03-21 20:20:29
【问题描述】：

假设我们有一个数据框，其中包含一组 3 个因变量和 6 个由分组变量标记的自变量。使用下面的示例代码生成此格式的示例：

library(tidyverse)
library(broom)
n  <- 15
df  <- data.frame(groupingvar= sample(letters[1:2], size = n, replace = TRUE),
                  y1 = rnorm(n,10,1), y2=rnorm(n,100,10), y3=rnorm(n,1000,100),
                  x1=  rnorm(n,10,1), x2=rnorm(n,10,1), x3=rnorm(n,10,1),
                  x4=rnorm(n,10,1), x5=rnorm(n,10,1), x6=rnorm(n,10,1))
df <- arrange(df,groupingvar)

如果我想在 x1 到 x6 的集合上回归 y1、y2、y3 中的每一个，我可以使用以下内容：

y <- as.matrix(select(df,y1:y3))
x <- as.matrix(select(df,x1:x6))
regs <-lm(y~x)
coeffs <- tidy(regs)
coeffs <- arrange(coeffs,response, term)

（通过使用 lm() 帮助中的以下行：“如果响应是矩阵，则线性模型通过最小二乘法分别拟合到矩阵的每一列。”）

但是，如果我需要先按分组变量进行分组，然后应用 lm 函数，那么我不太确定该怎么做。我尝试了以下方法，但它为两组产生了相同的系数集。

regs2 <- df %>% group_by(groupingvar) %>%
  do(fit2 = lm(as.matrix(select(df,y1:y3)) ~ as.matrix(select(df,x1:x6))))
coeffs2 <- tidy(regs2,fit2)
coeffs2 <- arrange(coeffs2,groupingvar, response)

【问题讨论】：

"然后应用 lm 函数" -> 你试过用lapply()吗？
我不知道如何正确使用它。我尝试使用元素“y1~x1+x2+...+x6”、“y2~x1+x2+...+x6”、“y3~x1+x2+...+x6”和试图将此列表传递给 lm()，但我认为我在正确的语法上绊倒了。
apply、sapply、lapply 等系列对您的理解绝对至关重要。必须学习。有无限的资源可以比这里的任何答案更好地教你。请参阅 Hadley 的 Advanced R（可在线获取）或 bookdown 库中的许多示例
谢谢参考，我会研究的。您介意分享一下您将如何使用 lapply 处理此示例的两个步骤（即 group_by() 步骤和迭代 lm() 的步骤）吗？

标签： r dplyr tidy broom

【解决方案1】：

在data.table 中，您可以通过groupingvar 和结果变量melt（重塑长-将结果变量堆叠在一列而不是存储在三列中）和lm：

library(data.table)
setDT(df)

#alternatively, set id.vars = c('groupingvar', paste0('x', 1:6)), etc.
longDT = melt(df, id.vars = grep('y', names(df), invert = TRUE))

#this helper function basically splits a named vector into
#  its two components
coefsplit = function(reg) {
  beta = coef(reg)
  list(var = names(beta), coef = beta)
}

#I personally wouldn't assign longDT, I'd just chain this onto
#  the output of melt;
longDT[ , coefsplit(lm(value ~ ., data = .SD)), by = .(groupingvar, variable)]
#     groupingvar variable         var          coef
#  1:           a       y1 (Intercept) -3.595564e+03
#  2:           a       y1          x1 -3.796627e+01
#  3:           a       y1          x2 -1.557268e+02
#  4:           a       y1          x3  2.862738e+02
#  5:           a       y1          x4  1.579548e+02
# ...
# 38:           b       y3          x2  2.136253e+01
# 39:           b       y3          x3 -3.810176e+01
# 40:           b       y3          x4  4.187719e+01
# 41:           b       y3          x5 -2.586184e+02
# 42:           b       y3          x6  1.181879e+02
#     groupingvar variable         var          coef

【讨论】：

谢谢。我想我遵循你的回答。我能够在 R 中运行它并且它可以工作。一些注意事项/问题：（1）我还尝试使用 linest() 函数在 Excel 中复制结果，结果仅与第一个分组变量匹配，而对于第二个分组变量它们有显着不同。 R和Excel的最小二乘算法有区别吗？ (2) 我同意我不应该分配 longDT，因为行数会显着增长，但是你将如何管道化的输出，以便每个系数也正确地标记为截距、x1 等？
结果应该是一样的（大概都是用(X'X)^(-1)X'Y）。您能否具体说明哪个组返回了错误的估计值，以及 R 与 Excel 的输出是什么？
让我们设置种子，例如使用 set.seed(100)。那么对于分组变量a，y1的截距为0.2256。这与 excel 的 linet() 函数的输出相匹配。但是，对于分组变量 b，y1 的截距为 -3.3755，而 excel 的输出为 6.8732。
我发现了问题：n=15为b分组变量产生了6行，与自变量个数相同。如果我将 n=20 与 set.seed(100) 一起使用，那么对于 R 和 Excel 组 a 和 b 的 y1 的截距分别为 1.05464 和 5.75549。我会在线程中接受您的解决方案。谢谢！
啊，对。好吧，我猜 R 会自动删除变量。 Excel 可能以不同的方式做同样的事情，因此存在分歧。

【解决方案2】：

我还找到了一种使用 cbind() 实现此目的的方法，如下所示：

library(tidyverse)
library(broom)
n  <- 20
df4  <- data.frame(groupingvar= sample(1:2, size = n, replace = TRUE),
                   y1 = rnorm(n,10,1), y2=rnorm(n,100,10), y3=rnorm(n,1000,100),
                   x1=  rnorm(n,10,1), x2=rnorm(n,10,1), x3=rnorm(n,10,1),
                   x4=rnorm(n,10,1), x5=rnorm(n,10,1), x6=rnorm(n,10,1))
df4 <- arrange(df4,groupingvar)

regs <- df4 %>% group_by(groupingvar) %>%
  do(fit = lm(cbind(y1,y2,y3) ~ . -groupingvar, data = .))
coeffs <- tidy(regs, fit)

【讨论】：