使用R语言进行单因素方差分析

[VIctoryRoad]

注意：笔者不是统计学专业出身，只是近期由于工作上的一些需要粗略学了一点 R 语言的统计分析。本文也算是个人的一个学习笔记，内容难免有疏漏、错误之处，仅供参考，欢迎大家批评指正、

1 准备数据

将数据和分组依据存入如下二维表中：

grp	y
…	…

其中，

• grp 为分组依据。
• y 是要比较的数据。

这里我们讨论的是单因素方差分析，即只有一个分组依据。或者说，只希望研究某一个因素是否会对某一指标产生影响。

2 分析方法选择

应用不同的方差分析方法，要求数据满足不同的条件。具体请参考下表：

是否服从正态分布？	各组方差是否相等？	可选用的方法
✓	✓	aov() 单因素方差分析
✓	✕	oneway.test()
✕	✕	kruskal.test() 非参数方差分析

3 函数用法

3.1 读取数据

首先我们将 1 中准备好的数据保存为一个 .csv 表格文件。比如 mydata.csv。

然后我们将这些数据保存在变量 mydata 之中。

mydata <- read.csv('mydata.csv')

然后，我们需要将其中的 grp 列的数据类型转换为 “因子”（frame），以便后续进行一些检验。

mydata$grp <- as.factor(mydata$grp)

3.2 执行单因素方差分析

执行如下命令：

aov.manu <- aov(y ~ grp, data=mydata)
summary(aov.manu)

上述命令还有一种简化的写法：
aov(y ~ grp, data=mydata) |> summary()

可得到如下形式的输出：

            Df Sum Sq Mean Sq F value  Pr(>F)
grp          2    520  260.00   9.176 0.00382 **
Residuals   12    340   28.33
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

可以看到，本案例中计算得到的 p 值为 0.00382，小于 0.01，表明各组之间存在显著差异。

为了解各组两两之间的差异有多大，我们需要开展多重比较。多重比较需要使用 multcomp 包。

执行多重比较的命令如下：

library(multcomp, quietly=TRUE)
aov.manu <- aov(y ~ grp, data=d.manu)
glht(aov.manu, linfct = mcp(grp = "Tukey")) |> summary()

可得到如下形式的结果：

The following object is masked from ‘package:MASS’:

    geyser


	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = y ~ grp, data = mydata)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)   
B - A == 0    4.000      3.367   1.188   0.4820   
C - A == 0  -10.000      3.367  -2.970   0.0291 * 
C - B == 0  -14.000      3.367  -4.159   0.0035 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

可以看到，在 0.05 水平下，A、B 两组之间没有显著差异，而 C 与 A、B 之间均存在显著差异。

3.3 执行 oneway.test()

执行如下命令：

oneway.test(y ~ grp, data=mydata)

可得到如下形式的输出：

	One-way analysis of means (not assuming equal variances)

data:  y and grp
F = 8.1976, num df = 2.0000, denom df = 7.9914, p-value = 0.01159

可以看到，本案例中计算得到的 p 值为 0.01159，表明各组数据之间的差异在 0.05 的水平下存在显著性。

对于服从方差不齐、但服从正态分布的样本，也有多重比较方法。可通过如下命令来进行：

with(d.manu, pairwise.t.test(y, grp, pool.sd=FALSE, p.adjust.method="holm"))

可得到如下形式的输出：

	Pairwise comparisons using t tests with non-pooled SD 

data:  y and grp 

  A     B    
B 0.258 -    
C 0.039 0.010

P value adjustment method: holm

矩阵中的数字即为某两组之间的显著性水平 p。

3.4 执行非参数方差分析

执行如下命令：

kruskal.test(y ~ grp, data=mydata)

可得到如下形式的输出：

	Kruskal-Wallis rank sum test

data:  y by grp
Kruskal-Wallis chi-squared = 7.7677, df = 2, p-value = 0.02057

可以看到，本案例中计算得到的 p 值为 0.02057，表明各组数据之间的差异在 0.05 的水平下存在显著性。

在通过非参数方差分析确认各组数据存在统计学意义上的差异后，我们还可以进行 Dunn Test，来比较各组两两之间是否有差异。

执行 Dunn Test 可借助 FSA 程序包。测试命令如下：

library(FSA)
dunnTest(y ~ grp, data = mydata, method = "bh")

其中 method 参数用于指定 p 值的校正方法。包括三种：

• bonferroni
• holm
• bh［控制错误发现率（FDR）］

执行上述命令后，可得到如下形式的输出：

Dunn (1964) Kruskal-Wallis multiple comparison
  p-values adjusted with the Benjamini-Hochberg method.

  Comparison         Z     P.unadj      P.adj
1      A - B -1.062559 0.287981854 0.28798185
2      A - C  1.700095 0.089113081 0.13366962
3      B - C  2.762654 0.005733348 0.01720004

其中 P.adj 是采用选定的方法校正后的 p 值。

参考资料

• 36 R方差分析|R语言教程
  https://www.math.pku.edu.cn/teachers/lidf/docs/Rbook/html/_Rbook/stat-aov.html

附：数据检验

S1 正态性检验

参考资料：https://www.jianshu.com/p/eee4889b6729

方法 1：直方图

执行如下命令：

library(ggplot2)
ggplot(mydata, aes(x = y)) + geom_histogram(aes(color=grp))

其中， mydata 是存放我们的实验数据的变量，y 是我们要检验的数据。

方法 2：密度图

服从正态分布的量，其概率密度分布图应该接近钟形曲线。执行如下命令：

library(ggplot2)
ggplot(mydata, aes(x=y)) + geom_density(aes(color=grp))

其中， mydata 是存放我们的实验数据的变量，y 是我们要检验的数据，grp 是分组依据。上述命令可绘制各组变量的概率密度分布图。

注：

1. + geom_density(aes(color=grp)) 中的 aes(color=grp) 表示用不同颜色的曲线区分不同组的概率密度分布情况。
2. 如果希望使用不同颜色填充不同组的概率密度曲线，请将上述参数改为 + geom_density(aes(fill=grp))。
3. 为了美观起见，还可以使用 alpha= 参数设置透明度。比如要实现半透明效果，可使用 + geom_histogram(aes(fill=grp), alpha=0.4)

方法 3：Q-Q 图

在 Q-Q 图中，数据点越接近参考线，代表数据的实际分布情况与理想分布（这里指的是我们检验的正态分布）越接近。

注意：该方法需要使用的函数在 R 程序包 car 中。该包需要手动安装。在 Ubuntu、Debian 系统上，安装该包前，需要先安装 cmake 软件包。
sudo apt install cmake
然后，在 R shell 中执行如下命令，安装 car 程序包：
install.packages('car')

执行如下命令：

library(car)
qqPlot(mydata$y)

其中， mydata 是存放我们的实验数据的变量，y 是我们要检验的数据。

方法 4：Shapiro-Wilk 检验

执行如下命令：

shapiro.test(mydata$y)

可得到如下形式的输出：

	Shapiro-Wilk normality test

data:  mydata$y
W = 0.93143, p-value = 0.2866

可以看到 p 值 > 0.05，表示这组数据的分布与正态分布没有显著差异，换言之可以认为这组数据服从正态分布。

S2 方差齐性检验

参考资料：https://www.jianshu.com/p/94eb0e5dadf3

aov() 函数要求被检验的各组数据方差相等（即方差齐）。在 R 中，可通过如下方法检验各组数据的方差是否相等。

方法 1：barrlett.test()

执行如下命令：

bartlett.test(mydata$y, mydata$grp)

其中，

• mydata 是之前我们导入的二维表形式的数据。
  • y 是数据所在的列。
  • grp 是分组依据所在的列。

可以得到如下形式的输出：

	Bartlett test of homogeneity of variances

data:  mydata$y and mydata$grp
Bartlett's K-squared = 0.024476, df = 2, p-value = 0.9878

其中，p 值大于 0.05，表示方差齐性。

方法 2：leveneTest()

该方法不仅适用于正态数据，也适用于非正态数据。

执行如下命令：

library(car) #加载 car 程序包
leveneTest(mydata$y, mydaya$grp)

其中，

• mydata 是之前我们导入的二维表形式的数据。
  • y 是数据所在的列。
  • grp 是分组依据所在的列。

可得到如下形式的输出：

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  2  0.0292 0.9713
      12
Warning message:
In leveneTest.default(mydata$y, mydata$grp) : mydata$grp coerced to factor.

其中，p 值大于 0.05，表示方差齐性。

[tiger (小虎王)]

学习 SPSS 课程的回忆开始攻击我