r - 在两个变量列表上进行循环迭代，以实现 R 中的多重回归

Question

我想在 R 中编写一个循环，以使用一个因变量和两个自变量列表(所有连续变量)运行多重回归。该模型是可加的，循环应通过迭代两个变量列表来运行，以便它从第一个列表中获取第一列+从第二个列表中获取第一列，然后对两个列表中的第二列进行相同的操作，依此类推。问题是我无法让它正确地迭代列表，而是我的循环运行了比应有的模型更多的模型。

我在这里描述的数据框只是一个子集，我实际上必须运行 3772 次(我正在研究 RNA-seq 转录本表达)。

我的数据框称为 dry，包含 22 个变量(列)和 87 个观察值(行)。第 1 列包含 genotypeID，第 2:11 列包含一组要迭代的自变量，第 12:21 列包含第二组要迭代的自变量，第 23 列包含我的因变量，称为 FITNESS_DRY。结构如下所示:

str(dry)
'data.frame':   87 obs. of  22 variables:
$ geneID     : Factor w/ 87 levels "e10","e101","e102",..: 12 15 17 24 25    30 35 36 38 39 ...
$ RDPI_T1    : num  1.671 -0.983 -0.776 -0.345 0.313 ...
$ RDPI_T2    : num  -0.976 -0.774 -0.532 -1.137 1.602 ...
$ RDPI_T3    : num  -0.197 -0.324 0.805 -0.701 -0.566 ...
$ RDPI_T4    : num  0.289 -0.92 1.117 -1.214 -0.447 ...
$ RDPI_T5    : num  -0.671 1.963 NA -1.024 -0.295 ...
$ RDPI_T6    : num  2.606 -1.116 -0.383 -0.893 0.119 ...
$ RDPI_T7    : num  -0.843 -0.229 -0.297 0.504 -0.712 ...
$ RDPI_T8    : num  -0.227 NA NA -0.816 -0.761 ...
$ RDPI_T9    : num  0.754 -1.304 1.867 -0.514 -1.377 ...
$ RDPI_T10   : num  1.1352 -0.1028 -0.69 2.0242 -0.0925 ...
$ DRY_T1     : num  0.6636 -0.64508 -0.24643 -1.43231 -0.00855 ...
$ DRY_T2     : num  1.008 0.823 -0.658 -0.148 0.272 ...
$ DRY_T3     : num  -0.518 -0.357 1.294 0.408 0.771 ...
$ DRY_T4     : num  0.0723 0.2834 0.5198 1.6527 0.4259 ...
$ DRY_T5     : num  0.1831 1.9984 NA 0.0923 0.1232 ...
$ DRY_T6     : num  -1.55 0.366 0.692 0.902 -0.993 ...
$ DRY_T7     : num  -2.483 -0.334 -1.077 -1.537 0.393 ...
$ DRY_T8     : num  0.396 NA NA -0.146 -0.468 ...
$ DRY_T9     : num  -0.694 0.353 2.384 0.665 0.937 ...
$ DRY_T10    : num  -1.24 -1.57 -1.36 -3.88 -1.4 ...
$ FITNESS_DRY: num  1.301 3.365 0.458 0.346 1.983 ...

目标是运行 10 个多重回归，如下所示:

lm1<-lm(FITNESS_DRY~DRY_T1+RDPI_T1)
lm2<-lm(FITNESS_DRY~DRY_T2+RDPI_T2)

依此类推，遍历两个列表的所有十列这在索引方面相当于以下内容

lm1<-lm(FITNESS_DRY~dry[,12]+dry[,2])
lm1<-lm(FITNESS_DRY~dry[,12]+dry[,2])

等等

然后，我的循环应计算每个模型的摘要，并将所有 pvalue(lm 摘要的第 4 列)合并到输出对象中。

我首先定义了我的变量列表

var_list<-list(
var1=dry[,12:21],
var2=dry[,2:11]
)

这是我尝试过的循环，但无法正常工作:

lm.test1<-name<-vector()
for (i in 12:length(var_list$var1)){
    for (j in 2:length(var_list$var2)){
lm.tmp<-lm(FITNESS_DRY~dry[,i]+dry[,j], na.action=na.omit, data=dry)
sum.tmp<-summary(lm.tmp)
lm.test1<-rbind(lm.test1,sum.tmp$coefficients[,4]) }
}

循环返回此错误消息:

Warning message:
In rbind(lm.test6, sum.tmp$coefficients[, 4]) :
number of columns of result is not a multiple of vector length (arg 2)

我可以调用对象“lm.test1”，但该对象有 27 行，而不是我想要的 10 行，因此迭代在这里无法正常工作。有人可以帮忙吗？另外，如果我可以将每个变量列表的列名称添加到摘要中，那就太好了。我尝试对每个变量列表使用它，但没有成功:

name<-append(name, as.character(colnames(var_list$var1))

有什么想法吗？预先感谢您的帮助!

更新1:有关完整数据集的更多信息:我的实际数据仍将包含第一列“geneID”，然后我有 3772 个列名称 DRY_T1...DRY_T3772，然后另外 3772 个列名称 RDPI_T1...RDPI_T3772，最后是我的因变量“FITNESS_DRY”。我仍然想运行所有附加模型:

lm1<-lm(FITNESS_DRY~DRY_T1+RDPI_T1)
lm2<-lm(FITNESS_DRY~DRY_T2+RDPI_T2)
lm3772<-lm(FITNESS_DRY~DRY_T3772+RDPI_T3772)

我模拟了这样的数据集:

set.seed(2)
dat3 = as.data.frame(replicate(7544, runif(20)))
names(dat3) = paste0(rep(c("DRY_T","RDPI_T"),each=3772), 1:3772)
dat3 = cbind(dat3, FITNESS_DRY=runif(20))

然后我运行 for 循环:

models = list()
for(i in 1:3772) {
vars = names(dat3)[grepl(paste0(i,"$"), names(dat3))]
models2[[as.character(i)]] = lm(paste("FITNESS_DRY ~ ", paste(vars, collapse=" 
+")),
                                 data = dat3)
}

这在数据模拟上工作得很好，但是当我在以完全相同的方式设置的真实数据集上尝试它时，它不起作用。该循环可能在处理两位或更多位数的数字时出现问题。我收到此错误消息:

Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
0 (non-NA) cases

更新 2:确实，该模型在处理两位或更多数字的数字时存在问题。为了看看原始版本中出了什么问题，我使用了以下命令:(我的数据集称为“dry2”):

names(dry2)[grepl("2$", names(dry2))]

这将返回所有数字包含“2”的 DRY_T 和 RDPI_T 变量，而不是仅一对 DRY_T 和 RDPI_T。

为了解决这个问题，这个新代码可以工作:

models = list()

for(i in 1:3772) {
vars = names(dry2)[names(dry2) %in% paste0(c("DRY_T", "RDPI_T"), i)]
models[[as.character(i)]] = lm(paste("FITNESS_DRY ~ ", paste(vars, collapse=" +   ")),
data = dry2)
}

Answer 1

有多种方法可以设置迭代模型公式。这是一种方法，我们演示了使用 purrr 包中的 for 循环或 map 进行迭代。然后我们使用 broom 包中的 tidy 来获取系数和 p 值。

library(tidyverse)
library(broom)

# Fake data
set.seed(2)
dat = as.data.frame(replicate(20, runif(20)))
names(dat) = paste0(rep(c("DRY_T","RDPI_T"),each=10), 0:9)
dat = cbind(dat, FITNESS_DRY=runif(20))

# Generate list of models

# Using for loop
models = list()

for(i in 0:9) {

  # Get the two column names to use for this iteration of the model
  vars = names(dat)[grepl(paste0(i,"$"), names(dat))]

  # Fit the model and add results to the output list
  models[[as.character(i)]] = lm(paste("FITNESS_DRY ~ ", paste(vars, collapse=" + ")),
                                 data = dat)
}

# Same idea using purrr::map to iterate
models = map(0:9 %>% set_names(), 
             ~ {
               vars = names(dat)[grepl(paste0(.x,"$"), names(dat))]
               form = paste("FITNESS_DRY ~ ", paste(vars, collapse=" + "))
               lm(form, data = dat)
             })

# Check first two models
models[1:2]
#> $`0`
#> 
#> Call:
#> lm(formula = form, data = dat)
#> 
#> Coefficients:
#> (Intercept)       DRY_T0      RDPI_T0  
#>      0.4543       0.3025      -0.1624  
#> 
#> 
#> $`1`
#> 
#> Call:
#> lm(formula = form, data = dat)
#> 
#> Coefficients:
#> (Intercept)       DRY_T1      RDPI_T1  
#>     0.64511     -0.33293      0.06698

# Get coefficients and p-values for each model in a single data frame
results = map_df(models, tidy, .id="run_number")

results
#> # A tibble: 30 x 6
#>    run_number term        estimate std.error statistic p.value
#>    <chr>      <chr>          <dbl>     <dbl>     <dbl>   <dbl>
#>  1 0          (Intercept)   0.454      0.153     2.96  0.00872
#>  2 0          DRY_T0        0.303      0.197     1.53  0.143  
#>  3 0          RDPI_T0      -0.162      0.186    -0.873 0.395  
#>  4 1          (Intercept)   0.645      0.185     3.49  0.00279
#>  5 1          DRY_T1       -0.333      0.204    -1.63  0.122  
#>  6 1          RDPI_T1       0.0670     0.236     0.284 0.780  
#>  7 2          (Intercept)   0.290      0.147     1.97  0.0650 
#>  8 2          DRY_T2        0.270      0.176     1.53  0.144  
#>  9 2          RDPI_T2       0.180      0.185     0.972 0.345  
#> 10 3          (Intercept)   0.273      0.187     1.46  0.162  
#> # … with 20 more rows

^{由reprex package于2019-06-28创建(v0.2.1)}

如果不需要保存模型对象，则可以只返回系数和 p 值的数据框:

results = map_df(0:9 %>% set_names(), 
            ~ {
              vars = names(dat)[grepl(paste0(.x,"$"), names(dat))]
              form = paste("FITNESS_DRY ~ ", paste(vars, collapse=" + "))
              tidy(lm(form, data = dat))
            }, .id="run_number")

更新:为了回答您的评论，如果您将 0:9 的所有实例替换为 1:10(抱歉，没有'请注意，您的列后缀是从 1:10 而不是 0:9)，以及带有 dry2 (或您的任何名称)的 dat (我的假数据)的所有实例用于您的数据框)，只要列名称与您在问题中使用的列名称相同，代码就会使用您的数据运行。如果您使用不同的列名称，则需要调整代码，方法是对新名称进行硬编码，或者创建一个可以接受您用于模型的任何列名称的函数。生成。

解释一下代码的作用:首先，我们需要获取要在模型的每次迭代中使用的列的名称。例如，在 for 循环版本中:

vars = names(dry2)[grepl(paste0(i,"$"), names(dry2))]

例如，当 i=2 时，解析为:

vars = names(dry2)[grepl("2$", names(dry2))]
vars

[1] "RDPI_T2" "DRY_T2"

所以这些是我们要用来生成回归公式的两列。 "2$" 是一个正则表达式(正则表达式是一种字符串匹配语言)，意思是:匹配 names(dry2) 中以数字“2”结尾的值。

为了创建我们的公式，我们这样做:

paste(vars, collapse=" + ")

[1] "RDPI_T2 + DRY_T2"

form = paste("FITNESS_DRY ~ ", paste(vars, collapse=" + "))
form

[1] "FITNESS_DRY ~  RDPI_T2 + DRY_T2"

现在我们有了在 lm 中使用的回归公式。

每次迭代(使用for或map，或者在@RomanLuštrik的建议中，mapply)都会生成连续的模型。

更新 2: 正如我在评论中指出的，我意识到正则表达式 paste(i, "$") 将失败(通过匹配多个每种类型的自变量列)当最终数字超过一位数时。因此，请尝试这样做(对于 map 版本也类似):

models = list()

for(i in 1:3772) {

  # Get the two column names to use for this iteration of the model
  vars = names(dry2)[names(dry2) %in% paste0(c("DRY_T", "RDPI_T"), i)]

  # Fit the model and add results to the output list
  models[[as.character(i)]] = lm(paste("FITNESS_DRY ~ ", paste(vars, collapse=" + ")),
                                 data = dry2)
}

要查看原始版本中出现的问题，请运行例如 names(dry2)[grepl("2$", names(dry2))]