加速 R 中的 lmer 函数答案

【问题标题】：Speed up lmer function in R加速 R 中的 lmer 函数
【发布时间】：2015-11-17 14:17:52
【问题描述】：

我想分享一些我在尝试使用lme4 包改善R 中线性混合效应模型的模型拟合时间时的一些想法。

数据集大小：数据集大约由 400.000 行和 32 列组成。遗憾的是，无法共享有关数据性质的任何信息。

假设和检查：假设响应变量来自正态分布。在模型拟合过程之前，使用相关表和 R 中提供的alias 函数测试变量的共线性和多重共线性。

为了帮助收敛，对连续变量进行了缩放。

模型结构：模型方程包含31个固定效应（包括截距）和30个随机效应（不包括截距）。对于具有 2700 个水平的特定因子变量，随机效应是随机的。协方差结构是方差分量，因为它假设随机效应之间存在独立性。

模型方程示例：

lmer(Response ~ 1 + Var1 + Var2 + ... + Var30 + (Var1-1| Group) + (Var2-1| Group) + ... + (Var30-1| Group), data=data, REML=TRUE)

模型已成功拟合，但是大约需要 3.1 小时才能提供结果。 SAS 中的相同模型需要几秒钟。网上有关于如何使用非线性优化算法nloptwrap 并关闭优化完成后执行的耗时导数计算来减少时间的文献calc.derivs = FALSE：

https://cran.r-project.org/web/packages/lme4/vignettes/lmerperf.html

时间减少了 78%。

问题：是否有任何其他替代方法可以通过相应地定义lmer 参数输入来减少模型拟合时间？ R 和 SAS 在模型拟合时间方面存在很大差异。

欢迎提出任何建议。

【问题讨论】：

我的第一步是对数据集（行）进行二次抽样并进行拟合。
嗨@Roman Luštrik。我需要使用整个数据集来获得系数估计值。
我认为lmer 和 SAS proc 混合中使用的算法是不同的，尽管您必须让 @benbolker 确认。您还可以检查source code
几秒钟与 0.75 小时（在您完成优化之后）确实非常令人惊讶。 SAS 和lmer 结果是否相互一致，即你真的适合同一个模型吗？所有Varx 变量都是数字的吗？您能否提供一个可重现的示例（我知道您的数据是机密的，但您真正需要做的就是发布模拟代码，生成与您的数据大约相同的结构......）
嗨，本博克。结果完全相同，变量也是数字。我将在一周后回来，所以我很乐意分享大致相同的数据结构。感谢您的宝贵时间。

标签： r performance lme4 mixed-models

【解决方案1】：

lmer() 通过针对随机效应的协方差矩阵中的参数优化配置的对数似然或配置的 REML 标准来确定参数估计。在您的示例中，将有 31 个这样的参数，对应于 31 个项中每个项的随机效应的标准差。这种规模的受限优化需要时间。

SAS PROC MIXED 可能具有特定的优化方法或更复杂的方法来确定初始估计值。 SAS 是一个闭源系统意味着我们不知道他们在做什么。

顺便说一下，你可以把随机效应写成(1+Var1+Var2+...+Var30||Group)

【讨论】：

对于这些类型的模型，可以更快地评估 lmer 中的异形对数似然评估。当我们编写 lme4 时，我们在每次迭代中求解了固定效应系数和随机效应的模式。这涉及对具有 n 行和长度为 n 的向量的矩阵进行操作，其中 n 是观察数。最近我意识到这没有必要。我有修改后算法的Julia 实现，但 R 中还没有。
我是 stackoverflow 的新手，显然错过了一些微妙之处。我通常认为自己是 Doug Bates，而不是 user1864481
亲爱的道格拉斯。 E. Bates，你的回答令人惊讶。我设法在Julia 中安装了相同的模型。结果表明，在拟合特定模型时，Julia 比R 快大约 74 倍。这是性能上的巨大差异。此外，结果是一致的。有没有计划在不久的将来在R lme4 进行同样的开发？谢谢您的帮助。康斯坦丁诺斯妈妈
我目前有两个 Julia 包用于拟合线性混合效应模型 - MixedModels 中的旧算法和 ReTerms 中的新算法。该计划是将 ReTerms 折叠到 MixedModels 中并记录新算法。我目前不打算在lme4 中实现新算法。我认为，一些关键部分很难在 R/Rcpp 中实现。
这很清楚，也为我打开了一个全新的任务；是否使用Julia 而不是R 以拟合广义线性模型或线性混合效应模型。非常感谢您的宝贵帮助。我认为这次对话会让很多人重新考虑使用哪种统计语言来拟合混合模型。

【解决方案2】：

我们已经在 R 包 Rfast 中实现了假设复合对称的随机截距回归。命令是 rint.reg。它比相应的 lme4 函数快 30 倍以上。我不知道这是否有帮助，但以防万一。

https://cran.r-project.org/web/packages/Rfast/index.html

【讨论】：

【解决方案3】：

如果您使用glmer 而不是lmer，则有一个参数nAGQ。我发现设置nAGQ=0 显着减少了拟合相当复杂模型所需的时间（13 个固定效应，一个具有不同截距和斜率的随机效应，300k 行）。这基本上告诉glmer 对 GLMM 使用不太精确的参数估计形式。详情请参阅?glmer，或this 发帖。

【讨论】：

感谢您的意见，我一定会记住的！