先载入需要的R包。
library(lmerTest)
library(emmeans)
来看一个简单的例子，这里是用了lmerTest包里的名为ham的data。我们假设产品的种类（分类变量，共4种产品）会影响顾客对产品的喜好（连续变量），具体而言，我们想知道4种产品之间是否有显著的区别。同时，虽然年龄和性别并非我们关心的变量，但因为我们假设这些变量也会影响喜好，因此一并纳入到模型中。最后，顾客被认为是一种随机因子。构建模型如下。
fm <- lmer(Informed.liking ~ Product + Gender + Age + (1|Consumer) , data=ham)
summary(fm)
固定效应的结果如下。
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 6.013956 0.728086 84.304883 8.260 1.79e-12 ***
Product2 -0.703704 0.232456 564.000000 -3.027 0.00258 **
Product3 0.283951 0.232456 564.000000 1.222 0.22240
Product4 0.117284 0.232456 564.000000 0.505 0.61408
Gender2 -0.251114 0.268044 78.000000 -0.937 0.35173
Age -0.001729 0.014086 78.000000 -0.123 0.90262
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
因为代码默认将Product1作为基线，所以模型的结果告诉我们Product2对顾客喜好的影响显著区别于Product1（beta = -0.70, SE=0.23, p = 0.003）。Product3、4以此类推。顺便一说，可以用以下代码在跑模型之前指定基线，例如设置Product4为基线。
ham$Product <- relevel(as.factor(ham$Product), ref=4)
可以用emmeans包进行多重比较，代码和结果如下。可以发现这里给出了模型中的系数、标准误和自由度。唯一的区别就是t值和p值，这是因为emmeans默认使用了Tukey方法进行校正，该方法在组间样本量不一致的被试间设计中具有较强的检验力。
> EMM1 <- emmeans(fm, ~ Product)
> pairs(EMM1)
contrast estimate SE df t.ratio p.value
1 - 2 0.704 0.232 564 3.027 0.0137
1 - 3 -0.284 0.232 564 -1.222 0.6134
1 - 4 -0.117 0.232 564 -0.505 0.9580
2 - 3 -0.988 0.232 564 -4.249 0.0001
2 - 4 -0.821 0.232 564 -3.532 0.0025
3 - 4 0.167 0.232 564 0.717 0.8903
Results are averaged over the levels of: Gender
Degrees-of-freedom method: kenward-roger
P value adjustment: tukey method for comparing a family of 4 estimates
我自己（基于快要遗忘的统计学知识）的理解。模型关注的是每个条件是否显著区别于基线，例如Product2对顾客喜好的影响是否区别于1，又例如Product3是否区别于1，所以可以被认为是做了多个独立的检验，不涉及多重比较。而因为事后比较（即，多重比较）的目的是同时考虑多个对比的区别，由于一类错误（Type one error，也就是false positive）会随着检验次数的增加而上升，所以需要通过校正来排除这种情况。不过，是否需要做事后比较，还是说通过其它方法在模型中得到更多的对比的结果（例如，设置不同的基线，或对变量进行编码），大概还是要根据研究问题和具体的情境来做决定。
如果指定不使用任何校正手段（adjust=NULL），则结果和模型完全一致，如下。这里仅仅是作为解释，实际的数据分析中是需要做校正的。
> EMM1 <- emmeans(fm, ~ Product)
> pairs(EMM1, adjust=NULL)
contrast estimate SE df t.ratio p.value
1 - 2 0.704 0.232 564 3.027 0.0026
1 - 3 -0.284 0.232 564 -1.222 0.2224
1 - 4 -0.117 0.232 564 -0.505 0.6141
2 - 3 -0.988 0.232 564 -4.249 <.0001
2 - 4 -0.821 0.232 564 -3.532 0.0004
3 - 4 0.167 0.232 564 0.717 0.4737
Results are averaged over the levels of: Gender
Degrees-of-freedom method: kenward-roger
也可以用Bonferroni方法进行校正（adjust="bonferroni"），如下。一般而言，该方法在被试内设计中的检验力是高于Tukey方法的（这是因为Tukey方法的方差齐性假设在被试内设计中常常难以实现）。
> EMM1 <- emmeans(fm, ~ Product)
> pairs(EMM1, adjust="bonferroni")
contrast estimate SE df t.ratio p.value
1 - 2 0.704 0.232 564 3.027 0.0155
1 - 3 -0.284 0.232 564 -1.222 1.0000
1 - 4 -0.117 0.232 564 -0.505 1.0000
2 - 3 -0.988 0.232 564 -4.249 0.0002
2 - 4 -0.821 0.232 564 -3.532 0.0027
3 - 4 0.167 0.232 564 0.717 1.0000
Results are averaged over the levels of: Gender
Degrees-of-freedom method: kenward-roger
P value adjustment: bonferroni method for 6 tests
或者FDR方法（adjust="fdr"），如下。
> EMM1 <- emmeans(fm, ~ Product)
> pairs(EMM1, adjust="fdr")
contrast estimate SE df t.ratio p.value
1 - 2 0.704 0.232 564 3.027 0.0052
1 - 3 -0.284 0.232 564 -1.222 0.3336
1 - 4 -0.117 0.232 564 -0.505 0.6141
2 - 3 -0.988 0.232 564 -4.249 0.0002
2 - 4 -0.821 0.232 564 -3.532 0.0013
3 - 4 0.167 0.232 564 0.717 0.5684
Results are averaged over the levels of: Gender
Degrees-of-freedom method: kenward-roger
P value adjustment: fdr method for 6 tests
详细的解释以及更复杂的分析，可参考emmeans的文档[1][2]。
References
Quick start guide for emmeans
Comparisons and contrasts in emmeans
Correction for multiple testing in Multiple regression analysis
Multiple group comparisons in a linear mixed-effect model
Multinomial logistic regression with same DV but different baseline, how to correct p value
Maxwell, S. E., & Delaney, H. D. (2017). Designing experiments and analyzing data: a model comparison perspective, 3rd. Taylor and Francis, New York.
2023/06/28更新
修改了讨论多重比较的段落。
2023/06/29更新
添加了对Tukey和Bonfferoni方法的简单说明。
