如何有效地实现 dplyr do call for lmer 函数？答案

【问题标题】：How to efficiently implement dplyr do call for lmer function?如何有效地实现 dplyr do call for lmer 函数？
【发布时间】：2017-08-22 17:56:48
【问题描述】：

我有一个包含 ~400000 行的数据集，我正在尝试使用 R 中的 dplyr do 调用来提取 lme4 混合模型方差分量。函数是：

myfunc <- function(dat) {
    if (sum(!is.na(dat$value)) > 840) {  # >70% data present 
           v = data.frame(VarCorr(lmer(value ~ 0 + (1|gid) + (1|trial:rep) + (1|trial:rep:block), data=dat)))
           data.frame(a=round(v[1,4]/(v[1,4]+(v[4,4]/2)),2), b=round(v[1,4],2), c=round(v[4,4],2), n_obs=nrow(dat), na_obs=sum(is.na(dat$value))) 
    } else { 
        data.frame(a=NA, b=NA, c=NA, n_obs= nrow(dat), na_obs=sum(is.na(dat$value)))
    }
}

在按四个分组变量对数据进行分组后，使用dplyr do 调用调用此函数。最后的dplyr 调用是：

system.time(out <- tst %>% group_by(iyear,ilocation,trait_id,date) %>% 
          do(myfunc(.)))

现在，当此代码在 11000 行的较小测试数据帧上运行时，大约需要 25 秒。但是在完整的 443K 行上运行它大约需要 8-9 小时才能完成，这非常慢。很明显，有一部分代码会降低性能，但我似乎无法弄清楚是lmer 部分还是dplyr 导致了性能下降。我感觉函数处理矢量化操作的方式有问题，但不确定。我尝试在函数调用之外初始化“out”矩阵，但并没有提高性能。
不幸的是，我没有更小的可重现数据集可供分享。但想听听您对如何使这段代码更高效的想法。

【问题讨论】：

速度慢的是lmer部分，它适合复杂的模型。 dplyr 所做的只是给它正确的数据。您有一些低效率的问题，每次迭代都会增加几毫秒（为什么要让v 成为data.frame?? 将其保留为矩阵），但这可以忽略不计。您正在拟合的某些模型可能信号较弱并且需要很长时间才能收敛。看看?lmerControl，也许您可以增加容差并减少最大迭代次数以加快速度。
在您的函数中使用打印语句来跟踪您正在进行的迭代，或在lmer() 中设置verbose = 1 或verbose = 2 以查看步骤级诊断。
如果您的个人数据集的结构足够相似，您可能能够做一些更优雅的事情，但需要更多的工作/更深入地了解 @ 的内部结构987654339@（例如，请参阅?refit ...如果您通过并行化有足够的解决方案，我可能不会打扰。

标签： r dplyr lme4

【解决方案1】：

解决方案： parallel 包中的 mclapply 函数来救援。正如@gregor 正确指出的那样，可能是lmer 部分正在减慢速度。最后我最终并行化了函数调用：

myfunc <- function(i) {
     dat = tst[tst$comb==unique(tst$comb)[i],]  #comb is concatenated iyear,ilocation....columns
     if (sum(!is.na(dat$value)) > 840) {  # >70% data present per column
         v = data.frame(VarCorr(lmer(value ~ 0 + rand_factor + nested_random_factor), data=dat)))
         data.frame(trait=unique(tst$comb)[i], a=round(v[1,4])/5, b=round(v[1,4],2), c=round(v[4,4],2), n_obs=nrow(dat), na_obs=sum(is.na(dat$value))) 
     } else {
          data.frame(trait=unique(tst$comb)[i], a=NA, b=NA, c=NA, n_obs= nrow(dat), na_obs=sum(is.na(dat$value))) 
     }
}

#initialize an empty matrix
out <- matrix(NA,length(unique(tst$comb)),6)

## apply function in parallel. output is list
n_cores = detectCores() - 2
system.time(my.h2 <- mclapply(1:length(unique(tst$comb)),FUN = myfunc, mc.cores = n_cores))

一台十二核 unix 机器大约需要 2 分钟才能完成。

【讨论】：