最近,我一直在尝试将许多随机效果模型拟合到相对较大的数据集。假设在多达25个时间点观察到约50,000人(或更多)。由于样本量如此之大,我们包含了许多我们要调整的预测因素-可能有50种左右的固定效果。我使用R中的lme4::glmer
将模型拟合为二进制结果,并对每个主题随机截取。我无法详细说明数据,但是我使用的glmer
命令的基本格式为:
fit <- glmer(outcome ~ treatment + study_quarter + dd_quarter + (1|id),
family = "binomial", data = dat)
其中
study_quarter
和dd_quarter
都是因子,每个因子约有20个水平。当我尝试在R中拟合此模型时,它会运行约12-15小时,并返回无法收敛的错误。我做了很多故障排除(例如,遵循these准则),但没有任何改善。而且收敛甚至还没有结束(最大梯度大约在5-10之间,而我认为收敛标准是0.001)。
然后,我尝试使用melogit命令在Stata中拟合模型。该模型在2分钟内拟合完毕,没有收敛问题。相应的Stata命令是
melogit outcome treatment i.study_quarter i.dd_quarter || id:
是什么赋予了? Stata只是具有更好的拟合算法,还是针对大型模型和大型数据集进行了更好的优化?令人惊讶的是运行时间有多么不同。
最佳答案
使用可选参数glmer
可以使nAGQ=0L
适合更快。您有许多固定效果参数(study_quarter
和dd_quarter
的每个20级,总共产生28个对比度),默认优化方法(对应于nAGQ=1L
)将所有这些系数放入常规的非线性优化调用中。使用nAGQ=0L
,可以在更快的惩罚性迭代重加权最小二乘(PIRLS)算法中优化这些系数。从某种意义上说,缺省值通常会提供更好的估计值,因为该估计值的偏差较小,但是该差值通常很小,并且时差很大。
我将这些算法的差异写成 Jupyter
笔记本 nAGQ.ipynb
。该写操作使用 MixedModels
包而不是 Julia
来使用 lme4
包,但是方法相似。 (我是lme4
的作者之一和MixedModels
的作者。)
如果您打算进行很多GLMM拟合,我会考虑在Julia
和MixedModels
中进行。即使使用R
中的所有复杂代码,它通常也比lme4
快得多。
关于r - lme4::glmer vs.Stata的melogit命令,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/44677487/