之前看到不少人思维导图，我也来凑凑热闹。用幕布做了一个，以做为这篇的笔记的总提纲

什么是方差分析

方差分析分析(Analysis of Variance)，简写为ANOVA，不仅是一种方法，更是一种分析思路，是变异分解的思路。这种思路不仅可以用于多组均值差异的比较，也可以用于其他统计学方法中，比如线性回归、Logistic回归中也有方差分析 ——白话说统计

方差分析的基本思想

将数据的总变异分解为来源于不同因素的相应变异，从而明确各个变异因素在总变异中所占的重要程度

方差分析一般基于两类误差：随机误差和系统误差

随机误差：是指同一因素下，样本各观察值之间的差异，这种差异可以看做随机因素所引起的

系统误差：是指不同因素下，样本各观察值之间的差异，这种差异可能是由于抽样的随机性所引起的，也可能是因素所造成的(也就是系统性因素所造成的)

方差分析实质比较的是以上两类误差，这误差可以用组内(within groups)/组间(between groups)离差平方和表示

考虑到离差平方和会随着样本数增加而增大，所以将离差平方和转变为方差来表示

进而将其中的误差方差作为和其他因素方差比较的标准，从而推断总变异是由误差引起的还是由因素导致的

但是在方差分析中，我们是将所有样本响应变量的方差称为Total Sum of Squares(SST)，也叫总离差平方和，全部观测值与总平均值的离差平方和，反映全部观测值的离散情况

由因素不同水平间差异引起的(可以由模型中因素解释的部分)方差称为Model Sum of Squares(SSM)，也叫组间离差平方和，各组观测值的平均值与总平均值的离差平方和，反映各组样本均值之间的差异程度，包括随机误差和系统误差

由抽本过程本身所引起的部分方差称为Error Sum of Squares，(SSE)，也叫组内离差平方和，各组观测值与其组平均值的离差平方和，反映组内各观测值的离散情况，也反映了随机误差的大小

\[

SST=SSM+SSE

\]

如果我们想衡量上述SSM和SSE中哪个占显著比例，可以通过衡量上述两部分比例大小的统计量F

从上述离差平方和的公式(翻书哈)可看出，其大小跟观测值的数目有关，因此为了消除观测值数目对其的影响，需要将其平均，也就是转化为方差(离差平方和除以对应自由度)；SST的自由度为n-1， 其中n为全部观测值的个数，SSM的自由度为k-1， 其中k为因子水平的个数，SSE的自由度为n-k

那么统计量F的计算方式如下：

\[

F=\frac{SSM/(k-1)}{SSE/(n-k)}

\]

根据模型的自由度(k-1)以及误差自由度(n-k)，可以确定一个F分布；再由上述公式计算出的F0进一步确定P值；再根据显著性水平来判断是否能拒绝原假设

方差分析的前提假设：

样本数据独立

每组数据的总体服从正态分布

每组数据方差齐性

样本数据独立与否比较好判断，根据实验设计来即可确定，保证每组数据之间无关联即可

正态性检验：

使用Shapiro-Wilk正态检验方法来检验样本是否符合正态分布

使用Q-Q图来检验正态性

方差齐性检验：

Bartlett检验，适用于数据服从正态分布；而当数据非正态时则容易导致假阳性

Levene检验，在非正态数据表现较为稳健，对正态不敏感

Fligner-Killeen检验，非参数检验，完全不依赖已知分布

离群点检验：

由于方差分析对离群点很敏感，所以需要对数据检测是否有离群点，最常用(可能也较为好使的)是以图形展示，如箱线图、散点图等形式：

aq.plot()作图

car包的outlierTest()检测

在实际应用，特别是生物统计中，往往并不要求数据严格服从正态分布(因为方差分析对正态不是很敏感，特别是对单因素方差分析)；如果数据近似服从正态分布(当样本数足够多的时候，一般认为观测的分布是服从正态的)，或者对数据进行一定的变化，使其接近于正态分布

如果数据不符合方差齐性，可以参看方差分析为何要进行方差齐次检验？中提到的Handbook of Biological Statistics，作者认为在每组样本数平衡的情况下，单因素方差分析对方差齐性不敏感，除非标准差相差3倍以上并且样本数小于10个，不然假阳性还是比较低的。当然，如果样本数不平衡以及标准差相差较大的情况下，可以考虑使用Welch's anova，R语言的oneway.test()函数以将其封装进去了

如果正态分布也不满足，方差齐性也不满足，可以使用非参的秩和检验kruskal.test()，但是有时针对单因素方差分析而言，kruskal.test()并不好使，还不如使用Welch's anova

单因素方差分析

单因素方差分析是方差分析中最容易理解的一种方法了，其主要为了检验在一个因素的多个水平下各组均值的差异

以R中的数据集PlantGrowth为例，来将上面所说的内容都实现下

先初步展示下数据，其包含了植物在空白(ctrl)以及两个处理下(trt1，trt2)的植物的重量数据

data

>head(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

再看下散点图分布

library("ggpubr")

ggline(data, x = "group", y = "weight",

add = c("mean_se", "jitter"),

order = c("ctrl", "trt1", "trt2"),

ylab = "Weight", xlab = "Treatment")

oneway

接着正态检验(对3个水平下的每组数据都做一次正态检验)，可看出P值均大于0.05，都满足正态分布

group_data

>unlist(lapply(group_data, function(x){

shapiro.test(x)$p.value

}))

# ctrl trt1 trt2

0.7474734 0.4519440 0.5642519

去其中ctrl组做个Q-Q图

library(car)

qqPlot(group_data[[1]])

Q-Qplot

方差齐性检验，使用car包的leveneTest()，可看出，P值大于0.05，满足方差齐性

> leveneTest(weight~group, data = data)

Levene's Test for Homogeneity of Variance (center = median)

Df F value Pr(>F)

group 2 1.1192 0.3412

27

查看是否有离群点

> outlierTest(lm(weight~group, data=data))

No Studentized residuals with Bonferonni p < 0.05

Largest |rstudent|:

rstudent unadjusted p-value Bonferonni p

17 2.537341 0.017509 0.52526

最后使用aov函数进行单因素方差分析，结果显示P值小于0.05，说明不同处理对植物的重量有显著影响

> summary(aov(weight~group, data = data))

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

pvalue F)"][1]

> pvalue

[1] 0.01590996

Df为自由度，Sum Sq为总方差和，Mean Sq为平均方差和，F value为F统计量，Pr(>F)为P值

如果方差不齐性(并且每组样本数也不平衡)，那么可以使用oneway.test()，Welch's anova相当于根据每组方差加了个不同的权重调整了F统计量

> oneway.test(weight~group, data = data, var.equal = F)

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

data: weight and group

F = 5.181, num df = 2.000, denom df = 17.128, p-value = 0.01739

> oneway.test(weight~group, data = data, var.equal = F)$p.value

[1] 0.01739282

单因素方差分析中还有一种是单因素协方差分析，简单的提一下，以multcomp包中litter数据集为例，该数据集将怀孕小鼠分为4组，每组接受不同剂量的药物处理(dose值)，dose为自变量，生出的幼崽体重(weight)为因变量，怀孕时间(gesttime)为协变量，那么用aov()函数协方差分析如下：

> summary(aov(weight~gesttime+dose,data = litter))

Df Sum Sq Mean Sq F value Pr(>F)

gesttime 1 134.3 134.30 8.049 0.00597 **

dose 3 137.1 45.71 2.739 0.04988 *

Residuals 69 1151.3 16.69

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

这里aov的函数顺序有点要求，需要将协方差写在前面，接着是自变量，然后是双因素的交互项等等(如果是双因素协方差分析的话。。。)，从上结果可看出：gesttime与weight相关，在控制gesttime后，dose与weight相关

获得去除协变量效应后的组均值：

> effect("dose", aov(weight~gesttime+dose,data = litter))

dose effect

dose

0 5 50 500

32.35367 28.87672 30.56614 29.33460

协方差分析和方差分析一样，都需要正态性和同方差性假设，前者还假定回归斜率相同，所以需要检验下回归斜率的同质性，以本数据集为例，原假设则是斜率相同，其意思就是每组通过怀孕时间来预测出生体重的回归斜率都相同

> summary(aov(weight~gesttime*dose,data = litter))

Df Sum Sq Mean Sq F value Pr(>F)

gesttime 1 134.3 134.30 8.289 0.00537 **

dose 3 137.1 45.71 2.821 0.04556 *

gesttime:dose 3 81.9 27.29 1.684 0.17889

Residuals 66 1069.4 16.20

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

通过上述结果可看出：gesttime和dose交互作用不显著，因此原假设成立；因为如果交互作用显著，则说明怀孕时间与出生体重的关系依赖于药物剂量dose，那么需要使用非参协方差分析函数(不需要假设回归斜率同质性)，如sm包的sm.ancova()函数

用HH包可视化下结果，ancova()函数可以绘制因变量、协变量和因子之间的关系图

> ancova(weight~gesttime+dose,data = litter)

Analysis of Variance Table

Response: weight

Df Sum Sq Mean Sq F value Pr(>F)

gesttime 1 134.30 134.304 8.0493 0.005971 **

dose 3 137.12 45.708 2.7394 0.049883 *

Residuals 69 1151.27 16.685

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

oneway_ANCOVA

从上图可看出：用怀孕时间来预测出生体重的回归线相互平行(因为之间检验的了回归斜率的同质性，因此只是截距项不同)，随着怀孕时间增加，幼崽出生体重也会增加；还可以看到dose与weight的关系，0剂量组截距项最大，5剂量组截距项最小

PS. 若用anvova(weight~gesttime*dose),生成的图形将允许斜率和截距项依据组别而发生变化，这对可视化那些违背回归斜率同质性的实例非常有用

双因素方差分析

如果不考虑双因素之间的交互作用的话，那么双因素方差分析和单因素方差分析没啥区别

以一个R数据集为例来看看如何分析双因素的交互作用(有重复)

data

> str(data)

'data.frame': 60 obs. of 3 variables:

$ len : num 4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...

$ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...

$ dose: num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...

可看出ToothGrowth数据集的dose列并没有设为因子，而后续双因素方差分析其中一个因素就是dose，因此需将其设为因子(不然结果是错误的！)

data$dose

粗略用点线图展示下

library("ggpubr")

ggline(data, x = "dose", y = "len", color = "supp",

add = c("mean_se", "dotplot"),

palette = c("#00AFBB", "#E7B800"))

twoway

接着就是用aov()来双因素方差分析

> summary(aov(len ~ supp * dose, data = data))

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

从上述结果看跟单因素差别在于其还有了个交互作用的结果(主要其原假设中就是假设没交互作用)，各个因素的自由度也是其水平数减1，supp:dose对应的行则是其交互作用的结果

HH包可视化展示

interaction2wt(len~supp*dose)

twoway_jiaohu

这张图分别展示了任意顺序的因素和其交互作用

资源：

还有一个大神总结的方差分析资料方差分析，可以另存为下来本地慢慢看

参考资料：