如何对 lmer 模型进行事后测试？

机器算法验证 r lme4-nlme 事后
2022-02-01 22:18:19

这是我的数据框：

Group   <- c("G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G1","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G2","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3","G3")
Subject <- c("S1","S2","S3","S4","S5","S6","S7","S8","S9","S10","S11","S12","S13","S14","S15","S1","S2","S3","S4","S5","S6","S7","S8","S9","S10","S11","S12","S13","S14","S15","S1","S2","S3","S4","S5","S6","S7","S8","S9","S10","S11","S12","S13","S14","S15")
Value   <- c(9.832217741,13.62390117,13.19671612,14.68552076,9.26683366,11.67886655,14.65083473,12.20969772,11.58494621,13.58474896,12.49053635,10.28208078,12.21945867,12.58276212,15.42648969,9.466436017,11.46582655,10.78725485,10.66159358,10.86701127,12.97863424,12.85276916,8.672953949,10.44587257,13.62135205,13.64038394,12.45778874,8.655142642,10.65925259,13.18336949,11.96595556,13.5552118,11.8337142,14.01763101,11.37502161,14.14801305,13.21640866,9.141392359,11.65848845,14.20350364,14.1829714,11.26202565,11.98431285,13.77216009,11.57303893)

data <- data.frame(Group, Subject, Value)

然后我运行一个线性混合效应模型来比较 3 组在“值”上的差异，其中“主题”是随机因素：

library(lme4)
library(lmerTest)
model <- lmer (Value~Group + (1|Subject), data = data)
summary(model)

结果是：

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept) 12.48771    0.42892 31.54000  29.114   <2e-16 ***
GroupG2     -1.12666    0.46702 28.00000  -2.412   0.0226 *  
GroupG3      0.03828    0.46702 28.00000   0.082   0.9353

但是，如何比较 Group2 和 Group3？学术文章中的约定是什么？

3个回答

您可以使用emmeans::emmeans()或lmerTest::difflsmeans()， 或multcomp::glht()。

我更喜欢emmeans（以前lsmeans）。

library(emmeans)
emmeans(model, list(pairwise ~ Group), adjust = "tukey")

下一个选项是difflsmeans。注意difflsmeans不能纠正多重比较，并且默认使用 Satterthwaite 方法来计算自由度，而不是 emmeans 默认使用的Kenward-Roger 方法，因此最好明确指定您喜欢的方法。

library(lmerTest)
difflsmeans(model, test.effs = "Group", ddf="Kenward-Roger")

Hack-R 在这个问题的另一个答案中描述了该multcomp::glht()方法。

此外，您可以通过加载lmerTest然后使用anova.

library(lmerTest)
lmerTest::anova(model)

为了清楚起见，您打算对每个主题评估 3 次价值，对吗？看起来 Group 是“学科内”，而不是“学科间”。

拟合lmer模型后，您可以对模型对象执行 ANOVA、MANOVA 和多重比较过程，如下所示：

library(multcomp)
summary(glht(model, linfct = mcp(Group = "Tukey")), test = adjusted("holm"))

   Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = Value ~ Group + (1 | Subject), data = data)

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)  
G2 - G1 == 0 -1.12666    0.46702  -2.412   0.0378 *
G3 - G1 == 0  0.03828    0.46702   0.082   0.9347  
G3 - G2 == 0  1.16495    0.46702   2.494   0.0378 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- holm method)

至于学术论文中的惯例，这将因领域、期刊和特定主题而有很大差异。因此，对于这种情况，只需查看相关文章并查看它们的作用。

为什么不使用模型中的拟合值在您的组之间使用 holm 或 bonferroni 校正进行成对 t.test，因为您看到您的 group2 在您的线性模型中变化很大？然后，您可以从数据中对所有 3 组进行比较。

在这种情况下，你可以写：

PT <- pairwise.t.test(fitted.Values,Group, method=bonferroni)

其它你可能感兴趣的问题
上一篇如何使用主成分分析白化数据？ 下一篇机器学习或深度学习算法能否用于“改进”MCMC 技术的采样过程？