用 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 单因素方差分析表

实际计算时常采用方差分析表，如表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 不同伤寒杆菌对小白鼠存活天数影响的方差分析表

可以认为不同伤寒杆菌对小白鼠存活天数有极显著影响。

上述最后一行命令执行的结果为p=0.0038，实际判定时，此值更能说明显著性程度。

当各组样本容量相等时，Matlab自带了单因素方差分析函数anova1，并且自动返回类似的方差分析表,调用格式为anova1(x),这里x为矩阵，按列分组，第一列为第一组数据，要求每组数据相同。

例5.2 某班学生共分别住在四个宿舍，某次英语水平考试成绩如下：

表5-3 各宿舍学生英语成绩表

问不同宿舍学生英语水平有无显著差异？

解 利用复制粘贴的办法输入矩阵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组有极显著的差异。