r - glmer 的事后检验

Question

我正在使用广义线性混合模型(glmer、lme4-package)用 R 分析我的二项式数据集。我想使用 Tukey 的事后测试(glht，multcomp-package)对某个固定效果(“声音”)进行成对比较。

大部分工作正常，但我的一个固定效果变量(“SoundC”)根本没有变化(“1”的 96 次和“0”的 0 次)，而且 Tukey 的测试似乎无法处理。与此“SoundC”的所有成对比较给出的 p 值为 1.000，而有些则明显显着。

作为验证，我将 96 个“1”中的一个更改为“0”，然后我再次获得了正常的 p 值，并且在我预期的地方出现了显着差异，而在我手动更改后差异实际上变得更小了。

有人有解决方案吗？如果不是，是否可以使用我修改后的数据集的结果并报告我的手动更改？

可重现的例子:

Response <- c(1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,1,1,0,
              0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,
              1,1,0,1,1,0,1,1,1,1,0,0,1,1,0,1,1,0,1,1,0,1,1,0,1)    
Data <- data.frame(Sound=rep(paste0('Sound',c('A','B','C')),22),
                   Response,
                   Individual=rep(rep(c('A','B'),2),rep(c(18,15),2)))


# Visual
boxplot(Response ~ Sound,Data)

# Mixed model
library (lme4)
model10 <- glmer(Response~Sound + (1|Individual), Data, family=binomial)

# Post-hoc analysis
library (multcomp)
summary(glht(model10, mcp(Sound="Tukey")))

Answer 1

这接近于 CrossValidated题;您肯定会看到完全分离，其中将您的响应完美地划分为 0 对 1 结果。这导致 (1) 参数的无限值(由于计算缺陷，它们仅被列为非无限)和 (2) Wald 标准误差的疯狂/无用值和相应的 $p$ 值(这就是您在这里看到)。给出讨论和解决方案here , here , 和 here ，但我将在下面进行更多说明。

暂时成为统计上的提示:无论如何，您真的不应该尝试仅使用 3 个级别来拟合随机效应(参见例如 http://glmm.wikidot.com/faq )......

Firth 校正逻辑回归:

library(logistf)
L1 <- logistf(Response~Sound*Individual,data=Data,
        contrasts.arg=list(Sound="contr.treatment",
         Individual="contr.sum"))

                                 coef se(coef)            p
(Intercept)              3.218876e+00 1.501111 2.051613e-04 
SoundSoundB             -4.653960e+00 1.670282 1.736123e-05 
SoundSoundC             -1.753527e-15 2.122891 1.000000e+00 
IndividualB             -1.995100e+00 1.680103 1.516838e-01 
SoundSoundB:IndividualB  3.856625e-01 2.379919 8.657348e-01 
SoundSoundC:IndividualB  1.820747e+00 2.716770 4.824847e-01

标准误差和 p 值现在是合理的(A 与 C 比较的 p 值是 1，因为实际上没有区别......)

具有弱先验的混合贝叶斯模型:

library(blme)
model20 <- bglmer(Response~Sound + (1|Individual), Data, family=binomial,
                  fixef.prior = normal(cov = diag(9,3)))

##              Estimate Std. Error    z value     Pr(>|z|)
## (Intercept)  1.711485   2.233667  0.7662221 4.435441e-01
## SoundSoundB -5.088002   1.248969 -4.0737620 4.625976e-05
## SoundSoundC  2.453988   1.701674  1.4421024 1.492735e-01

规范书 diag(9,3)固定效应方差-协方差矩阵产生

$$
\剩下(
\begin{数组}{ccc}
9 & 0 & 0\
0 & 9 & 0\
0 & 0 & 9
\end{数组}
\对)
$$

换句话说，3 指定矩阵的维度(等于固定效应参数的数量)，而 9 指定方差——这对应于标准差 3 或约 $\pm 的 95% 范围6$，对于logit-scaled 响应来说是相当大/弱/无信息量的。

这些是 大致一致(模型非常不同)

library(multcomp)
summary(glht(model20, mcp(Sound="Tukey")))

##                     Estimate Std. Error z value Pr(>|z|)    
## SoundB - SoundA == 0   -5.088      1.249  -4.074 0.000124 ***
## SoundC - SoundA == 0    2.454      1.702   1.442 0.309216    
## SoundC - SoundB == 0    7.542      1.997   3.776 0.000397 ***

正如我上面所说，我会 不是 无论如何，在这种情况下推荐混合模型......