用于回归lm(y~x)的R for循环

Question

例子：

df <- data.frame(A=1:5, B=2:6, C=3:7,D=4:8,E=5:9,F=6:10)

我希望能够使用回归循环lm(y,x)。其中，y代表第1列和第2列，x代表其余的列。

我的想法是：

lmf <- function (y,x) {
                         f <- lm(y ~ x, data=df)
                         cbind(summary(f)$r.squared,summary(f)$coefficients)                  
                        }
 for(y in 1:3)
  {
    R<- apply(df[,3:6], 2, lmf(y,x)); R
  }

错误：在model.frame.default（formula = y〜x，data = df，drop.unused.levels = TRUE）中出现错误：变量长度不同（发现'x'）

我给出的示例非常小，但我的数据有50列用于y和300列用于x。

我的要求是相同的：lm（df $1〜df $3，data = df）; lm（df $1〜df $4，data = df），[...]，lm（df $2〜df $3，data = df）...但以自动方式执行。此外，我想提取结果$coefficients和$r.squared。

- Giffredo

---

你离正确的方法还有很远的距离。你没有正确使用apply，也没有使用正确的函数参数，也没有正确使用for循环迭代器。从简单的开始。 - Señor O

尝试一下这个 lapply(2:6, function(x) lm(df[ , 1] ~ df[, x], data = df ) )，虽然子集有点棘手（它是一个嵌套列表），但它可以完成任务。 - SabDeM

感谢回复，代码lapply(2:6, function(x) lm(df[ , 1] ~ df[, x], data = df ) )给我提供了使用固定y（在这种情况下为列df [，1]）和不同x的回归函数。我已经达到了这个目标... 我想要更多的东西，添加使用不同y的可能性。 - Giffredo

Answer 1

我有一个使用dplyr、tidyr和broom包的替代版本。思路是指定你想要作为Y和X处理的变量。基于这些Y和X集合创建2个不同的数据集。然后重塑数据集，以便能够将每个Y与一个X组合。最后，对于每个组合运行线性回归，并将模型输出保存为数据集。

# Check whether package name is installed...
check_package <- function(package_name) {
  if (!(package_name %in% rownames(installed.packages()))) {
    install.packages(package_name, dependencies = TRUE)
  }
}

check_package("broom")
check_package("dplyr")
check_package("tidyr")

library(dplyr)
library(broom)
library(tidyr)

# example dataset (picking 4 columns)
dt <- data.frame(mtcars) %>% select(mpg, disp, cyl, wt)

# specify which columns we want as y (dependent) and x (independent)
ynames <- c("disp","mpg")
xnames <- c("cyl","wt")

# create and reshape datasets
dt1 <- dt[,ynames]
dt1 <- gather(dt1,y,yvalue)

dt2 <- dt[,xnames]
dt2 <- gather(dt2, x, xvalue)



dt1 %>% 
  group_by(y) %>%                       # group by dependent variable
  do(data.frame(.,dt2)) %>%             # combine each y with all x
  group_by(y,x)%>%                      # get combinations of y and x to regress
  do(tidy(lm(yvalue~xvalue, data=.)))   # return lm output as dataframe


#      y   x        term    estimate  std.error statistic      p.value
# 1 disp cyl (Intercept) -156.608976 35.1805064 -4.451584 1.090157e-04
# 2 disp cyl      xvalue   62.598925  5.4693168 11.445474 1.802838e-12
# 3 disp  wt (Intercept) -131.148416 35.7165961 -3.671918 9.325668e-04
# 4 disp  wt      xvalue  112.478138 10.6353299 10.575896 1.222320e-11
# 5  mpg cyl (Intercept)   37.884576  2.0738436 18.267808 8.369155e-18
# 6  mpg cyl      xvalue   -2.875790  0.3224089 -8.919699 6.112687e-10
# 7  mpg  wt (Intercept)   37.285126  1.8776273 19.857575 8.241799e-19
# 8  mpg  wt      xvalue   -5.344472  0.5591010 -9.559044 1.293959e-10

- AntoniosK

---

看起来非常不错！！！但是我的R无法识别gather函数...也许是因为这个问题："> library(dplyr)附加包：‘dplyr’以下对象来自‘package:stats’：filter以下对象来自‘package:base’：intersect，setdiff，setequal，union" - Giffredo

Gather 函数是从 tidyr 包中的。你不应该遇到问题。你装了这个包吗？还有其他的方法来执行这一步骤。你可以使用其他重塑命令。如果你想要另一个选择，请告诉我。 - AntoniosK

@Giffredo，我更新了解决方案，使其安装包及其依赖项。 - Stereo

我有一个几乎完全相同的问题，只是我想在原始数据框中为一堆不同的组执行此操作。例如，在mtcars中每个gear的回归。我应该提出一个新问题，但如果很容易添加，group_by(gear)将无法工作。@Stereo 有什么想法？ - LauraR

@LauraR 如果您能提出一个新问题并说明您想要获取哪种组合，那将会更容易些。 - AntoniosK

Answer 2

我举个例子，以鸢尾花数据集中的数值型变量为例，但你可以根据需要更改为任何数据集。

我会根据我喜欢的名称构建公式，而不是使用数字来索引您感兴趣的列。

我建议，

 result <- sapply(names(iris)[1 : 4], 
   function(x) { 
     lapply(names(iris)[1 : 4], 
            function(y) {
              if (x != y) {
                model <- lm(as.formula(paste0(y, "~", x)), iris) 
                return(list(x = x, 
                            y = y, 
                            r.squared = summary(model)$r.squared, 
                            coefficients =  summary(model)$coefficients))
              }
              })
            })


 result
 ## Sepal.Length Sepal.Width Petal.Length Petal.Width
 ## [1,] NULL         List,4      List,4       List,4     
 ## [2,] List,4       NULL        List,4       List,4     
 ## [3,] List,4       List,4      NULL         List,4     
 ## [4,] List,4       List,4      List,4       NULL       

 result[1, 2]
 ## $Sepal.Width
 ## $Sepal.Width$x
 ## [1] "Sepal.Width"
 ## 
 ## $Sepal.Width$y
 ## [1] "Sepal.Length"
 ## 
 ## $Sepal.Width$r.squared
 ## [1] 0.01382265
 ## 
 ## $Sepal.Width$coefficients
 ## Estimate Std. Error   t value     Pr(>|t|)
 ## (Intercept)  6.5262226  0.4788963 13.627631 6.469702e-28
 ## Sepal.Width -0.2233611  0.1550809 -1.440287 1.518983e-01

或者，您可以将结果存储在列表中，并编写一个单独的函数来遍历该列表以创建仅包含您感兴趣的信息的矩阵。

- Stereo

---

这个很好用！我调整了返回部分以适应我的需求，现在没问题了！谢谢。 - Giffredo