使用 ForEach 逐步执行 1000 次回归的列

Question

首先是一些数据。为协变量和我感兴趣的回归结果制作一个数据框，为解释变量制作一个数据框。

我正在做的是单步执行 lm(outcome ~ mycovs + ith column of betas)，对于这个例子，收集残差。

set.seed(123) # for repeatability
mycovs = data.frame(outcome = rnorm(100,20,5),
                    race = rep(c("white","black","hispanic","other"),25),
                    income = rep(c("high","low"),50), 
                    age = rnorm(100,30,3))
betas = data.frame(replicate(10000,rnorm(100,50,6)/100))

要对 betas 中的每个变量执行此操作，我编写了以下代码：

get_resids <- function(x){
  mydata = cbind(mycovs,x)
  cpg = names(mydata)[ncol(mydata)]
  as.vector(resid(lm(formula(paste("outcome ~ as.factor(race) + as.factor(income) + age + ", cpg )),
                     data = mydata)))
  }
head(get_resids(betas[1]))
[1] -1.8525090 -0.7299173  6.4941289  0.5357159 -0.1771154  7.7554550

然后我可以使用 do.call(lapply()) 为我的 betas 数据框中的每个 10,000 个变量生成这些残差的矩阵，如下所示。

system.time(
myresids <- do.call(cbind, lapply(betas, get_resids))
)
   user  system elapsed 
  20.63    0.06   20.76 
> dim(myresids)
[1]   100 10000
> myresids[1:5,1:10]
             X1         X2          X3         X4          X5          X6         X7          X8         X9        X10
[1,] -1.8525090 -3.2651298 -3.54352587 -3.2962217 -2.95237520 -2.52995146 -3.0971490 -3.07625585 -2.8306409 -2.6454698
[2,] -0.7299173 -1.7982698 -2.54966496 -1.8009449 -1.60265484 -0.35825398 -1.6771846 -1.55455681 -1.2834764 -1.0941130
[3,]  6.4941289  6.6330879  5.88252329  7.1254892  6.88332171  7.79059098  6.9549380  6.84726299  6.9756743  6.3790811
[4,]  0.5357159 -0.0629098  0.06064112  0.3261975 -0.05377268 -0.04489599  0.1968423  0.02764062  0.2472463 -0.6944623
[5,] -0.1771154  0.1974865  0.56104333 -0.1188214  0.40202835  1.37694954  0.2904445  0.22634565  1.0650977  0.3231615

不错。我正在进行 10,000 次回归并将所有回归的残差存储在一个矩阵中，这需要 20 多秒的时间。请注意，这是一个单线程操作，按顺序执行 10,000 次回归。

好吧，这些暴露实际上是遗传 CpG 甲基化分数，我有大约一百万个要做，所以我想使用 foreach() 和 doParallel 来进行多线程处理，但我一直无法做到搞清楚。

这是我试过的。我首先将 betas 矩阵分解为 4 个命名数据帧，每个部分有 1/4 列：

mylist <- list(b1 = betas[1:2500], b2 = betas[2501:5000], b3 = betas[5001:7500], b4 = betas[7501:10000])
names(mylist); length(mylist)
[1] "b1" "b2" "b3" "b4"
[1] 4

然后我尝试按如下方式实现doParallel：

myresids_par <- foreach(i = 1:length(mylist), .combine = "cbind") %dopar% {
    do.call(cbind, lapply(mylist[i], get_resids))
  }
stopCluster(cl)

但我得到的是以下内容；只有 4 组残差如下，我不确定它做了什么：

> dim(myresids_par)
[1] 100   4
> head(myresids_par)
             b1         b2         b3          b4
[1,] -1.1051559 -3.2815443 -4.0951682 -2.97181934
[2,] -1.7884883 -1.5842009 -2.2403507 -1.48095064
[3,]  6.0211664  6.8417766  7.0208282  6.93438155
[4,] -0.4692244  0.1247481  0.9653631 -0.08206986
[5,] -0.1857339  0.2945526  1.8936715  0.30034781
[6,]  8.7706564  7.9744631  8.5240021  8.05232223

Answer 1

这里的问题是 mylist[i] 正在访问长度为 1 的 sub-list（不是存储在列表的 i-th 元素中的数据框；您需要 mylist[[i]] 代替）。

所以你可以使用：

myresids_par <- foreach(i = 1:length(mylist), .combine = "cbind") %dopar% {
  do.call(cbind, lapply(mylist[[i]], get_resids))
}

或更好，只需使用：

myresids_par <- foreach(i = seq_along(mylist), .combine = "c") %dopar% {
  lapply(mylist[[i]], get_resids)
}

如果你想要一个矩阵，然后使用 do.call(cbind, myresids_par)，或者如果你想要一个数据框，则只使用 as.data.frame(myresids_par)。

PS：注意这里的lapply是有效的，因为数据框也是一个列表。如果列表中有矩阵，则需要使用 apply(MAT, 2, FUN).

使用 ForEach 逐步执行 1000 次回归的列

Using ForEach to Step Through columns for 1000's of regressions

parallel-processing

foreach

r

doparallel