在 R 中跨数据框的几列运行多重线性回归

Question

我的数据集结构如下： enter image description here

我想使用 V1、V2...等运行线性回归模型和方差分析。作为自变量，g 列作为每种情况下的因变量（即 lm(V1 ~ g)、lm(V2 ~ g) 等）。这很简单，只是这些线性回归需要在 pair 列中按级别分组，例如，我的输出包含对 1.1 的所有行的 lm(V1 ~ g) 和所有行的 lm(V1 ~ g)对 1.201 等

我已经尝试了多种使用 for 循环、lapply 和 data.table 包的方法，但没有什么能提供我想要的输出。谁能告诉我解决这个问题的最佳方法？

编辑： 我的完整数据集在 pair 列和 100 V 列（V1...V100）中有 7056 个不同的对。我对这个问题的最新尝试：

df$pair <- as.factor(df$pair)
out <- list()
for (i in 3:ncol(df)){
    out[[i]] <- lapply(levels(df$pair), function(x) {
    data.frame(df=x, g = coef(summary(lm(df[,i]~ df$g, data=df[df$pair==x,])),row.names=NULL))})
    }

Answer 1

让我们在这里获得一些tidyverse 的力量，以及broom，并放弃所有这些循环......

首先我会做一个虚拟表：

df <- data.frame(
  g = runif(50), 
  pair = sample(x = c("A", "B", "C"), size = 50, replace = TRUE), 
  V1 = runif(50), 
  V2 = runif(50), 
  V3 = runif(50), 
  V4 = runif(50), 
  V5 = runif(50),
  stringsAsFactors = FALSE
)

这大概是您的数据结构的样子。现在进入代码的核心：

library(tidyverse)
library(broom)

df %>% 
  as_tibble %>% 
  gather(key = "column", value = "value", V1:V5) %>%       # first set the data in long format
  nest(g, value) %>%                                       # now nest the dependent and independent factors
  mutate(model = map(data, ~lm(g ~ value, data = .))) %>%  # fit the model using purrr
  mutate(tidy_model = map(model, tidy)) %>%                # clean the model output with broom
  select(-data, -model) %>%                                # remove the "untidy" parts
  unnest()                                                 # get it back in a recognizable data frame

这给了我们以下信息：

# A tibble: 30 x 7
   pair  column term        estimate std.error statistic  p.value
   <chr> <chr>  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
 1 C     V1     (Intercept)  0.470       0.142    3.31   0.00561 
 2 C     V1     value        0.125       0.265    0.472  0.645   
 3 B     V1     (Intercept)  0.489       0.142    3.45   0.00359 
 4 B     V1     value       -0.0438      0.289   -0.151  0.882   
 5 A     V1     (Intercept)  0.515       0.111    4.63   0.000279
 6 A     V1     value       -0.00569     0.249   -0.0229 0.982   
 7 C     V2     (Intercept)  0.367       0.147    2.50   0.0265  
 8 C     V2     value        0.377       0.300    1.26   0.231   
 9 B     V2     (Intercept)  0.462       0.179    2.59   0.0206  
10 B     V2     value        0.0175      0.322    0.0545 0.957   
# … with 20 more rows

是的，这很漂亮！请注意，我使用了lm(g ~ value) 而不是lm(value ~ g)，因为这是您的文字描述所暗示的。

Answer 2

使用tidyverse 包过滤您的数据框：

library(tidyverse)

lm(V1~g, data=filter(yourData, pair==1.1))
lm(V2~g, data=filter(yourData, pair==1.201))

这可确保您为每个回归模型消除不包含所需 pair 值的行。您可能可以创建一个循环来执行此操作，但我认为手动继续过滤 pair 值更容易。如果你真的想使用循环，这里有一个相当简单的方法：

for (i in levels(yourData$pair)) {
  if (i==1.1) {
    mod1 <- lm(V1~g, data=filter(yourData, pair==i))
  }

  if (i==1.201) {
    mod2 <- lm(V2~g, data=filter(yourData, pair==i))
  }
}

但这仍然是手动循环通过pair 的级别。我必须查看您的整个数据集才能自动执行循环过程。

另外，如果g 列包含您的依赖值，则调用应该是lm(g~V1)、lm(g~V2) 等。它不应该是lm(V1~g)。