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

转载作者：行者123 更新时间：2023-12-02 20:54:59

我使用 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)输出。我应该吗？

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

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

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

最佳答案

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

统计建模的一条经验法则(来自 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>

关于r - 混合效应逻辑回归 : different results with MASS and lme4，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/40813711/

文章推荐： android - 如何在android BLE中写入 “0xFF”作为特征值？

r - 错误: package or namespace load failed for ‘raster’ in inDL(x, as.逻辑(本地)，as.逻辑(现在)，...):
当我尝试加载库 Raster 时，我收到如下错误: 错误:inDL(x, as.logic(local), as.logic(now), ...) 中的“raster”的包或命名空间加载失败:无法加载
r - 错误: package or namespace load failed for ‘raster’ in inDL(x, as.逻辑(本地)，as.逻辑(现在)，...):
当我尝试加载库 Raster 时，我收到如下错误: 错误:inDL(x, as.logic(local), as.logic(now), ...) 中的“raster”的包或命名空间加载失败:无法加载
PowerShell 逻辑
望着help section about_Comparison_Operators of PowerShell我是这样理解的: PS C:\> $false,$false -eq $true PS C
JavaScript 逻辑
我刚刚修改了旧代码，现在似乎没有任何效果。请您指导我哪里出错了。一些不起作用的事情是: 以前，焦点始终停留在屏幕上唯一的输入字段上。 (现在不行了)，代码中的 if else 条件也不起作用。 On
Javascript持久化>逻辑
请帮我找到一个使用普通 'ol javascript 的解决方案(我无法使用外部框架)。此外，CSS :hover 选择器不适用于现实世界的实现。注册事件发生的事情设置所有调用最后注册事件数组项。
平衡负载和确定总线路线的算法/逻辑
我想创建一个软件来为残障 child 交通规划公交路线(及其最佳载客量)。这些总线具有以下规范: m 个座位(最多 7 个 - 因为有司机和助理) o 轮椅“座位”(最多 4 个) 固定的最大负载量
excel - VBA日期+逻辑
有人能帮我吗？似乎我的 for 逻辑根本不起作用，因为它一直在上午 12:00 返回我的开始时间这是我的代码 Sub forlogic() Dim i As Single Dim t
powerbi - 如何为来自同一数据集的属性在两个切片器之间设置 OR 逻辑？
我正在尝试设置 OR两个切片器过滤器之间的逻辑。两个切片器来自相同的数据集。以下是更多详细信息: 我的源表: 带切片器的视觉效果: 我的目标是，如果我从切片器 1 和切片器 2 中选择任何值，我的视觉
C位 bool 逻辑
我有以下 C 语句: int res = x & (x ^ y); 有没有办法做同样的事情，但每次只使用一次x和y？例如: x | (~x & y) == x | y 最佳答案是的，通过扩展 xo
Azure 逻辑 - 错误网关
我正在创建 Azure 逻辑应用程序以将新的 Sharepoint 文件添加到 Azure Blob。 Sharepoint 由我的公司运行，我使用我的凭据登录来为逻辑应用程序创建 Sharepoin
boolean 逻辑 - 简化产品的总和
我有一个问题要求为给定函数合成最简单的乘积表达式总和。基本上，如果 AB == CD，则函数为 1，否则为 0，结果如下: (!A && !B && !C && !D) || (!A && B &&
c - 逻辑 NOT (!) 运算符不适用于按位语句
我正在尝试确定是否可以在不溢出的情况下计算两个 32 位整数的总和，同时仅使用某些按位运算符和其他运算符。因此，如果整数 x 和 y 可以相加而不会溢出，则以下代码应返回 1，否则返回 0。 ((((
refactoring - 如何简化复杂的业务 "IF"逻辑？
处理乍一看需要许多嵌套 if 语句的复杂业务逻辑的好方法是什么？例子: 折扣券。可能: 1a) 超值折扣 1b) 百分比折扣 2a) 正常折扣 2b) 累进折扣 3a) 需要访问优惠券 3b) 不需
iphone - 逻辑 - 需要找到重叠的数字范围
假设我有一个“numbers”对象数组，其中包含“startNo”整数和“endNo”整数。数组中可以有多个“数字”，我想获取一个包含修改对象的新数组，该数组仅具有不重叠的范围。例如:如果数组有:
javascript - 从数组中查找事件索引(逻辑)
我在这个问题上遇到了困难。我正在使用 JavaScript。我有一个文本区域，用于检测 @ 输入并将其位置存储在数组中。 var input = "@a @b @c" //textarea var
SQL IN 子句中包含 AND 逻辑
默认 IN 使用 OR 基本逻辑。有没有办法在范围内使用 AND 基本逻辑。例如下面的查询 SELECT ItemId,CategoryID FROM ItemCategories WHERE Ca
JavaScript 逻辑/行为错误
我想在您将鼠标悬停在网站图像上时添加叠加层。我在这里实现了这个，它工作正常http://jsfiddle.net/stujLbjh/ 这是js代码: var divs = document.query
c++ - 逻辑/算术移位越少位越快吗？
这个问题在这里已经有了答案: Which is faster: x>2 是否比 x>>31 快？换句话说，sar x, 2 是否比 sar x, 31 快？我做了一些简单的测试，他们似乎有相同的速度
grails - Grails条件查询检查 `OR`逻辑
我有grails criteriaQuery，我在这里再次检查OR逻辑，就像这样一个状态变量: or { eq("status", Status.ONE) eq("status",
grails - Grails条件查询检查 `OR`逻辑
我有grails criteriaQuery，我在这里再次检查OR逻辑，就像这样一个状态变量: or { eq("status", Status.ONE) eq("status",

行者123

个人简介

我是一名优秀的程序员,十分优秀！

作者热门文章

滴滴打车优惠券免费领取

全站热门文章

首页

博学

6Ren·AI

商城

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

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

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

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