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

本节书摘来自华章计算机《数学建模：基于r》一书中的第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,<code>`</code>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"),

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方差齐性检验.&gt;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的数据作单因素方差分析.&gt;oneway.test(x<code>`</code>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

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 &lt;- 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 &lt;- 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(&gt;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 &lt;- data.frame(

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(&gt;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的均值&gt; attach(tree); tapply(y, a, mean)得到a1a2a3

19.5523.5517.75所以选择第2种树对生长有利.下面计算因素b(位置)的均值&gt; 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<code>`</code>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,