glmer() 以一个因子水平作为基线收敛，但在将基线重新调整到另一个水平时无法收敛答案

【问题标题】：glmer() converges with one factor level as baseline, but fails to converge when re-level baseline to the other levelglmer() 以一个因子水平作为基线收敛，但在将基线重新调整到另一个水平时无法收敛
【发布时间】：2020-08-10 01:21:25
【问题描述】：

背景：我正在拟合具有两个固定效应和两个交叉随机效应的逻辑混合模型。构成固定效应的两个变量都是二元分类变量（即没有数字预测变量）。

数据来自眼动仪，随着时间的推移，每个受试者都有数千个 0 和 1。我想考虑的一个模型是逻辑混合模型，它考虑了个人内部 0 和 1 的相关性，以及个人在给定时间查看的图像（在整个实验中，他们查看了六张不同的图像， 12 种可能）。（顺便说一句，我也打算用增长曲线分析来分析这些数据，但也很想尝试一下，暂时忽略了时间进程。）

问题：我使用两个二元预测变量的虚拟编码运行模型，每个预测变量都有默认的基线水平。模型收敛。根据我已经制作的数据图，这也是有道理的：估计在预期的方向上，固定效应的重要性与预期的一样，等等。

然后，我想在基线的 other 级别使用预测变量之一拟合重新参数化的模型，以获得相关的 p 值（参见，例如：https://stats.stackexchange.com/questions/449723/calculating-simple-effects-from-standard-regression-output)。我使用relevel() 来执行此操作，然后再次拟合完全相同的模型。然而，这一次出现了收敛警告。我又试了一次（因为为什么不），但又失败了。我想不出基线水平会影响收敛（或缺乏收敛）的原因。

数据结构：这是一个小的虚拟数据集，显示了数据的一般结构（大大比实际数据集减少了大约 5000 个观察值（0 和 1） ) 每个人每张图像，总共超过 100 万行（实际数据集中有大约 60 个主题）。注意 Time_stamp 不在此处的模型中，因为我没有考虑此分析中的时间过程。 AOI_look 中的“AOI”代表“感兴趣区域”，这是眼动追踪数据中的常用术语。

structure(list(Subject = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), Image = c("A", "A", 
"A", "B", "B", "B", "C", "C", "C", "D", "D", "D", "E", "E", "E", 
"F", "F", "F", "G", "G", "G", "E", "E", "E", "D", "D", "D", "A", 
"A", "A", "J", "J", "J", "L", "L", "L"), Time_stamp = c(0L, 1L, 
10L, 0L, 8L, 16L, 0L, 7L, 16L, 0L, 1L, 10L, 0L, 8L, 16L, 0L, 
7L, 16L, 0L, 8L, 16L, 0L, 7L, 16L, 0L, 1L, 9L, 0L, 8L, 16L, 0L, 
8L, 16L, 0L, 1L, 9L), Intervention = c("Pre", "Pre", "Pre", "Pre", 
"Pre", "Pre", "Pre", "Pre", "Pre", "Post", "Post", "Post", "Post", 
"Post", "Post", "Post", "Post", "Post", "Pre", "Pre", "Pre", 
"Pre", "Pre", "Pre", "Pre", "Pre", "Pre", "Post", "Post", "Post", 
"Post", "Post", "Post", "Post", "Post", "Post"), Sex = c("Female", 
"Female", "Female", "Female", "Female", "Female", "Female", "Female", 
"Female", "Female", "Female", "Female", "Female", "Female", "Female", 
"Female", "Female", "Female", "Male", "Male", "Male", "Male", 
"Male", "Male", "Male", "Male", "Male", "Male", "Male", "Male", 
"Male", "Male", "Male", "Male", "Male", "Male"), AOI_look = c(1L, 
1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 
1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 
1L, 1L, 1L)), class = "data.frame", row.names = c(NA, -36L))

模型 1：干预 = 发布为基线：没有收敛警告

logit_model01 <- glmer(AOI_look ~ factor(Sex)*factor(Intervention) + 
                                  (1 | Image) + 
                                  (1 | Subject), 
                       data = data01, 
                       family="binomial")

模型 1 的输出（收敛）：

summary(logit_model01)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod]
Family: binomial  ( logit )
Formula: AOI_look ~ factor(Sex) * factor(Intervention) + (1 | Image) +  
         (1 | Subject)
Data: data01

      AIC       BIC    logLik  deviance  df.resid 
1519023.6 1519095.5 -759505.8 1519011.6   1182007 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6200 -0.9733  0.5504  0.8444  2.2317 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 0.2320   0.4816  
 Image    (Intercept) 0.1318   0.3631  
Number of obs: 1182013, groups:  Subject, 59; Image, 9

Fixed effects:
                                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)                               0.052133   0.036172   1.441     0.15    
factor(Sex)Male                           0.187736   0.039967   4.697 2.64e-06 ***
factor(Intervention)Pre                   0.257262   0.006136  41.924  < 2e-16 ***
factor(Sex)Male:factor(Intervention)Pre  -0.220240   0.007932 -27.766  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
              (Intr) fc(S)M fc(I)P
fctr(Sex)M    -0.006               
factr(Int)Pr   0.011  0.012        
f(S)M:(I)P    -0.013 -0.018  -0.747

模型 2：与模型 1 相同，但干预 = Pre 在运行前设置为基线：获得收敛警告：

logit_model02 <- glmer(AOI_look ~ factor(Sex)*factor(Intervention) + 
                                  (1 | Image) + 
                                  (1 | Subject), 
                       data = data01, 
                       family="binomial")

收敛警告：

Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge with max|grad| = 0.00210935 (tol = 0.002, component 1)

虽然有警告，但我仍然可以使用 summary() 获得 Model 2 的输出：

summary(logit_model02)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod]
Family: binomial  ( logit )
Formula: AOI_look ~ factor(Sex) * factor(Intervention) + (1 | Image) +  
         (1 | Subject)
Data: data01

      AIC       BIC    logLik  deviance  df.resid 
1519023.6 1519095.5 -759505.8 1519011.6   1182007 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6200 -0.9733  0.5504  0.8444  2.2317 

Random effects:
 Groups    Name        Variance Std.Dev.
 Subject   (Intercept) 0.2320   0.4816  
 Image     (Intercept) 0.1318   0.3631  
Number of obs: 1182013, groups:  Subject, 59; Image, 9

Fixed effects:
                                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)                               0.309400   0.038106   8.120 4.68e-16 ***
factor(Sex)Male                          -0.032522   0.042520  -0.765    0.444    
factor(Intervention)Post                 -0.257262   0.006158 -41.774  < 2e-16 ***
factor(Sex)Male:factor(Intervention)Post  0.220234   0.007926  27.787  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
              (Intr) fc(S)M fc(I)P
fctr(Sex)M    -0.107               
fctr(Int)Pst   0.009  0.024        
f(S)M:(I)P    -0.013 -0.025  -0.748
convergence code: 0
Model failed to converge with max|grad| = 0.00210935 (tol = 0.002, component 1)

为了清晰/完整，在 Ben Bolker 的答案下方的 cmets 中添加了信息：

模型 3：与模型 1 相同，但基于聚合数据（0 和 1 按主题和图片组合）：

结果与模型 1 的结果不相同（拟合非聚合的 0/1 数据），即使模型 1 没有收敛警告。即，所有点估计都相似，但并非所有 SE -subjects) 和截距明显不同：

logit_model03 <- glmer(cbind(Successes, Failures) ~ factor(Sex)*factor(Intervention) +
                                                    (1 | Image) + 
                                                    (1 | Subject), 
                       data = data01_agg, 
                       family="binomial")

summary(logit_model03)

Fixed effects:
                                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)                              0.052151   0.150564   0.346    0.729    
factor(Sex)Male                          0.187734   0.125563   1.495    0.135    
factor(Intervention)Pre                  0.257261   0.006227  41.314   <2e-16 ***
factor(Sex)Male:factor(Intervention)Pre -0.220239   0.008061 -27.322   <2e-16 ***

为了查看原始数据集的绝对大小（超过 167 万行）是否会导致问题，我从原始数据集中随机抽取了 15000 个观察值样本，并使用非聚合 (0/ 1）和汇总（成功、失败）数据。（我对 Intervention 的两个基线级别都这样做了，但我只是展示了 Intervention = 在此处发布。）：

模型 4：15000 个 obs 的随机样本。来自实际数据集，非汇总 (0/1) 结果，干预 = 作为基线发布：

logit_model04 <- glmer(AOI_look ~ factor(Sex)*factor(Intervention) + 
                                  (1 | Image) + 
                                  (1 | Subject), 
                       data = smaller_raw_data, 
                       family="binomial")

summary(logit_model04)

Fixed effects:
                                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)                             -0.03212    0.15564  -0.206 0.836505    
factor(Sex)Male                          0.28384    0.14121   2.010 0.044429 *  
factor(Intervention)Pre                  0.22727    0.06449   3.524 0.000425 ***
factor(Sex)Male:factor(Intervention)Pre -0.25061    0.08410  -2.980 0.002884 **

模型 5：15000 个 obs 的随机样本。来自实际数据集，汇总（成功、失败）结果，干预 = 发布为基线：

logit_model05 <- glmer(cbind(Successes, Failures) ~ factor(Sex)*factor(Intervention) + 
                                                    (1 | Image) + 
                                                    (1 | Subject), 
                       data = smaller_agg_data, 
                       family="binomial")

summary(logit_model05)

Fixed effects:
                                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)                             -0.03211    0.15569  -0.206 0.836582    
factor(Sex)Male                          0.28383    0.14126   2.009 0.044500 *  
factor(Intervention)Pre                  0.22727    0.06450   3.524 0.000425 ***
factor(Sex)Male:factor(Intervention)Pre -0.25061    0.08412  -2.979 0.002890 **

如您所见，模型 4 和 5 的固定效应的点估计值和 SE 至少在小数点后 5 位是相同的。我还以 Intervention = Pre 作为基线重复了这一点，结果在非聚合模型和聚合模型之间的相似性方面是等效的（正如我所预期的那样）。

因此，使用非聚合或聚合数据在理论上应该无关紧要，至少对于大小合理（即不太大）的数据集而言。这让我想知道当数据集“太大”时模型是否存在不稳定性。来自?lme4::convergence，有理由相信这可能是这种情况：“探索性分析表明（1）对于大型数据集（观察数大于大约 1e5），Hessian 的幼稚估计可能会失败。”所以也许这就是我在上面的模型 1-3 中观察到的根源。

在我的情况下，我可以选择聚合数据，因为我的两个预测变量都是二进制分类的，但是对于不允许聚合的数字预测变量，这不是在拟合之前减少计算负载的选项型号。

【问题讨论】：

1/2 @Ben Bolker - 这与您刚刚回答的问题 (stackoverflow.com/questions/63360751/…) 直接相关。实际上，我写这个问题是为了“支持”这个问题，以防我对模型的编码方式有问题。因为，对于我的另一个问题，您认为我使用的代码是合理的，我想知道您是否对为什么会发生这种情况有任何倾向。
2/2 我将尝试添加(1|Subject:Trial) 作为随机效果，看看是否可以解决它。注意这个是和stackoverflow.com/questions/63360751/…完全一样的场景，我刚才在这里提到它是眼动追踪数据，所以结果被称为AOI_look而不是Binary_outcome。再次感谢您的帮助。（另外，这个问题是否也应该迁移到 CrossValidated 中，还是在这里可以？对我来说，这个问题更清楚地是关于编码。）
我想我可能会做一些事情：当我在聚合而不是原始数据上运行模型时，我得到的参数估计值与上面收敛的模型非常相似（直到小数点后第 5 位或第 6 位）（第一个，以Post 作为基线），但 Sex 的 p 值非常不同（0.135 与上面给出的 2.64e-06），如以及截距项。这让我想知道 lmer() 是否由于原始数据的大小（约 160 万行）而“超时”（或者在上面的后一个模型中存在收敛问题），从而给出不准确的推断。正在处理它...将报告回来。
目前正在查看rstudio-pubs-static.s3.amazonaws.com/…，但认为我需要尽可能增加内存并尝试让模型运行更长时间。

标签： r logistic-regression lme4 mixed-models convergence

【解决方案1】：

不要将收敛检查作为全部和最终目标

虽然确实理论上正确的是，改变随机项的顺序、重新调整因子、居中/缩放/等调整根本不应该改变模型拟合的结果，但确实如此稍微改变线性代数，因此可以使模型诊断超过一个报告不收敛的阈值。要记住的主要事情是模型收敛的警告需要小心谨慎：它们并不一定意味着模型有任何问题，只是你应该加倍-查看。事实上，如果您仔细查看了汇总输出（例如，比较对数似然和随机效应方差，这应该几乎相同）和预测（同上），那么您可以得出结论，收敛警告并不重要.这是我检查一些简单预测的代码：

cc1 <- c(0.052133, 0.187736, 0.257262, -0.220240)
cc2 <- c(0.309400,-0.032522,-0.257262,0.220234)

df2 <- expand.grid(Sex=c("Female","Male"),Intervention=c("Pre","Post"))
df1 <- transform(df2,Intervention=relevel(Intervention,"Post"))


X1 <- model.matrix(~Sex*Intervention, data=df1)
X2 <- model.matrix(~Sex*Intervention, data=df2)
cbind(drop(X1 %*% cc1), drop(X2 %*% cc2))
##       [,1]     [,2]
## 1 0.309395 0.309400
## 2 0.276891 0.276878
## 3 0.052133 0.052138
## 4 0.239869 0.239850

all.equal(drop(X1 %*% cc1), drop(X2 %*% cc2))
## [1] "Mean relative difference: 4.78203e-05"

这些在4e-05 的水平上有所不同，因此除非您期望并需要预测到该准确度水平，否则我不会担心。（鉴于您在下面的 cmets 中所说的关于标准误差计算的差异，我也可以尝试merTools::predictInterval() 来寻找不确定性的差异。）另请参阅：?lme4::convergence、?lme4::troubleshooting、?lme4::allFit

不要费心重新找平和改装

由于所有水平调整都可以通过简单的线性代数来完成，因此有多种方便的下游工具可让您根据需要进行比较而无需重新调整和重新调整。这更快，更方便！我将从emmeans 包开始； effects 包也可能有用。

【讨论】：

评论不用于扩展讨论；这个对话是moved to chat。

【解决方案2】：

将一个因素重新参数化以具有不同的参考水平应该不会影响模型的收敛性，也不会影响其结果的解释。

首先处理结果，Intervention 是一个两级因素。因此，效果或reordering 它使Pre 成为参考水平（显然是明智之举）是改变相应参数估计的符号。（因此，在您的结果中，0.257262 更改为 -0.257262。）没有其他内容应该更改：p 值、标准错误等都保持不变。不仅用于干预，还用于模型中的所有其他术语。

不仅事情会发生变化，而且模型无法收敛，这表明您的模型存在问题。

logit_model <- glmer(AOI_look ~ factor(Sex)*factor(Intervention) + 
                                (1 | Image) + 
                                (1 | Subject)

我认为这里可能存在几个问题。 Intervention 不应该是主题内效果吗？（您不能随机化受试者内的干预顺序。）Image 不应该是受试者之间的效果吗？（显示给患者 1 的图像“A”与显示给患者 2 的图像“A”相同。）可能还包括患者体内Image 效应。

拟合纵向混合效应模型是一项复杂的任务。我认为您可能需要就您的具体案例采取详细的专家建议，以了解哪种模型适合您的研究以及如何适应它。模型收敛的事实并不意味着它是正确的模型。

【讨论】：

1/9 @Limey 感谢您的 cmets。请注意，我并不是因为它收敛而假设它是正确的模型，而是因为我考虑了很长时间应该包含的每个变量，考虑了两个相关级别（主题和试验），并且模型给出了合理的估计和p 值基于数据在图中向我显示的内容（我在下面讨论了许多 cmets，因为响应很长）。
2/9 您对干预的看法是正确的——它在主体内。除了为参与者和试验指定随机效应外，我不/不知道如何进一步解释这一点。你有什么建议吗？
3/9 是的，trail A 总是 Trial A，例如，但 Trial 本身并不是一个感兴趣的变量，这就是为什么我没有将它作为固定效应拟合，只是随机的。我只是想说明这样一个事实，即当给定的参与者查看某个图像时，他或她可能对该图像比另一个图像或多或少感兴趣，因此由于图像可能存在额外的聚类相关性，超出相关性已经存在于个人体内。关于如何更好/适当地建模/说明这一点的其他想法？
4/9 这是我对您所写内容的其余回复，该回复继续跨越多个 cmets：由于交互作用，更改基线的影响超过该参数的估计符号.正如你所说，对于干预，量级将保持不变，只有符号会改变。但是，所有其他估计都发生了变化（为此，我更新了我的帖子以显示以 Intervention = Pre 作为基线的重新参数化模型的输出）。这很重要，因为我想比较两个二元预测变量的不同组合的估计 ln(odds)。
5/9 这就是为什么我什至想要更改基线 - 尽管我可以从第一个代数模型中获得我可能想要的所有估计 ln(odds)，我想要干预前后的 Sex 的 p 值。例如，使用第一个参数化（“Post”作为基线），我可以得出以下结论：对于 Intervention = Post（基线），男性和女性之间的 ln(odds) 存在显着差异，根据 p - 与性相关的值 (2.64e-06