循环回归模型项的组合

Question

我正在运行以

的形式进行回归

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

代替最后一项 z1，我想遍历一组不同的变量，z1 到 z10，运行进行回归对于每个以它作为最后一项。例如。第二运行我想使用

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

第 3 运行:

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

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

Answer 1

有了这个虚拟数据：

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 调用中。

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

Answer 2

虽然 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 3

这是使用 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 值占据同一列。 use_series() 从 magrittr package return the list of lm objects instead of a data.frame. If you load the 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))

Answer 4

根据您的最终目标，可以更快拟合基本模型，用 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 以重新拟合模型。

循环回归模型项的组合

Looping over combinations of regression model terms

r

apply

lm