假设我们有一个响应变量和一个解释变量（5 个水平）。
con.data <- data.frame(x = c(rnorm(75), c(rnorm(75)+5)),
category = rep(c("A","B1","B2","C1","C2"), each=30))
如果我们做方差分析，然后是经典的TukeyHSD post-hoc ...
m1 <- aov(con.data$x ~ con.data$category); summary(m1)
TukeyHSD(m1)
diff lwr upr p adj
B1-A 0.1877877 -0.9109837 1.286559 0.9897357
B2-A 2.2050543 1.1062829 3.303826 0.0000013
C1-A 4.5951737 3.4964022 5.693945 0.0000000
C2-A 4.7790488 3.6802773 5.877820 0.0000000
B2-B1 2.0172666 0.9184952 3.116038 0.0000117
C1-B1 4.4073860 3.3086145 5.506157 0.0000000
C2-B1 4.5912611 3.4924896 5.690033 0.0000000
C1-B2 2.3901193 1.2913479 3.488891 0.0000001
C2-B2 2.5739944 1.4752230 3.672766 0.0000000
C2-C1 0.1838751 -0.9148963 1.282647 0.9905238
...我们获得所有可能的组合。
如果您只想进行您感兴趣的比较，请使用CONTRASTS。
在 R 中有几个默认的对比矩阵（概述）：
治疗，赫尔默特，总和，多项式......但你可以创建自己的。
首先 -确定您感兴趣的比较。
然后创建一个对比矩阵：
contrasts(con.data$category) <- cbind(c(1,-1/4,-1/4,-1/4,-1/4),
c(0,-1/2,-1/2,1/2,1/2),
c(0,0,0,1/2,-1/2),
c(0,-1/2,1/2,0,0))
查看此表以供参考：
contrasts(con.data$category)
[,1] [,2] [,3] [,4]
A 1.00 0.0 0.0 0.0
B1 -0.25 -0.5 0.0 -0.5
B2 -0.25 -0.5 0.0 0.5
C1 -0.25 0.5 0.5 0.0
C2 -0.25 0.5 -0.5 0.0
请注意，每列的总和等于 0。
如果一个因子有 5 个水平，则由于自由度的原因，只能进行 4 次比较。
在第一列中，您将“A”的平均值与所有其他类别的平均值进行比较。
在第二列中，您仅比较类别 "B" 和 "C" (B1+B2 vs. C1+C2)。
在第三列中，您仅在“C”类别（C1 与 C2）内进行比较。
在第四列中，您仅在“B”类别（B1 与 B2）内进行比较。
要查看结果，请使用创建的对比度矩阵重新制作 ANOVA 。
summary(lm(con.data$x, con.data$category))
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.5910 0.1258 20.599 < 2e-16 ***
con.data$category1 -2.3534 0.2516 -9.355 < 2e-16 ***
con.data$category2 3.4907 0.2813 12.411 < 2e-16 ***
con.data$category3 -0.1839 0.3978 -0.462 0.645
con.data$category4 2.0173 0.3978 5.072 1.19e-06 ***
...并且表中的每一行对应于所做的每个比较（对比矩阵中的列）。（即 con.data$category1 很重要，因此“A”的平均值与所有其他组的平均值之间存在显着差异......等）
简而言之：
尝试制作一个仅包含您感兴趣的比较的对比矩阵。使用上面的示例应该不难。
然而 ！！！
我不会立即对数据使用事后（或对比）。这就像在模型“走”之前教模型“跑”。所以首先我会创建一个包含所有变量的模型。随后，根据边际规则删除所有非显着变量（或其交互作用）。此过程（减少）将确定您是否需要进行所有所需的比较。
所以尝试制作完整的模型：
m1 <- glm(measurement ~ my.profile*my.contrast*my.disease, data = my.data)
anova(m1, test = "F")
例如，如果因素“疾病”不显着（单独或相互作用 - 它不应该包含在 post-hoc.
假设m1模型的结果如下所示：
my.profile 0.00012 **
my.contrast 0.00231 *
my.disease 0.07690 .
my.profile:my.contrast 0.26159
my.profile:my.disease 0.07709 .
my.contrast:my.disease 0.21256
my.profile:my.contrast:my.disease 0.23319
使用边际规则：
更新号 1：
您可以看到三重交互不重要 - 让我们删除它
m2 <- update(m1, ~.-my.profile:my.contrast:my.disease)
现在，显示更新模型的方差分析：
anova(m2, test = "F")
my.profile 0.00012 **
my.contrast 0.00231 *
my.disease 0.07690 .
my.profile:my.contrast 0.26159
my.profile:my.disease 0.07709 .
my.contrast:my.disease 0.21256
更新号 2：
您可以看到双重交互不显着 - 让我们删除它们（从具有最高 p 值的双重交互开始）
m3 <- update(m2, ~.-my.profile:my.contrast)
anova(m3, test = "F")
my.profile 0.00012 **
my.contrast 0.00231 *
my.disease 0.07690 .
my.profile:my.disease 0.07709 .
my.contrast:my.disease 0.21256
更新号 3：
去掉另一个双重交互
m4 <- update(m3, ~.-my.contrast:my.disease)
anova(m4, test = "F")
my.profile 0.00012 **
my.contrast 0.00231 *
my.disease 0.07690 .
my.profile:my.disease 0.07709 .
更新号 4：
去掉最后的双重交互
m5 <- update(m4, ~.-my.profile:my.disease)
anova(m5, test = "F")
my.profile 0.00012 **
my.contrast 0.00231 *
my.disease 0.07690 .
更新号 5：
去除不显着因素
m6 <- update(m5, ~.-my.disease)
anova(m6, test = "F")
my.profile 0.00012 **
my.contrast 0.00231 *
模型 m6 是我们的最终模型。
不幸的是，很明显进行比较（Y.NEG.NO、X.NEG.NO 和其他）是没有意义的，因为三重和双重交互也是不显着的。从 TukeyHSD 中选择所需的行是不正确的（即使这样的事后会显示出显着的差异！！！）。相信我，这样的方法在同行评审过程中很难辩护。因此，您只能在配置文件（X vs. Y）和对比度（NO vs. YES）中进行比较。疾病因素不显着。不要难过——即使是不重要的结果也是结果。
边际规则是奥卡姆剃刀的实际应用（另见 Crawley, MJ Statistics: An Introduction Using R. 2nd ed. Wiley, 2014, Ch. 10 Multiple Regression, p. 195）。一个好的模型总是最简单的——并且解释了数据变化的最大部分。
您可以将此结果作为“完整模型”
（包含所有因素及其相互作用）或“最小充分模型”（MAM，仅包含显着影响）发布在文章中。我更愿意将这两个版本都包含在手稿中，让审稿人决定更喜欢哪一个。
当方差分析结果不显着时，重点不是在事后检验中寻找 p 值。