格洛默的事后测试 [英] Post-hoc test for glmer

问题描述

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

I'm analysing my binomial dataset with R using a generalized linear mixed model (glmer, lme4-package). I wanted to make the pairwise comparisons of a certain fixed effect ("Sound") using a Tukey's post-hoc test (glht, multcomp-package).

大多数它工作正常，但是我的一个固定效果变量("SoundC")根本没有方差(96乘以"1"，零乘以"0")，看来Tukey的测试不能处理.所有与此"SoundC"进行的成对比较，其p值均为1.000，而其中的一些值显然很有意义.

Most of it is working fine, but one of my fixed effect variables ("SoundC") has no variance at all (96 times a "1" and zero times a "0") and it seems that the Tukey's test cannot handle that. All pairwise comparisons with this "SoundC" give a p-value of 1.000 whereas some are clearly significant.

作为验证，我将96个"1"中的一个更改为"0"，此后，我又得到了正常的p值，并且在我期望的位置有明显的差异，而在我使用手册后，差异实际上变小了改变.

As a validation I changed one of the 96 "1"'s to a "0" and after that I got normal p-values again and significant differences where I expected them, whereas the difference had actually become smaller after my manual change.

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

Does anybody have a solution? If not, is it fine to use the results of my modified dataset and report my manual change?

可复制的示例:

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")))

推荐答案

这是一个 CrossValidated 问题的验证；您肯定会看到完全分离，其中您的响应可以完美地划分为0对1的结果.这导致(1)参数的无限值(由于计算不完善，它们仅被列为非无限值)，以及(2)Wald标准错误的疯狂/无用值和相应的$ p $值(这就是您所需要的)在这里看到).在此处，此处和

This is verging on a CrossValidated question; you are definitely seeing complete separation, where there is a perfect division of your response into 0 vs 1 results. This leads to (1) infinite values of the parameters (they're only listed as non-infinite due to computational imperfections) and (2) crazy/useless values of the Wald standard errors and corresponding $p$ values (which is what you're seeing here). Discussion and solutions are given here, here, and here, but I'll illustrate a little more below.

暂时不做统计:您真的不应该尝试只用3个级别来拟合随机效果(请参见例如 http://glmm.wikidot.com/faq )...

To be a statistical grouch for a moment: you really shouldn't be trying to fit a random effect with only 3 levels anyway (see e.g. http://glmm.wikidot.com/faq) ...

经胎儿校正的逻辑回归:

Firth-corrected logistic regression:

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，因为实际上没有区别...)

Standard errors and p-values are now reasonable (p-value for the A vs C comparison is 1 because there is literally no difference ...)

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

Mixed Bayesian model with weak priors:

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)产生

The specification diag(9,3) of the fixed-effect variance-covariance matrix produces

$$ \剩下( \ begin {array} {ccc} 9& 0& 0 \ 0& 9& 0 \ 0& 0& 9 \ end {array} \正确的) $$

$$ \left( \begin{array}{ccc} 9 & 0 & 0 \ 0 & 9 & 0 \ 0 & 0 & 9 \end{array} \right) $$

换句话说，3表示矩阵的维数(等于固定效果参数的数量)，9表示方差-这对应于标准偏差3或95％的范围$ \ pm 6 $，对于logit缩放的响应来说，这是很大/弱/无意义的.

In other words, the 3 specifies the dimension of the matrix (equal to the number of fixed-effect parameters), and the 9 specifies the variance -- this corresponds to a standard devation of 3 or a 95% range of about $\pm 6$, which is quite large/weak/uninformative for logit-scaled responses.

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

These are roughly consistent (the model is very different)

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 ***

如上所述，在这种情况下，我将不推荐混合模型...

As I said above, I would not recommend a mixed model in this case anyway ...

这篇关于格洛默的事后测试的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！