第十四讲 R-单因素方差分析1

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

前面第八到第十三讲，我们介绍了单样本与总体的比较，独立两样本的比较，以及配对两样本的比较的统计学方法及R实现。根据数据是否满足正态分布／方差齐性，又分为参数检验和非参数检验跟大家进行了介绍。

从今天开始，我们将开始进入多组样本相比较的学习。

单因素方差分析
单因素方差分析（one-way analysis of variance, one-way ANOVA），也称单向方差分析（one-factor analysis of variance）,是对两独立样本t检验的延伸（可参考第十讲 R-两独立样本t检验），它用于在存在两个以上的相互独立的组之间进行比较。
在单因素方差分析中，数据基于一个分组变量（也称为因子变量）被分为几组，而要对几组之间有无统计学差异进行比较。

方差分析假设：

• 零假设：不同组的均值相同

• 备择假设：至少一个组样本均值与其他均值不相等。

请注意

如果只有两个组，则可以使用两独立样本t检验（第十讲 R-两独立样本t检验）。在这种情况下，方差分析F检验和t检验是等效的。

1. ANOVA的假设检验条件

方差分析的使用，需要满足以下几个条件：

• 1. 观察值之间在所研究的因子水平，是相互独立且随机从总体中获得的。

• 2. 每个因子水平的数据（每一组）均呈正态分布。

• 3. 每个组的样本的方差之间没有差异 。（Leven’s 检验可用于检查这一点。）

2. 单因素方差分析原理

假设有3组（A，B，C）进行比较：

1. 整个实验多组间的总变异来自两部分：组内变异和组间变异。

2. 组内变异：计算每组数据内各个样本间的差异，这称为组内的方差（S2within）或残差（residual variance）。

3. 组间变异：计算每个组（样本均值）之间的方差，如下所示：

• 计算每组的平均值

• 计算样本均值之间的方差（S2between）

产生F统计量(比率)：S2within/S2between 。

注意当比率<1时，表示比较的样本平均值之间没有显著差异。较高的比率（比率 >1）表示组间至少存在一组均值与其他组均值间的差异很大。

3. 可视化数据并在R中计算单因素方差分析

3.1 将数据导入R

在这里，我们将使用名为PlantGrowth的内置R数据集。它包含对照组和两种不同处理条件下获得的植物重量（共三组）。

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

3.2 检查你的数据

为了了解数据的外观

# 显示前六行
head(my_data)

输出结果

weight group
1   4.17  ctrl
2   5.58  ctrl
3   5.18  ctrl
4   6.11  ctrl
5   4.50  ctrl
6   4.61  ctrl

# 显示每组有多少样本
table(my_data$group)

输出结果

ctrl trt1 trt2
10   10   10

在R术语中，“group”列称为因子(factor)，而不同类别（“ ctrl”，“ trt1”，“ trt2”）称为因子水平。因子默认按字母顺序排列。

# 显示因子水平
levels(my_data$group)

[1] "ctrl" "trt1" "trt2"

如果水平不是自动按照正确的顺序排列，可以按以下顺序重新排列：

my_data$group <- ordered(my_data$group,
levels = c("ctrl", "trt1", "trt2"))

使用dplyr包可以按组计算信息统计量（平均值和标准差）。

• 按组计算摘要统计信息-计数，均值，标准差：

library(dplyr)
group_by(my_data, group) %>%
summarise(
count = n(),
mean = mean(weight, na.rm = TRUE),
sd = sd(weight, na.rm = TRUE)
)

输出结果

# A tibble: 3 x 4
group count  mean    sd
1 ctrl     10  5.03 0.583
2 trt1     10  4.66 0.794
3 trt2     10  5.53 0.443

3.3 可视化数据

3.3.1 使用箱形图可视化数据

library("ggpubr")
ggboxplot(my_data, x = "group", y = "weight",
color = "group", palette = c("red", "blue", "black"),
order = c("ctrl", "trt1", "trt2"),
ylab = "Weight", xlab = "Treatment")

# 箱形图也可以用如下代码
boxplot(weight ~ group, data = my_data,
xlab = "Treatment", ylab = "Weight",
frame = FALSE, col = c("red", "blue", "black"))

3.3.2 使用点图可视化数据

# 添加均数和标准差: mean_se
# (其他选项包括: mean_sd, mean_ci, median_iqr, ....)
library("ggpubr")
ggline(my_data, x = "group", y = "weight",
add = c("mean_se", "jitter"), order = c("ctrl", "trt1", "trt2"),
ylab = "Weight", xlab = "Treatment")

# 显示均数和95%可信区间：mean_ci
ggline(my_data, x = "group", y = "weight",
add = c("mean_ci", "jitter"), order = c("ctrl", "trt1", "trt2"),
ylab = "Weight", xlab = "Treatment")

4. 计算单因素方差分析

我们想知道在3种实验条件下植物的平均重量之间是否存在显着差异。

R函数aov()可以回答这个问题。

函数summary.aov()用于总结方差分析的结果。

# 计算方差分析
res.aov <- aov(weight ~ group, data = my_data)

# 总结分析结果
summary(res.aov)

输出结果

Df Sum Sq Mean Sq F value Pr(>F)
group        2  3.766  1.8832   4.846 0.0159 *
Residuals   27 10.492  0.3886
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

输出包括列F值（4.846）和P（> F）对应于测试的p值 （0.016）。由于p值小于显着性水平0.05，我们可以得出结论，三组中至少存在一组均值与其他组不同。            

在R输出结果以“ *”突出显示的组之间存在显著差异。

在发现至少存在一组与其他组之间存在差异以后，我们想知道具体是哪两组间存在差异。同时，如开头说到，方差分析需要满足几个条件，他们的假设检验条件的检验方法是什么？如果数据不满足假设条件，要怎么办？下一讲，我们会进行详细解说。

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

快扫二维码撩客服，

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

让我们共同进步！

↓↓

- END -

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