lme4::glmer 对比 Stata 的 melogit 命令
lme4::glmer vs. Stata's melogit command
最近我一直在尝试将大量随机效应模型拟合到相对较大的数据集。假设在最多 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 小时，并且 returns 出现无法收敛的错误。我做了很多故障排除（例如，遵循 these 指南），但没有任何改进。最后收敛甚至不接近（最大梯度在 5-10 左右，而我认为收敛标准是 0.001）。
然后我尝试使用 melogit 命令在 Stata 中拟合模型。模型拟合不到 2 分钟，没有收敛问题。对应的Stata命令为
melogit outcome treatment i.study_quarter i.dd_quarter || id:
什么给了？ Stata 只是有一种更好的拟合算法，还是一种针对大型模型和大型数据集进行了更好优化的算法？ 运行 时代的差异真是令人惊讶。
您确定 Stata 正在读取整个文件吗？
http://www.stata.com/manuals13/rlimits.pdf
我问的原因是因为在我看来你已经对 5 万人进行了 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
使用可选参数 nAGQ=0L，glmer 拟合可能会快得多。您有许多固定效应参数（study_quarter 和 dd_quarter 各有 20 个水平，总共产生 28 个对比）并且默认优化方法（对应于 nAGQ=1L）将所有这些系数进入一般的非线性优化调用。使用 nAGQ=0L，这些系数在更快的惩罚迭代重新加权最小二乘 (PIRLS) 算法中得到优化。默认值通常会提供更好的估计，因为估计的偏差较低，但差异通常很小，时间差异很大。
我将这些算法的差异写成 Jupyter notebook nAGQ.ipynb. That writeup uses the MixedModels package for Julia instead of lme4，但方法相似。 （我是lme4的作者之一，也是MixedModels的作者。）
如果您要进行大量 GLMM 拟合，我会考虑在 Julia 和 MixedModels 中进行。它通常比 R 快得多，即使 lme4.
中的所有复杂代码也是如此
最近我一直在尝试将大量随机效应模型拟合到相对较大的数据集。假设在最多 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 小时，并且 returns 出现无法收敛的错误。我做了很多故障排除（例如，遵循 these 指南），但没有任何改进。最后收敛甚至不接近（最大梯度在 5-10 左右，而我认为收敛标准是 0.001）。
然后我尝试使用 melogit 命令在 Stata 中拟合模型。模型拟合不到 2 分钟，没有收敛问题。对应的Stata命令为
melogit outcome treatment i.study_quarter i.dd_quarter || id:
什么给了？ Stata 只是有一种更好的拟合算法，还是一种针对大型模型和大型数据集进行了更好优化的算法？ 运行 时代的差异真是令人惊讶。
您确定 Stata 正在读取整个文件吗？
http://www.stata.com/manuals13/rlimits.pdf
我问的原因是因为在我看来你已经对 5 万人进行了 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
使用可选参数 nAGQ=0L，glmer 拟合可能会快得多。您有许多固定效应参数（study_quarter 和 dd_quarter 各有 20 个水平，总共产生 28 个对比）并且默认优化方法（对应于 nAGQ=1L）将所有这些系数进入一般的非线性优化调用。使用 nAGQ=0L，这些系数在更快的惩罚迭代重新加权最小二乘 (PIRLS) 算法中得到优化。默认值通常会提供更好的估计，因为估计的偏差较低，但差异通常很小，时间差异很大。
我将这些算法的差异写成 Jupyter notebook nAGQ.ipynb. That writeup uses the MixedModels package for Julia instead of lme4，但方法相似。 （我是lme4的作者之一，也是MixedModels的作者。）
如果您要进行大量 GLMM 拟合，我会考虑在 Julia 和 MixedModels 中进行。它通常比 R 快得多，即使 lme4.