lme4::glmer与Stata的melogit命令比较

Question

最近我一直在尝试对相对较大的数据集拟合很多随机效应模型，比如说大约有50,000人（或更多）在最多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左右，我想）。

然后，我尝试使用 Stata 中的 melogit 命令拟合模型。该模型在不到2分钟的时间内拟合完成，没有出现收敛问题。对应的Stata命令是：

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

是什么使得Stata表现更好？是拥有更好的拟合算法，还是针对大型模型和大型数据集进行了更好的优化呢？两者的运行时间差异如此之大确实让人惊讶。

- Jonathan Gellar

---

这里是相反情况的一个例子 - R快速处理了看似相同但在Stata中无法估计的SEM模型。 - radek

我不确定，但是在glmer中二项分布的默认选项是probit而不是logit吗？也许你可以添加family = binomial(link = "logit")然后再试试看？ - Erdne Htábrob

@radek - 感谢您的分享，但我特别指的是混合效应建模，而不是结构方程模型。我知道有很多情况下 R "胜过" Stata。 - Jonathan Gellar

@EndreBorbáth - 不，默認值是logit。 - Jonathan Gellar

我和你描述的经历完全一样 - 在R中使用二项式glmer时，需要运行数天并出现收敛警告 - 而在Stata中使用melogit只需几分钟即可完成拟合，并获得与相同规范下几乎相同的结果。 - Tom Wagstaff

Answer 1

< p >使用可选参数nAGQ=0L，glmer拟合速度可能会更快。您有许多固定效应参数（study_quarter和dd_quarter各有20个级别，共生成28个对比度），默认的优化方法（对应nAGQ=1L）将所有这些系数放入一般非线性优化调用中。使用nAGQ=0L，这些系数在更快的惩罚迭代加权最小二乘（PIRLS）算法内进行优化。默认情况通常会以偏差较低的意义提供更好的估计值，但差异通常非常小，时间差异巨大。

我有一篇关于这些算法差异的写作，它是一个 Jupyter 笔记本 nAGQ.ipynb。那篇文章使用了 Julia 的 MixedModels 包而不是 lme4，但方法类似。（我是 lme4 的作者和 MixedModels 的作者之一。）

如果你要进行大量的 GLMM 拟合，我建议你使用 Julia 和 MixedModels 进行拟合。即使在 lme4 中有所有复杂的代码，它通常比 R 快得多。

- Douglas Bates

---

我尝试将nAGQ设置为0，模型在大约5分钟内适配完成，没有任何收敛问题。虽然比Stata慢一点，但已经足够可接受了。 - Jonathan Gellar

我想知道，如果lmer首先使用PIRLS（对应于nAGQ=0）拟合模型以获得更好的起始值，然后再使用拉普拉斯近似（nAGQ=1）来优化模型，那么默认nAGQ=1L的模型是否会更快地收敛？也许这已经在内部完成了，但是如果是这样的话，我很惊讶在运行15个小时后，该模型仍然没有接近收敛（最大梯度约为5-10）。 - Jonathan Gellar

一个初始的 nAGQ=0L 拟合用于创建 nAGQ=1L 的起始估计值。它有所帮助，但不像预期的那样显著。问题在于非线性优化算法中要优化的参数数量。您还可以尝试将优化器从 nloptr 包中的默认设置更改为 NLOPT_LN_BOBYQA。 - Douglas Bates

一个后续问题：当你说估计值对于nAGQ=1L比PIRLS更好时，你是指固定效应估计值、随机效应估计值还是两者都是？如果我感兴趣的参数都是固定效应参数，那么使用PIRLS近似方法是否比如果我对随机效应感兴趣更“安全”？ - Jonathan Gellar

“更好”仅仅指的是Laplace逼近偏差的较低值。也就是说，估计结果更接近于最大似然估计值。 - Douglas Bates

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

- nsring

---

源代码已经过时了两个版本。请参考Stata 15文档：http://www.stata.com/manuals/rlimits.pdf。无论如何，在任何严肃的Stata中，100万个观测值（行，在您的术语中）本身并不是一个问题。 - Nick Cox

我还没有尝试过nlme，我同意它值得一试。但我怀疑“有点更快”是否能够让15小时（无收敛）和2分钟（有收敛）之间产生巨大的差异。我想拟合算法可能会有很多共同点，因为作者之间有重叠，但这只是一种假设。 - Jonathan Gellar

并且确认@NickCox所说的，模型输出表明它正在使用数据集中的所有观测值，因此我认为这不是一个问题。 - Jonathan Gellar

@nsring - 关于您第二次编辑中的链接，此主题关注使用REML拟合的lmer拟合和使用ML拟合的lme拟合之间的差异。一旦两个模型都使用REML进行拟合，这两种方法给出了几乎相同的结果。由于我的结果是二项式的，我正在使用glmer，它使用ML（请参见http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#reml-for-glmms）。Stata也使用ML，因此这不能是它们拟合之间的区别。 - Jonathan Gellar