R语言学习笔记7_方差分析

七、方差分析

方差分析的主要工作就是将观测数据的总变异按照变异原因的不同分解为因子效应与试验误差，并对其做出数量分析，比较各种原因在总变异中所占的重要程度，以此作为进一步统计推断的依据。
前提条件：独立、正态、方差齐性

7.1 单因子方差分析

7.1.1 数学模型

设试验仅有一个因子（因素）A，其有r个水平A₁,A₂,…,A_r，现在水平A_i下进行n_i次独立观测，得到观测数据X_ij，j=1,2,…,n_i，i=1,2,…,r
SS_T 总离差平方和（总变差）：X_ij与总平均值之差的平方和，描绘所有数据的离散程度。
SS_E 误差平方和（组内平方和）：对固定的水平i，观测值X_ij之间的差异大小的度量。
SS_A 因素A的效应平方和（组间平方和）：因子A各水平下的样本均值和总平均值的平方和。

aov(formula, data = NULL, projections = FALSE, qr = TRUE, contrasts = NULL, ...)

formula 是方差分析的公式，在单因素方差分析中表示为 x~A, data 是数据框。

例： 对除杂方法进行选择。在实验中选用5种不同的除杂方法，每种方法做4次试验，即重复4次，如下表所示。

除杂方法Ai	除杂量Xij	平均量
A1	25.6 22.2 28.0 29.8	26.4
A2	24.4 30.0 29.0 27.5	27.7
A3	25.0 27.7 23.0 32.2	27.0
A4	28.8 28.0 31.5 25.9	28.6
A5	20.6 21.2 22.0 21.2	21.3

> x<- c(25.6, 22.2, 28.0, 29.8,24.4 ,30.0 ,29.0, 27.5,25.0 ,27.7, 23.0 ,32.2,28.8, 28.0, 31.5 ,25.9,20.6, 21.2, 22.0 ,21.2)
> A<-factor(rep(1:5,each=4))
> miscellany<-data.frame(x,A)
> aov.mis<-aov(x~A,data=miscellany)
> summary(aov.mis)
            Df Sum Sq Mean Sq F value Pr(>F)  
A            4    132    33.0    4.31  0.016 *
Residuals   15    115     7.7                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Df 表示自由度，Sum Sq 表示平方和，Mean Sq表示均方和，F value 表示F检验统计量的值，Pr(>F)表示检验p的值， A 就是因素， Residuals 为残差。

可以看出F=4.31>F_0.05(5-1,20-5)=3.06，或p=0.0016<0.05说明有理由拒绝原假设，即认为五种除杂方法的差异显著。

方差分析表

据上述结果可以填写下面的方差分析表：

方差来源	自由度df	平方和SS	均方和MS	F比	p值
因素A	4	132	33.0	4.31	0.016
残差E	15	115	7.7
总和T	19	247

> plot(miscellany$x~miscellany$A)

从上图也可以看出5种除杂方法产生的除杂量有显著差异，特别是第5种与前4种，而方法1和3，方法2和4的差异不明显。

7.1.2 均值的多重比较

之前在进行方差分析后发现各效应的均值之间有显著差异，只能知道有某些均值彼此不同，但不知道具体是哪些均值不同。

多重 t 检验方法

当多次重复使用该方法时，会增大犯第一类错误的概率，从而使得结论不一定可靠。所以在进行较多次重复比较时，要对p值进行调整：

p.adjust(p, method = p.adjust.methods, n = length(p))

p是p值构成的向量，method是修正方法。
当比较次数较多时，Bonferroni方法的效果较好。
得到多重比较的p值：

pairwise.t.test(x, g, p.adjust.method = p.adjust.methods, pool.sd = !paired, paired = FALSE, alternative = c("two.sided", "less", "greater"), ...)

x是响应变量构成的向量，g是分组向量（因子）, p.adjust.methods是上面提到的调整p值的方法，p.adjust.methods=none表示不做任何调整。默认值按Holm方法进行调整。

例：对上例作均值的多重比较，进一步检验。
下面用三种方法进行多重比较：

1）不对p作出调整

> pairwise.t.test(x,A,p.adjust.method = 'none')

	Pairwise comparisons using t tests with pooled SD 

data:  x and A 

  1     2     3     4    
2 0.509 -     -     -    
3 0.773 0.707 -     -    
4 0.289 0.679 0.434 -    
5 0.019 0.005 0.010 0.002

P value adjustment method: none 

2）按默认值”holm“对p值进行调整

> pairwise.t.test(x,A)

	Pairwise comparisons using t tests with pooled SD 

data:  x and A 

  1    2    3    4   
2 1.00 -    -    -   
3 1.00 1.00 -    -   
4 1.00 1.00 1.00 -   
5 0.13 0.04 0.08 0.02

P value adjustment method: holm 

3）按”bonferroni“对p值进行调整

> pairwise.t.test(x,A,p.adjust.method