數理統計與Matlab: 第5章 方差分析

5.1 單因素方差分析

5.1.1 方差分析的基本概念

在實際問題中，人們常常需要在不同的條件下對所研究的對象進行對比試驗，從而得到若干組數據（樣本）。方差分析就是一種分析、處理多組實驗數據間均值差異的顯著性的統計方法。其主要任務是，通過對數據的分析處理，搞清楚各實驗條件對實驗結果的影響，以便更有效地指導實踐，提高經濟效益或者科研水平。

在統計中，人們稱受控制的條件為因素，因素所處的狀態稱為水平。

如果只讓一個因素變動，取該因素的多個不同水平進行試驗，而其他因素保持不變，稱該試驗為單因素試驗。例如小麥種植產量，只考慮"品種"這一因素，研究4個不同品種產量的差異，其它諸如施肥方案、灌溉方案等因素保持一致，就是一個4水平單因素試驗。

如果同時考慮兩個因素，例如4個小麥品種在3種不同施肥方案下的產量，就是一個雙因素試驗。

對於組實驗數據，我們假定都來自正態總體，並且具有相同的方差（稱為方差齊性），要檢驗這相互獨立的個正態總體

均值間有無差異，即：

H₀：； H₁：諸不全相同

前面我們講過兩正態總體均值的假設檢驗，有T檢驗的方法。自然有一個想法，對於，分別檢驗是否成立，若所有搭配均不拒絕，則接受H₀，只要有一種搭配拒絕原假設認為，那就拒絕H₀，看起來也不算麻煩。不妨稱上述想法為"兩兩T檢驗法"。

回憶前面內容，設為來自正態總體的簡單隨機樣本，為來自正態總體的簡單隨機樣本，且兩樣本獨立。為比較兩個總體的期望，提出如下原假設：

H₀：

當H₀成立時，檢驗統計量

我們給出函數t2test.m，解決上述計算問題。

function T=t2test(x,y)

m=length(x);

n=length(y);

vx=var(x);

vy=var(y);

a=(mean(x)-mean(y));

b=m+n-2;

c=(m-1)*vx+(n-1)*vy;

d=sqrt(m*n/(m+n));

T=a*d*sqrt(b/c);

以下給出m=10,n=10,且兩總體皆服從標准正態分布的情形下，萬次模擬的拒絕頻率。以下命令文件保存為PnT2.m

N=10000;

m=10; n=10;

alpha=0.05;

t0=tinv(1-alpha/2,m+n-2);

P=0;

for k=1:N

x=randn(1,10); y=randn(1,10);

T=t2test(x,y);

if abs(T)>t0

P=P+1;

end

end

P=P/N

執行上述程序，發現每次頻率都在0.05附近，說明上述兩個正態總體均值的T檢驗的確是水平為的檢驗。

我們設想有8組數據，客觀上都是來自標准正態分布，沒有差異，每組樣本容量都是10。現在用前述"兩兩T檢驗法"進行檢驗，下述程序計算出了萬次模擬中拒絕的頻率。

N=10000;

n=10; r=8;

alpha=0.05;

t0=tinv(1-alpha/2,n+n-2);

P=0;

for k=1:N

x=randn(8,10);

E=mean(x,2);

[EE,I]=sort(E);

X=x(I,:);

T=t2test(X(1,:),X(8,:));

if abs(T)>t0

P=P+1;

end

end

P=P/N

上述程序模擬發現，拒絕頻率大約在0.45左右，嚴重偏離0.05，說明依照"兩兩T檢驗"犯第一類錯誤的概率嚴重增大，判定結果很不可靠。

對於8組數據，兩兩比較共種組合，若每種組合接受原假設的概率為0.95，則28種組合都接受原假設的概率大致估計為，拒絕概率大致估計為0.76。由於相關性，拒絕概率沒有達到0.76，但0.45也相當大了。

為了避免上述問題的出現，1923年，波蘭數學家R.A.Fisher提出了方差分析(Analysis of Variance簡稱ANOVA) 法，可以同時判定多組數據均值間差異的顯著性檢驗問題。其檢驗統計量在H₀成立時服從F分布，這里F分布就是以Fisher姓氏的第一個字母命名的。

5.1.2 單因素方差分析的計算

設有組數據，表示因素A的個水平，每組有個觀測值。我們已知實際結果具有以下結構：

(; )

表示水平Ai下的理論均值，為實驗誤差，諸相互獨立且服從正態分布。

為了看出因素A個水平影響的大小，將進行分解，令

，

表示水平Ai對試驗結果的影響，稱為Ai的水平效應。顯然

這時數據有如下結構：

(; ) (5-1)

於是，我們需要進行的假設檢驗為：

H₀：； H₁：諸不全為零 (5-2)

 

記

，() ，

(5-3)

稱為總離差平方和，它反映了樣本觀測值之間的總的變異程度。以下我們將分解為兩部分，以便區別水平效應與隨機誤差的影響。

其中

記

， (5-4)

稱為組內平方和，它反映了每組的組內隨機誤差。稱為組間平方和，反應的是組與組之間的差異。上述推導說明，總離差平方可以分解為

(5-5)

一個自然的想法是：如果在總離差平方和中，所占比例很大，則拒絕原假設，認為客觀上存在水平效應。

在H₀成立時容易計算

； (5-6)

因此，當H₀成立時，有

， (5-7)

(5-8)

對於自由度，求臨界值，當時拒絕H₀即可。

表5-1 單因素方差分析表

方差來源	平方和	自由度	均方	F值	臨界值	顯著性
組間	SA
誤差	SE
總和	ST

實際計算時常采用方差分析表，如表5-1所示。當時，稱為不顯著，即認為各組均值之間沒有顯著差異，在顯著性一欄不做任何標記。當時，稱為較顯著，即認為各組均值之間有較顯著差異，在顯著性一欄用(*)標記。當時，稱為顯著，即認為各組均值之間有顯著差異，在顯著性一欄用*標記。當時，稱為極顯著，即認為各組均值之間有極顯著差異，在顯著性一欄用**標記。

上述傳統的方差計算表，在計算機普及后稍有變動，表中最后兩列可以變動為直接計算H₀成立時F分布大於此F值的概率，是否顯著一看自明。

例5.1 為了研究三種不同傷寒桿菌對於小白鼠存活天數的影響，分三組實驗，實驗數據如下：

A1: 2, 4, 3, 2, 4, 7, 7, 2, 5, 4

A2: 5, 6, 8, 5, 10,7, 12,6, 6

A3: 7, 11,6, 6, 7, 9, 5, 10,6, 3,10

試檢驗不同傷寒桿菌對於小白鼠存活天數有無顯著影響？

解 原假設H₀沒有顯著差異。以下利用Matlab計算，並將計算結果匯總到表5-2中。

x1=[2, 4, 3, 2, 4, 7, 7, 2, 5, 4];

x2=[5, 6, 8, 5, 10,7, 12,6, 6];

x3=[7, 11,6, 6, 7, 9, 5, 10,6, 3,10];

n1=length(x1);

n2=length(x2);

n3=length(x3);

X1=sum(x1), mx1=mean(x1),

X2=sum(x2), mx2=mean(x2),

X3=sum(x3), mx3=mean(x3),

n=n1+n2+n3, X=X1+X2+X3, mx=X/n,

SE=(x1-mx1)*(x1-mx1)'+(x2-mx2)*(x2-mx2)'+(x3-mx3)*(x3-mx3)',

SA=n1*(mx1-mx)^2+n2*(mx2-mx)^2+n3*(mx3-mx)^2,

F=(SA/2)/(SE/27),

finv(0.9,2,27),

finv(0.95,2,27),

finv(0.99,2,27),

p=1-fcdf(F,2,27),

表5-2 不同傷寒桿菌對小白鼠存活天數影響的方差分析表

方差來源	平方和	自由度	均方	F值	臨界值	顯著性
組間	70.4293	2	35.2146	6.9030	2.5106	**
誤差	137.7374	27	5.1014		3.3541
總和	208.1667	29			5.4881

可以認為不同傷寒桿菌對小白鼠存活天數有極顯著影響。

上述最后一行命令執行的結果為p=0.0038，實際判定時，此值更能說明顯著性程度。

當各組樣本容量相等時，Matlab自帶了單因素方差分析函數anova1，並且自動返回類似的方差分析表,調用格式為anova1(x),這里x為矩陣，按列分組，第一列為第一組數據，要求每組數據相同。

例5.2 某班學生共分別住在四個宿舍，某次英語水平考試成績如下：

表5-3 各宿舍學生英語成績表

宿舍一	86	74	78	76	73	82
宿舍二	59	76	63	66	75	65
宿舍三	79	71	82	66	87	96
宿舍四	86	91	95	77	88	85

問不同宿舍學生英語水平有無顯著差異？

解 利用復制粘貼的辦法輸入矩陣x，由於anova1要求按列分組，故使用x=x'命令，然后再執行命令

p=anova1(x),

返回值為p =0.002，故差異極顯著。Matlab同時返回了兩個圖形。

圖5-1 Matlab中anova1返回的方差分析表

 

圖5-2 Matlab中anova1返回的各組數據圖

圖5-2中各線含義依次為：最小值、1/4分位數、中位數、3/4分位數、最大值。

5.1.3 單因素方差分析的多重比較

經過方差分析之后，如果拒絕原假設，認為各組之間的均值有顯著差異，那么，這個判斷是對整體而言的，並不是說每兩個不同的組之間均值都存在顯著差異。那么，如何確定哪兩個組之間有顯著差異、無顯著差異呢？這就要對每種搭配做一對一的比較，即多重比較。Matlab提供了multcompare函數用於進行多重比較，為了使用這個函數，在用anova1做方差分析的時候要使用如下的三個輸出的調用格式

[P,ANOVATAB,STATS] = anova1(x)

在例5.2的數據輸入之后，假定x已經經過轉置使得每組數據按列排列，執行上述命令后，則除了返回上述兩個圖形外，還返回

P =

0.0020

ANOVATAB =

'Source' 'SS' 'df' 'MS' 'F' 'Prob>F'

'Columns' [1.1963e+003] [ 3] [398.7778] [7.0768] [0.0020]

'Error' [1.1270e+003] [20] [ 56.3500] [] []

'Total' [2.3233e+003] [23] [] [] []

STATS =

gnames: [4x1 char]

n: [6 6 6 6]

source: 'anova1'

means: [78.1667 67.3333 80.1667 87]

df: 20

s: 7.5067

其中最后一個輸出結果STATS用於多重比較函數調用。繼續執行命令

COMPARISON = multcompare(STATS,'alpha',0.1)

得到輸出結果

 

COMPARISON =

1.0000 2.0000 0.2251 10.8333 21.4415

1.0000 3.0000 -12.6082 -2.0000 8.6082

1.0000 4.0000 -19.4415 -8.8333 1.7749

2.0000 3.0000 -23.4415 -12.8333 -2.2251

2.0000 4.0000 -30.2749 -19.6667 -9.0585

3.0000 4.0000 -17.4415 -6.8333 3.7749

上述結果中，逐行顯示一對一比較的結果。第一列與第二列是參與比較的組號，第三列與第五列是均值差的置信區間，置信水平由輸入變量alpha=0.1確定。第四列為兩組樣本均值的差。若第三列置信下限與第五列置信上限正負符號相同，則差異顯著。由於alpha=0.1，我們可以發現第1，2組差異較顯著，第2，3組差異較顯著，第2，4組差異較顯著，其它對比差異不顯著。

再修改上述輸入變量alpha的值，重新計算

COMPARISON = multcompare(STATS,'alpha',0.05)

輸出結果為

COMPARISON =

1.0000 2.0000 -1.2972 10.8333 22.9639

1.0000 3.0000 -14.1305 -2.0000 10.1305

1.0000 4.0000 -20.9639 -8.8333 3.2972

2.0000 3.0000 -24.9639 -12.8333 -0.7028

2.0000 4.0000 -31.7972 -19.6667 -7.5361

3.0000 4.0000 -18.9639 -6.8333 5.2972

只有第2，3組，第2，4組有顯著差異。繼續執行

COMPARISON = multcompare(STATS,'alpha',0.01)

輸出結果為

COMPARISON =

1.0000 2.0000 -4.5448 10.8333 26.2115

1.0000 3.0000 -17.3781 -2.0000 13.3781

1.0000 4.0000 -24.2115 -8.8333 6.5448

2.0000 3.0000 -28.2115 -12.8333 2.5448

2.0000 4.0000 -35.0448 -19.6667 -4.2885

3.0000 4.0000 -22.2115 -6.8333 8.5448

發現只有第2，4組有極顯著的差異。

5.2 雙因素方差分析

5.2.1 有重復實驗的雙因素方差分析

如果實驗同時受到兩個因素的共同影響，因素有個水平，因素有個水平，一共有種搭配方案，每種搭配有做個實驗數據，實驗數據表依照表5-4所示排列。

首先引進記號

在的三個下標中，若有的下標為點，則表示對該下標求和，例如

， ，

相應的平均值記為

， ，

依此類推。

將上述諸看成是隨機變量，並且假設滿足以下線性模型

， （5-9）

，，

表5-4 雙因素有重復試驗數據表

		因素
因 素

其中是各種搭配下的總平均數的理論值，是因素的第個水平的效應，是因素的第個水平的效應，是與的交互作用，即搭配效應。是相互獨立的正態隨機變量，均值為零，方差為，並且諸相互獨立。滿足

， （5-10）

對於上述有重復試驗的雙因素實驗，方差分析檢驗如下三個假設：

H₀₁： （5-11）

H₀₂： （5-12）

H₀₃：， （5-13）

類似單因素方差分析，記總離差平方和為

（5-14）

則有如下分解式

自由度為 （5-15）

其中

， 自由度為 （5-16）

， 自由度為 （5-17）

， 自由度為 （5-18）

， 自由度為 （5-19）

由此可得檢驗統計量，當H₀₁成立時，

（5-20）

時，拒絕H₀₁，認為因素作用顯著。

當H₀₂成立時，

（5-21）

時，拒絕H₀₂，認為因素作用顯著。

當H₀₃成立時，

（5-22）

時，拒絕H₀₃，認為交互作用顯著。

方差計算可以類似地用方差計算表進行。

Matlab自帶的函數anova2用於處理雙因素方差分析，調用格式為

[P,TABLE,STATS] = anova2(x,n)

其中P返回概率值，TABLE返回方差計算表，STATS返回的信息用於多重比較。輸入變量n表示每種搭配下樣本容量，記號同前述公式。x為數據矩陣，為na行b列，格式與表5-4完全一致。

例5.3 楊樹一年中生長高度受兩種因素影響，A：施肥方案，B：深翻方案。對於4種施肥方案及3種深翻方案，共12種搭配，各實驗3株，實驗結果如表5-5所示。試問施肥方案、深翻方案、兩者的交互作用對於苗高有無顯著影響？

表5-5 楊樹苗增高實驗數據表

	B1	B2	B3
A1	52 43 39	41 47 53	49 38 42
A2	48 37 29	50 41 30	36 48 47
A3	34 42 38	36 39 44	37 40 32
A4	45 58 42	44 46 60	43 56 41

 

解 利用復制粘貼的辦法對矩陣x賦值，注意到anova2的要求，重復數據按列排列，故對x進行轉置，x=x'。轉置之后x為9行4列矩陣，每列表示4種施肥方案的數據，行對應的是3種深翻方案。執行Matlab命令：

[P,TABLE,STATS] = anova2(x,3)

計算結果為：

P =

0.0259 0.7504 0.9540

TABLE =

'Source' 'SS' 'df' 'MS' 'F' 'Prob>F'

'Columns' [ 562.0833] [ 3] [187.3611] [3.6838] [0.0259]

'Rows' [ 29.5556] [ 2] [ 14.7778] [0.2906] [0.7504]

'Interaction' [ 76.6667] [ 6] [ 12.7778] [0.2512] [0.9540]

'Error' [1.2207e+003] [24] [ 50.8611] [] []

'Total' [1.8890e+003] [35] [] [] []

STATS =

source: 'anova2'

sigmasq: 50.8611

colmeans: [44.8889 40.6667 38 48.3333]

coln: 9

rowmeans: [42.2500 44.2500 42.4167]

rown: 12

inter: 1

pval: 0.9540

df: 24

同時返回圖形，圖形中顯示了方差分析表。諸列（Columns）表示的是4種施肥方案，諸行（Rows）表示的是3種深翻方案，從上述方差分析表中可以看出，施肥方案有顯著影響，深翻方案無顯著影響，兩因素間無交互作用。

COMPARISON = multcompare(STATS,'alpha',0.05)

結果顯示

Note: Your model includes an interaction term. A test of main

effects can be difficult to interpret when the model includes

interactions.

COMPARISON =

1.0000 2.0000 -5.0520 4.2222 13.4964

1.0000 3.0000 -2.3853 6.8889 16.1631

1.0000 4.0000 -12.7187 -3.4444 5.8298

2.0000 3.0000 -6.6075 2.6667 11.9409

2.0000 4.0000 -16.9409 -7.6667 1.6075

3.0000 4.0000 -19.6075 -10.3333 -1.0591

多重比較的結果發現，第3，4種施肥方案差異顯著。Matlab返回的注釋說明，對於有交互項的情形，上述多重比較僅供參考。

5.2.2 無重復實驗的雙因素方差分析

要把交互作用與隨機誤差區別開，就必須對每種搭配進行重復試驗。如果兩個因素間確無交互作用，線性模型可以簡化為：

， （5-23）

此時檢驗H₀₁與H₀₂即可。

利用anova2仍可進行方差分析，其原理仍是平方和分解與檢驗。不再一一羅列公式。

例5.4 某養豬場進行豬增重試驗，選擇4個品種的豬和3種飼料，共12中搭配方案，每種飼養一頭，三個月后增重數據如表5-6所示。試研究品種與飼料對於豬增重的影響。

解 復制粘貼輸入數據矩陣x，執行

[P,TABLE,STATS] = anova2(x)

表5-6 豬增重數據表

	飼料1	飼料2	飼料3
品種1	51	53	52
品種2	56	57	58
品種3	45	49	47
品種4	42	44	43

計算結果為：

P =

0.0156 0.0000

TABLE =

'Source' 'SS' 'df' 'MS' 'F' 'Prob>F'

'Columns' [ 10.5000] [ 2] [ 5.2500] [ 9] [ 0.0156]

'Rows' [332.2500] [ 3] [110.7500] [189.8571] [2.4683e-006]

'Error' [ 3.5000] [ 6] [ 0.5833] [] []

'Total' [346.2500] [11] [] [] []

STATS =

source: 'anova2'

sigmasq: 0.5833

colmeans: [48.5000 50.7500 50]

coln: 4

rowmeans: [52 57 47 43]

rown: 3

inter: 0

pval: NaN

df: 6

可以看出，諸列（Columns）間差異顯著，說明飼料作用顯著；諸行（Rows）間差異極顯著，說明在這4個豬的品種間，增重差異極顯著。

對於無重復實驗可以可靠地進行多重比較。繼續計算：

COMPARISON = multcompare(STATS,'alpha',0.05, 'estimate','column')

結果為

COMPARISON =

1.0000 2.0000 -3.9071 -2.2500 -0.5929

1.0000 3.0000 -3.1571 -1.5000 0.1571

2.0000 3.0000 -0.9071 0.7500 2.4071

發現第1，2種飼料差異顯著。

COMPARISON = multcompare(STATS,'alpha',0.05, 'estimate','row')

結果為

COMPARISON =

1.0000 2.0000 -7.1588 -5.0000 -2.8412

1.0000 3.0000 2.8412 5.0000 7.1588

1.0000 4.0000 6.8412 9.0000 11.1588

2.0000 3.0000 7.8412 10.0000 12.1588

2.0000 4.0000 11.8412 14.0000 16.1588

3.0000 4.0000 1.8412 4.0000 6.1588

4個豬的品種，每種搭配對比，增重差異都顯著。

COMPARISON = multcompare(STATS,'alpha',0.01, 'estimate','row')

結果為

COMPARISON =

1.0000 2.0000 -8.1014 -5.0000 -1.8986

1.0000 3.0000 1.8986 5.0000 8.1014

1.0000 4.0000 5.8986 9.0000 12.1014

2.0000 3.0000 6.8986 10.0000 13.1014

2.0000 4.0000 10.8986 14.0000 17.1014

3.0000 4.0000 0.8986 4.0000 7.1014

4個豬的品種，每種搭配對比，增重差異都極顯著。