R 中的逻辑混合模型收敛警告，但 Stata 中没有 [已关闭]

Question

我正在使用带有

lme4

的单变量逻辑混合模型来研究血浆参数 (M) 对协变量疾病结果的影响：

formula <- Disease ~ Age + Sex + BMI + M + (1 | Family)

样本量根据疾病状况信息的可用性而有所不同，但在 1000-3000 范围内（组数约为总样本量的 2/3）。

我使用的型号是：

glmer(formula, data = df, family = binomial, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))

对于

Disease

和

M

的某些组合，我收到收敛警告。对于某些疾病，我只会收到一些血浆参数的警告，而对于其他疾病，我会收到几乎一半的警告。

在遵循关于 lme4 收敛警告：故障排除指南的建议后，这些仍然存在（缩放变量、增加容差、尝试不同的求解器、从之前的拟合重新启动）。

很少有男性患有某些疾病的实例（

table(df$Sex, df$Disease)

），但是，我没有收到完美/准完美的分离警告。

问题：

1. 如果我的同事使用相同的数据集和模型公式在统计中运行分析，他们不会收到任何警告。怎么会这样？

2. 我不想出于实际目的放弃使用 R。我还能做些什么来解决 R 中的这些问题吗？

Answer 1

此信息可能是从某处重复的，但是：

• lme4

  中的收敛容差过于敏感。维护者很清楚这一点，并考虑过完全放弃这些特定的测试，但是如果没有大量的工作，很难确定您不会错过收敛警告是合理的重要案例。收敛警告暗示您的分析中的某些内容可能在数值上不稳定，而不是问题的明确迹象。
• “我的同事如何在 Stata 中运行这个程序却没有收到任何收敛警告？” Stata 正在拟合一个以相同方式定义的模型，但实现完全不同（请参阅参考手册了解更多详细信息）

• 对于小簇大小（基于“组数是总样本量的 2/3”，这意味着每组的平均观察数仅为 1.5 左右），您可能应该使用高斯-埃尔米特求积而不是拉普拉斯无论如何近似（即将

  nAGQ

  设置为大于 1 的值，例如 10 或 20）
• 我还能做些什么来解决 R 中的这些问题吗？
  • 还有其他用于混合逻辑回归的包（请参阅混合模型任务视图）；特别是，我会尝试

    GLMMadaptive

    和

    glmmTMB

    （虽然后者不能做GHQ，请参阅上一点......）
  • 对于顽固的收敛警告，我的一般建议是运行

    allFit()

    （参见下面的示例）；重点不是要查看是否有一个特定的优化器在运行时没有收敛警告，而是要查看所有优化器的结果（或至少那些实现类似对数似然的优化器）对于您的应用程序来说是否足够接近。

这是之前的SO问题的示例，已更新。

df <- read.csv("http://www.math.mcmaster.ca/bolker/misc/SO23478792_dat.csv")

library(lme4)
library(broom.mixed)
library(dplyr)
library(ggplot2); theme_set(theme_bw())

df <- transform(df, SUR.ID = factor(SUR.ID),
                replicate = factor(replicate),
                Unit = factor(seq(nrow(df))))
m1 <- glmer(cbind(ValidDetections, FalseDetections) ~ tm:Area + tm:c.distance +
            c.distance:Area + c.tm.depth:Area +
            c.receiver.depth:Area + c.temp:Area +
            c.wind:Area +
            c.tm.depth + c.receiver.depth +
            c.temp +c.wind + tm + c.distance + Area +
            replicate +
            (1|SUR.ID) + (1|Day) + (1|Unit) ,
            data = df, family = binomial)


## rescale
numcols <- grep("^c\\.",names(df))
dfs <- df
dfs[,numcols] <- scale(dfs[,numcols])
m4 <- update(m1, data=dfs)

aa <- allFit(m4)

使用

broom.mixed

中的功能：

glance(aa) |> select(optimizer, NLL_rel) |> arrange(NLL_rel)

这表明所有优化器都给出了非常相似的整体拟合，除了两个失败的优化器（对数似然尺度上 0.016 的差异很小）：

  optimizer                        NLL_rel
  <chr>                              <dbl>
1 nlminbwrap                      0       
2 bobyqa                          1.32e-10
3 optimx.L-BFGS-B                 4.94e- 5
4 nloptwrap.NLOPT_LN_BOBYQA       3.10e- 3
5 Nelder_Mead                     1.63e- 2
6 nmkbw                         Inf       
7 nloptwrap.NLOPT_LN_NELDERMEAD Inf       

查看固定效应估计，± 2 个标准误差（假设您对固定效应感兴趣；您还可以研究随机效应 [方差/协方差/标准差/相关性] 参数估计）

tt <- tidy(aa, effects = "fixed") |> filter(term != "(Intercept)") |> select(optimizer, term, estimate, std.error)
ggplot(tt, aes(x = estimate, y = term, colour = optimizer)) +
    geom_pointrange(aes(xmin=estimate-2*std.error, xmax = estimate+2*std.error),
                    position = position_dodge(width = 0.2)) +
    facet_wrap(~term, scales = "free") +
    theme(axis.text.y = element_blank())

相对于不确定性的大小，优化器之间的差异很小。或者你可以看看绝对差异：

ggplot(tt, aes(x = estimate, y = term, colour = optimizer)) +
    geom_point(position = position_dodge(width = 0.2)) +
    facet_wrap(~term, scales = "free") +
    theme(axis.text.y = element_blank())

您必须自己决定这些差异的大小对于您的问题是否重要。 （在这种特殊情况下，我怀疑还发生了一些完全分离问题，可能是因为示例中使用了一小部分数据。）