方差分析是小编工作中最常用到的统计分析方法，以往对于这个分析方法总是懵懵懂懂，在学习《R语言实战》的过程中发现只有对统计知识的原理学深学透，才能在运用R进行分析时做到知其然并知其所以然。因此，本文通过参考论著，百度搜索相关知识点等方式对方差分析的原理进行了梳理，且通过案例一步一步进行最原始的计算，再在R中运行，对比结果以理解原理。
方差分析的假设
原假设H0：不同组的均值相等。
备择假设H1：至少一个样本均值与其他均值不相等。
方差分析的3个前提假定条件
独立性：每个样本的个体相互独立且随机从总体中获得。
正态性：每个组的值均服从正态分布。
方差齐性：每个组的样本的方差之间没有差异。
方差分析的基本思路
方差分析的基本思路是将数据波动（变异）分解为若干部分，除了有一部分代表随机误差，其余每个部分的变异分别代表了某个影响因素的作用(包括交互作用形成的因素)。通过比较因素所致的变异与随机误差的大小，借助F分布和F统计量做出推断：该因素对因变量的影响是否显著存在。F统计量=组间方差/组内方差。
当F统计量<1时，表示比较的样本平均值之间没有显著差异，F统计量>1时，表示组间至少存在一组均值与其他组均值间的差异很大。
在判断是否显著时，常常给定显著性水平，
分布对应的临界值为
，当
时，拒绝
。
方差分析推导公式
根据方差分析原理，总误差可分解为组间误差和组内误差，即：
：总误差，总离差平方和，其值越大，表示测定指标值之间的差异越大；
：组间误差，组间离差平方和，各个水平下样本平均值与数据总平均差异的平方和，反映的是试验因素的水平理论平均值不同而带来的影响；
：组内误差，组内离差平方和，各个水平下，样本观察值与样本均值差异的平方和。
式中，是因素的平均平方和，
是对应的自由度，
为试验因素设定的水平数。
式中，是随机误差的平均平方和，
是对应的自由度，
为样本数，也就是观测值个数，
为因素设定的水平数。
F统计量的计算。
当 时，接受原假设
，否则接受备择假设
。
是给定的显著水平。
方差分析推导
以数据集df为例。 df数据集中nitrogen包含两个水平，N1和N2，试验中测定变量v1，在两个水平下分别测定12次，现在我们来检验试验测定的变量v1在N1和N2条件下的均值是否差异显著？
df <- read.table(file = "D:/Documents/R wd/df.csv", header = T, sep = ",") # 数据导入。
df # 查看数据。
## year nitrogen variety block v1 v2 v3 v4 v5
## 1 2020 N1 a 1 1.26 2.14 0.4 5.0 3.25
## 2 2020 N1 a 2 1.20 2.90 0.1 5.3 1.27
## 3 2020 N1 a 3 1.30 3.00 0.3 5.6 2.24
## 4 2020 N1 b 1 1.08 1.72 1.8 2.8 1.00
## 5 2020 N1 b 2 1.05 1.65 1.7 2.5 3.12
## 6 2020 N1 b 3 1.15 1.35 1.5 3.1 4.57
## 7 2020 N2 a 1 1.32 3.78 1.6 6.0 5.85
## 8 2020 N2 a 2 1.28 4.32 1.4 6.1 6.48
## 9 2020 N2 a 3 1.35 3.95 1.3 6.2 7.21
## 10 2020 N2 b 1 1.33 3.47 2.8 4.1 6.56
## 11 2020 N2 b 2 1.28 2.72 2.4 4.3 8.43
## 12 2020 N2 b 3 1.30 3.90 2.2 4.5 7.55
## 13 2021 N1 a 1 1.19 3.61 0.8 6.0 3.11
## 14 2021 N1 a 2 1.21 3.29 0.5 5.7 2.54
## 15 2021 N1 a 3 1.24 3.26 0.7 5.6 1.28
## 16 2021 N1 b 1 1.09 2.71 1.8 4.0 3.24
## 17 2021 N1 b 2 1.28 2.32 1.6 4.2 1.27
## 18 2021 N1 b 3 1.35 1.95 1.3 4.3 1.15
## 19 2021 N2 a 1 1.45 4.35 1.8 7.2 5.74
## 20 2021 N2 a 2 1.40 3.80 1.2 7.0 6.85
## 21 2021 N2 a 3 1.37 4.23 1.6 6.8 7.42
## 22 2021 N2 b 1 1.28 2.72 2.4 5.1 8.20
## 23 2021 N2 b 2 1.15 3.35 2.5 5.5 5.70
## 24 2021 N2 b 3 1.24 3.46 2.7 4.9 6.00
研究对象信息表
| nitrogen | v1( |
合计 |
平均 |
|---|---|---|---|
| 1.26 1.20 1.30 1.08 1.05 1.15 1.19 1.21 1.24 1.09 1.28 1.35 | 14.40 | 1.20 | |
| 1.32 1.28 1.35 1.33 1.28 1.30 1.45 1.40 1.37 1.28 1.15 1.24 | 15.75 | 1.31 | |
| 合计 |
下面按照方差分析的思路进行一步一步的推导。
1、提出假设
H0：v1在N1和N2条件下的均值相等。
H1：v1在N1和N2条件下的均值不相等。
2、计算各组平均值
df7 <- df[order(df$nitrogen),] # 先对df以nitrogen进行排序。
df7 # 返回排序后结果。
## year nitrogen variety block v1 v2 v3 v4 v5
## 1 2020 N1 a 1 1.26 2.14 0.4 5.0 3.25
## 2 2020 N1 a 2 1.20 2.90 0.1 5.3 1.27
## 3 2020 N1 a 3 1.30 3.00 0.3 5.6 2.24
## 4 2020 N1 b 1 1.08 1.72 1.8 2.8 1.00
## 5 2020 N1 b 2 1.05 1.65 1.7 2.5 3.12
## 6 2020 N1 b 3 1.15 1.35 1.5 3.1 4.57
## 13 2021 N1 a 1 1.19 3.61 0.8 6.0 3.11
## 14 2021 N1 a 2 1.21 3.29 0.5 5.7 2.54
## 15 2021 N1 a 3 1.24 3.26 0.7 5.6 1.28
## 16 2021 N1 b 1 1.09 2.71 1.8 4.0 3.24
## 17 2021 N1 b 2 1.28 2.32 1.6 4.2 1.27
## 18 2021 N1 b 3 1.35 1.95 1.3 4.3 1.15
## 7 2020 N2 a 1 1.32 3.78 1.6 6.0 5.85
## 8 2020 N2 a 2 1.28 4.32 1.4 6.1 6.48
## 9 2020 N2 a 3 1.35 3.95 1.3 6.2 7.21
## 10 2020 N2 b 1 1.33 3.47 2.8 4.1 6.56
## 11 2020 N2 b 2 1.28 2.72 2.4 4.3 8.43
## 12 2020 N2 b 3 1.30 3.90 2.2 4.5 7.55
## 19 2021 N2 a 1 1.45 4.35 1.8 7.2 5.74
## 20 2021 N2 a 2 1.40 3.80 1.2 7.0 6.85
## 21 2021 N2 a 3 1.37 4.23 1.6 6.8 7.42
## 22 2021 N2 b 1 1.28 2.72 2.4 5.1 8.20
## 23 2021 N2 b 2 1.15 3.35 2.5 5.5 5.70
## 24 2021 N2 b 3 1.24 3.46 2.7 4.9 6.00
xT <- mean(df7$v1) # v1总平均值。
xT # 返回结果。
## [1] 1.25625
xA <- mean(df7$v1[1:12]) # N1条件下v1的平均值。
xA # 返回结果。
## [1] 1.2
xE <- mean(df7$v1[13:24]) # N2条件下v1的平均值。
xE # 返回结果。
## [1] 1.3125
3、计算各误差平方和
SST <- sum((df7$v1-xT)^2) # 总平方和SST。
SST # 返回结果。
## [1] 0.2383625
SSA <- 12*(xA-xT)^2+12*(xE-xT)^2 # 组间平方和。
SSA # 返回结果。
## [1] 0.0759375
SSE <- sum((df7$v1[1:12]-xA)^2)+sum((df7$v1[13:24]-xE)^2)# 组内平方和。
SSE # 返回结果。
## [1] 0.162425
4、计算F统计量
SST的自由度为，
为观测值的个数；
SSA的自由度为，
为因素水平的个数；
SSE的自由度为。
MSA <- SSA/(2-1) # 组间均方误差。
MSA # 返回结果。
## [1] 0.0759375
MSE <- SSE/(24-2) # 组内均方误差。
MSE # 返回结果。
## [1] 0.007382955
F <- MSA/MSE # 计算统计量F。
F # 返回F值结果。
## [1] 10.28552
经查表，，本例中计算得10.28552，大于临界值，因此拒绝原假设，接受备择假设，说明v1在N1和N2条件下的平均值是有显著差异的。
5、p值计算
library(stats) # 调用stats包，这个包应该是集成在R基础包中的吧。
1-pf(10.29,1,22) # pf中第一项填计算得的F值，第2项为组间自由度，第3项为组内差异的自由度。
## [1] 0.004056791
可以看出，p值小于0.05，说明v1在N1和N2条件下的平均值间差异显著。
下面按照《R语言实战》中的操作步骤再来分析一下，对比一下结果。
9.3 单因素方差分析
单因素方差分析中，你感兴趣的是比较分类因子定义的两个或多个组别中的因变量均值。
对于完全随机设计试验且处理数大于2时可以用单因素方差分析（等于2 时用t检验）。
df # 显示数据df。
## year nitrogen variety block v1 v2 v3 v4 v5
## 1 2020 N1 a 1 1.26 2.14 0.4 5.0 3.25
## 2 2020 N1 a 2 1.20 2.90 0.1 5.3 1.27
## 3 2020 N1 a 3 1.30 3.00 0.3 5.6 2.24
## 4 2020 N1 b 1 1.08 1.72 1.8 2.8 1.00
## 5 2020 N1 b 2 1.05 1.65 1.7 2.5 3.12
## 6 2020 N1 b 3 1.15 1.35 1.5 3.1 4.57
## 7 2020 N2 a 1 1.32 3.78 1.6 6.0 5.85
## 8 2020 N2 a 2 1.28 4.32 1.4 6.1 6.48
## 9 2020 N2 a 3 1.35 3.95 1.3 6.2 7.21
## 10 2020 N2 b 1 1.33 3.47 2.8 4.1 6.56
## 11 2020 N2 b 2 1.28 2.72 2.4 4.3 8.43
## 12 2020 N2 b 3 1.30 3.90 2.2 4.5 7.55
## 13 2021 N1 a 1 1.19 3.61 0.8 6.0 3.11
## 14 2021 N1 a 2 1.21 3.29 0.5 5.7 2.54
## 15 2021 N1 a 3 1.24 3.26 0.7 5.6 1.28
## 16 2021 N1 b 1 1.09 2.71 1.8 4.0 3.24
## 17 2021 N1 b 2 1.28 2.32 1.6 4.2 1.27
## 18 2021 N1 b 3 1.35 1.95 1.3 4.3 1.15
## 19 2021 N2 a 1 1.45 4.35 1.8 7.2 5.74
## 20 2021 N2 a 2 1.40 3.80 1.2 7.0 6.85
## 21 2021 N2 a 3 1.37 4.23 1.6 6.8 7.42
## 22 2021 N2 b 1 1.28 2.72 2.4 5.1 8.20
## 23 2021 N2 b 2 1.15 3.35 2.5 5.5 5.70
## 24 2021 N2 b 3 1.24 3.46 2.7 4.9 6.00
str(df) # 查看数据结构。
## 'data.frame': 24 obs. of 9 variables:
## $ year : int 2020 2020 2020 2020 2020 2020 2020 2020 2020 2020 ...
## $ nitrogen: chr "N1" "N1" "N1" "N1" ...
## $ variety : chr "a" "a" "a" "b" ...
## $ block : int 1 2 3 1 2 3 1 2 3 1 ...
## $ v1 : num 1.26 1.2 1.3 1.08 1.05 1.15 1.32 1.28 1.35 1.33 ...
## $ v2 : num 2.14 2.9 3 1.72 1.65 1.35 3.78 4.32 3.95 3.47 ...
## $ v3 : num 0.4 0.1 0.3 1.8 1.7 1.5 1.6 1.4 1.3 2.8 ...
## $ v4 : num 5 5.3 5.6 2.8 2.5 3.1 6 6.1 6.2 4.1 ...
## $ v5 : num 3.25 1.27 2.24 1 3.12 4.57 5.85 6.48 7.21 6.56 ...
df$nitrogen <- as.factor(df$nitrogen) # 设置nitrogen为因子。
table(df$nitrogen) # 显示nitrogen各组样本大小。
##
## N1 N2
## 12 12
aggregate(df$v1, by=list(df$nitrogen), FUN=mean) # nitrogen各组的均值。
## Group.1 x
## 1 N1 1.2000
## 2 N2 1.3125
aggregate(df$v1, by=list(df$nitrogen), FUN=sd) # nitrogen各组的标准差。
## Group.1 x
## 1 N1 0.09332251
## 2 N2 0.07782556
fit10 <- aov(df$v1 ~ df$nitrogen) # 检验组间差异。
summary(fit10) # 显示结果。
## Df Sum Sq Mean Sq F value Pr(>F)
## df$nitrogen 1 0.07594 0.07594 10.29 0.00406 **
## Residuals 22 0.16243 0.00738
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
这里可以看到，前面小编一步一步核算的结果，只用一行代码就得到了，df为自由度，组间自由度为1，组内自由度为22；组间离均差平方和Sum Sq中df-nitrogen为0.07594，对应前面单独核算的SSA，组内离均差平方和Sum Sq中Residuals为0.16243，对应前面单独核算的SSE；组间均平方和Mean Sq中df-nitrogen为0.07594，对应前面单独核算的MSA，组内均平方和Mean Sq中Residuals为0.00738，对应前面单独核算的MSE，F value为10.29，对应前面单独核算的F；Pr(>F)为0.00406，对应前面单独核算的1-pf(10.29,1,22) 结果。
gplots包中的plotmeans()可以用来绘制带有置信区间的组均值图形。
library(gplots) # 调用gplots包。
plotmeans(df$v1 ~ df$nitrogen, xlab = "nitrogen", ylab = "v1", main = "Mean Plot\nwith 95% CI") # 绘制各组均值及置信区间。
9.3.1 多重比较
多重比较（multiple comparisons）是指方差分析后对各样本平均数间是否有显著差异的假设检验的统称。方差分析只能判断各总体平均数间是否有差异，多重比较可用来进一步确定哪两个平均数间有差异，哪两个平均数间没有差异。比较方法有N－K（Newman－Keuls）检验、邓肯（DunCan）检验、图基（Tukey）检验邓尼特（Dunnett）检验、最小显著差检验及谢费（Scheffé）检验等它们的理论依据和应用条件都有所不同。
多重比较的方法很多，根据试验设计的目的不同有不同的应用。
若试验设计之初，便明确要比较某几个组均数间是否有差异，称为事前比较。常用的事前比较方法有LSD、Bonferroni和Dunnett法。
若研究目的是方差分析有统计学差异后，想知道哪些组间的均数有差异，便是事后比较。事后比较的常用方法有SNK、Turkey、Scheffe 和Bonferroni法。
LSD检验适用于在专业上有特殊意义的样本均数间的比较，是在设计之初，就已明确要比较某几个组均数间是否有差异。
Bonferroni检验用途最广，几乎可用于任何多重比较的情形。一般认为Bonferroni法是最为保守的。
Dunnett-t检验由multcomp包中glht()函数实现，适用于k-1个试验组与一个对照组均数差异的多重比较。
SNK-q检验适用于多个样本均数两两之间的全面比较，与LSD-t检验相似，可能存在假阳性。
Tukey法较LSD法保守，即较LSD不易发现显著差异。Tukey法要求比较的样本容量相差不大，一般用于样本容量相同的组之间均数的比较。
Scheffe检验在各组样本数相等或不等均可以使用，但是以各组样本数不相等时使用较多。
Duncan法的全称为Duncan's new multiple range test (MRT)，也称为新复极差法。该方法是对SNK法的修正，但是提高了一类错误概率，降低了二类错误的概率，通常用于农业研究。
注意TukeyHSD()函数与本章使用的HH包存在兼容性问题：若载入HH包，TukeyHSD()函数将会失效。使用detach("package::HH")将它从搜寻路径中删除，然后再调用TukeyHSD()。
TukeyHSD(fit10) # 成对组间比较。
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = df$v1 ~ df$nitrogen)
##
## $`df$nitrogen`
## diff lwr upr p adj
## N2-N1 0.1125 0.03975191 0.1852481 0.0040635
par(las = 2) # 扭转轴标签。
plot(TukeyHSD(fit10), xlim=c(0,0.2)) # 成对比较图形。
图形中置信区间包含0的疗法说明差异不显著（p>0.05）。
multcomp包中的glht()函数提供了多重均值比较更为全面的方法，既适用于线性模型，也适用于广义线性模型。
df # 显示数据。
## year nitrogen variety block v1 v2 v3 v4 v5
## 1 2020 N1 a 1 1.26 2.14 0.4 5.0 3.25
## 2 2020 N1 a 2 1.20 2.90 0.1 5.3 1.27
## 3 2020 N1 a 3 1.30 3.00 0.3 5.6 2.24
## 4 2020 N1 b 1 1.08 1.72 1.8 2.8 1.00
## 5 2020 N1 b 2 1.05 1.65 1.7 2.5 3.12
## 6 2020 N1 b 3 1.15 1.35 1.5 3.1 4.57
## 7 2020 N2 a 1 1.32 3.78 1.6 6.0 5.85
## 8 2020 N2 a 2 1.28 4.32 1.4 6.1 6.48
## 9 2020 N2 a 3 1.35 3.95 1.3 6.2 7.21
## 10 2020 N2 b 1 1.33 3.47 2.8 4.1 6.56
## 11 2020 N2 b 2 1.28 2.72 2.4 4.3 8.43
## 12 2020 N2 b 3 1.30 3.90 2.2 4.5 7.55
## 13 2021 N1 a 1 1.19 3.61 0.8 6.0 3.11
## 14 2021 N1 a 2 1.21 3.29 0.5 5.7 2.54
## 15 2021 N1 a 3 1.24 3.26 0.7 5.6 1.28
## 16 2021 N1 b 1 1.09 2.71 1.8 4.0 3.24
## 17 2021 N1 b 2 1.28 2.32 1.6 4.2 1.27
## 18 2021 N1 b 3 1.35 1.95 1.3 4.3 1.15
## 19 2021 N2 a 1 1.45 4.35 1.8 7.2 5.74
## 20 2021 N2 a 2 1.40 3.80 1.2 7.0 6.85
## 21 2021 N2 a 3 1.37 4.23 1.6 6.8 7.42
## 22 2021 N2 b 1 1.28 2.72 2.4 5.1 8.20
## 23 2021 N2 b 2 1.15 3.35 2.5 5.5 5.70
## 24 2021 N2 b 3 1.24 3.46 2.7 4.9 6.00
attach(df)
fit10 <- aov(v1 ~ nitrogen) # 差异显著性检验。
library(multcomp) # 调用multcomp包。
par(mar=c(5,4,6,2)) # 扩展图形边界。
tuk <- glht(fit10, linfct=mcp(nitrogen="Tukey")) # 多重比较。
plot(cld(tuk, level=0.05), col = "red")
有相同字母的组（用箱线图表示）说明均值差异不显著。
结果解读：本研究中N1和N2标注分别为a和b，无相同字母，说明差异显著。
9.3.2 评估检验的假设条件
我们对于结果的信心依赖于做统计检验时数据满足假设条件的程度。单因素方差分析中，我们假设因变量服从正态分布，各组方差相等。
正态性检验
library(car) # 调用car包。
qqPlot(lm(v1 ~ nitrogen), data = df, simulate=TRUE, main="Q-Q Plot", labels=FALSE) # 绘制QQ图判断数据正态性。
数据落在95%的置信区间范围内，说明满足正态性假设。
结果解读：所有数据点均落在置信区间内，说明数据符合正态分布。
方差齐性检验
bartlett.test(v1 ~ nitrogen, data = df) # 方差齐性检验。
##
## Bartlett test of homogeneity of variances
##
## data: v1 by nitrogen
## Bartlett's K-squared = 0.34507, df = 1, p-value = 0.5569
结果解读：p值为0.56，说明N1和N2之间的方差并没有显著不同。
离群点检测
library(car) # 调用car包。
outlierTest(fit10) # 离群点检测。
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 23 -2.127769 0.045366 NA
结果解读：（当p>1时将产生NA），说明本研究中无离群点。
实际应用中，还可以用lm进行差异分析。
summary(lm(v1~nitrogen,data = df)) # 用lm函数进行方差分析。
##
## Call:
## lm(formula = v1 ~ nitrogen, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16250 -0.03687 0.00375 0.05813 0.15000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.20000 0.02480 48.379 < 2e-16 ***
## nitrogenN2 0.11250 0.03508 3.207 0.00406 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08592 on 22 degrees of freedom
## Multiple R-squared: 0.3186, Adjusted R-squared: 0.2876
## F-statistic: 10.29 on 1 and 22 DF, p-value: 0.004063
可以看出最后的p值结果是一致的
因为本例中只有两个处理，小编认为也可以用t检验。
t.test(v1~nitrogen,data = df) # 用t检验检测结果。
##
## Welch Two Sample t-test
##
## data: v1 by nitrogen
## t = -3.2071, df = 21.312, p-value = 0.004178
## alternative hypothesis: true difference in means between group N1 and group N2 is not equal to 0
## 95 percent confidence interval:
## -0.18538442 -0.03961558
## sample estimates:
## mean in group N1 mean in group N2
## 1.2000 1.3125
虽然p值不同，但结果都是差异显著
参考资料：
- 《R语言实战》（中文版），人民邮电出版社，2013.
- 如何理解和使用方差分析，https://zhuanlan.zhihu.com/p/49479042
- 第十四讲 R-单因素方差分析1，https://mp.weixin.qq.com/s?__biz=MzU1Mzc3OTIwNg==&mid=2247488583&idx=2&sn=de31ac2282e9204f8b9500c8071f276a&chksm=fbecfd2ecc9b743815227ace281f424ca97812cc89f9103420e29dbae9117b8405904ddbab3b&scene=21#wechat_redirect
- R语言统计与绘图：组间差异的多重比较，https://www.sci666.com.cn/58293.html
- 方差分析之多重比较，https://zhuanlan.zhihu.com/p/44880434
