拟合线性混合模型到一个非常大的数据集中

Question

拟合线性混合模型到一个非常大的数据集中

rparallel-processingbigdatalme4mixed-models

9

我想使用 lme4::lmer 在以下格式的 6000 万个观测值上运行混合模型；所有预测/依赖变量都是分类变量（因子），除了连续依赖变量 tc；patient 是随机拦截项的分组变量。我有 64 位 R 和 16Gb RAM，并且从中央服务器工作。RStudio 是最新的服务器版本。

model <- lmer(tc~sex+age+lho+atc+(1|patient),
              data=master,REML=TRUE)

lho sex tc      age         atc patient
18  M   16.61   45-54       H   628143
7   F   10.52   12-15       G   2013855
30  M   92.73   35-44       N   2657693
19  M   24.92   70-74       G   2420965
12  F   17.44   65-69       A   2833610
31  F   7.03    75 and over A   1090322
3   F   28.59   70-74       A   2718649
29  F   4.09    75 and over C   384578
16  F   67.22   65-69       R   1579355
23  F   7.7     70-74       C   896374

我遇到了一个错误信息：无法分配25.5Gb的向量。我在服务器上被分配了40Gb，目前已经使用了25Gb，所以我需要额外10Gb的空间。但我不认为我能获得额外的空间。

我对并行处理一窍不通，除了现在正在使用其中四个核心之外一无所知。有人能为这个模型提供并行代码或其他解决方法吗？

- steve

尝试使用Doug Bates的MixedModels包/库来进行Julia编程？ - Ben Bolker

首先，“分配向量”意味着RAM分配，因此服务器磁盘内存无关紧要。接下来，你真的需要花时间学习诸如“parallel”、“snow”和“bigdata”等软件包。否则，即使有人为你编写了工具，你也不会理解如何修改它或找到速度瓶颈。 - Carl Witthoft

谢谢Ben，我会看一下的。 - steve

3

时间不允许我等待，卡尔，这就是为什么我要发布这个问题的原因。 - steve

看一下你的设计矩阵。我猜它非常庞大。你需要认真考虑你要拟合的模型是否适当。例如，你是否真的需要每个“atc”水平？你能否将“atc”建模为随机效应？或者可能聚合atc级别？此外，你已经检查过你的数字变量在R中是否是“numeric”吗？无论如何，并行化不会解决内存问题。相反，由于每个CPU都需要RAM来完成其任务，因此并行化会加剧问题。 - Roland

谢谢Roland。我在atc上不会失去任何东西。没有想到RAM的问题。 - steve

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ben Bolker · Accepted Answer

如卡尔·威特福特所指出的那样，R中的标准并行化工具使用共享内存模型，因此它们会使情况变得更糟（它们的主要目的是通过使用多个处理器加速计算密集型作业）。

在短期内，您可以通过将分类固定效应预测器（`age`，`atc`）视为随机效应，但强制它们的方差很大来节省一些内存。我不知道这是否足以拯救您；它将大大压缩固定效应模型矩阵，但模型框架仍将与模型对象一起存储/复制...

dd1 <- read.table(header=TRUE,
text="lho sex tc      age         atc patient
18  M   16.61   45-54       H   628143
7   F   10.52   12-15       G   2013855
30  M   92.73   35-44       N   2657693
19  M   24.92   70-74       G   2420965
12  F   17.44   65-69       A   2833610
31  F   7.03    75_and_over A   1090322
3   F   28.59   70-74       A   2718649
29  F   4.09    75_and_over C   384578
16  F   67.22   65-69       R   1579355
23  F   7.7     70-74       C   896374")
n <- 1e5
set.seed(101)
dd2 <- with(dd1,
      data.frame(tc=rnorm(n,mean=mean(tc),sd=sd(tc)),
                 lho=round(runif(n,min=min(lho),max=max(lho))),
                 sex=sample(levels(sex),size=n,replace=TRUE),
                 age=sample(levels(age),size=n,replace=TRUE),
                 atc=sample(levels(atc),size=n,replace=TRUE),
                 patient=sample(1:1000,size=n,replace=TRUE)))
library("lme4")
m1 <- lmer(tc~sex+(1|lho)+(1|age)+(1|atc)+(1|patient),
           data=dd2,REML=TRUE)

随机效应会自动按照水平数量从大到小排序。根据 ?modular 帮助页面中的描述，进行如下操作：

lmod <- lFormula(tc~sex+(1|lho)+(1|age)+(1|atc)+(1|patient),
                  data=dd2,REML=TRUE)
names(lmod$reTrms$cnms)  ## ordering
devfun <- do.call(mkLmerDevfun, lmod)
wrapfun <- function(tt,bigsd=1000) {
    devfun(c(tt,rep(bigsd,3)))
}
wrapfun(1)
opt <- optim(fn=wrapfun,par=1,method="Brent",lower=0,upper=1000)
opt$fval <- opt$value  ## rename/copy
res <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr=lmod$fr)
res

您可以忽略分类项的报告差异，使用ranef()来恢复它们（未收缩）的估计值。

从长远来看，解决这个问题的正确方法可能是使用分布式内存模型并行化。换句话说，您需要将数据分成不同的块分配到不同的服务器上；使用?modular中描述的机制设置一个似然函数（实际上是一个REML准则函数），该函数将给出作为参数的子集的REML准则函数；然后运行一个中央优化器，通过将一组参数提交给每个服务器、检索每个服务器的值并相加来评估REML准则函数。我看到实现这个有两个问题：(1)我实际上不知道如何在R中实现分布式内存计算（基于这份介绍文档，似乎Rmpi/doMPI包可能是正确的选择）；(2)在默认情况下，lmer的实现方式是剖析掉固定效应参数，而不是显式地将其作为参数向量的一部分。