循环回归模型项的组合

Question

我正在以表格形式进行回归

reg=lm(y ~ x1+x2+x3+z1,data=mydata)

在最后一个术语 z1 的位置上，我想遍历一组不同的变量，z1 到 z10，对每个变量进行回归，并将其作为最后一个术语。例如。在第二次运行中我想使用

reg=lm(y ~ x1+x2+x3+z2,data=mydata)

第三次运行：

reg=lm(y ~ x1+x2+x3+z3,data=mydata)

如何通过循环遍历 z 变量列表来自动执行此操作？

Answer 1

虽然 Sam 提供的解决方案很有效，而且是一个很好的解决方案，但我个人更愿意稍微不同。他的回答已经被接受了，所以为了完整起见，我只是发布这个。

dat1 <- data.frame(y = rpois(100, 5),
                   x1 = runif(100),
                   x2 = runif(100),
                   x3 = runif(100),
                   z1 = runif(100),
                   z2 = runif(100))

lapply(colnames(dat1)[5:6],
       function(x, d) lm(as.formula(paste("y ~ x1 + x2 + x3", x, sep = " + ")), data = d),
       d = dat1)

不是循环遍历数据框的实际列，而是仅循环遍历名称字符串。这提供了一些速度改进，因为在迭代之间复制的东西更少。

library(microbenchmark)

microbenchmark({ lapply(<what I wrote above>) })
# Unit: milliseconds
# expr
# {lapply(colnames(dat1)[5:6], function(x, d) lm(as.formula(paste("y ~ x1 + x2 + x3", x, sep = "+")), data = d), d = dat1)}
#       min       lq     mean   median       uq      max neval
#  4.014237 4.148117 4.323387 4.220189 4.281995 5.898811   100

microbenchmark({ lapply(<other answer>) })
# Unit: milliseconds
# expr
# {lapply(dat1[, 5:6], function(x) lm(dat1$y ~ dat1$x1 + dat1$x2 + dat1$x3 + x))}
#       min       lq     mean   median       uq    max neval
#  4.391494 4.505056 5.186972 4.598301 4.698818 51.573   100

对于这个玩具示例，差异相当小，但随着观察和预测变量数量的增加，差异可能会变得更加明显。

Answer 2

根据您的最终目标，拟合基本模型、使用 add1 更新它并提取您想要的 F-test/AIC 可以快得多：

> basemodel <- lm(y~x1+x2+x3, dat1)
> 
> add1(object=basemodel, grep("z\\d", names(dat1), value=TRUE), test="F")
Single term additions

Model:
y ~ x1 + x2 + x3
       Df Sum of Sq    RSS    AIC F value Pr(>F)
<none>              477.34 164.31               
z1      1    0.0768 477.26 166.29  0.0153 0.9019
z2      1    5.1937 472.15 165.21  1.0450 0.3093

另请参阅?update 以改装模型。

Answer 3

有了这个虚拟数据：

dat1 <- data.frame(y = rpois(100,5),
x1 = runif(100),
x2 = runif(100),
x3 = runif(100),
z1 = runif(100),
z2 = runif(100)
)

您可以通过这种方式获取两个 lm 对象的列表：

 lapply(dat1[5:6], function(x) lm(dat1$y ~ dat1$x1 + dat1$x2 + dat1$x3 + x))

遍历这两列并将它们作为参数替换为lm 调用。

正如亚历克斯在下面所说，最好通过公式传递名称，而不是像我在这里所做的那样传递实际的数据列。

Answer 4

这是使用 dplyr / tidyr 系列中的包的另一种方法。它将数据重组为长格式，然后使用 dplyr 包中的group_by() 而不是lapply()：

library(dplyr)
library(tidyr)
library(magrittr) # for use_series ()
dat1 %>%
  gather(varname, z, z1:z2) %>% # convert data to long form
  group_by(varname) %>%
  do(model = lm(y ~ x1 + x2 + x3 + z, data = .)) %>%
  use_series(model)

这会将数据转换为使用gather 的长格式，其中z 值占据同一列。 magrittr 包中的 use_series() 返回 lm 对象列表，而不是 data.frame。如果你加载broom包，你可以在这段代码中提取模型系数：

library(broom)
dat1 %>%
  gather(varname, z, z1:z2) %>%
  group_by(varname) %>%
  do(model = lm(y ~ x1 + x2 + x3 + z, data = .)) %>%
  glance(model) # or tidy(model)

Source: local data frame [2 x 12]
Groups: varname

  varname  r.squared adj.r.squared    sigma statistic   p.value df    logLik      AIC      BIC deviance df.residual
1      z1 0.06606736    0.02674388 2.075924  1.680099 0.1609905  5 -212.3698 436.7396 452.3707 409.3987          95
2      z2 0.06518852    0.02582804 2.076900  1.656192 0.1666479  5 -212.4168 436.8337 452.4647 409.7840          95

数据：

dat1 <- data.frame(y = rpois(100, 5), x1 = runif(100),
                   x2 = runif(100), x3 = runif(100),
                   z1 = runif(100), z2 = runif(100))