r - 混合效应逻辑回归 : different results with MASS and lme4

Question

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

我的模型如下:

DV:synt_num(0 与 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)输出。我应该吗？

我正在阅读贝茨的观察

(https://cran.r-project.org/web/packages/lme4/vignettes/Theory.pdf)。

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

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

非常感谢所有勇敢(和善良)阅读这篇文章直到最后的人

Answer 1

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

统计建模的一条经验法则(来自 Harrell 的回归建模策略一书)是，您需要为每个要估计的参数提供 10-20 个观测值的有效样本量。对于二进制数据，“有效样本量”为 min(0 响应的数量，1 响应的数量)，因此在最好的情况下，它将是观测值数量的一半(如果大多数观测值是 0 或大多数观测值是 1，则该值会更少) ，大约。你的情况是 60。你有 9 个参数(8 个固定效应加上随机效应方差)，所以这正在插入它。也许您最感兴趣的是一组减少的交互？

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

• 您的数据每个簇的有效样本量非常小:也就是说，这些是二进制数据，每个受试者平均少于两个观察值(总共 120 个观测值，67 个受试者)。在这种情况下，近似似然中令人讨厌的积分的方式最为重要:一般来说，PQL 最差，其次是拉普拉斯近似(即 glmer 默认值，其中 nAGQ= 1)，然后是具有更多正交点的高斯-埃尔米特求积(glmer 和 nAGQ>1；我会尝试 nAGQ=10 >).
• 一般来说，对于 GLM(M)，响应是无单位的(对数或对数标度)，大的效应量(例如 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>