《数学建模：基于R》一一2.2　方差分析

本节书摘来自华章计算机《数学建模：基于R》一书中的第2章，第2.2节,作者：薛　毅 更多章节内容可以访问云栖社区“华章计算机”公众号查看。

2.2　方差分析

方差分析是分析试验数据的一种方法.对于抽样得到的试验数据,由于观测条件不同(同一因素不同水平或不同因素的各个水平)会引起试验结果有所不同;另一方面,由于各种随机因素的干扰,实验结果也会有所不同.由观测条件不同所引起的实验结果的差异是系统的,而随机因素引起的差异是偶然的.
方差分析的目的在于从实验数据中分析出各个因素的影响以及各个因素间的交互影响,以确定各个因素作用的大小,从而将由于观测条件不同而引起实验结果的不同与由于随机因素而引起实验结果的差异以数量的形式区别开来,以确定在试验中有没有系统的因素在起作用.
2.2.1　单因素方差分析
在一个试验中,影响试验结果的因素有很多,如果其他因素能控制在一定的范围之内,为研究方便,可以认为仅有一个因素在变化,只分析单个因素(也称为因子)对试验的影响,这就是单因素方差分析.
在试验中,为了评价试验的性质需要进行多次测量,将测量结果称为指标.因素所处的状态或所取的等级称为水平.
1.数学模型
设因素A有r个水平A1,A2,…,Ar,每个水平Ai进行ni次独立观测,将水平Ai下的试验结果xi1,xi2,…,xini看成来自第i个正态总体Xi~N(μi,σ2)的样本观测值,其中μi,σ2均未知,并且每个总体Xi都相互独立.考虑线性统计模型xij=μi+εij,i=1,2,…,r,　j=1,2,…,ni
εij~N(0,σ2)且相互独立（2.17）其中μi为第i个总体的均值,εij为相应的试验误差.
令αi=μi－μ,其中μ=1n∑ri=1niμi,n=∑ri=1ni,则模型(2.17)可以等价写成xij=μ+αi+εij,i=1,2,…,r,　j=1,2,…,ni
εij~N(0,σ2)且相互独立（2.18）其中μ表示总和的均值,αi表示水平Ai对指标的效应,称模型(2.18)为单因素方差分析的数学模型,它是一种线性模型.
检验因素A的r个水平是否有显著差异归结为检验这r个总体的均值是否相同,即检验假设H0:α1=α2=…=αr=0,　H1:α1,α2,…,αr不全为零（2.19）如果H0被拒绝,则说明因素A的各水平的效应之间有显著的差异；否则差异不显著.
为了导出H0的检验统计量，方差分析法建立在平方和分解和自由度分解的基础上,考虑统计量ST=∑ri=1∑nij=1(xij－x)2,　x=1n∑ri=1∑nij=1xij称ST为总离差平方和(或称为总变差),它是所有数据xij与总平均值x差的平方和,描绘了所有观测数据的离散程度.经计算可以证明如下的平方和分解公式ST=SE+SA（2.20）其中SE=∑ri=1∑nij=1(xij－xi·)2,　xi·=1ni∑nij=1xij,
SA=∑ri=1∑nij=1(xi·－x)2=∑ri=1ni(xi·－x)2SE表示随机误差的影响.这是因为对于固定的i来讲,观测值xi1,xi2,…,xini是来自同一个正态总体N(μi,σ2)的样本.因此,它们之间的差异是由随机误差所导致的.而∑nij=1(xij－xi·)2是这ni个数据的变动平方和,正是它们差异大小的度量.将r组这样的变动平方和相加,就得到了SE,通常称SE为误差平方和或组内平方和.
SA表示在Ai水平下的样本均值与总平均值之间的差异之和,它反映了r个总体均值之间的差异,因为xi·是第i个总体的样本均值,是μi的估计,因此,r个总体均值μ1,μ2,…,μr之间的差异越大,样本均值x1·,x2·,…,xr·之间的差异也就越大.平方和∑ri=1ni(xi·－x)2正是这种差异大小的度量,这里ni反映了第i个总体样本大小在平方和SA中的作用，称SA为因素A的效应平方和或组间平方和.
式(2.20)表明,总平方和ST可按其来源分解成两部分,一部分是误差平方和SE,它是由随机误差引起的.另一部分是因素A的平方和SA,是由因素A的各水平的差异引起的.
由模型假设(2.19),经过统计分析得到E(SE)=(n－r)σ2,即SE/(n－r)为σ2的一个无偏估计且SEσ2~χ2(n－r)如果原假设H0成立,则有E(SA)=(r－1)σ2,即SA/(r－1)也是σ2的无偏估计,且SAσ2~χ2(r－1)并且SA与SE相互独立.因此,当H0成立时,有F=SA/(r－1)SE/(n－r)~F(r－1,n－r)（2.21）于是F(也称为F比)可以作为H0的检验统计量.通过计算P值的方法来决定是接受还是拒绝原假设H0.当P值小于α时拒绝原假设H0;当P值大于α时,则无法拒绝原假设,所以应接受原假设H0.通常将计算结果列成表2.5的形式,称为方差分析表.
表2.5　单因素方差分析表方差来源自由度平方和均方F比P值因素Ar－1SAMSA=SAr－1F=MSAMSEp误差n－rSEMSE=SEn－r总和n－1ST
2.计算
在R中,使用aov()函数完成方差分析,其使用格式为　　aov(formula, data = NULL, projections = FALSE, qr = TRUE,
　　　　contrasts = NULL, ...)参数formula为公式,形如X ~ A.
data为数据框,描述数据的指标、因素和相应的水平,默认值为NULL.
与lm()函数一样,该函数需要用summary()函数提取计算结果(方差分析表).
由模型(2.18),模型(2.24)和模型(2.25)可知,方差分析模型本质上是线性模型的一种,所以可以用lm()函数计算,用anova()函数给出方差分析表.
例2.8　利用4种不同配方的材料A1,A2,A3和A4生产出来的元件,测得其使用寿命如表2.6所示.那么4种不同配方下元件的使用寿命是否有显著差异呢？
表2.6　元件寿命数据材料使用寿命A11600161016501680170017001780A215001640140017001750A316401550160016201640160017401800A4151015201530157016401600
解　　该问题有一个因素(材料)和4个水平(4种不同的配方),属于单因素问题.用数据框的格式输入数据,aov()函数计算,summary()函数提取信息,程序(程序名：exam0208.R)如下：lamp <- data.frame(
　X = c(1600, 1610, 1650, 1680, 1700, 1700, 1780,
1500, 1640, 1400, 1700, 1750, 1640, 1550,
1600, 1620, 1640, 1600, 1740, 1800, 1510,
1520, 1530, 1570, 1640, 1600),
　　 A = factor(rep(1:4, c(7, 5, 8, 6)))
)
lamp.aov <- aov(X ~ A, data = lamp)
summary(lamp.aov)在程序中,数据输入采用数据框结构,其中X为数据,A为对应的因子,factor为因子函数,将变量转化为因子.在aov()函数中,公式X ~ A表示作单因素方差分析,用summary()函数提取方差分析表,其结果如下：　　DfSum SqMean SqF valuePr(>F)
A349212164042.1660.121
Residuals221666227574上述数据与方差分析表2.5中的内容相对应,其中Df表示自由度,Sum Sq表示平方和,Mean Sq表示均方,F value表示F值,即F比.Pr(>F)表示P值,A就是因素A,Residuals是残差,即误差.
从P值(0.121>0.05)可以看出,没有充分的理由拒绝H0,也就是说,4种配方生产出的元件的平均寿命无显著差异.
例2.9　小白鼠在接种了三种不同菌型的伤寒杆菌后的存活天数如表2.7所示.判断小白鼠被注射三种菌型后的平均存活天数有无显著差异？
表2.7　白鼠试验数据菌型存活天数124324772254256851071212663711667955106310
解　这是一个单因素3水平问题,使用lm()函数计算,anova()函数列出方差分析表.程序(程序名: exam0209.R)如下：mouse <- data.frame(
　　X = c( 2,　4, 3,　2,　4, 7, 7, 2,　2, 5, 4, 5, 6, 8,
　　　　　5,10, 7,12,12, 6, 6, 7,11, 6, 6, 7, 9, 5,
　　　　　5,10, 6,　3, 10),
　　A = factor(rep(1:3, c(11, 10, 12)))
)
mouse.lm <- lm(X ~ A, data = mouse)
anova(mouse.lm)计算结果如下：Analysis of Variance Table
Response: X
　　DfSum SqMean SqF valuePr(>F)　　
A294.25647.1288.48370.001202 **
Residuals30166.6535.555在计算结果中,P值远小于0.01,因此,拒绝原假设,即认为小白鼠在接种三种不同菌型的伤寒杆菌后的存活天数有显著差异.
2.2.2　多重均值检验
在单因素方差分析中,如果F检验的结论是拒绝H0,则说明因素A的r个水平效应有显著的差异,也就是说r个均值之间有显著差异.但这并不意味着所有均值之间都存在着显著差异,这时还需要对每一对μi和μj作一对一的比较,即多重比较.
1.多重t检验
通常采用多重t检验方法进行多重比较,这种方法本质上就是针对每组数据进行t检验,只不过估计方差时利用的是全体数据,因而自由度变大.具体地说,要比较第i个总体与第j个总体的均值是否相同,即检验H0:μi=μj,　H1:μi≠μj,　i≠j,`javascript
i,j=1,2,…,r在R中,用pairwise.t.test()函数完成多重t检验,其使用格式为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.method为P值的调整方法,默认值为holm(Holm方法),关于调整方法,将在稍后介绍.
pool.sd为开关变量,表示是否使用公共的标准差.
paired为逻辑变量,表示是否打算作成对数据的t检验,默认值为FALSE.
alternative为备择假设,取"two.sided"(默认值)表示双侧检验;取"less"表示单侧小于,"greater"表示单侧大于.
...为附加参数.
2. P值的调整
多重t检验方法的优点是使用方便,但在均值的多重检验中,如果因素的水平较多,而检验又是同时进行的,则多重t检验会增大犯第一类错误的概率,所得到的“有显著差异”的结论不一定可靠.
为了克服多重t检验方法的缺点,统计学家们提出了许多有效的方法来调整P值,目前已有的修正方法有:Holm修正(1979),Hochberg修正(1988),Hommel修正(1988),Bonferroni修正,Benjamini & Hochberg修正(1995), Benjamini & Yekutieli修正(2001).
在pairwise.t.test()函数中,p.adjust.methods(调整方法)对应的取值分别为:"holm"(默认值), "hochberg"， "hommel"， "bonferroni"， "BH"， "BY"，或者使用同义名"fdr"， "none"表示不作修正.
在R中, p.adjust()函数是专门用于调整P值的函数,其使用格式为p.adjust(p, method = p.adjust.methods, n = length(p))参数p为P值构成的数值向量,允许在向量中使用缺失值(NA).
method为P值调整的方法,p.adjust.methods的取值分别为"holm"(默认值), "hochberg", "hommel","bonferroni","BH","BY"和"fdr"之一.
n为正整数,至少为p的维数.
例2.10（续例2.9）　由于在例2.9中F检验的结论是拒绝H0,应进一步检验H0:μi=μj,　H1:μi≠μj,　i,j=1,2,3,　i≠j.解　首先计算各个因子间的均值,再用pairwise.t.test()作多重t检验.程序和计算结果如下：>attach(mouse)
>tapply(X, A, mean)
123
3.8181827.7000007.083333
>pairwise.t.test(X, A)
　　 Pairwise comparisons using t tests with pooled SD
data:　X and A 
　1　　　2
2 0.0021 -
3 0.0048 0.5458
P value adjustment method: holm程序中的attach()函数是连接函数,连接数据框或列表,这样可以直接使用数据框或列表中的元素.tapply()函数是应用函数,计算各水平的均值.
各水平的样本均值分别为:3.82,7.70和7.08.多重t检验得到的结果是:μ1与μ2,μ1与μ3有显著差异,μ2与μ3没有显著差异,即小白鼠所接种的三种不同菌型的伤寒杆菌中,第一种与后两种使得小白鼠的平均存活天数有显著差异.而后两种的差异不显著.
这里虽然没有指定P值的调整方法,但程序使用默认值Holm修正.如果不打算调整P值,将参数调整为"none",也可以选择其他的修正方法.通过计算会发现,无论采用何种修正P值的方法,调整后P值都会增大.因此,在一定程度上会克服多重t检验方法的缺点.
还可以用plot()函数,其命令为> plot(X ~ A)绘出各水平的箱线图(图形略),这样可以直观地看出,哪些水平之间有显著差异.
**2.2.3　进一步讨论**
1.正态性检验
对于单因素方差分析模型xij=μi+εij,i=1,2,…,r,　j=1,2,…,ni
εij~N(0,σ2i)且相互独立（2.22）需要假设模型的误差项εij服从正态分布,并且满足σ21=σ22=…=σ2r=σ2（2.23）关于正态性检验,可以使用1.4节介绍的分布的检验(如Shapiro-Wilk正态性检验)来检验xij是满足从正态分布.
例2.11　用Shapiro-Wilk正态性检验分析例2.9中的数据是否满足正态性条件.
解　用shapiro.test()函数作正态性检验,程序和计算结果如下：> source("ex```javascript
am0209.R")
>with(lamp, tapply(X, A, shapiro.test))
$'1'
　　Shapiro-Wilk normality test
data:X［［1L］］
W=0.8464,p-value=0.03828
$'2'
　　Shapiro-Wilk normality test
data:X［［2L］］ 
W = 0.8424, p-value = 0.04708
$'3'
　　Shapiro-Wilk normality test
data:X［［3L］］

W=0.9407,p-value=0.5068从计算结果来看,例2.9的第1组数据和第2组数据并不满足正态性要求.
2.方差的齐性检验
当式(2.23)成立时称为齐方差,关于齐方差的检验称为方差的齐性检验.在统计中,Bartlett检验是针对方差的齐性检验设计的,其原假设和备择假设为H0:σ21=σ22=…=σ2r，　H1:σ2i不全相同在R中,bartlett.test()函数完成Bartlett方差齐性检验,其使用格式有两种,一种是向量因子形式bartlett.test(x, g,...)参数x为数据构成的向量或列表.
g为因子构成的向量,当x为列表时,该项无效....为附加参数.
另一种是公式形式bartlett.test(formula, data, subset, na.action, ...)参数formula为形如lhs ~ rhs的公式,其中lhs表示数据,rhs表示数据对应的分组.
data为数据构成的数据框.subset为可选向量,表示选择样本的子集.
na.action为函数,表示处理缺失值(NA)的方法....为附加参数.
用bartlett.test()函数对例2.9的数据作Bartlett方差齐性检验.>bartlett.test(X ~ A, data = mouse)
　　Bartlett test of homogeneity of variances
data:　X by A
Bartlett's K-squared = 1.2068,df = 2, p-value = 0.5469数据满足方差齐性要求.
3.非齐性方差数据的方差分析
当数据只满足正态性,但不满足方差齐性的要求时,可用oneway.test()函数作方差分析,其使用格式为oneway.test(formula, data, subset, na.action,
var.equal = FALSE)参数formula为形如lhs ~ rhs的公式,其中lhs表示数据, rhs表示数据对应的分组.
data为数据构成的数据框.subset为可选向量,表示选择样本的子集.na.action为函数,表示处理缺失值(NA)的方法.
var.equal为逻辑变量,取TRUE表示在方差齐性的条件下计算;默认值为FALSE.
例如,用函数oneway.test()对例2.9的数据作单因素方差分析.>oneway.test(X`javascript
~ A, data = mouse)
　　One-way analysis of means (not assuming equal variances)
data:　X and A
F = 9.7869,num df = 2.000, denom df = 19.104,
p-value = 0.001185

oneway.test(X ~ A, data = mouse, var.equal = TRUE)
　　One-way analysis of means

data:　X and A

F = 8.4837,num df = 2, denom df = 30, p-value = 0.001202注意到, 在齐方差的假设下,其计算结果与例2.9的计算结果相同.
**2.2.4　秩检验**
1. Kruskal-Wallis秩和检验
如果需要分析的数据既不满足正态性要求,又不满足方差齐性要求,就不能用前面介绍的方法作方差分析,这就需要用到Kruskal-Wallis秩和检验.
在R中,kruskal.test()函数完成Kruskal-Wallis秩和检验,其使用格式有两种,一种是向量因子形式kruskal.test(x, g, ...)参数x为数据构成的向量或列表.
g为因子构成的向量,当x为列表时,该项无效....为附加参数.
另一种格式是公式形式kruskal.test(formula, data, subset, na.action, ...)参数formula为形如lhs ~ rhs的公式,其中lhs表示数据,rhs表示数据对应的分组.
data为数据构成的数据框.subset为可选向量,表示选择样本的子集.
na.action为函数,表示处理缺失值(NA)的方法....为附加参数.
例2.12　对例2.9中的数据作Kruskal-Wallis秩和检验.
解　虽然例2.9的数据满足方差齐性的要求,但并不满足正态性要求,严格地说,需要作秩检验.>kruskal.test(X ~ A, data = mouse)
　　Kruskal-Wallis rank sum test
data:　X by A 
Kruskal-Wallis chi-squared=12.0258, df=2, p-value=0.002447Kruskal-Wallis秩和检验的P值<0.05,说明三种不同菌型的伤寒杆菌相互之间有差异.
2.多重Wilcoxon秩和检验
如果Kruskal-Wallis秩和检验得到的结论是拒绝原假设,即各水平之间有显著差异,还需要进一步检验各水平之间谁与谁有差异.由于此时没有正态性和方差齐性的条件,所以需要作多重Wilcoxon秩和检验.
在R中,pairwise.wilcox.test()函数完成多重Wilcoxon秩和检验,其使用格式为pairwise.wilcox.test(x, g, 
　　p.adjust.method = p.adjust.methods,
　　paired = FALSE, ...)参数x为数据构成的向量.g为因子构成的向量或因子.
p.adjust.method为P值的调整方法,p.adjust.methods的取值分别为"holm"(默认值), "hochberg","hommel", "bonferroni","BH","BY","fdr"和"none"之一.
paired为逻辑变量,表示数据是否为成对数据....为附加参数.
例如,继续对例2.9中的数据作多重Wilcoxon秩和检验.>with(mouse, pairwis```javascript
e.wilcox.test(X, A))
　　Pairwise comparisons using Wilcoxon rank sum test
data:　X and A 
　1　　　2
2 0.0076 -
3 0.0086 0.6878

P value adjustment method: holm多重Wilcoxon秩和检验说明:第1种与第2种有显著差异,第1种与第3种有显著差异,而第2种与第3种无显著差异,这个结论与多重t检验是相同的.P值调整方法使用的是默认值,所以采用的是Holm修正.注意:由于数据打结(即有相同的秩),程序会给出警告.
2.2.5　双因素方差分析
所谓双因素,就是考虑两个因素——因素A和因素B,其中因素A有r个水平A1,A2,…,Ar,因素B有s个水平B1,B2,…,Bs.
1.不考虑交互效应
双因素方差分析分两种情况,一种是不考虑交互效应,每组条件下只取一个样本.假定xij~N(μij,σ2)(i=1,2,…,r,j=1,2,…s)且各xij相互独立,数据可以分解为
xij=μ+αi+βj+εij,　i=1,2,…,r,　j=1,2,…,s
εij~N(0,σ2)且相互独立(2.24)其中μ=1rs∑ri=1∑sj=1μij为总平均,αi为因素A的第i个水平的效应,βj为因素B的第j个水平的效应.
在线性模型(2.24)下,方差分析的主要任务是,系统分析因素A和因素B对试验指标影响的大小,因此,在给定显著性水平α下,提出如下统计假设:
对于因素A,“因素A对试验指标影响是否显著”等价于检验H01:α1=α2=…=αr=0,　H11:α1,α2,…,αr不全为零对于因素B,“因素B对试验指标影响是否显著”等价于检验H02:β1=β2=…=βs=0,　H12:β1,β2,…,βs不全为零双因素方差分析与单因素方差分析的统计原理基本相同,也是基于平方和分解公式ST=SE+SA+SB其中ST=∑ri=1∑sj=1(xij-x)2,　x=1rs∑ri=1∑sj=1xij
SA=s∑ri=1(xi·-x)2,　xi·=1s∑sj=1xij,　i=1,2,…,r
SB=r∑sj=1(x·j-x)2,　x·j=1r∑ri=1xij,　j=1,2,…,s
SE=∑ri=1∑sj=1(xij-xi·-x·j+x)2这里ST为总离差平方和,SE为误差平方和,SA是由因素A的不同水平所引起的离差平方和(称为因素A的平方和).类似地,SB称为因素B的平方和.可以证明当H01成立时SA/σ2~χ2(r－1)且与SE相互独立,而SE/σ2~χ2((r－1)(s－1))于是当H01成立时FA=SA/(r－1)SE/［(r－1)(s－1)］~F(r－1,(r－1)(s－1))类似地,当H02成立时,FB=SB/(s－1)SE/［(r－1)(s－1)］~F　(s－1,(r－1)(s－1))分别以FA,FB作为H01,H02的检验统计量,将计算结果列成方差分析表,如表2.8所示.
表2.8　双因素方差分析表方差来源平方和自由度均方F比P值因素ASAr－1MSA=SAr－1FA=MSAMSEpA因素BSBs－1MSB=SBs－1FB=MSBMSEpB误差SE(r－1)(s－1)MSE=SE(r－1)(s－1)总和STrs－1
双因素方差分析的计算与单因素相同,仍然用到aov()函数与summary()函数,或者是lm()函数与anova()函数.
例2.13　在一个农业试验中,考虑4种不同的种子品种A1,A2,A3和A4,三种不同的施肥方法B1,B2和B3,得到产量数据如表2.9所示.试分析种子与施肥对产量有无显著影响.
表2.9　农业试验数据（单位:kg）B1B2B3A1325292316A2317310318A3310320318A4330370365
解　这是一个双因素试验,因素A(种子)有4个水平,因素B(施肥)有3个水平.由于每组条件下只取一个样本,所以不考虑交互效应.程序(程序名:exam0213.R)如下：agriculture <- data.frame(
　 Y = c(325, 292, 316, 317, 310, 318,
　　　　 310, 320, 318, 330, 370, 365),
　 A = gl(4, 3),
　 B = gl(3, 1, 12)
)
agriculture.aov <- aov(Y ~ A + B, data = agriculture)
summary(agriculture.aov)在程序中,gl()生成因子,因素A为4水平,每个水平重复3次,因素B为3水平,每个水平只有1次,需要对应12个变量,所以长度为12.Y ~ A + B表示考虑双因素.计算得到　　　　　　Df Sum Sq Mean Sq F value Pr(>F)　
A　　　　　　3　 3824　1274.8　 5.226 0.0413 *
B　　　　　　2　　163　　81.3　 0.333 0.7291
Residuals　　6　1463　 243.9根据P值说明不同品种(因素A)对产量有显著影响,而没有充分理由说明施肥方法(因素B)对产量有显著的影响.
事实上,在应用模型(2.24)时,遵循一种假定,即因素A和因素B对指标的效应是可以叠加的,而且认为因素A的各水平效应的比较与因素B在什么水平无关.这里并没有考虑因素A和因素B的各种水平组合(Ai,Bj)的不同给产量带来的影响,而这种影响在许多实际工作中是应该给予足够的重视的,这种影响被称为交互效应.这就导出下面所要讨论的问题.
2.考虑交互效应
在考虑交互效应的情况下,每组条件下要取多个样本.假定xijk~N(μij,σ2),　i=1,2,…,r;　j=1,2,…,s;　k=1,2,…,t各xijk相互独立,所以数据可以分解为xijk=μ+αi+βj+δij+εijk,
εijk~N(0,σ2),且相互独立，
i=1,2,…,r,　j=1,2,…,s,　k=1,2,…,t(2.25)其中μ=1rs∑ri=1∑sj=1μij,αi为因素A的第i个水平的效应,βj为因素B的第j个水平的效应,δij表示Ai和Bj的交互效应.此时判断因素A、因素B,以及交互效应的影响是否显著等价于检验下列假设:H01:　α1=α2=…=αr=0,　H11:α1,α2,…,αr不全为零
H02:　β1=β2=…=βs=0,　H12:β1,β2,…,βs不全为零
H03:　δij=0,　H13:δij不全为零，　i=1,2,…,r,　j=1,2,…,s在这种情况下,方差分析法与前面类似,有下列计算公式:ST=SE+SA+SB+SA×B其中ST=∑ri=1∑sj=1∑tk=1(xijk－x)2,x=1rst∑ri=1∑sj=1∑tk=1xijk
SE=∑ri=1∑sj=1∑tk=1(xijk－xij·)2
xij·=1t∑tk=1xijk,　i=1,2,…,r,　j=1,2,…,s
SA=st∑ri=1(xi··－x)2,　xi··=1st∑Sj=1∑tk=1xijk,　i=1,2,…,r
SB=rt∑Sj=1(x·j·－x)2,　x·j·=1rt∑ri=1∑tk=1xijk,　j=1,2,…,s
SA×B=t∑ri=1∑Sj=1(xij·－xi··－x·j·+x)2
这里ST为总离差平方和,SE为误差平方和,SA为因素A的平方和,SB为因素B的平方和,SA×B为交互效应平方和.可以证明：当原假设成立时,FA=SA/(r－1)SE/［rs(t－1)］~F(r－1,rs(t－1))
FB=SB/(s－1)SE/［rs(t－1)］~F(s－1,rs(t－1))
FA×B=SA×B/［(r－1)(s－1)］SE/［rs(t－1)］~F((r－1)(s－1),rs(t－1))将FA,FB,FA×B作为检验统计量,构造方差分析表(见表2.10).

例2.14　研究树种与地理位置对松树生长的影响,对4个地区的三种同龄松树的直径进行测量得到数据如表2.11所示.A1,A2,A3表示三个不同树种,B1,B2,B3,B4表示4个不同地区.对每一种水平组合,进行了5次测量,对此试验结果进行方差分析.

解　这是一个双因素问题,并且每组条件下取5组数据,因此需要考虑两因素间的交互效应.程序(程序名:exam0214.R)如下：tree <- data.frame(

Y = c(23, 25, 21, 14, 15, 20, 17, 11, 26, 21,
　　　　　16, 19, 13, 16, 24, 20, 21, 18, 27, 24,
　　　　　28, 30, 19, 17, 22, 26, 24, 21, 25, 26,
　　　　　19, 18, 19, 20, 25, 26, 26, 28, 29, 23,
　　　　　18, 15, 23, 18, 10, 21, 25, 12, 12, 22,
　　　　　19, 23, 22, 14, 13, 22, 13, 12, 22, 19),
　　A = gl(3, 20, 60, labels = paste("A", 1:3, sep = "")),
　　B = gl(4, 5, 60, labels = paste("B", 1:4, sep = ""))
)
tree.aov <- aov(Y ~ A+B+A:B, data=tree)

summary(tree.aov)在程序中,paste()为粘连函数,表示因子为A1, A2, A3和B1, B2, B3, B4.A:B表示交互效应.使用Y ~ A*B与Y ~ A+B+A:B的意义是等价的.计算结果为　　　　　　Df Sum Sq Mean Sq F value　 Pr(>F)
A2352.5176.278.9590.000494 *
B387.529.171.4830.231077
A:B671.711.960.6080.722890
Residuals48944.419.68可见在显著性水平α=0.05下,树种(因素A)效应是高度显著的,而位置(因素B)效应和交互效应并不显著.
得到结果后,如何使用它呢?一种简单的方法是计算各因素的均值.由于树种(因素A)效应是高度显著的,也就是说,选什么树种对树的生长很重要,因此,要选那些生长粗壮的树种.计算因素A的均值> attach(tree); tapply(Y, A, mean)得到A1A2A3
19.5523.5517.75所以选择第2种树对生长有利.下面计算因素B(位置)的均值> tapply(Y, B, mean)得到B1B2B3B4
19.8666720.6000018.6666722.00000是否选择位置4最有利呢?不必了.因为计算结果表明:关于位置效应并不显著,也就是说,所受到的影响是随机的.因此,选择成本较低的位置种树就可以了.
本题关于交互效应更不显著,因此,没有必要计算交互效应的均值.如果需要计算其均值,可用命令matrix(tapply(Y, A:B, mean), nr = 3, nc = 4, byrow = T,
　　　 dimnames=list(levels(A), levels(B)))得到B1B2B3B4
A119.619.017.622.0
A223.224.420.226.4
A316.818.418.217.6如果交互效应是显著的,则可根据上述结果选择最优的方案.
3.交互效应图
除了用P值来判断两因素是否有交互效应外,还可以通过图形来判断.在R中,interaction.plot()函数就是为这种需求设计的,其使用格式为interaction.plot(x.fa`javascript
ctor, trace.factor,
　　response, fun = mean,
　　type = c("l", "p", "b", "o", "c"), legend = TRUE,
　　trace.label = deparse(substitute(trace.factor)),
　　fixed = FALSE,
　　xlab = deparse(substitute(x.factor)),
　　ylab = ylabel,
　　ylim = range(cells, na.rm = TRUE),
　　lty = nc:1, col = 1, pch = c(1:9, 0, letters),
　　xpd = NULL, leg.bg = par("bg"), leg.bty = "n",
　　xtick = FALSE, xaxt = par("xaxt"), axes = TRUE,

　　...)参数x.factor为画在X轴上的因子,trace.factor为另一个因子,它构成水平的迹.
response为影响变量.fun为计算函数,默认值为计算均值.type为绘图的类型,它与plot()函数中的意义相同.
legend为逻辑变量,表示是否包含图例,默认值为TRUE. 
trace.label为图例的总标记.
fixed为逻辑变量,表示图例是否按trace.factor的水平排序.
xlab, ylab, ylim, lty,col, pch等参数与plot()函数的意义相同,这里就不一一介绍了.
绘出例2.14的交互效应图,命令为> interaction.plot(A, B, Y, lwd = 2, col = 2:5)其图形如图2.5所示.如果图中折线接近平行,则说明交互效应不显著;如果折线相互交叉,则说明两因素的交互效应明显.从图2.5的结果可以得到,因素A和因素B的交互效应不显著,这与方差分析的结果是相同的.
图2.5　交互效应图