方差分析（Analysis of variance, ANOVA）是统计学中的经典方法，1923年由英国统计学家R. A. Fisher提出，其实质是关于观测值变异原因的数量分析，主要用途是研究外界因素或试验条件对观测结果影响的显著性。方差分析方法就是根据平方和的加和性原理，将总平方和分解为不同影响因素产生的平方和以及代表重复观测值之间的随机变化的误差平方和，并利用假设检验（F-检验）手段分析不同影响因素对总平方和贡献的显著性程度，进而判断这些因素是否对观测结果产生明显影响。
方差分析的前提假设有：采样随机性，总体中每个个体被采集的机会均等；个体独立性，每个样本的采集不受上一个样本的影响；分布正态性，无需多言，数理统计中最重要的假定。对于非正态分布的总体，可以通过数据变换或改用非参数检验方法；方差同质性，要求各个样本（同一因素不同水平的一组重复观测值）的方差大小没有显著差异；方差加和性，即一组观测值的总平方和可以分解成不同影响因素造成的平方和之和。
对于完全随机设计试验且处理数大于2时可以用单因素方差分析（等于2 时用t检验）。单因素方差分析（One-way ANOVA）应用举例：
随机在A、B、C三地抽取家蝇，测量翅膀长度，一共五十个样本，数字特征分别为μ1、μ2和μ3。问题：三地家蝇翅膀长度是否有差异？H0假设：μ1=μ2=μ3，即三地家蝇翅膀长度无显著差异；Ha假设：μ1，μ2，μ3不完全相等，即三地家蝇翅膀长度至少有一个与其他样地有显著差异，α=0.05。R代码和结果：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
| # 输入数据site1 <- c(45, 44, 43, 47, 48, 44, 46, 44, 40, 45, 42, 40, 43, 46, 47, 45, 46, 45, 43, 44)site2 <- c(45, 48, 47, 43, 46, 47, 48, 46, 43, 49, 46, 43, 47, 46, 47, 46, 45, 46, 44, 45, 46, 44, 43, 42, 45)site3 <- c(47, 48, 45, 46, 46, 44, 45, 48, 49, 50, 49, 48, 47, 44, 45, 46, 45, 43, 44, 45, 46, 43, 42)fly.survey <- data.frame(length = c(site1, site2, site3), site = factor(c(rep("1", 20), rep("2", 25), rep("3", 23))))# 检查数据options(digits = 3) # default value = 7tapply(fly.survey$length, fly.survey$site, mean)## 1 2 3## 44.4 45.5 45.9tapply(fly.survey$length, fly.survey$site, var)## 1 2 3## 4.56 3.26 4.48boxplot(length ~ site, data = fly.survey, xlab = "Sites", ylab = "Length")# Bartlett Test方差齐性检验（参数）bartlett.test(length ~ site, data = fly.survey)#### Bartlett test of homogeneity of variances#### data: length by site## Bartlett's K-squared = 0.764, df = 2, p-value = 0.6825# 单因子方差分析One Way ANOVAfit <- aov(length ~ site, data = fly.survey)summary(fit)## Df Sum Sq Mean Sq F value Pr(>F)## site 2 26.3 13.15 3.24 0.045 *## Residuals 65 263.4 4.05## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(fit) |
单因子方差分析结果显示F value = 3.24 ，Pr(>F) = 0.045，因此拒绝H0假设，即认为三地家蝇翅膀长度在统计学上有显著差异。
两因素方差分析（Two-way ANOVA）可用于随机区组实验设计，用来分析两个因素的不同水平对结果是否有显著影响，以及两因素之间是否存在交互效应。应用举例：
有一个关于检验毒品强弱的试验，给48只老鼠注射I, II, III三种毒药(因素A)，同时有A, B, C, D 4种治疗方案(因素B)，这样的试验在每一种因素组合下都重复四次测试老鼠的存活时间（年）。试分析毒药和治疗方案以及它们的交互作用对老鼠存活时间有无显著影响。
H01假设：α1=α2=α3=0，即因素A毒药种类对老鼠存活时间没有影响；Ha1假设：α1，α2和α3不全相等，即因素A毒药种类对老鼠存活时间有影响。
H02假设：β1=β2=β3=β4=0，即因素B治疗方案对老鼠存活时间没有影响；Ha2假设：β1，β2，β3和β4不全相等，即因素B治疗方案对老鼠存活时间有影响。
H03假设：δ11=δ12=δ13=δ14=δ21=δ22=δ23=δ24=δ31=δ32=δ33=δ34=0，即A因素和B因素没有交互作用；Ha3假设：δ11，δ12，…，δ34不全相等，即A因素和B因素有交互作用。α=0.05。
R代码和结果：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
| # 输入数据rats <- data.frame(time = c(0.39, 0.45, 0.46, 0.43, 0.72, 0.8, 0.78, 0.92, 0.53, 0.45, 0.63, 0.56, 0.65, 0.71, 0.66, 0.62, 0.36, 0.29, 0.23, 0.33, 0.42, 0.61, 0.79, 1.04, 0.94, 0.85, 0.81, 0.6, 0.76, 0.82, 0.71, 0.65, 0.22, 0.31, 0.28, 0.33, 0.65, 0.67, 0.69, 0.78, 0.43, 0.45, 0.54, 0.48, 0.38, 0.36, 0.41, 0.33), toxicant = gl(3, 16, 48, labels = c("I", "II", "III")), cure = gl(4, 4, 48, labels = c("A", "B", "C", "D")))# 毒素和治疗方案两因素各自效应分析op <- par(mfrow = c(1, 2))plot(time ~ toxicant + cure, data = rats)# Bartlett Test检查方差齐性bartlett.test(time ~ toxicant, data = rats)#### Bartlett test of homogeneity of variances#### data: time by toxicant## Bartlett's K-squared = 4.085, df = 2, p-value = 0.1297bartlett.test(time ~ cure, data = rats)#### Bartlett test of homogeneity of variances#### data: time by cure## Bartlett's K-squared = 6.464, df = 3, p-value = 0.09111# 方差分析rats.aov <- aov(time ~ toxicant * cure, data = rats)summary(rats.aov)## Df Sum Sq Mean Sq F value Pr(>F)## toxicant 2 0.304 0.152 14.87 2.0e-05 ***## cure 3 0.998 0.333 32.47 2.4e-10 ***## toxicant:cure 6 0.307 0.051 4.99 0.00082 ***## Residuals 36 0.369 0.010## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 |
由p值可以判断，拒绝H01和H02假设，即因素A毒素和因素B治疗方案对存活时间有显著影响；同时也拒绝H03假设，即认为因素A和因素B交互作用对老鼠存活时间有显著影响。
参考资料：
[1]陶澍编著.应用数理统计方法[M].北京市：中国环境科学出版社.1994.
[2]明道绪主编.生物统计附试验设计[M].北京市：中国农业出版社.2008.
[3]汤银才主编.R语言与统计分析[M].北京市：高等教育出版社.2008.
用R做t检验
研究者对研究对象总体按一定方法抽样，为了尽可能准确地根据这一样本的参数估计总体，验证某些假设，这就需要用到假设检验，而t检验就是假设检验中最常用的一种检验方法，可以通过R软件非常方便地实现。
第一种情况是单样本t检验（one sample t-test），目的是检验一个正态分布的总体均值是否满足零假设的已知均值。例如：已知A地家蝇翅膀长度均值47mm，为研究B地家蝇翅膀长度和A地有无差异，随机抽取B地100只家蝇，测量其翅膀长度。注意，研究对象是B地所有的家蝇，这个数目是远大于100的，如果直接把这100只的数据平均然后比较是不准确的，统计学上用T检验解决了这个问题。已知B地家蝇翅膀长度符合正态分布，样本符合随机原则。零假设H0：μ=47，即B地家蝇翅膀长度与A地一致；备择假设Ha：μ≠47，即B地家蝇翅膀长度与A地不一致，显著性水平α为0.05。R下的做法是：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
| #输入数据wl <- c(36,37,38,38,39,39,40,40,40,40,41,41,41,41,41,41,42,42,42,42,42,42,42,43,43,43,43,43,43,43,43,44,44,44,44,44,44,44,44,44,45,45,45,45,45,45,45,45,45,45,46,46,46,46,46,46,46,46,46,46,47,47,47,47,47,47,47,47,47,48,48,48,48,48,48,48,48,49,49,49,49,49,49,49,50,50,50,50,50,50,51,51,51,51,52,52,53,53,54,55)#t检验t.test(wl, mu = 47, var.equal = FALSE)#mu表示平均值，var.equal = FALSE表示方差未知One Sample t-testdata: wlt = -3.8269, df = 99, p-value = 0.0002274alternative hypothesis: true mean is not equal to 4795 percent confidence interval:44.72226 46.27774sample estimates:mean of x45.5 |
从上述结果可以看出，基于100个样本的p-value = 0.0002274 < 0.05，因此有充足的理由拒绝H0假设，即B地家蝇翅膀长度与A地不一致，95%置信区间为（44.72226,46.27774），B地家蝇翅膀长度比A地短。
第二种情况是两样本t检验（two-sample t-test），目的是检验两个正态分布的总体均值之间是否有差异（是否相等）。例如：为验证1和2两地家蝇翅膀长度是否有差异，随机在1和2两地抽取家蝇，测量翅膀长度，一共40个样本，数字特征分别为μ1和μ2。零假设H0：μ1=μ2，即两地家蝇翅膀长度无显著差异；备择假设Ha：μ1≠μ2，即两地家蝇翅膀长度有显著差异，显著性水平α为0.05。R下的做法是：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
| # 输入、整理数据lines <- "ID length site\n1 45 1\n2 44 1\n3 45 2\n4 43 1\n5 48 2\n6 47 2\n7 47 1\n8 48 1\n9 43 2\n10 44 1\n11 46 2\n12 46 1\n13 44 1\n14 47 2\n15 40 1\n16 48 2\n17 45 1\n18 46 2\n19 42 1\n20 43 2\n21 40 1\n22 43 1\n23 49 2\n24 46 1\n25 47 1\n26 46 2\n27 43 2\n28 45 1\n29 46 1\n30 47 2\n31 46 2\n32 47 2\n33 45 1\n34 46 2\n35 43 1\n36 45 2\n37 44 1\n38 46 2\n39 44 2\n40 45 2"wl <- read.table(con <- textConnection(lines), header = TRUE)close(con)x = wl[wl$site == 1, 2] # 样地1数据y = wl[wl$site == 2, 2] # 样地2数据# 检查前提条件shapiro.test(x) # 样地1数据正态性检验Shapiro-Wilk normality testdata: xW = 0.9522, p-value = 0.4017#P值为0.4017，说明没有足够的证据拒绝样地1数据来自正态分布> shapiro.test(y) # 样地2数据正态性检验Shapiro-Wilk normality testdata: yW = 0.9397, p-value = 0.2363#P值为0.2363，说明没有足够的证据拒绝样地2数据来自正态分布# Bartlett Test方差齐性检验（参数）bartlett.test(length ~ site, data = wl)Bartlett test of homogeneity of variancesdata: length by siteBartlett's K-squared = 0.9775, df = 1, p-value = 0.3228#t检验t.test(x, y)Welch Two Sample t-testdata: x and yt = -2.4616, df = 36.141, p-value = 0.01874alternative hypothesis: true difference in means is not equal to 095 percent confidence interval:-2.7356655 -0.2643345sample estimates:mean of x mean of y44.35 45.85 |
基于1和2个地点t检验的p-value = 0.01874 < 0.05，因此拒绝H0假设，即认为两地地家蝇翅膀长度有显著差异。
第三种情况是配对t检验或称重复测量t检验（ "paired" or "repeated measures" t-test），用于检验同一统计量的两次测量值之间的差异是否为零。例如测量病人治疗前和治疗后的某种指标水平，检验其差异是否显著；或者要比较同一河流上下游两地污染情况，同一指标在上下游各测量一次，共有N个指标，这都要用到配对T检验。R中配对t检验的代码为t.test(x, y, paired=TRUE)，x和y为两组变量。懒得找数据举具体例子了。
