第十六讲主要解说了双向方差分析的概念和前期数据展现及假设条件验证；第十七讲我们将介绍具体针对平衡设计的实验和不平衡设计的实验的方差分析的R实现。
首先，我们来回忆一下平衡设计和不平衡设计的随机区组设计的概念。
- 当单元内的样本大小相等时，我们具有所谓的平衡设计。在这种情况下，可以应用标准的双向方差分析。比如，5个区组中每个区组内的处理组小鼠和对照组中小鼠都是10只。
- 当自变量每个级别内的样本量不相同时（不平衡设计的情况），应采用不同的方差分析方法。比如，区组间处理组和对照组小鼠数量不一样多。
1将数据导入R
在这里，我们将使用名为ToothGrowth的内置R数据集。它包含一项评估维生素C对豚鼠牙齿生长的影响的研究数据。实验在60只豚鼠上进行，其中每只豚鼠通过两种递送方法（橙汁，OJ，或抗坏血酸，VC）分别接受三种剂量水平的维生素C量（0.5、1和2 mg /天， VC）。实验者测量了牙齿生长的长度，数据示例如下所示。
# 导入R内自带的ToothGrowth数据集
library(datasets)
data(ToothGrowth)
# 将数据存储在变量my_data中
my_data <- ToothGrowth2 R分析
2.1 双向方差分析：针对平衡设计
研究牙齿的长度是否取决于供应方式和剂量。
R函数aov()可以用来回答这个问题。函数summary.aov()用于总结方差分析。
res.aov2 <- aov(len ~ supp + dose, data = my_data)
summary(res.aov2)Df Sum Sq Mean Sq F value Pr(>F)
supp 1 205.4 205.4 14.02 0.000429 ***
dose 2 2426.4 1213.2 82.81 < 2e-16 ***
Residuals 56 820.4 14.7
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1输出包括列F值和PR（> F）对应于分析的p值。
从方差分析表中，我们可以得出结论，供应方式和剂量均具有统计学意义。剂量是最重要的因素变量（P值小于2e-16）。这些结果使我们相信，改变递送方式（供应方式）或维生素C的剂量将显着影响平均牙齿长度。
不能将上述拟合模型称为加性模型，因为这需要假设两个因素变量是相互独立的。如果您认为这两个变量可能相互作用以产生协同效应／交互作用，请按如下所示用星号（*）代替加号（+）。
# 双向方差分析（考虑交互作用）
# 以下两种代码方式性质一样
res.aov3 <- aov(len ~ supp * dose, data = my_data)
res.aov3 <- aov(len ~ supp + dose + supp:dose, data = my_data)
summary(res.aov3)Df Sum Sq Mean Sq F value Pr(>F)
supp 1 205.4 205.4 15.572 0.000231 ***
dose 2 2426.4 1213.2 92.000 < 2e-16 ***
supp:dose 2 108.3 54.2 4.107 0.021860 *
Residuals 54 712.1 13.2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1可以看出，两个主要作用（补给和剂量）及其相互作用均具有统计学意义。
请注意，在交互作用不显着的情况下，应使用加性模型对最后结果进行解释。
2.1.1解释结果
从方差分析结果中，您可以基于p值和0.05的显着性水平得出以下结论：
- supp的p值为0.000429（显着），这表明供应方式与牙齿长度相关。
- 剂量的p值<2e-16（显着），表明剂量水平与牙齿长度有关。
- supp * dose之间相互作用的p值为0.02（显着），这表明剂量与牙齿长度之间的关系还依赖于补充方式。
2.1.2计算数据的统计量
- 使用dplyr R软件包按组计算均值和SD：
require("dplyr")
group_by(my_data, supp, dose) %>%summarise(count = n(),mean = mean(len, na.rm = TRUE),sd = sd(len, na.rm = TRUE))
Source: local data frame [6 x 5]
Groups: supp [?]supp dose count mean sd(fctr) (fctr) (int) (dbl) (dbl)
1 OJ D0.5 10 13.23 4.459709
2 OJ D1 10 22.70 3.910953
3 OJ D2 10 26.06 2.655058
4 VC D0.5 10 7.98 2.746634
5 VC D1 10 16.77 2.515309
6 VC D2 10 26.14 4.797731- 也可以如下使用function model.tables()：
model.tables(res.aov3, type="means", se = TRUE)2.2 各组均值之间的多重成对比较
在方差分析检验中，显着的p值表示某些组均值不同，但我们不知道哪些组对不同。
可以执行多个成对比较，以确定特定组对之间的平均差异是否具有统计显着性。
2.2.1 Tukey多重成对比较
由于方差分析检验很重要，我们可以计算Tukey HSD（Tukey Honest Significant Differences，R函数：TukeyHSD()），以在各组均值之间进行多次成对比较。函数TukeyHD()将拟合的方差分析作为参数。
我们不需要对“ supp”变量进行测试，因为它只有两个级别，方差分析测试已经证明这两个级别有显着差异。因此，Tukey HSD测试将仅针对因素变量“ dose”进行。
TukeyHSD(res.aov3, which = "dose")Tukey multiple comparisons of means95% family-wise confidence level
Fit: aov(formula = len ~ supp + dose + supp:dose, data = my_data)
$dosediff lwr upr p adj
D1-D0.5 9.130 6.362488 11.897512 0.0e+00
D2-D0.5 15.495 12.727488 18.262512 0.0e+00
D2-D1 6.365 3.597488 9.132512 2.7e-06- diff：两组平均值之间的差异
- lwr，upr：置信区间的上下端点为95％（默认值）
- p adj：调整后的多个比较的p值。
从输出中可以看出，所有成对比较均具有显着性，调整后的p值<0.05。
2.2.2 使用multcomp软件包进行多次比较
可以使用函数glht()[在multcomp包中]对方差分析执行多个比较。glht代表一般线性假设检验。简化格式如下：
glht(model, lincft)
- model：拟合模型，例如aov()返回的对象。
- lincft()：要测试的线性假设内容，它为通过mcp()函数输出的结果。
使用glht()执行多个成对比较：
library(multcomp)
summary(glht(res.aov2, linfct = mcp(dose = "Tukey")))Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: aov(formula = len ~ supp + dose, data = my_data)
Linear Hypotheses:Estimate Std. Error t value Pr(>|t|)
D1 - D0.5 == 0 9.130 1.210 7.543 <1e-05 ***
D2 - D0.5 == 0 15.495 1.210 12.802 <1e-05 ***
D2 - D1 == 0 6.365 1.210 5.259 <1e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)2.2.3 成对t检验
pairwise.t.test()函数还可以用于计算组级别之间的成对比较，并对多次检验进行矫正。
pairwise.t.test(my_data$len, my_data$dose,p.adjust.method = "BH")Pairwise comparisons using t tests with pooled SD
data: my_data$len and my_data$dose D0.5 D1
D1 1.0e-08 -
D2 4.4e-16 1.4e-05
P value adjustment method: BH 3 双向方差分析：针对不平衡设计
不平衡的设计存在于各个单元间研究对象数量不相等。
在不平衡设计中运行方差分析的方法有三种不同的方法。它们被称为I型，II型和III型平方和。为简单起见，请注意，推荐的方法是III型平方和。
平衡设计时，这三种方法给出的结果相同。但是，当设计不平衡时，它们不会给出相同的结果。
car包中的函数Anova()可用于计算非平衡设计的双向方差分析测试。
首先将软件包安装在您的计算机上。在R中，键入install.packages（‘car’）。
library(car)
my_anova <- aov(len ~ supp * dose, data = my_data)
Anova(my_anova, type = "III")
Anova Table (Type III tests)
Response: lenSum Sq Df F value Pr(>F)
(Intercept) 1750.33 1 132.730 3.603e-16 ***
supp 137.81 1 10.450 0.002092 **
dose 885.26 2 33.565 3.363e-10 ***
supp:dose 108.32 2 4.107 0.021860 *
Residuals 712.11 54
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1如果您觉得我说的对您有帮助，请点赞让我感到您的支持，您的支持是我写作最大的动力~
ijournal：高颜值的期刊检索网站，助您快速找到理想目标期刊（weixin小程序也上线了哦）
投必得：全专业中英文论文润色编辑助力您的论文快速发表，点击了解业务详情
专栏传送门：
投必得科研软件安装使用手册；投必得：SCI期刊介绍与选择；投必得，教你写论文；投必得统计分析大讲堂；投必得科研生活解忧杂货店