lme4::glmer与Stata的melogit命令

Question

最近，我一直在尝试将许多随机效应模型拟合到相对较大的数据集中。假设大约有5万人(或更多)在25个时间点观察到。有了这么大的样本数量，我们就包括了很多我们正在调整的预测因子--也许是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个小时，并返回了一个错误，它无法收敛。我做了一堆故障排除(例如，遵循这些准则)，但没有改进。最终收敛甚至不接近(最大梯度在5-10左右，而收敛准则是0.001，我认为)。

然后，我尝试使用melogit命令在Stata中拟合模型。该模型适用于2分钟以下，不存在收敛问题。对应的Stata命令是

melogit outcome treatment i.study_quarter i.dd_quarter || id:

怎么回事？Stata只是有一个更好的拟合算法，还是对大型模型和大型数据集进行了更好的优化？运行时间有多么不同，实在令人惊讶。

Answer 1

使用可选参数glmer nAGQ=0L，fit可能要快得多。您有许多固定效果参数(每个study_quarter和dd_quarter生成总共28个对比的20个级别)，默认的优化方法(对应于nAGQ=1L)将所有这些系数放入一般的非线性优化调用中。使用nAGQ=0L，这些系数在更快的惩罚迭代重加权最小二乘(PIRLS)算法中得到优化。缺省值通常会提供一个更好的估计值，因为估计值的偏差较小，但差异通常很小，而且时间差很大。

我把这些算法的差异写成了一个Jupyter笔记本nAGQ.ipynb。该写包使用MixedModels包用于Julia，而不是lme4，但方法类似。(我是lme4的作者之一，也是MixedModels的作者。)

如果你要做很多GLMM适合，我会考虑在Julia与MixedModels。它通常比R快得多，即使在lme4中有这么多复杂的代码。

Answer 2

你确定Stata在读整份文件吗？

http://www.stata.com/manuals13/rlimits.pdf

我之所以问这个问题，是因为在我看来，你已经对50k人进行了25次观察(1,250,000行)；根据你使用的Stata版本，你可能会得到截断的结果。

编辑，因为它不是一个文件长度问题，您是否尝试过其他包的混合效果，如nlme？我怀疑非线性混合效应模型会更快地接受你的数据。

编辑这个资源可能比关于不同模型的任何东西都更有用：https://stats.stackexchange.com/questions/173813/r-mixed-models-lme-lmer-or-both-which-one-is-relevant-for-my-data-and-why