Jon*_*lar 8 r lme4 stata mixed-models
最近我一直试图将大量随机效应模型适用于相对较大的数据集.假设在最多25个时间点观察到大约50,000人(或更多).如此大的样本量,我们包含了许多我们正在调整的预测因子 - 可能有50个左右的固定效应.我lme4::glmer在R中使用模型拟合二元结果,每个主题都有随机截距.我不能详细介绍数据,但glmer我使用的命令的基本格式是:
fit <- glmer(outcome ~ treatment + study_quarter + dd_quarter + (1|id),
family = "binomial", data = dat)
其中两个study_quarter和dd_quarter与各约20个级别的因素.
当我尝试在R中使用此模型时,它会运行大约12-15个小时,并返回一个无法收敛的错误.我做了一堆故障排除(例如,遵循这些指导原则),没有任何改进.并且最终收敛甚至不接近(最大梯度在5-10附近,而收敛标准是0.001我认为).
然后我尝试使用melogit命令在Stata中拟合模型.该模型适合在2分钟内完成,没有收敛问题.相应的Stata命令是
melogit outcome treatment i.study_quarter i.dd_quarter || id:
是什么赋予了？Stata是否只有更好的拟合算法,或者更好地针对大型模型和大型数据集进行了优化？令人惊讶的是,运行时间有多么不同.
glmer使用可选参数,拟合可能会快得多nAGQ=0L.您有许多固定效应参数(每个参数有20个级别study_quarter并且dd_quarter总共产生28个对比度),默认优化方法(对应于nAGQ=1L)将所有这些系数放入一般非线性优化调用中.利用nAGQ=0L这些系数在快得多的惩罚迭代重加权最小二乘(PIRLS)算法中进行优化.默认值通常会提供更好的估计,即估计值的偏差较低,但差异通常非常小,时间差异很大.
我把这些算法的差异写成Jupyter笔记本nAGQ.ipynb.该writeup使用MixedModels包Julia而不是lme4方法,但方法是相似的.(我是其中一位作者lme4和作者MixedModels.)
如果你打算做很多GLMM的适合我会考虑这样做Julia有MixedModels.R即使使用了所有复杂的代码,它通常也要快得多lme4.