我正在使用不同的 R 包（effects, ggeffects, emmeans, lmer）来计算线性混合模型中边际均值的置信区间。我的问题是，effects与其他方法相比，该软件包产生的 CI 更小。这是一个例子：
library(effects)
library(ggeffects)
library(emmeans)
library(lmerTest)
library(ggplot2)
ham$Age_scale <- scale(ham$Age)
# fit model without the intercept
fit <- lmer(Informed.liking ~ -1 + Product + Age_scale +
(1|Consumer/Product), data=ham)
型号总结：
Random effects:
Groups Name Variance Std.Dev.
Product:Consumer (Intercept) 3.1464 1.7738
Consumer (Intercept) 0.3908 0.6252
Residual 1.6991 1.3035
Number of obs: 648, groups: Product:Consumer, 324; Consumer, 81
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
Product1 5.809e+00 2.327e-01 3.110e+02 24.960 <2e-16 ***
Product2 5.105e+00 2.327e-01 3.110e+02 21.936 <2e-16 ***
Product3 6.093e+00 2.327e-01 3.110e+02 26.180 <2e-16 ***
Product4 5.926e+00 2.327e-01 3.110e+02 25.464 <2e-16 ***
Age_scale 9.621e-03 1.311e-01 7.900e+01 0.073 0.942
然后，我使用不同的方法来计算 的边际均值的 CI Product：
term = 'Product'
## ggpredict
c0 <- as.data.frame(ggpredict(fit, terms = term))
c0 <- c0[,4:5]
## confint.merMod
c1 <- confint(fit, method='profile')
c1 <- c1[4:7,]
## confint.merMod
c2 <- confint(fit,method='Wald')
c2 <- c2[4:7,]
## confint.merMod
c3 <- confint(fit,method='boot')
c3 <- c3[4:7,]
## effect
c4 <- with(effect(term,fit),cbind(lower,upper))
## emmeans,'kenward-roger'
c5 <- with(summary(emmeans(fit,spec=term)),cbind(lower.CL,upper.CL))
我将所有 CI 放在一起：
tmpf <- function(method,val) {
data.frame(method=method,
v=LETTERS[1:4],
setNames(as.data.frame(tail(val,4)),
c("lwr","upr")))
}
allCI <- rbind(tmpf("ggpredict",c0),
tmpf('profile',c1),
tmpf('wald',c2),
tmpf('boot',c3),
tmpf("effects",c4),
tmpf("emmeans",c5))
ggplot(allCI,aes(v,ymin=lwr,ymax=upr,colour=method))+
geom_linerange(position=position_dodge(width=0.8),size=1) + theme_bw()
包产生的置信区间effects比其他的要小得多。我注意到之前有人问过类似的问题，但它的答案表明由 生成的 CIeffects与其他人非常相似。我想知道这里发生了什么。我错过了什么？
更新：@Daniel 运行完全相同的代码，发现效果包没有偏差。ggeffects包0.8.0 版中的ggeffect()和ggemmeans()也产生了非常相似的 CI。
但是，对于效果（4.0.1 版）和ggeffects （0.7.0 版） ，我仍然得到较小的 CI，置信水平为 0.95。
## effect
eff = effect(term,fit)
c4 <- with(eff,cbind(lower,upper))
## emmeans,'kenward-roger'
c5 <- with(summary(emmeans(fit,spec=term)),cbind(lower.CL,upper.CL))
## ggeffect
c6 <- ggeffect(fit, terms = term)
c6 <- c6[,4:5]
packageVersion('ggeffects')
‘0.7.0’
packageVersion('effects')
‘4.0.1’
eff$confidence.level
0.95
然而，在升级到最新版本后，所有方法都产生了相似的 CI。
packageVersion('ggeffects')
‘0.8.0’
packageVersion('effects')
‘4.1.0’
## effect
eff = effect(term,fit)
c4 <- with(eff,cbind(lower,upper))
## emmeans,'kenward-roger'
c5 <- with(summary(emmeans(fit,spec=term)),cbind(lower.CL,upper.CL))
## ggeffect
c6 <- ggeffect(fit, terms = term)
c6 <- c6[,4:5]
## ggemmeans (only the 0.8.0 version supports this)
c7 <- ggemmeans(fit, terms = term)
c7 <- c7[,4:5]
