第十七讲 R-双向方差分析2

在“R与生物统计专题”中，我们会从介绍R的基本知识展开到生物统计原理及其在R中的实现。以从浅入深，层层递进的形式在投必得医学公众号更新。

第十六讲 R-双向方差分析1主要解说了双向方差分析的概念和前期数据展现及假设条件验证；

第十七讲我们将介绍具体针对平衡设计的实验和不平衡设计实验的方差分析R实现。

首先，我们来回忆一下平衡设计和不平衡设计的随机区组设计的概念。

平衡设计和不平衡设计

• 当单元内的样本大小相等时，我们称为所谓的平衡设计。在这种情况下，可以应用标准的双向方差分析。比如，5个区组中每个区组内的处理组小鼠和对照组中小鼠都是10只。

• 当自变量每个级别内的样本量不相同时（不平衡设计的情况），应采用不同的方差分析方法，后文有具体说明。比如，区组间处理组和对照组小鼠数量不一样多。

1. 将数据导入R

在这里，我们将使用名为ToothGrowth的内置R数据集。它包含一项评估维生素C对豚鼠牙齿生长的影响的研究数据。

实验在60只豚鼠上进行，其中每只豚鼠通过两种维生素C补充供应方式（橙汁，OJ，或抗坏血酸，VC），分别接受三种剂量水平的维生素C量（0.5、1和2 mg /天）。实验者测量了牙齿生长的长度，数据示例如下所示。

# 导入R内自带的ToothGrowth数据集
library(datasets)
data(ToothGrowth)
# 将数据存储在变量my_data中
my_data <- ToothGrowth

2 双向方差分析：针对平衡设计

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值和对应p值（Pr）

从方差分析表中，我们可以得出结论，供应方式和剂量均具有统计学意义。剂量是最重要的因素变量（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.2 解释结果

从方差分析结果中，您可以基于p值和0.05的显着性水平得出以下结论：

• supp的p值为0.000429（显着），这表明供应方式与牙齿长度相关。

• 剂量的p值<2e-16（显着），表明剂量水平与牙齿长度有关。

• supp * dose之间相互作用的p值为0.02（显着），这表明剂量与牙齿长度之间的关系依赖于补充方式的选择。

2.3 计算数据的统计量

• 使用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.4 各组均值之间的多重成对比较

在方差分析检验中，显着的p值表示某些组均值不同，但我们不知道哪些组对不同。

可以执行多组间两两成对比较，以确定特定两组之间的平均差异是否具有统计学差异。

2.4.1 Tukey多重成对比较

由于方差分析检验很重要，我们可以计算Tukey HSD（Tukey Honest Significant Differences，R函数：TukeyHSD()），以在各组均值之间进行多次两两比较。函数TukeyHD()将拟合的方差分析（res.aov3）作为参数。

我们不需要对“ supp”变量进行测试，因为它只有两个级别，方差分析测试已经证明这两个级别有显着差异。因此，Tukey HSD测试将仅针对因素变量“ dose”进行。

TukeyHSD(res.aov3, which = "dose")

结果输出

Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = len ~ supp + dose + supp:dose, data = my_data)
$dose
diff       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.4.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.4.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 双向方差分析：针对不平衡设计

3.1 数据分析

不平衡的设计存在于各个单元间研究对象数量不相等的情况。

在不平衡设计中，运行方差分析有三种不同的方法。它们被称为I型，II型和III型平方和。为简单起见，请注意，推荐的方法是III型平方和。

注意

平衡设计时，这三种方法给出的结果相同。但是，当设计不平衡时，它们不会给出相同的结果。

car包中的函数Anova()可用于计算非平衡设计的双向方差分析测试。

首先将软件包安装在您的计算机上。

在R中，键入

#安装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: len
Sum 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

对结果的解释与前面类似。该你解释啦～

好了，本期讲解就先到这里。小伙伴们赶紧试起来吧。

在之后的更新中，我们会进一步为您介绍R的入门，以及常用生物统计方法和R实现。欢迎关注，投必得医学手把手带您走入R和生物统计的世界。

提前打个预告，下一期我们将学习“多元方差分析”。

当然啦，R语言的掌握是在长期训练中慢慢积累的。一个人学习太累，不妨加入“R与统计交流群”，和数百位硕博一起学习。

快扫二维码撩客服，

带你进入投必得医学交流群，

让我们共同进步！

↓↓

- END -

长按二维码关注「投必得医学」，更多科研干货在等你！