【问题标题】：mixed-effects logistic regression: different results with MASS and lme4混合效应逻辑回归：MASS 和 lme4 的不同结果
【发布时间】：2017-04-10 08:52:14
【问题描述】：

我使用 MASS 和 lme4 运行混合效应逻辑回归，但我得到不同的结果，我想知道是否（以及在哪里）出了问题。

我的模型如下：

DV: synt_num (0 vs 1)

IV： voc1_tvl：连续

wo：二分法（“vs”和“vo”）

动词：二分法（“c”和“t”）

大量使用

如果我输入以下公式：

synt1 = glmmPQL (synt_num ~ voc1_tvl*verb*wo, random = ~1|participants, data=opz, family="binomial")

我得到以下输出

Linear mixed-effects model fit by maximum likelihood
 Data: ops 
  AIC BIC logLik
   NA  NA     NA

Random effects:
 Formula: ~1 | participants
        (Intercept)    Residual
StdDev:    21.50518 0.003233762

Variance function:
 Structure: fixed weights
 Formula: ~invwt 
Fixed effects: synt_num ~ voc1_tvl * verb * wo 
                                Value Std.Error DF    t-value p-value
(Intercept)                 -294.2627  22.41770 65 -13.126356  0.0000
voc1_tvl                       5.0748   0.27293 65  18.593945  0.0000
verb[T.t]                    152.6929  36.62349 47   4.169260  0.0001
wo[T.vs]                     327.9534  36.62150 47   8.955216  0.0000
voc1_tvl:verb[T.t]            -3.3897   0.45180 47  -7.502784  0.0000
voc1_tvl:wo[T.vs]             -5.3135   0.45178 47 -11.761381  0.0000
verb[T.t]:wo[T.vs]          -252.0190  73.23996 47  -3.441004  0.0012
voc1_tvl:verb[T.t]:wo[T.vs]    4.3519   0.90346 47   4.816905  0.0000
 Correlation: 
                            (Intr) vc1_tv vr[T.] w[T.v] vc1_tvl:v[T.]
voc1_tvl                    -0.992                                   
verb[T.t]                   -0.588  0.601                            
wo[T.vs]                    -0.588  0.601  1.000                     
voc1_tvl:verb[T.t]           0.588 -0.601 -1.000 -1.000              
voc1_tvl:wo[T.vs]            0.588 -0.601 -1.000 -1.000  1.000       
verb[T.t]:wo[T.vs]           0.588 -0.601 -1.000 -1.000  1.000       
voc1_tvl:verb[T.t]:wo[T.vs] -0.588  0.601  1.000  1.000 -1.000       
                            vc1_tvl:w[T.] v[T.]:
voc1_tvl                                        
verb[T.t]                                       
wo[T.vs]                                        
voc1_tvl:verb[T.t]                              
voc1_tvl:wo[T.vs]                               
verb[T.t]:wo[T.vs]           1.000              
voc1_tvl:verb[T.t]:wo[T.vs] -1.000        -1.000

Standardized Within-Group Residuals:
          Min            Q1           Med            Q3           Max 
-3.8945233603 -0.0025293449 -0.0001866073  0.0018900747  3.2648586389 

Number of Observations: 120
Number of Groups: 67

现在，有人告诉我，如果 SE 超过估计值的一半，则该模型不可靠。然而，情况似乎并非如此。也就是说，我认为输出没有任何问题。如果我将数据绘制成图表，这些结果是有道理的（而且它也符合预测）。

但是，当我使用 lme4 时.....

使用 LME4

如果我使用 lme4 和以下代码

>synt = glmer (synt_num ~ voc1_tvl*verb*wo + (1|participants), data=opz, family="binomial")

我收到以下警告

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge with max|grad| = 0.00757738 (tol = 0.001, component 1)
Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?

现在，为什么我会收到关于 lme4 而不是 MASS 的警告？？

据我了解，最后一个警告可以通过将我拥有的唯一连续变量 (voc1_tvl) 居中来解决。我错了吗？

其他警告呢？

如果我强制输出最后一个模型，我会得到与之前的 MASS 输出完全不同的东西

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: synt_num ~ voc1_tvl * verb * wo + (1 | participants)
   Data: opz

     AIC      BIC   logLik deviance df.resid 
   117.8    142.9    -49.9     99.8      111 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8743 -0.3513  0.2158  0.3564  2.5211 

Random effects:
 Groups       Name        Variance Std.Dev.
 participants (Intercept) 1.074    1.036   
Number of obs: 120, groups:  participants, 67

Fixed effects:
                             Estimate Std. Error z value Pr(>|z|)  
(Intercept)                 -9.857319   4.856734  -2.030   0.0424 *
voc1_tvl                     0.153749   0.070039   2.195   0.0281 *
verb[T.t]                   -1.321681   6.377127  -0.207   0.8358  
wo[T.vs]                    11.089804   6.181132   1.794   0.0728 .
voc1_tvl:verb[T.t]          -0.004359   0.082863  -0.053   0.9580  
voc1_tvl:wo[T.vs]           -0.135696   0.081618  -1.663   0.0964 .
verb[T.t]:wo[T.vs]          -6.648285   8.266898  -0.804   0.4213  
voc1_tvl:verb[T.t]:wo[T.vs]  0.068410   0.105191   0.650   0.5155  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
              (Intr) vc1_tv vr[T.] w[T.v] vc1_tvl:v[T.] vc1_tvl:w[T.] v[T.]:
voc1_tvl      -0.986                                                        
verb[T.t]     -0.528  0.493                                                 
wo[T.vs]      -0.834  0.829  0.431                                          
vc1_tvl:v[T.]  0.594 -0.579 -0.985 -0.490                                   
vc1_tvl:w[T.]  0.838 -0.848 -0.445 -0.986  0.518                            
vr[T.]:[T.]    0.583 -0.566 -0.794 -0.719  0.802         0.710              
v1_:[T.]:[T   -0.596  0.591  0.802  0.724 -0.828        -0.737        -0.986
convergence code: 0
Model failed to converge with max|grad| = 0.00757738 (tol = 0.001, component 1)
Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?

现在，可能是我的模型数据太少，事实上，如果我在没有交互的情况下进入模型，我不会收到警告和以下输出

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: synt_num ~ voc1_tvl + verb + wo + (1 | participants)
   Data: opz

     AIC      BIC   logLik deviance df.resid 
   118.7    132.7    -54.4    108.7      115 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.2222 -0.4314  0.2455  0.4756  2.7784 

Random effects:
 Groups       Name        Variance Std.Dev.
 participants (Intercept) 0.4567   0.6758  
Number of obs: 120, groups:  participants, 67

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.96657    1.63767  -2.422 0.015432 *  
voc1_tvl     0.08047    0.02380   3.381 0.000721 ***
verb[T.t]   -2.13885    0.65960  -3.243 0.001184 ** 
wo[T.vs]    -0.18627    0.50623  -0.368 0.712901    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
          (Intr) vc1_tv vr[T.]
voc1_tvl  -0.933              
verb[T.t]  0.271 -0.531       
wo[T.vs]  -0.168 -0.011  0.051

因此，正如 Ben 在下面建议的那样，我可能应该追求一组更有限/更少的交互。如果只使用两个变量进入模型，但有一个交互如下：

> synt = glmer (synt_num ~ voc1_tvl*verb + (1|participants), data=opz, family="binomial")

我再次收到以下警告：

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?

我得到的输出表明动词不重要

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: synt_num ~ voc1_tvl * verb + (1 | participants)
   Data: opz

     AIC      BIC   logLik deviance df.resid 
   118.5    132.5    -54.3    108.5      115 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.6510 -0.4008  0.2483  0.4341  3.0910 

Random effects:
 Groups       Name        Variance Std.Dev.
 participants (Intercept) 0.5154   0.7179  
Number of obs: 120, groups:  participants, 67

Fixed effects:
                   Estimate Std. Error z value Pr(>|z|)  
(Intercept)        -3.24172    2.16991  -1.494   0.1352  
voc1_tvl            0.06912    0.03060   2.259   0.0239 *
verb[T.t]          -3.80171    3.10388  -1.225   0.2206  
voc1_tvl:verb[T.t]  0.02172    0.03888   0.559   0.5764  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) vc1_tv vr[T.]
voc1_tvl    -0.968              
verb[T.t]   -0.594  0.518       
vc1_tv:[T.]  0.653 -0.619 -0.977
convergence code: 0
Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables

为什么当我只输入主效应时它会告诉我动词是显着的，而当我输入交互作用并且只输入两个预测变量时却不是？

这很可疑，因为如果我使用回归树，动词会被拾取并且很重要。

好像我用lme4的时候，模型不能进行交互

大问题

我可以信任第一个 (MASS) 输出吗？下面的友好答案表明不，我不能

我应该担心 lme4 说什么吗？

在第二种情况下，我的数据集是否超出了希望，因为我正在调查的是交互的作用？

本说：

您可能会因为完全分离而遇到问题（即问题在于您的数据，而不是拟合过程），但如果您处理这种可能性，然后尝试使用其他优化器拟合 glmer 模型 - 并得到类似的结果 - 然后您可以忽略警告。

我尝试使用以下公式使用“nloptwrap” (http://angrystatistician.blogspot.co.uk/2015/10/mixed-models-in-r-bigger-faster-stronger.html)

> synt.model = synt_num ~ voc1_tvl*verb*wo + (1|participants)
> fit <- glmer(synt.model, data=opz, verbose=TRUE, family=poisson(link=log), nAGQ=0, control=glmerControl(optimizer = "nloptwrap"))

但是，我收到以下警告

Warning in optwrap(optimizer, devfun, start, rho$lower, control = control,  :
  'verbose' not yet passed to optimizer 'nloptwrap'; consider fixing optwrap()

因此我摆脱了冗长的部分并输入了

> fit <- glmer(synt.model, data=opz, family=poisson(link=log), nAGQ=0, control=glmerControl(optimizer = "nloptwrap"))

和

> summary(fit)

 ( log )
Formula: synt_num ~ voc1_tvl * verb * wo + (1 | participants)
   Data: opz
Control: glmerControl(optimizer = "nloptwrap")

     AIC      BIC   logLik deviance df.resid 
   234.0    259.1   -108.0    216.0      111 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-0.96391 -0.45210  0.08698  0.26229  2.25782 

Random effects:
 Groups       Name        Variance Std.Dev.
 participants (Intercept) 0        0       
Number of obs: 120, groups:  participants, 67

Fixed effects:
                            Estimate Std. Error z value Pr(>|z|)
(Intercept)                 -2.29833    1.57388  -1.460    0.144
voc1_tvl                     0.02457    0.01779   1.381    0.167
verb[T.t]                   -1.33202    2.46242  -0.541    0.589
wo[T.vs]                     2.05098    1.89449   1.083    0.279
voc1_tvl:verb[T.t]           0.01354    0.02807   0.482    0.630
voc1_tvl:wo[T.vs]           -0.02277    0.02207  -1.032    0.302
verb[T.t]:wo[T.vs]          -2.22007    3.31623  -0.670    0.503
voc1_tvl:verb[T.t]:wo[T.vs]  0.02114    0.03799   0.556    0.578

Correlation of Fixed Effects:
              (Intr) vc1_tv vr[T.] w[T.v] vc1_tvl:v[T.] vc1_tvl:w[T.] v[T.]:
voc1_tvl      -0.992                                                        
verb[T.t]     -0.639  0.634                                                 
wo[T.vs]      -0.831  0.824  0.531                                          
vc1_tvl:v[T.]  0.629 -0.634 -0.992 -0.522                                   
vc1_tvl:w[T.]  0.799 -0.806 -0.511 -0.988  0.511                            
vr[T.]:[T.]    0.475 -0.471 -0.743 -0.571  0.737         0.564              
v1_:[T.]:[T   -0.464  0.468  0.733  0.574 -0.739        -0.581        -0.991

输出与我之前的输出不同，这里没有什么重要的（原文如此）。因此，我应该担心第一个（lme4）输出。我应该吗？

我会尝试不同的优化器，看看结果是否不同。

但是，如果结果在达到显着性的方式与原始模型和使用 nloptwrap 优化的模型确实不同时不同，这是否意味着即使 lme4 也不可靠？

非常感谢任何勇敢（和善良）到最后阅读这篇文章的人

【问题讨论】：

您是否尝试过按照警告的建议重新调整变量？您的预测变量似乎在非常不同的范围内。
嗨，约翰，您的意思是使连续变量居中吗？或者重新缩放是否意味着其他/不同？ - 谢谢
缩放...就像将以毫米为单位的距离转换为米。这就是它要求你做的事情。它可能只是将一个变量乘以 100 或除以它。

标签： r lme4 mixed-models

【解决方案1】：

为什么整个分析可能有点粗略

统计建模的一个经验法则（来自 Harrell 的 Regression Modeling Strategies 一书）是，对于每个要估计的参数，您需要 10-20 个观测值的有效样本量。对于二进制数据，“有效样本量”为最小值（0 个响应的数量，1 个响应的数量），因此在最好的情况下，它将是您观察次数的一半（如果大多数观察为 0 或大多数观察为 1，则更少） , 大约你的情况是60。你有 9 个参数（8 个固定效应加上随机效应方差），所以这是推动它。也许您最感兴趣的互动减少了？

为什么我不信任 MASS::glmmPQL 结果

您的数据每个集群的有效样本量非常小：也就是说，这些是二元数据，平均每个受试者少于两个观察值（总共 120 个 obs，67 个受试者）。这正是在可能性中逼近讨厌的积分的方式最重要的情况：通常，PQL 最差，其次是拉普拉斯逼近（即glmer 默认，nAGQ=1），其次是 Gauss-Hermite具有更多正交点的正交（glmer 和 nAGQ>1；我会尝试 nAGQ=10）。
一般而言，对于 GLM(M)，其响应为无单位（log 或 logit 标度），大效应量（例如 abs(beta)>10）通常是可疑的，表明数据集中完全分离 (例如见here) 或其他一些问题。
固定效应之间的相关性显示出许多完美的负 (-1.0) 和正 (+1.0) 值...表明模型存在问题

为什么 lme4::glmer 结果没有那么糟糕

先验我们预计glmer（默认为拉普拉斯近似）比 PQL 做得好一点，如上所述
参数估计和相关性没有那么极端
固定效应之间的参数估计和相关性没有那么差
?lme4::convergence 有一些关于警告和故障排除的信息。根据我的经验，0.007 的梯度还不错，可能是误报。由于完全分离，您可能会遇到问题（即问题在于您的数据，而不是拟合过程），但如果您处理这种可能性，然后尝试将 glmer 模型与其他优化器拟合 - 并获得类似的结果 - 然后您可以忽略警告。

更新

我不会使用nAGQ=0，这是错误的方向（您可能应该使用nAGQ>1 而不是默认的nAGQ=1）
如果你得到不同的结果使用不同的优化器，使用nAGQ>1，那么我会担心。你可以试试

source(system.file("utils", "allFit.R", package="lme4"))
allFit(original_fit)

在原始拟合上尝试所有不同的可用优化器。一般来说，具有最低 AIC/负对数似然的拟合是你最能相信的……） - 那时我可能会说是的，你的模型对于数据来说太复杂了，所以你需要做一些更简单的事情 - 但无论如何，我建议查看您的数据，看看是否有任何令人惊讶的现象（例如，个别数据点、参与者或因素组合级别的极端值/异常值） p>

【讨论】：

嗨，本，谢谢您的回答。我按下了有用的答案箭头。我仍然有一些疑问。我通过编辑我的问题来表达我的想法，因为评论太长了。