雙因素一元方差分析

（1.）雙因素一元方差分析的MATLAB實現

MATLAB統計工具箱中提供了anova2函數，用來做雙因素一元方差分析，其調用格式如下：

<1>p=anova2(X,reps)

根據樣本觀測值均值X進行均衡實驗的雙因素一元方差分析。X的每一列對應因素A的一個水平，每行對應因素B的一個水平，X還應滿足方差分析的基本假定。reps表示因素A和B下的每一個水平組合下重複實驗的次數。

anova2函數檢驗矩陣X的各列是否具有相同的均值，即檢驗因素A對實驗指标的影響是否顯著，原假設為：

H0A：X的各列具有相同的均值（或因素A對實驗指标的影響不顯著）

anova2函數還檢驗矩陣X的各行是否具有相同的均值，即檢驗因素B對實驗指标的影響是否顯著，原假設為：

H0B：X的各行有相同的均值（或因素B對實驗指标的影響不顯著）

若參數reps的取值大于1（默認值為1），anova2函數還檢驗因素A和因素B的交互作用是否顯著，原假設為：

H0AB: A和B的交互作用不顯著

anova2函數返回檢驗的p值，若參數reps的取值等于1，則p是一個包含2個元素的行向量；若參數是reps的取值大于1，則p是一個包含3個元素的行向量，其元素粉筆是與H0A，H0B，H0AB對應的檢驗的p值。當檢驗的p值小于或等于給定的顯著性水平時，應拒絕原假設。

anova2函數還生成1個圖形，用來顯示一個标準的雙因素一元方差分析表。方差分析表把數據之間的差異分為三部分（當reps=1時）或四部分（當reps=2時）：

<2>p=anova2(X,reps,displayopt)

通過displayopt參數指定是否顯示帶有标準雙因素一元方差分析表的圖形窗口，當displayopt參數設置為‘on’（默認情況）時，顯示方差分析表；當displayopt參數設定為‘off’時，不顯示方差分析表。

<3>[p,table]=anova(.....)

返回元胞數組形式的方差分析表table（包含列标簽和行标簽）。

<4>[p,table,stats]=anova2(......)

返回一個結構體變量stats，用于進行後續的多重比較。

（2）例：為了研究肥料使用量對水稻産量的影響，某研究所做了氮（因素A）、磷（因素B）兩種肥料施用量的二因素試驗。氮肥用量設低、中、高三個水平，分布使用N1，N2和N3表示；磷肥用量設低、高2個水平，分别用P1，P2表示。供3x2=6個處理，重複4次，随機區組設計，測得水稻産量如下表

根據上表中的數據，不考慮區組因素，分析氮、磷兩種肥料的施用量對水稻産量是否有顯著性影響，并分析交互作用是否顯著。取顯著性水平=0.05；

注意：這裡不需要進行正态性和方差性齊次性檢驗，因素數據少，在數據比較少的情況下正态檢驗的結果是不可靠的，即使不滿足方差分析的假定，方差分析的結果通常也是比較穩定的。

雙因素一元方差分析首先要把數據矩陣處理一下，要把矩陣裝換成每一列對應因素A的一個水平，每行對應因素B的一個水平。本例中，每一列對應一個A因素（氮）水平，每一行對應一個B因素（磷）水平；反過來也可以。

先轉置為

把第2列，第4列，第6列接到第1列，第3列，第5列下面

提出公共項

%定義一個矩陣，輸入原始數據

yield=[38 29 36 40

45 42 37 43

58 46 52 51

67 70 65 71

62 64 61 70

58 63 71 69];

yield=yield'; %矩陣轉置

%将數據矩陣yield轉換成8行3列的矩陣，列對應因素A（氮），行對應因素B（磷）

yield=[yield(:,[1,3,5]);yield(:,[2,4,6])];

%定義元胞數組，以元胞數組形式顯示轉換後的數據

top={'因素','N1','N2','N3'};

left={'P1';'P1';'P1';'P1';'P2';'P2';'P2';'P2'};

%顯示數據

[top;left,num2cell(yield)]

%調用anova2函數作雙因素方差分析，返回檢驗的p值向量，方差分析表，結構體标量stats

[p,table,stats]=anova2(yield,4)

ans =

'因素' 'N1' 'N2' 'N3'

'P1' [38] [58] [62]

'P1' [29] [46] [64]

'P1' [36] [52] [61]

'P1' [40] [51] [70]

'P2' [45] [67] [58]

'P2' [42] [70] [63]

'P2' [37] [65] [71]

'P2' [43] [71] [69]

p =

0.0000 0.0004 0.0080

table =

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

'Columns' [ 3067] [ 2] [1.5335e 03] [78.3064] [1.3145e-09]

'Rows' [368.1667] [ 1] [ 368.1667] [18.8000] [3.9813e-04]

'Interaction' [250.3333] [ 2] [ 125.1667] [ 6.3915] [ 0.0080]

'Error' [352.5000] [18] [ 19.5833] [] []

'Total' [ 4038] [23] [] [] []

stats =

source: 'anova2'

sigmasq: 19.5833

colmeans: [38.7500 60 64.7500]

coln: 8

rowmeans: [50.5833 58.4167]

rown: 12

inter: 1

pval: 0.0080

df: 18

因素A、因素B以及他們的交互作用對應的檢驗p值均小于給定的顯著性水平0.05，所以可以認為氮、磷兩種肥料的施用量對水稻的産量均有顯著性影響，并且他們之間的交互作用也是非常顯著的。由于氮、磷兩種肥料的用量對水稻的産量均有非常顯著的影響，可以作進一步分析，例如進行多重分析，找出因素A、B在哪種水平的組合下水稻的平均産量最高。

（3）多重比較

下面調用multcompare函數，把anova2函數返回的結構體變量stats作為它的輸入，進行多重比較。

%對列（因素A）進行多重比較

[c_A,m_A]=multcompare(stats,'estimate','column')

%對行（因素B）進行多重比較

[c_B,m_B]=multcompare(stats,'estimate','row')

Note: Your model includes an interaction term that is significant at the level

you specified. Testing main effects under these conditions is questionable.

c_A =

1.0000 2.0000 -26.8971 -21.2500 -15.6029 0.0000

1.0000 3.0000 -31.6471 -26.0000 -20.3529 0.0000

2.0000 3.0000 -10.3971 -4.7500 0.8971 0.1084

m_A =

38.7500 1.5646

60.0000 1.5646

64.7500 1.5646

Note: Your model includes an interaction term that is significant at the level

you specified. Testing main effects under these conditions is questionable.

c_B =

1.0000 2.0000 -11.6289 -7.8333 -4.0378 0.0004

m_B =

50.5833 1.2775

58.4167 1.2775

由上面結果可以看出，若單獨考慮A因素，它的第1個水平與後兩個水平差異顯著，它的第2個水平與第3個水平差異不顯著，并且當A去第3個水平（N3）時，水稻産量的均值達到最大（64.75）；如果單獨考慮B因素，它的兩個水平差異顯著，在因素A，B的水平組合N3P2下，水稻的平均産量達到最大值，然而這确實錯誤的，因為A,B之間存在着非常顯著的交互作用，在這種情況下對主效應進行檢驗可能存在問題，這時應該對因素A、B的每種水平組合進行多重比較，找出所要的水平組合。這就是下一節的内容。