将线性混合模型拟合到非常大的数据集

Question

我想对以下格式的 60M 观测值运行混合模型（使用

lme4::lmer

）；除了连续因变量之外，所有预测变量/因变量都是分类变量

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

我收到

cannot allocate a vector of 25.5Gb

错误。我在服务器上分配了 40Gb，而我正在使用 25Gb，所以我想这意味着我还需要 10Gb 左右。我认为我无法分配任何额外的空间。

除了我目前正在使用四个核心之一之外，我对并行处理一无所知。任何人都可以建议该模型的并行代码，或者可能是不同的修复程序？

Answer 1

正如 Carl Witthoft 所指出的，R 中的标准并行化工具通过为每个工作人员制作数据副本来工作（我最初称其为 共享内存模型，但我认为这不准确），因此运行它们在具有单个内存块的系统上（与例如“简单的工作站网络”（SNOW）相反，其中每个进程位于不同的机器上），它们将使事情变得更糟而不是更好（它们的主要目的是加速计算密集使用多个处理器的作业）。

短期内，您可以通过将分类固定效应预测变量（

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 中实现分布式内存计算（基于 

this 介绍文档，似乎 Rmpi/doMPI 包可能是正确的方法）； (2) 在实现 lmer

 的默认方式中，固定效应参数被分析出来，而不是明确地成为参数向量的一部分。