matlab多元线性回归分析（MATLAB多因素一元方差分析）

多因素一元方差分析

（1）多因素一元方差分析的MATLAB实现

MATLAB工具工具箱中提供了anovan函数，用来根据样本观测值向量进行均衡或非均衡试验的多因素一元方差分析，检验多个（N个）因素的主效应或交互式效应是否显著，调用格式如下：

<1>p=anovan(y,group)

根据样本观测值向量y进行均衡或非均衡试验的多因素一元方差分析，检验多个因素的主效应是否显著。输入参数group是一个元胞数组，它的每一个元胞对应一个因素，是该因素的水平列表，与y等长，用来标记y中的每个观测所对应的因素的水平。

anovan函数还生成1个图形，用来显示一个标准的多因素一元方差分析表。

<2>p=anovan(y,group,param1,val1,param2,val2,......)

通过指定一个或多个成对出现的参数名与参数值来控制多因素一元方差分析。可用的参数名与参数值如下表

参数名	参数值	说明
‘alpha’	（0,1）内的标量	指定置信水平
‘continuous’	下标向量	用来指明哪些分组变量被作为连续变量，而不是离散的分类变量
‘display’	‘on’‘off’	用来指定是否显示方差分析表
‘model’	‘linear’，‘interaction’‘full’	用来指定所用模型的类型，‘linear’（默认）只对N个主效应进行检验，不考虑交互效应，‘interaction’对主效应和两个因素的交互效应减小检验‘full’对N个主效应和全部的交互效应进行检验。也可以通过0和1的矩阵自定义效应项
‘nested’	由0和1构成的矩阵M	指定分组变量之间的嵌套关系
‘random’	下标向量	用来指明哪些分组变量是随机的。
‘sstype’	1,2或3	指定平方和的类型，默认值为3
‘varnames’	字符矩阵或字符串元胞数组	指定分组变量的名称。当没有指定时，默认用‘X1’，‘X2’，…’XN’做为名称

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

返回元胞数组形式的方差分析表table。

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

返回一个结构体变量stats，用于进行后续的多重比较，。当某因素对实验指标的影响显著时，在后续的分析中，可用调用multcompare函数，把stats作为其输入，进行多重比较。

处理区组1234N1P138293640N1P245423743N2P158465251N2P267706571N3P162646170N3P258637169

调用anovan函数进行分析，在不考虑区组因素的情况下，分析氮、磷两种肥料的施用量对水稻的产量是否有显著影响，并分析交互作用是否显著，然后找出在因素A，B的哪种水平组合下水稻的平均产量高，显著性水平为0.05.

%定义一个矩阵，输入原始数据

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按列拉长成24行1列的向量

yield=yield(:);

%定义因素A(氮)的水平列表向量

A=strcat({'N'},num2str([ones(8,1);2*ones(8,1);3*ones(8,1)]));

%定义因素B(磷)的水平列表向量

B=strcat({'P'},num2str([ones(4,1);2*ones(4,1)]));

B=[B;B;B];

%将因素A、B的水平列表向量与yield向量放在一起构成一个元胞数组，以元胞数组形式显示出来

[A,B,num2cell(yield)]

%指定因素名称，A表示氮肥施用量，B表示磷肥施用量

varnames={'A','B'};

%调用anovan函数作双因素一元方差分析，返回主效应A、B和交互效应AB所对应的p值向量

%还返回方差分析表table，结构体变量stats，标识模型效应项的矩阵term

[p,table,stats,term]=anovan(yield,{A,B},'model','full','varnames',varnames)

ans =

'N1' 'P1' [38]

'N1' 'P1' [29]

'N1' 'P1' [36]

'N1' 'P1' [40]

'N1' 'P2' [45]

'N1' 'P2' [42]

'N1' 'P2' [37]

'N1' 'P2' [43]

'N2' 'P1' [58]

'N2' 'P1' [46]

'N2' 'P1' [52]

'N2' 'P1' [51]

'N2' 'P2' [67]

'N2' 'P2' [70]

'N2' 'P2' [65]

'N2' 'P2' [71]

'N3' 'P1' [62]

'N3' 'P1' [64]

'N3' 'P1' [61]

'N3' 'P1' [70]

'N3' 'P2' [58]

'N3' 'P2' [63]

'N3' 'P2' [71]

'N3' 'P2' [69]

p =

0.0000

0.0004

0.0080

table =

'Source' 'Sum Sq.' 'd.f.' 'Singular?' 'Mean Sq.' 'F' 'Prob>F'

'A' [3.0670e 03] [ 2] [ 0] [1.5335e 03] [78.3064] [1.3145e-09]

'B' [ 368.1667] [ 1] [ 0] [ 368.1667] [18.8000] [3.9813e-04]

'A*B' [ 250.3333] [ 2] [ 0] [ 125.1667] [ 6.3915] [ 0.0080]

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

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

stats =

source: 'anovan'

resid: [24x1 double]

coeffs: [12x1 double]

Rtr: [6x6 double]

rowbasis: [6x12 double]

dfe: 18

mse: 19.5833

nullproject: [12x6 double]

terms: [3x2 double]

nlevels: [2x1 double]

continuous: [0 0]

vmeans: [2x1 double]

termcols: [4x1 double]

coeffnames: {12x1 cell}

vars: [12x2 double]

varnames: {2x1 cell}

grpnames: {2x1 cell}

vnested: []

ems: []

denom: []

dfdenom: []

msdenom: []

varest: []

varci: []

txtdenom: []

txtems: []

rtnames: []

term =

1 0

0 1

1 1

返回的向量p是主效应A，B和交互效应AB所对应的检验的p值，table是元胞数组的方差分析表，stats是一个结构体变量，可用于后续的分析中（如多重比较），矩阵term的3行分布表示了3个效应项：主效应项A，主效应项B和交互式效应项AB，还生成了一个方差分析表。

返回的结果与调用anova2函数得到的结果一样，因素A，因素B以及他们的交互作用所对应的p值均小于给定的显著想水平0.05，所以可以认为氮、磷两种肥料的施用量对水稻的产量均有非常显著的影响，并且他们之间的交互作用也是非常显著的。

（3）多重比较。

调用multcompare函数，把anovan函数返回的结构体变量stats作为它的输入，对因素A、B的每种水平的组合进行多重比较，找出在因素A、B的哪种水平组合下水稻的平均产量最高。

%多重比较

[c,m,h,gnames]=multcompare(stats,'dimension',[1 2])

%查看各处理的均值

[gnames,num2cell(m)]

c =

1.0000 2.0000 -25.9446 -16.0000 -6.0554 0.0009

1.0000 3.0000 -38.4446 -28.5000 -18.5554 0.0000

1.0000 4.0000 -15.9446 -6.0000 3.9446 0.4236

1.0000 5.0000 -42.4446 -32.5000 -22.5554 0.0000

1.0000 6.0000 -39.4446 -29.5000 -19.5554 0.0000

2.0000 3.0000 -22.4446 -12.5000 -2.5554 0.0093

2.0000 4.0000 0.0554 10.0000 19.9446 0.0483

2.0000 5.0000 -26.4446 -16.5000 -6.5554 0.0006

2.0000 6.0000 -23.4446 -13.5000 -3.5554 0.0047

3.0000 4.0000 12.5554 22.5000 32.4446 0.0000

3.0000 5.0000 -13.9446 -4.0000 5.9446 0.7926

3.0000 6.0000 -10.9446 -1.0000 8.9446 0.9995

4.0000 5.0000 -36.4446 -26.5000 -16.5554 0.0000

4.0000 6.0000 -33.4446 -23.5000 -13.5554 0.0000

5.0000 6.0000 -6.9446 3.0000 12.9446 0.9251

m =

35.7500 2.2127

51.7500 2.2127

64.2500 2.2127

41.7500 2.2127

68.2500 2.2127

65.2500 2.2127

h =

2

gnames =

'A=N1,B=P1'

'A=N2,B=P1'

'A=N3,B=P1'

'A=N1,B=P2'

'A=N2,B=P2'

'A=N3,B=P2'

ans =

'A=N1,B=P1' [35.7500] [2.2127]

'A=N2,B=P1' [51.7500] [2.2127]

'A=N3,B=P1' [64.2500] [2.2127]

'A=N1,B=P2' [41.7500] [2.2127]

'A=N2,B=P2' [68.2500] [2.2127]

'A=N3,B=P2' [65.2500] [2.2127]

返回的矩阵c是6种处理（N1P1，N2P1，N3P1，N1P2，N2P2，N3P2）间多重比较的结果矩阵，每一行的前两列是进行比较的两个处理的编号，第4列是两个处理的均值之差，第3列是两个处理均值差的95%置信下限，第5列是两个处理均值差的95%置信下限，当两个处理均值差的95%置信区间不包含0时，说明在显著性水平0.05下，这两个处理均值间差异是显著的。m矩阵列出了6种处理的平均值，很明显第5个处理（即N2P2）的平均值最大，由矩阵c或交互式多重比较的图中可以得到，处理5与处理3,6差异不显著，所以可以认为第3个和第6个处理也是可以的，所以，综上，可以在处理3,5,6中做出选择，即N3P1，N2P2，N3P2。

,