《數學模組化：基于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,