R：线性混合效应模型的事后比较

Question

我运行对我的线性混合效应模型的事后比较有疑问。我将尝试用一个快速构建的不完美示例来解释它：

这是我的示例数据：

Variable<-as.factor(rep(c(1,2,3),5))
Random<-rep(c(1,2,2),5)
Result<-rnorm(15,mean=10,sd=2)

Data<-as.data.frame(cbind(Variable,Random,Result))

我的模型中实际上包含了几个固定和随机效应，但这足以说明我的问题：

library(lme4)
LME=lmer(Result~Variable+(1|Random))
summary(LME)

查看固定效应输出，与截距相比，我仅了解变量不同水平的重要性

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)   9.5104     1.3685 12.0000   6.949 1.54e-05 ***
Variable2     0.9155     1.9354 12.0000   0.473    0.645    
Variable3     1.7386     1.9354 12.0000   0.898    0.387

但是，我现在想比较变量级别 1 和级别 2 以及变量级别 2 和级别 3，所以我尝试了以下操作：

library(multcomp)
summary(glht(LME, linfct=c("Variable2-Variable1=0","Variable3-Variable2=0")))

给我留下这个错误：

Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'object' in selecting a method for function 'summary': multcomp:::chrlinfct2matrix: variable(s) ‘Variable1’ not found

如果我排除变量级别 1 并且只查看 2 与 3 的比较，代码工作正常：

summary(glht(LME, linfct=c("Variable3-Variable2=0")))

     Simultaneous Tests for General Linear Hypotheses

Fit: lmer(formula = Result ~ Variable + (1 | Random))

Linear Hypotheses:
                           Estimate Std. Error z value Pr(>|z|)
Variable3 - Variable2 == 0   0.8231     1.6694   0.493    0.622
(Adjusted p values reported -- single-step method)

我也可以运行 linfct函数与Tukey对比：

summary(glht(LME, linfct= mcp(Variable="Tukey")),test=adjusted("none"))

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = Result ~ Variable + (1 | Random))

Linear Hypotheses:
           Estimate Std. Error z value Pr(>|z|)
2 - 1 == 0   0.9155     1.9354   0.473    0.636
3 - 1 == 0   1.7386     1.9354   0.898    0.369
3 - 2 == 0   0.8231     1.6694   0.493    0.622
(Adjusted p values reported -- none method)

鉴于我对 3 与 1 的比较不感兴趣，我只会使用其他 2 个 p 值并逐步调整它们，但这并不是我真正想要的解决方案。我的数据结果不仅仅是这里显示的两个比较，所以带有 Tukey 对比的选项会让我有很多我并不真正感兴趣的比较。

有没有办法从 LME 得到 Variable1？在固定效果中，它作为拦截包含在内，将 Variable1 替换为 Intercept 或我能想到的任何组合都没有成功。还是通常有更好的方法来实现我正在寻找的比较？

如有任何帮助，我们将不胜感激！

Answer 1

你真的已经得到了你想要的。在这种情况下，由于Variable=1是参考组，它的系数固定为0，方差为0。所以，Variable1=Variable2=0是否真的只是Variable2=0的测试。同样 Variable3。您可以从以下两段代码产生相同的输出这一事实看出这一点：

summary(glht(LME, linfct=c("Variable2=0","Variable3=0", "Variable3-Variable2=0")))

# Simultaneous Tests for General Linear Hypotheses
# Fit: lmer(formula = Result ~ Variable + (1 | Random))
# Linear Hypotheses:
#                            Estimate Std. Error z value Pr(>|z|)
# Variable2 == 0              -0.6524     2.0145  -0.324    0.942
# Variable3 == 0              -2.0845     2.0145  -1.035    0.545
# Variable3 - Variable2 == 0  -1.4321     1.1199  -1.279    0.396
# (Adjusted p values reported -- single-step method)

summary(glht(LME, linfct=mcp(Variable="Tukey")))

# Simultaneous Tests for General Linear Hypotheses
# Multiple Comparisons of Means: Tukey Contrasts
# Fit: lmer(formula = Result ~ Variable + (1 | Random))
# Linear Hypotheses:
#            Estimate Std. Error z value Pr(>|z|)
# 2 - 1 == 0  -0.6524     2.0145  -0.324    0.942
# 3 - 1 == 0  -2.0845     2.0145  -1.035    0.545
# 3 - 2 == 0  -1.4321     1.1199  -1.279    0.396
# (Adjusted p values reported -- single-step method)

因此，如果您只想与参考进行调整后的比较，您可以这样做：

summary(glht(LME, linfct=c("Variable2=0","Variable3=0")))

# Simultaneous Tests for General Linear Hypotheses
# Fit: lmer(formula = Result ~ Variable + (1 | Random))
# Linear Hypotheses:
#                Estimate Std. Error z value Pr(>|z|)
# Variable2 == 0  -0.6524     2.0145  -0.324    0.889
# Variable3 == 0  -2.0845     2.0145  -1.035    0.404
# (Adjusted p values reported -- single-step method)

R：线性混合效应模型的事后比较

R: Posthoc comparison for linear mixed effects model

r

lme4

posthoc