matlab正态分布样本方差（MATLAB双因素一元方差分析）

X的每一列对应因素A的一个水平，每行对应因素B的一个水平，X还应满足方差分析的基本假定。reps表示因素A和B下的每一个水平组合下重复实验的次数。当检验的p值小于或等于给定的显著性水平时，应拒绝原假设。anova2函数还生成1个图形，用来显示一个标准的双因素一元方差分析表。氮肥用量设低、中、高三个水平，分布使用N1，N2和N3表示；磷肥用量设低、高2个水平，分别用P1，P2表示。

双因素一元方差分析

（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的每种水平组合进行多重比较，找出所要的水平组合。这就是下一节的内容。