第3章方差分析课程

3.1 单因素方差分析的数学模型

首先让我们看两个例子：

例3.1 设甲、乙、丙三块麦田的基本苗数（按面积大小抽取样本点数）得表3.1，问三块地的基本苗数是否有显著差别？

表3.1三块麦田的基本苗数 甲 21 29 24 22 25 30 27 26 . .

例3.2为了研究淬火温度和等温温度对铣刀硬度的影响三

种不同淬火温度和三种等温温度淬火，测得铣刀平均硬度如表3.2，检验淬火温度及等温温度是否对硬度y有显著影响。

表3.2 淬火温度与等温温度对硬度的影响 淬火温度 等温温度 乙 20 25 25 23 29 31 24 26 20 21 丙 24 22 28 25 21 26 . . . . B1 66 65 B2 66 68 67 B3 68 67 68 A1 A2 A3 这两个例子和以前回归分析的问题不同．首先，它们都只考察某种因素(地块、淬火温度、等温温度)在一系列试验中对产品某个指标(寿命、得率)的影响是否显著，而不要求建立回归方程；其次，这些因素可以不是定量的(如地块)，即或这些因素是定量(如淬火温度与等

温温度)但其数值也不作为回归中变量的观察值，而只是代表一种处理（试验方案）；最后，当因素确定后，可以作反复的试验。这两个例子和以前均值假设检验也不同，均值假设检验不考虑因素问题，而方差分析要考虑。

在许多科学研究中都遇到和这两个例子类似的问题。尤其是科学研究中常涉及许多因素，例如研究作物栽培时，要考虑播种期、品种、土质、施肥方式、灌溉方式对产量的影响；在化学反应中要观察原料成分、剂量、催化剂、温度、压力，搅拌速度等对得率的影响。这些因素中要选出影响大的，以进一步安排更细致的试验，而判断一个因素的影响“是否大”的主要方法就是方差分析。

我们所考察的。影响产品指标的因素(如产地，温度等)也称为因子，用大写字母A,B,C表示，例1有一个因子（地块），例2有2个因子（如淬火温度与等温温度）。因素所能处的状况，如甲、乙、丙；60，65，70，75，……，称为因子的水平，简称为水平。例1的因子有三个水平（甲、乙、丙），例2每个因子恰也有3个水平，水平常以A1,A2,B1,B2,...表示。 因子也可以看成是一种变量，其取值不是数，而是水平。例如“产地”是一个变量，它取的值是“北京”、“上海”、“南京”等。这种变量称为属性变量，定性变量或分类变量．本节只讨论一个因子，即一个分类变量的方差分析——单因子方差分析。

方差分析的目的在于找出自变量与因变量之间的线性关系，或自变量对因变量的实验效果。这种实验效果可分为：主效果、交互效果、镶嵌效果。

Qualitative Variable(自变量，又称变量、定性变量)，Classification Variable(分类变量，其数值多半是不连续的。Response Variable（因变量，又称反应变量，其数值则是连续的）

一般地，假设因素A有k个水平：A1,...Ak。第j个水平做实验nj次，得指标y，的nj个数据y1j,y2j,...ynjj。例3.1中n18;n210;n36。通常作如下假设：

（1） 同一个水平Aj下得到的观测值y1j,y2j,...,ynjj，是由于实验过程中各种偶然因素的干扰及测量误差所致，每次实验中这些偶然因素的总和称为实验误差，它们是方差

相同的零均值正态随机变量；

（2） 所有误差相互；

（3） 由于水平的不同，可能会给yij一个定量的确定性的影响，其大小是未知的。 于是我们建立单因子方差分析数学模型

yijiij,i1,2,...,nj j1,2,...,k （3.1）

~N(0,2),它们相互ij其中ij相互，ij~N(0,2)。(4.1)式称为单因素方差分析的数学模型。

判断这个因素的影响是否显著就是要检验假设：

令 nH0:1...k,H1:1...k不全相等

nj，yj(yji)/nj，yyji/n

容易证明yj是j的最小二乘估计 作方差分解

ST(yijy)(yijyjyjy)2

2j1i1j1i1knjknj(yijyj)(yjy)(yijyj)nj(yjy)2

222j1i1j1i1j1i1j1knjknjknjk并令

SSE(yijyj)j1i1knj2，SSAn(yjj1kjy)2

即:

ST(yijy)2-----总的误差平方和

j1i1knjSE(yijyj)2组内差，反映试验误差影响的大小。

j1i1knj

SAnj(yjy)2j1k------组间差，反映因素A的各个水平不同引起

的误差，若A的水平引起的误差显著时，SA 就比较大，反之就比较小。 则有：

STSASE

分别称为组内差和组间差。组内差又称为误差，用以估计实验误差影响的大小；组间差反映因素A的水平不同引起的系统差异。若A的水平不同引起的系统差异(即组间的差距)显著时，SSA就比较大；反之，当A引起的系统差异不显著时，SSA就比较小．而SSE主要是由试验误差引起的。SSA由k个平方之和形成，但有一个恒等式

nj1kj(yjy)0约

束，只有k-1个自由度；同理SSE有n-k个自由度。下列定理给出SSA/SSE的分布

定理： 对于所给的模型，若H0:1...k成立时，

则 FSA/(k1)~F(k1,nk)

SE/(nk)如果F的值超过临界值（通常取为F0.05），就否定H0；当F超过F0.01时，就称为高度

显著。当H0:1...k成立时，F的值不应太大，若F的值大于临界值时，就应否定H0，即认为1,...k间存在显著差异。 一般总用方差分析表表示计算结果，其形式为： 方差来源 因素A 误差 总和

平方和 自由度 K-1 N-K N-1 均方 F值 F 临界值 SA SA/(k1) F0.05 SE ST SE/(k1) 3.2 方差分析的计算

方差分析的计算一般都很复杂，可用SAS软件计算。SAS中有GLM，ANOVA和NESTED等过程可用方差分析。其中GLM过程和ANOVA最常用。

3.2.1 PROC ANOVA

PROC ANOVA 执行对均衡数据的方差分析，其应变量取连续值。所谓均衡数据是指对每个类变量的组合有同样多的响应变量观察值数。如果数据的非均衡的就要用PROC GLM,做方差分析。 ANOVA的语法为： PROC ANOVA DATA= SAS-data-set MANOVA MULTIPASS

OUTSTAT= SAS-data-set;

CLASS variables; /* required */ MODEL dependents=effects / options; /* required */ ABSORB variables; BY variable-list; FREQ variable;

MANOVA H= effects E= effect M= equations...

MNAMES= names PREFIX= name / options; MEANS effects / options;

REPEATED factorname levels(levelvalues) transformation<,...> / options;

TEST H= effects E= effect;

其中PROC ANOVA、CLASS variables和MODEL dependents=effects / options是必须的。 CLASS variables给出分类变量名，这些变量可以是数值的也可以是字符型的，而且CLASS必须在MODEL语句之前。 MODEL dependents=effects ，给出应变量和自变量。其后有若干选项。

3．2．2 PROC GLM

PROC GLM 用最小二乘法拟合线性模型。该过程可以用来做回归、方差分析、协方差分析等。PROC GLM是在广义线性模型的框架下分析数据。处理的数据可以是类变量、离散变量或连续变量。其只要功能是：

A、simple regression

B、 multiple regression C、 analysis of variance (ANOVA), especially for unbalanced data D、 analysis of covariance E、response-surface models F、weighted regression

G、 polynomial regression H、 partial correlation

I、 multivariate analysis of variance (MANOVA) J、repeated measures analysis of variance. 其语法为：

PROC GLM options ; CLASS variable-list;

MODEL dependents= independents / options; /* required */ ABSORB variable-list;

BY variable-list; FREQ variable; ID variable-list; WEIGHT variable; CONTRAST 'label' effect values... / options; ESTIMATE 'label' effect values... / options; LSMEANS effects / options;

MANOVA H= effects E= effect M= equations...

MNAMES= names PREFIX= name / options; MEANS effects / options;

OUTPUT OUT= SAS-data-set keywords= names... ; RANDOM effects / options;

REPEATED factorname levels (levelvalues) transformation<,...> / options;

TEST H= effects E= effect / options;

其中PROC GLM options 、CLASS variable-list、MODEL dependents= independents /是必须的。 MODEL dependents= independents 的选项有：

NOINT INTERCEPT NOUNI SOLUTION TOLERANCE E E1 E2 E3 E4

SS1 SS2 SS3 SS4 ALPHA= p CLM CLI P XPX INVERSE SINGULAR= value ZETA= value *|@||

上述选项可分为五大类：

第一类选项 与截距的界定有关，有两个选项：

（1） NOINT 要求GLM将截距的参数排除在模型之外 （2） INT 要求GLM印出截距的统计鉴定。 第二类选项 与报表的打印有关，有三个选项： （1） NOUNI （2） SOLUTION （3） TOLERANCE

第三类选项 与原假设的检验有关，有九个： E E1 E2 E3 E4 SS1 SS2 SS3 SS4 第四类选项 与控制计算过程有关： XPX INVERSE

第五类选项 可用来调整统计的精度：ZETA

GLM过程主要有四个语句：PROC GLM，CLASS，MODEL和LSMEANS语句

PROC GLM语句用以调用GLM过程，有许多选项，其中DATA＝……选项用以说明GLM过程所加工的数据集。

CLASS语句说明哪些变量是分类变量。方差分析中的因素都是分类变量，响应变量和协方差分析中的协变量不是分类变量。例如class x z;就指示计算机把因子x ,z作为分类变量。

MODEL语句中有等号，等号前是响应变量.

LSMEANS语句用以求待估参数的最小二乘估计。

3．2．3 计算实例

例3.1（续），可用下列SAS程序

data seed; input area $ y; cards;

甲 21 甲 29 甲 24 甲 22 甲 25 甲 30 甲 27 甲 26

乙 20 乙 25 乙 25 乙 23 乙 29 乙 31 乙 24 乙 26 乙 20 乙 21 丙 24 丙 22 丙 28 丙 25 丙 21 丙 26 ;

proc glm; class area; model y=area; run;

执行此程序后，得到的主要输出有

The GLM Procedure Dependent Variable: y

Sum of

Source DF Squares Mean Square F Value Pr > F Model 2 6.7666667 3.3833333 0.32 0.7314

Error 21 223.7333333 10.6539683 Corrected Total 23 230.5000000

R-Square Coeff Var Root MSE y Mean 0.029356 13.18805 3.2042 24.75000

Source DF Type I SS Mean Square F Value Pr > F area 2 6.76666667 3.38333333 0.32 0.7314 Source DF Type III SS Mean Square F Value Pr > F area 2 6.76666667 3.38333333 0.32 0.7314

其中DF列表示自由度，Som of Squares列表示平方和。从 Error行DF列可查得误差自由度为2，从Som of Squares列可查得误差（残差平方和）为223.7333333；由area 行DF列可查得组间查自由度为21，从Type III SS列可查得组间差为6.76666667，从F Value列查得F值为0.32；从 Pr > F查得自由度为(2 ，21)的 F分布随机变量大于0.32的概率为0.7314。它们构成表3.3。

表3.3 例3.1的方差分析表

方差来源 平方和 因素area 6.76666667 误差 总和

表3.3中自由度为(2 ，21)的 F分布随机变量大于0.32的概率为0.7314，大于0.05，所以地块因素（area）作用不显著，即由于地块不同，基本苗数没有差异。

例3.1的SAS程序也可简化为

data seed;/*建立数据库seed*/

do area=1 to 3;/*循环语句开始，area=1，2，3分别表示甲乙丙三地块*/ input y @@;/*读入y的值*/

output;/*将area 、y的值存入建立数据库seed*/ end;/*结束循环*/ cards; 21 20 24 29 25 22 24 25 28 22 23 25 25 29 21 30 31 26 27 24 . 26 26 . . 20 . . 21 . ;

proc glm; class area; model y=area; run;

223.7333333 230.5000000 自由度 2 21 23 均方 3.38333333 10.6539683 F值 0.32 F分布随机变量大 于1.14的概率 0.7314 例3．3 本文件（CLOVER）包含一个自变量及一个因变量。自变量是苜蓿的培养基，下分六种（即3DOK1， 3DOK4，3DOK5，。。。，COMPOS等，因变量是红色苜蓿内氮气的含量。这是一个平衡的实验设计，我们用ANOVAC程序执行单因素的方差分析，并比较各培养基组的平均氮气含量。

data clover;

input strain $ nitrogen @@; cards;

3dok1 19.4 3dok1 32.6 3dok1 27.0 3dok1 32.1 3dok1 33.0 3dok5 17.7 3dok5 24.8 3dok5 27.9 3dok5 25.2 3dok5 24.3 3dok4 17.0 3dok4 19.4 3dok4 9.1 3dok4 11.9 3dok4 15.8 3dok7 20.7 3dok7 21.0 3dok7 20.5 3dok7 18.8 3dok7 18.6 3dok13 14.3 3dok13 14.4 3dok13 11.8 3dok13 11.6 3dok13 14.2 compos 17.3 compos 19.4 compos 19.1 compos 16.9 compos 20.8 ;

proc anova; class strain;

model nitrogen=strain;

means strain / duncan waller; means strain / lsd tukey cldiff; run;

结果分析： 六组含氮量经F检验后证明不尽相同。3DOK5组与3DOK4，3DOK13两组也有显著的不同。3DOK7与3DOK5组的平均数近似，3DOK4与3DOK13两组的结果十分接近。其他的不能达成共识。

这个例子是均衡数据，既可以用ANOVA也可以用GLM实现，而前一个 例子是非均衡数据只能用GLM进行计算。

SAS的GLM过程采用采用哑变量方法，以便作方差分析和协方差分析。该过程的数学原理是建立一般线性模型（类似回归模型）：

yb0b1x1...bmxm （3.2）

并用最小二乘法求解。在方差分析问题中，引进的变量

即只取0或1的变量。GLMx1...xm都是示性变量，

过程对每一因子的每一水平，通过class语句产生1个示性变量。对于例4.1的计算，一般线性模型是

yb0b1x1b2x2b3x3

其中对于地块甲x1=1，否则为0；对于地块乙x2=1，否则为0；对于地块丙x3=1，否则为0。从而设计矩阵是

1 1 1 1 1 0 0 1 0 0 1 0 0 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 0 1 0 1 0 1 0 1 0 1 0 1 对于一般线性模型，由于设计矩阵不满秩，b0,b1,...bm的解不是唯一的，但残差平方和是唯一确定的，从而F检验可行。这从以下例3.1的一般线性模型计算可见。 令S(b0b1x1k1n(k)b2x2(k)b3x3(k)y(k))2。为使S最小，分别对b0,b1,b2,b3

偏导得①、②、③、④

k1nnb0b1x1(k)b2x2(k)b3x3(k)y(k)0 ①

k1n(b0b1x1(k)b2x2(k)b3x3(k)y(k))x1(k)0 ②

k1n(b0b1x1(k)b2x2(k)b3x3(k)y(k))x2(k)0 ③

k1(b0b1x1(k)b2x2(k)b3x3(k)y(k))x3(k)0 ④

注意到xi是是示性函数，①=②+③+④,所以b0,b1,b2,b3的解不是唯一的。由②、③、④分别可得

b0y1b1 ⑤ b0y2b2 ⑥

b0y3b3 ⑦

代入可得最小值Snj(yjy)2。

j13SAS-GLM过程的假设检验总要用到3型平方和(Type III SS)。3型平方和定义如下：设方差分析问题有m个因子，不包含某个因子一般线性模型的残差平方和减去包含全体因子一般线性模型的残差平方和之差称为该因子的3型平方和。例3.1中Type III SS计算步骤如下：先求不含地块因子的最小值min再求求不含地块因子的最小值

k1n(b0y)(yiyji)2=223.7333333

(k)2i1j13nimin(b0b1x1k124(k)b2x2ni(k)b3x3(k)y3(k)2)ni(yiy)2= 230.5000000

i13Type III SS=

i1j13(yiyji)ni(yiy)2=230.5-223.7333333=6.76666667

2i1参见复旦大学编，概率论，第二册，第一分册29页，定理4(Cochran)，37页定理5。

3.3 多因子方差分析

实际问题中往往要考虑多因素的影响，例如，检验药效时，要用多种剂量的药对不同品种动物试验观察其效果；化工反应要考虑反应温度、浓度、反应时间等对得率的影响。和单因子方差分析一样，需要检验这些因子的影响是否显著，这就是多因子方差分析问题。例3.2中的因子就是淬火温度及等温温度，是双因子分析问题。由于因素越多，公式越复杂，我们着重介绍两因素的方差分析，更多因素的方差分析的原理是一样的，用SAS计算的程序也是相同的。

3．3．1一般双因素无重复试验无交互作用的模型

设有两个因子A、B，各有p、q个水平，(Ai,Bj)时y的观测值为yij。假设：

（1） A的第I个水平使y增加i,i1pi0；

（2） B的第j个水平使y增加j,j1qj0 ；

（3） 试验误差是方差相等均值为零的正态随机变量； （4） 试验误差相互。 则得到模型为：

yijijiji1,...,p;j1,...,qpq 0,0,iji1j12~N(0,),且相互ij要检验因素A、B的影响是否显著，就相当于分别检验以下假设：

H101...p0，H11:1...p不全为零； H20:1...q0，H21:1...q不全为零。

令y.j(yi1pij)/p，yi.(yij)/q ，y(yij)/(pq)

j1i1j1qpq可以证明：yi.是i的最小二乘估计；y.j是j的最小二乘估计。作方差分解

ST(yijy)((yijy.jyi.y)(yi.y)(y.jy))2

2i1j1qi1j1pqjpq(yijy.jyi.y)q(yi.y)p(y.jy)2

2pp2qi1j1i1j1并令SSE(yi1j1pqijy.jyi.y)，SSAq(yi.y)，SSBp(y.jy)2

2p2qi1j1其中SSA称为A的组间差，它由p个平方和组成，有一等式

(yi1pi.y)0

约束，自由度是p一1，反映A引起的系统差异；SSB称B的组间差，反映B引起的系统

差异，自由度是q一1；SSE称为误差，反映各种随机因素引起的差异，由pq个平方和组成，包含(p+q-1) 个等式

(yj1qijy.jyi.y)0i1,...p

(yi1ppijy.jyi.y)0j1,...p1

(等式

(yi1iqy.qyi.y)0可用上述等式中前p个之和减去后q-1个之和而得)所以自

由度为(p-1(q-1)． 即有：

STi1,j1p(yi1qp,qijy)2,SAq(yi.y)2,

SBp(y.jy)2j1SEi1,j1(yp,qijyi.y.jy)2则有方差分解：STSASBSE

其中SA反映了A引起的系统误差，SB反映了B引起的系统误差，SE反映了其它随机误差。我们有定理：

定理3。2： 当H01为真时，有：

SA/(p1)FA~F(p1,(p1)(q1))

SE/(p1)(q1)当当H02为真时，有：

SB/(q1)FB~F(q1,(p1)(q1))

SE/(p1)(q1)可由FA,FB是否超过临界值，判断响应的假设是否成立。 上述结果可以类似地用方差分析表的形式给出。

例5.2（续），可用下列SAS程序 data quench; do eqt=1 to 3; do qut=1 to 3; input hard @@; output; end; end; cards; 66 68 66 68 67 65 67 68 ;

proc glm;

class eqt qut;

model hard=eqt qut; run;

执行此程序得到的主要输出有

Dependent Variable: hard

Sum of

Source DF Squares Mean Square F Value Pr > F Model 4 13.11111111 3.27777778 4.21 0.0962 Error 4 3.11111111 0.77777778 Corrected Total 8 16.22222222

R-Square Coeff Var Root MSE hard Mean 0.808219 1.325084 0.881917 66.55556

Source DF Type I SS Mean Square F Value Pr > F eqt 2 1.55555556 0.77777778 1.00 0.4444 qut 2 11.55555556 5.77777778 7.43 0.0450

Source DF Type III SS Mean Square F Value Pr > F eqt 2 1.55555556 0.77777778 1.00 0.4444 qut 2 11.55555556 5.77777778 7.43 0.0450

从 Error行DF列可查得误差自由度为4，从Som of Squares列可查得误差（残差平方和）为3.11111111。由eqt行DF列可查得等温温度组间差自由度为2，从Type III SS列可查得组间差为1.55555556，从F Value列查得等温温度F值为1.00；从 Pr > F查得自由度为(2 ，2)的 F分布随机变量大于1.0的概率为0.444 4。由qut行DF列可查得淬火温度组间差自由度为2，从Type III SS列可查得淬火温度组间差为11.55555556，从F Value列查得淬火温度F值为7.43；从 Pr > F查得自由度为(2 ，2)的 F分布随机变量大于7.43的概率为0.0450。

对于例3.2全部计算结果可以列成方差分析表——表35.4

表3.4 方差来源 因素A（等温温度） 因素B（淬火温度） 误差 总和 平方和 1.5556 11.5556 3.1110 16.2222 自由度 2 2 4 8 平均平方和 0.7778 5.7778 0.7778 1.0 7.43 F 临界值 F0.05=6.94 F0.05=6.94 由表3.4可见：不否定H10，否定H20从而B因素的作用是显著的，A因素的作用不显著，即淬火温度对铣刀硬度有显著影响，而等温温度(即闷火温度)没有显著影响．

模型式的优点是不需做重复试验，即对A,B不同水平，每一组合只需做一次试验，这大大节省人力、物力，但也有缺点，如例4.3所示。

对于多因子方差分析，GLM过程对每一因子的每一水平，通过class语句产生1个示性变量。

例3.2中有6个示性变量（3个示性变量A1,A2,A3表示等温温度，3个示性变量。一般线性模型是 B1,B2,B3表示沾火温度）

yb0b1A1b2A2b3A3b4B1b5B2b6B3

由于设计矩阵不满秩，b0.b1,b2,b3,b4,b5,b6解不唯一。但残差平方和是唯一确定的，从而F检验可行。

例3.2中Type III SS计算步骤如下：先求包含所有因子的最小值

min(b0b1A1k19(k)b2A2(k)b3A3(k)b4B1(k)b5B2(k)b6B3(k)y(k))2=3.1111

再求不含等温温度因子的最小值

min(b0b4B1k19(k)b5B2(k)b6B3(k)y(k))2=4.66667

等温温度因子的3型平方和就是

Type III SS=4.66667-3.11111=1.55556

然后求不含沾火温度因子的最小值

min(b0b1A1k19(k)b2A2(k)b3A3(k)y(k))2=14.66667

沾火温度因子的3型平方和就是

Type III SS=14.66667-3.11111=11.55556

3．3．2 双因素有重复试验有交互作用模型

例3.4对三种推进器（M1、M2、M3）和四种燃料（F1、F2、F3、F4）试验火箭 的射程，得表5.5。推进器M与燃料F的作用是否显著?

表3.5 火箭的射程

推进器 燃料 M1 58.2，52.6 49.1，42.8 60.1，58.3 75.8，71.5 M2 56.2，41.2 .1，50.5 70.9，73.2 58.2，51.0 M3 65.3，60.8 51.69，48.4 39.2，40.7 48.7，41.4 F1 F2 F3 F4

以F与M为因素按模型(4.4)式分析时，可采用SAS程序 data rockey; do f=1 to 4; do m=1 to 3; do rep=1 to 2; input r@@; output; end; end; end; cards;

58.2 52.6 56.2 41.2 65.3 60.8 49.1 42.8 .1 50.5 51.6 48.4 60.1 58.3 70.9 73.2 39.2 40.7 75.8 71.5 58.2 51.0 48.7 41.4 ;proc glm; class f m; model r=f m; run;

从所得结果可见推进器和燃料引起的射程差异不显著，是否这些推进器和燃料对射程的影响不大呢?并不是这样的，从表（4.5）中可见，M1和F4的组合射程很远；M1和F3的组合射程很近，M1不同水平的射程值作了平均后就不大不小了的缘故。因而模型(4.4)不能反映推进器和燃料的搭配对射程的影响。为了反映两种因素搭配后的影响，必须引入“交互作用”的概念。

设有两个因子A、B，各有p、q个水平，对A、B的每个组合(Ai,Bj)作k次重复试验（为了检验交互作用必须有重复实验）。y的观测值为yijt,i1,...,p;j1,...,q;t1,...,k。假设：

（1） 的第I个水平使y增加i,i1pi0；

（2） B的第j个水平使y增加j,j1qj0 ；

pq（3） (Ai,Bj)的交互作用使Y增加ij,i1ij0,ij0,i1,...,p;j1,...,q。

j1（4） 试验误差是方差相等均值为零的正态随机变量；

（5） 试验误差相互。 于是有数学模型：

yijtijijitji1,...,p;j1,...,q,t1,...,kqp

i0,j0,j1i1qpij0,ij0,i1,...,p;j1,...,qj1i12ijt~N(0,),且相互其中{i},{j}分别称为因子A、B的主效应；{ij}称为A与B的交互效应。当考虑交互作用时，方差分析问题就要检验

H10:1...p0， H11:1...p不全为零； H20:1...q0，H21:1...q不全为零； H30:11...pq0，H31:11...pq不全为零。

yijyi.

1k1q1pyt1qkijtyj1pi1ij,,

y.jyijp,q1yyijpqi1,j1kpqjSTt1(yi1j1qijkijky)2t1pk((yi1j1qpyij)(yijy.jyi.y)(yi.y)(y.jy))2

k(yijy.jyi.y)kq(yi.y)2

2i1j1qpi1kp(y.jy)2j1t1k(yi1j1ijpqijkyij)2

并令SSABk(yi1j1pqy.jyi.y)，SSAkq(yi.y)2，

2pi1SSBkp(y.jy)，SSE2j1t1qk2(yy)ijkij i1j1pq其中SSA称为A的组间差，反映A引起的发散程度，统计它们包含等式个数容易看出，自

由度是p一1；SSB称B的组间差，反映B引起的发散，自由度是q一1；SSAB反映AB交互作用效应引起的发散；自由度是(p-1)(q-1),SSE称为误差，反映各种随机因素引起的发散，自由度是pq(k-1)。

p,q,kSTi1,j1,t1p(yi1qijty)2,SAqk(yi.y)2,SBpk(y.jy)2j1

SABk

i1,j1p,q,k2(yyyy)iji..jp,qSE

i1,j1,t1(y

ijtyij)2则有方差分解：STSASBSABSE

其中SA反映了A引起的系统误差，SB反映了B引起的系统误差，SAB反映了A与B的交互作用，SE反映了其它随机误差。我们有定理： 定理： 当H01为真时，有：

SA/(p1)FA~F(p1,pq(k1))

SE/pq(k1)当当H02为真时，有：

SB/(q1)FB~F(q1,pq(k1))

SE/pq(k1)当H03为真时，有

FABSAB/(p1)(q1)~F((p1)(q1),pq(k1))

SE/pq(k1)可由FA,FB,FAB是否超过临界值，判断响应的假设是否成立。而且这3个检验是相互的。所以当FA,FB,FAB之一大于临界值时，A,B因素的主效应或交互作用效应显著。上述结果可以类似地用方差分析表的形式给出。

表5.8 方差分析

方差来源 因子A 平方和 SSA 自由度 p一1 平均平方和 SSA／(p一1) F SSA/(p1) SSE/pq(k1)SSB/(q1) SSE/pq(k1)SSAB/(p1)(q1) SSE/pq(k1) 因子B SSB q一1 SSB／(q一1) 交互作用 误差 总和 SSAB (p-1)(q-1) SSAB/(p-1(q-1) SSE TSS pq(k-1) Pq（k-1）

用SAS软件计算交互作用时，只需在有关因子中加‘*’号。例如MODEL z=a b c a*b; 指示计算机把a,b,c作为因子;考虑a,b的交互作用。

例3.4 的解：考虑交互作用的方差分析可用下述模型

yijtijijijt其中

i1,...4j1,2,3t1,2

iji0，j0，ij0i1,...4，

0j1,2,3，

ijt~N(0,2)且相互。

并用程序

data rockey;/*建立数据库rockey*/

do f=1 to 4;/*开始循环，燃料因子取水平1、2、3、4*/ do m=1 to 3; /*循环，发动机因子取水平1、2、3*/

do rep=1 to 2;/*每个组合重复2次，变量rep计算中不起作用*/ input r@@;/*读入变量r*/ output;

end; end; end; cards;

58.2 52.6 56.2 41.2 65.3 60.8 49.1 42.8 .1 50.5 51.6 48.4 60.1 58.3 70.9 73.2 39.2 40.7 75.8 71.5 58.2 51.0 48.7 41.4 ;

proc glm; class f m;

model r=f m f*m;/*对主效应f、m和交互作用效应 f*m*作方差分析/

lsmeans f m f*m; /*求主效应f、m和交互作用效应 f*m*的最小二乘估计/ run;

执行此程序后得到主要输出如下

General Linear Models Procedure

Dependent Variable: R

Sum of Mean

Source DF Squares Square F Value Pr > F Model 11 2401.348333 218.304394 11.06 0.0001 Error 12 236.950000 19.745833 Corrected Total 23 2638.298333

R-Square C.V. Root MSE R Mean 0.910188 8.0809 4.443628 .99167

从 Error行DF列可查得误差自由度为12，从Som of Squares列可查得误差（残差平方和）

为236.95，从 Pr > F的值为0.0001可见，一般线性模型的线性关系是高度显著的。

Dependent Variable: R

Source DF Type I SS Mean Square F Value Pr > F F 3 261.675000 87.225000 4.42 0.0260 M 2 370.980833 185.490417 9.39 0.0035 F*M 6 1768.692500 294.782083 14.93 0.0001

Source DF Type III SS Mean Square F Value Pr > F F 3 261.675000 87.225000 4.42 0.0260 M 2 370.980833 185.490417 9.39 0.0035 F*M 6 1768.692500 294.782083 14.93 0.0001

由F行DF列可查得燃料组间差自由度为3，从Type III SS列可查得燃料组间差为261.675000，从F Value列查得燃料F值为4.42；从 Pr > F查得发动机自由度为(3 ，12)的 F分布随机变量大于4.42的概率为0.0260。由M行DF列可查得发动机组间差自由度为2，从Type III SS列可查得发动机组间差为370.980833，从F Value列查得淬火温度F值为9.39；从 Pr > F列查得自由度为(2 ，12)的 F分布随机变量大于9.39的概率为0.0035。

由F*M行DF列可查得交互作用自由度为6，从Type III SS列可查得交互作用组间差为1768.692500，从F Value列查得交互作用F值为14.93；从 Pr > F查得自由度为(2 ，12)的 F分布随机变量大于14.93的概率为0.0001。

此表配合前表即可得例3.3方差分析表——表3.4

表3.4 例3.3方差分析表 方差来源 F(燃料) M(发动机) F*M(交互作用) SSE(误差) 总和 自由度 3 2 6 12 23 平方和 261.657000 370.980933 1768.692500 236.9500 2638.29833 平均平方和 87.225000 185.49017 294.782086 19.745833 F值 4.42 4.39 14.93· Pr>F 0.0260 0.0035 0.0001

从方差分析表可见：F(燃料)对应的概率在0.01与0.05之间，说明燃料的作用显著；M(发动机)及F*M(交互作用)对应的概率小于0.01，说明发动机与交互作用的影响高度显著。

Least Squares Means F R LSMEAN 1 55.7166667 2 49.4166667 3 57.0666667 4 57.7666667

上表给出不同燃料的射程均值（的最小二乘估计）。

Least Squares Means M R LSMEAN 1 58.5500000 2 56.9125000 3 49.5125000

上表给出不同发动机的射程均值（的最小二乘估计）。

Least Squares Means F M R LSMEAN 1 1 55.4000000 1 2 48.7000000 1 3 63.0500000 2 1 45.9500000 2 2 52.3000000 2 3 50.0000000 3 1 59.2000000 3 2 72.0500000 3 3 39.9500000 4 1 73.6500000 4 2 .6000000 4 3 45.0500000

上表给出不同燃料与发动机组合的射程均值（的最小二乘估计）。由表可见，最远的射程是由第四种燃料与第一种发动机组合形成的，平均射程为73．65公里。

注意：考虑交互作用时，不要写赋值语句．

3．3．3 PROC ANOVA的进一步讨论和实例

PROC ANOVA 前面已经做了一些介绍。下面进一步作一些讨论。 假设有一个平衡的实验设计，含有三个自变量（A、B、C），因变量以Y记，则此三个因子的主效应方差分析可以用下面的程序来执行： PROC ANOVA； CLASS A B C； MODEL Y=A B C；

交互效应的统计模型

以上述三个因子为例，其对应的主效应以及交互效应可用下列程序来执行： PROC ANOVA； CLASS A B C ； MODEL Y=A B C A*B B*C A*C A*B*C； 当实验设计含多个自变量时，交互效应会变得繁杂。此时可用“|”来简化。如上例的MODEL指令可简化为： MODEL Y=A|B|C；

它相当于 MODEL Y=A

B C A*B B*C A*C A*B*C。

例 3．5

title 'randomized complete block'; data rcb;

input block trtment $ yield worth; cards; 1 a 32.6 112 1 b 36.4 130 1 c 29.5 106 2 a 42.7 139 2 b 47.1 143 2 c 32.9 112 3 a 35.3 124 3 b 40.1 134 3 c 33.6 116 proc anova;

class block trtment;

model yield worth = block trtment; run;

例3．6

*-----------------split plot--------------------*

| b defines subplots within a*block whole plots.| | the whole plot effects must be tested with a | | test statement against block*a. the subplot | | effects can be tested against the residual. | *-----------------------------------------------* ; data split;

input block 1 a 2 b 3 response ; cards; 142 40.0 141 39.5 112 37.9 111 35.4 121 36.7 122 38.2 132 36.4 131 34.8 221 42.7 222 41.6

212 40.3 211 41.6 241 44.5 242 47.6 231 43.6 232 42.8

proc anova;

class block a b;

model response = block a block*a b a*b ; test h=a e=block*a; title 'split plot design'; run;

例3．7

data beets;

do harvest=1 to 2; do rep=1 to 6; do col=1 to 6;

input variety y @; output; end; end; end; cards;

3 19.1 6 18.3 5 19.6 1 18.6 2 18.2 4 18.5 6 18.1 2 19.5 4 17.6 3 18.7 1 18.7 5 19.9 1 18.1 5 20.2 6 18.5 4 20.1 3 18.6 2 19.2 2 19.1 3 18.8 1 18.7 5 20.2 4 18.6 6 18.5 4 17.5 1 18.1 2 18.7 6 18.2 5 20.4 3 18.5 5 17.7 4 17.8 3 17.4 2 17.0 6 17.6 1 17.6 3 16.2 6 17.0 5 18.1 1 16.6 2 17.7 4 16.3 6 16.0 2 15.3 4 16.0 3 17.1 1 16.5 5 17.6 1 16.5 5 18.1 6 16.7 4 16.2 3 16.7 2 17.3 2 17.5 3 16.0 1 16.4 5 18.0 4 16.6 6 16.1 4 15.7 1 16.1 2 16.7 6 16.3 5 17.8 3 16.2 5 18.3 4 16.6 3 16.4 2 17.6 6 17.1 1 16.5 ;

proc anova;

class col rep variety harvest;

model y= rep col variety rep*col*variety harvest harvest*rep harvest*variety;

test h=rep col variety e=rep*col*variety; test h=harvest e=harvest*rep; run;

例3．8

*------------strip-split plot design: winter barley-------------*

| | |fertilizer treatments are laid out in vertical strips which are| |then split into calcium effect sub-plots. soil type is then | |stripped across the split plot experiment. the whole experiment|

|is replicated three times. | |data from the notes of g. cox and a. rotti *---------------------------------------------------------------*;

data barley;

do rep=1 to 4;

do soil=1 to 3; * 1=d 2=h 3=p; do fertil=0 to 3; do ca=0,1;

input y @; output; end; end; end; end; cards;

4.91 4.63 4.76 5.04 5.38 6.21 5.60 5.08 4.94 3.98 4. 5.26 5.28 5.01 5.45 5.62 5.20 4.45 5.05 5.03 5.01 4.63 5.80 5.90 6.00 5.39 4.95 5.39 6.18 5.94 6.58 6.25 5.86 5.41 5. 5.41 5.28 6.67 6.65 5.94 5.45 5.12 4.73 4.62 5.06 5.75 6.39 5.62 4.96 5.63 5.47 5.31 6.18 6.31 5.95 6.14 5.71 5.37 6.21 5.83 6.28 6.55 6.39 5.57 4.60 4.90 4.88 4.73 5. 6.20 5.68 5.72 5.79 5.33 5.13 5.18 5.86 5.98 5.55 4.32 5.61 5.15 4.82 5.06 5.67 5. 5.19 4.46 5.13 4.90 4.88 5.18 5.45 5.80 5.12 4.42 ;

* note that since the model is completely specified and several * error terms are present, the test statement must be used to * obtain the proper test statistics. the top portion of the output

|

| | * should be ignored, since the residual error term is not * meaningful here;

proc anova;

class rep soil ca fertil; model y=rep

fertil fertil*rep

ca ca*fertil ca*rep(fertil) soil soil*rep

soil*fertil soil*rep*fertil soil*ca soil*fertil*ca soil*ca*rep(fertil);

test h=fertil e=fertil*rep; test h=ca ca*fertil e=ca*rep(fertil); test h=soil e=soil*rep; test h=soil*fertil e=soil*rep*fertil; test h=soil*ca

soil*fertil*ca e=soil*ca*rep(fertil); means fertil ca soil ca*fertil; title 'strip-split plot'; run;

例3．9 （单因素的方差分析） 本资料文件（PLANTS）。该实验的目的在于比较金鱼草在七种土壤（TYPE）里的生长的情形，每一种土壤中种三朋金鱼草（BLOCK）。这个实验是一个平衡的设计。

*---------------snapdragon experiment---------------* | as reported by stenstrom, 1940, an experiment was | | undertaken to investigate how snapdragons grew in | | various soils. each soil type was used in three |

| blocks. | *---------------------------------------------------*;

data plants;

input type $ @; do block=1 to 3;

input stemleng @; output; end; cards;

clarion 32.7 32.3 31.5

clinton 32.1 29.7 29.1 knox 35.7 35.9 33.1 o'neill 36.0 34.2 31.2 compost 31.8 28.0 29.2 wabash 38.2 37.8 31.9 webster 32.5 31.1 29.7 ;

proc glm;

class type block;

model stemleng=type block; run;

proc glm order=data; class type block;

model stemleng=type block / solution; means type / waller regwq;

*-type-order------------clrn-cltn-knox-onel-cpst-wbsh-wstr; contrast 'compost vs others' type -1 -1 -1 -1 6 -1 -1; contrast 'river soils vs.non' type -1 -1 -1 -1 0 5 -1,

type -1 4 -1 -1 0 0 -1; contrast 'glacial vs drift' type -1 0 1 1 0 0 -1; contrast 'clarion vs webster' type -1 0 0 0 0 0 1; contrast 'knox vs oneill' type 0 0 1 -1 0 0 0; run;

结果分析：

金鱼草的生长情形随七种土壤以及三种盆栽的不同而改变。（F=10。80，P=0。0002） WABASH的土壤最优，COMPOST的土壤最劣，KNOX、O‘NELL两种不分轩轾，CLARION、WEBSTER两种土壤也无显著的差异。

例3．10 （双因素的方差分析）

这是一个非平衡的实验设计。两个自变量DRUG、DISEASE。注意其中的MODEL指令。 *-------------------------------------------------------------*

| a two-way analysis of variance example using the data from: | | kutner, michael h., the american statistician, aug.74, p.98.| | original data source: afifi and azen(1972), statistical | | analysis: a computer-oriented approach. academic press, ny, |

| p.166. | *------------------------------------------------------------ *;

data a;

input drug disease @; do i=1 to 6;

input y @; output; end; cards;

1 1 42 44 36 13 19 22 1 2 33 . 26 . 33 21 1 3 31 -3 . 25 25 24 2 1 28 . 23 34 42 13 2 2 . 34 33 31 . 36 2 3 3 26 28 32 4 16 3 1 . . 1 29 . 19 3 2 . 11 9 7 1 -6 3 3 21 1 . 9 3 . 4 1 24 . 9 22 -2 15 4 2 27 12 12 -5 16 15 4 3 22 7 25 5 12 . ;

proc glm;

class drug disease;

model y=drug disease drug*disease /ss1 ss2 ss3 ss4 ; run;

结果表明：DRUG的效应达到显著水平，然而DISEASE或两者的交互作用均未达到显著水平

例3．11 （三因素的方差分析与平均数的比较）

本资料文件（ONE）这个实验探讨电击对萎缩肌肉的影响。 REP—实验重复的次数（一次或二次） TIME—电流通过的时间（一到四级时间单位） CURRENT—电流的强度（一到四级） NUMBER—每天受电击的次数（一到三次）

Y—肌肉经电击后的重量（因变量） 每一位受试者接受两次电击，但电流通过肌肉的时间，电流的强度与每天接受的次数因人而异。注意CONTRAST指令的写法。

第一个CONTRAST指令是针对TIME的主效果以及CURRENT在每一个TIME组内的简单主效应

第二个CONTRAST指令O则直接比较CURRENT的第一与第二组的平均数。

data one;

do rep=1 to 2; do time=1 to 4;

do current=1 to 4;

do number=1 to 3; input y @@;

output; end; end; end; end; cards;

72 74 69 61 61 65 62 65 70 85 76 61 67 52 62 60 55 59 65 67 72 60 57 66 72 72 43 43 63 66 72 56 75 92 57 56 78 60 63 58 61 79 68 73 86 71 46 74 58 60 52 71 71 53 65 66 44 58 57 55 51 62 61 79 60 78 82 53 50 61 56 57 56 56 56 71 56 58 69 46 55 56 55 57 66 62 59 58 88 ;

proc glm;

class rep current time number; model y=rep current|time|number; contrast 'time in current 3'

time 1 0 0 -1 current*time 0 0 0 0 0 0 0 0 1 0 0 -1, time 0 1 0 -1 current*time 0 0 0 0 0 0 0 0 0 1 0 -1, time 0 0 1 -1 current*time 0 0 0 0 0 0 0 0 0 0 1 -1; contrast 'curr 1 vs. curr 2' current 1 -1; run;

结果表明整个实验设计只有一个效果打到显著水平，即电流的主效果，其余的效果均未达到显著水平。

例3．12 多变量的方差分析 本资料文件（SKULL），旨在决定性别（SEX）对四种反应（LENGTH，BASILAR，ZYGOMAT，POSTORB）的效果。

*---------multivariate analysis of variance-------* | data from a. anderson, oregon state university. | | four different response variables are measured. | | the hypothesis that we want to test is that sex | | does not affect any of the four responses. | *-------------------------------------------------*;

data skull;

input sex $ length basilar zygomat postorb @@; cards;

m 60 4962 3286 1100 m 6252 4773 3239 1061 m 5772 4480 3200 1097

m 62 4806 3179 10 m 6622 5113 3365 1071 m 6656 5100 3326 1012 m 41 4918 3153 1061 m 6281 4821 3133 1071 m 6606 5060 3227 10 m 6573 4977 3392 1110 m 6563 5025 3234 1090 m 6552 5086 3292 1010 m 6535 4939 3261 1065 m 6573 4962 3320 1091 m 6537 4990 3309 1059 m 6302 4761 3204 1135 m 49 4921 3256 1068 m 81 4887 3233 1124 m 6368 4824 3258 1130 m 6372 4844 3306 1137 m 6592 5007 3284 1148 m 6229 4746 3257 1153 m 6391 4834 3244 1169 m 6560 4981 3341 1038 m 6787 5181 3334 1104 m 6384 4834 3195 10 m 6282 4757 3180 1179 m 6340 4791 3300 1110 m 6394 4879 3272 1241 m 6153 4557 3214 1039 m 6348 4886 3160 991 m 6534 4990 3310 1028 m 6509 4951 3282 1104 f 6287 4845 3218 996 f 6583 4992 3300 1107 f 6518 5023 3246 1035 f 32 4790 3249 1117 f 50 4888 3259 1060 f 6379 4844 3266 1115 f 24 4855 3322 1065 f 6615 5088 3280 1179 f 6760 5206 3337 1219 f 6521 5011 3208 9 f 16 48 3200 1001 f 6511 4910 3230 1100 f 60 4997 3320 1078 f 6780 5259 3358 1174 f 6336 4781 3165 1126 f 72 49 3125 1178 f 76 46 3148 1066 f 6276 4709 3150 1134 f 6693 5177 3236 1131 f 6328 4792 3214 1018 f 6661 5104 3395 1141 f 6266 4721 3257 1031 f 6660 5146 3374 1069 f 6624 5032 3384 11 f 6331 4819 3278 1008 f 6298 4683 3270 1150 ;

proc glm; class sex;

model length basilar zygomat postorb=sex; manova h=sex / printe printh;

title 'multivariate analysis of variance'; run;

结果表明他们的值经转换为F值后，大小均完全一致。所以男女在这四种反映上的结果不相上下。

例3．13

data mileage;

input mph mpg @@; cards;

20 15.4 30 20.2 40 25.7 50 26.2 50 26.6 50 27.4 55 . 60 24.8 ;

proc glm;

model mpg=mph mph*mph / p clm; output out=pp p=mpgpred r=resid; run;

proc plot data=pp;

plot mpg*mph='a' mpgpred*mph='p' / overlay; run;

例3．14 data dogs;

input drug $ depl $ hist0 hist1 hist3 hist5; lhist0=log(hist0); lhist1=log(hist1); lhist3=log(hist3); lhist5=log(hist5); cards;

morphine n .04 .20 .10 .08 morphine n .02 .06 .02 .02 morphine n .07 1.40 .48 .24 morphine n .17 .57 .35 .24 morphine y .10 .09 .13 .14 morphine y .12 .11 .10 . morphine y .07 .07 .06 .07 morphine y .05 .07 .06 .07 trimeth n .03 .62 .31 .22 trimeth n .03 1.05 .73 .60 trimeth n .07 .83 1.07 .80 trimeth n .09 3.13 2.06 1.23 trimeth y .10 .09 .09 .08 trimeth y .08 .09 .09 .10 trimeth y .13 .10 .12 .12 trimeth y .06 .05 .05 .05 ;

proc glm;

class drug depl;

model lhist0--lhist5=drug depl drug*depl/nouni; repeated time 4 (0 1 3 5) polynomial/short summary; run;

例3．15

data machine;

input machine person rating @@; cards;

1 1 52.0 1 2 51.8 1 2 52.8 1 3 60.0 1 4 51.1 1 5 50.9 1 5 51.8 1 5 51.4 1 6 46.4 1 6 44.8 2 1 .0 2 2 59.7 2 2 60.0 2 2 59.0 2 3 68.6 2 4 63.2 2 4 62.8 2 4 62.2 2 5 .8 2 5 65.0 2 6 44.2 2 6 43.0 3 1 67.5 3 1 67.2 3 1 66.9 3 2 61.7 3 2 62.3 3 3 70.8 3 3 70.6 3 3 71.0 3 4 66.2 3 4 .0 3 5 72.1 3 5 72.0 3 5 71.1 3 6 61.4 3 6 60.5 ;

1 4 52.3 1 6 49.2 2 3 65.8 2 6 43.7 3 2 61.5 3 4 .1 3 6 62.0 proc glm;

class machine person;

model rating = machine person machine*person; random person machine*person/ test; run;

3.4协方差分析

在方差分析的试验中，因素是可以控制的，每一个水平都可以人为地达到。但实际问题中有些因素是普通变量而且它的值是不能控制或难以控制的。例如研究某几种治疗方案对疾病治愈率的影响时，每个病人自身健康程度的指标是不能控制的。提炼出这些因素的影响并作方差分析的方法称为协方差分析。换句话说，回归分析、方差分析和协方差分析都研究相关关系，回归分析是以普通变量为自变量和因变量的数学模型；方差分析是以属性变量为自变量，普通变量为因变量的数学模型；协方差分析则是以普通变量和分类变量为自变量，以普通变量为因变量的数学模型。协方差分析中普通变量作为自变量时称为协变量。 协方差分析有两种计算方法，一种是通过方差分解实施的，另一种是通过一般线性模型实施的。希望深入钻研前一种的读者可参阅中国科学院数学研究所统计组编的《方差分析》和王松桂《线性模型的理论及应用》等书。我们仅对单因子单协变量情形概述其原理，总用SAS的GLM 过程计算与分析，GLM 过程总以一般线性模型为计算模型。 首先让我们看两个例子，

例5.4 比较三种猪饲料A1,A2,A3对猪催肥的效果，测得每头猪增加的重量y和初始重 量x数据如表4.10。问三种饲料对猪的催肥有无显著的不同?哪种饲料效果最好?初始重量与猪增加的重量有无明显的关系?

表5.10猪体重增加数据

x y x y x y x y x y x y x y x y A1 15 85 13 83 11 65 A3 22 24 91 20 83 12 76 12 80 16 91 14 84 17 90 18 94 A2 17 97 16 90 18 100 18 95 21 103 22 106 19 99 23 95 25 100 27 102 30 105 32 110

例5.6 考察生产合成纤维时，对纤维弹性y影响的2因素：收缩率A和拉伸倍数B。A,B各取4个水平，电流周波x是不可控制因素。生产记录如表4.12。分析周波，收缩率和拉伸倍数的影响是否显著。

表5.12

A1 |A2 |A3 |A4

x y x y x y x y x y x y x y x y B1 B2 49.0 71 49.2 73 49.8 73 49.8 75 49.9 76 49.8 73 49.7 75 49.8 73 49.5 72 49.3 73 49.9 76 49.8 74 50.2 79 50.1 77 49.4 73 49.4 72 49.7 75 49.5 73 50.1 78 50.0 77 49.7 74 50.0 75 49.5 70 49.6 71 49.9 77 49.7 75 49.6 74 49.3 74 49.5 74 49.2 73 49.0 69 48.9 69 B3 B4 例5.4中猪饲料是属性变量，取值A1,A2,A3；每头猪初始重量x是普通变量，它的 值是不可控制的，需要分析它与增加重量的关系，我们把它作为普通变量处理，即协变量。例4.4中收缩率A和拉伸倍数B是分类变量，电流周波x是普通变量，它的值是不可控制，需要分析它与增加重量的关系，我们把它作为普通变量处理，即协变量。。

不考虑协变量时，例5.4的求解不合理：问题中A1,A2,A3作为分类变量A的值。直接作方差分析，可以执行下列SAS程序。

data pig;/*建立数据库pig*/

input a $ @;/*读入变量a的值并转读别的变量*/

do i=1 to 8;/*读变量i的值1-8，并循环*/ input x @;/*读入变量x的值并转读别的变量*/ input y @;/*读入变量y的值*/ output;/*将读入的值存入数据库*/ end;/*结束循环*/ cards;

A1 15 85 13 83 11 65 12 76 12 80 16 91 14 84 17 90 A2 17 97 16 90 18 100 18 95 21 103 22 106 19 99 18 94 A3 22 24 91 20 83 23 95 25 100 27 102 30 105 32 110 ; proc glm; class a; model y=a ; lsmeans a; run;

执行此程序主要得到二表格如下 General Linear Models Procedure Dependent Variable: Y Source DF Type I SS Mean Square F Value Pr > F A 2 1317.583333 658.791667 11.17 0.0005 Source DF Type III SS Mean Square F Value Pr > F A 2 1317.583333 658.791667 11.17 0.0005 从上表Pr > F概率0.005可见三种饲料的作用有高度显著差异 General Linear Models Procedure Level of --------------Y-------------- A N Mean SD A1 8 81.7500000 8.34522960 A2 8 98.0000000 5.12695956 A3 8 96.8750000 8.99900788 从上表可见第二种饲料最好，第三种次之．第一种最差。但这结论有问题：吃A1饲料猪初始重量比较低(平均13.75)，吃A2饲料猪初始重量（平均18.628）稍高，吃A3饲料的猪初始重量(平均25.375)最重，初重量越大的猪长得越快，应当把初始重量的影响扣除才能看出三种饲料对催肥的效果。由散点图可见，对每一种饲料而言，x与y有明显的线性关系。

于是我们应当采用协方差分析模型。假设：(1)猪的初重x使增重y增加b0b1x；(2)Ai饲料增重i；(3)猪的增重受随机误差影响值ij，ij相互，零均值、方差相同。得模型

yijb0ibxijij再令

i1,2,3j1,2...8

ib0ii0，得数学模型

i1,2,3j1,2...8 （5.5）

yijbxijiij 其中

i0，ij相互。

对于例5.5也可以采用协方差分析模型：假设周波对弹性有线性影响，我们以x为协变量建立数学模型

yijkbxijkijkijki,j1,2,3,4k1,2 (4.6)

其中i,j分别是A和B的主效应，ij是交互作用效应，x是协变量，ijk相互，服从N(0，)。

模型(5.5)式和(5.6)式可以容易地推广，例如有3个因素，2个协变量时，可以考虑模型(也可以加上交互作用项)

2yijklijkb1x(1)ijklb2x(2)ijklijkl

协方差分析的公式推演比较复杂，我们仅讨论一个因素一个协变量的情况。 设A有m个水平，每个水平试验r次。考虑数学模型

yijbxijiij

i,j1,.m..j1,2,.r..

i0ij~N(0,2)其中ij相互。

首先介绍第一种方法 令bii选bi，b使 SS(yijbibxij)2 （4.6）

最小，这相当在不同水平分别作回归。对bi偏导得

(yijbibxij)0 （4.7）

令n=mr,yiyj1rij/r,xixij/r,yj1j1rryi1mij/n,则可将(4.7)式化为

biyibxi （5.8）

(5.6)式对b偏导得

(yijbibxij)xij0

将（4.8）代入得

[(yb*(yijyibxibxij)(xijxi)0

解得回归平方和及b的估计值

ijyi)(xijxi)](xijxi)2[(y，SSEijyi)(xijxi)]2(xijxi)2

代入(4.8)式得bi的估计值bi*yib*xi。

SSE就是协方差模型的误差平方和，它由n个平方构成，但因m+1个参数b*,bi*是由

yij,xij算出的，因而要扣去m+1个自由度，所以自由度是n—m一1。

再将所有数据yij,xij混在一起，作回归，即选择,使S(yijaxij)2最

小，得残差平方和SE(yijy)2[(yijyi)(xijxi)]2，其自由度为n-2。

(xijxi)2但SE中含有因素A的影响，故 SSA=SE-SSE就是提炼出的A因素的影响。于是SSA的自由度为(n-2)-(n-m-1)=m-1。由此。我们得到单因素单协变量的方差分析表4.10

表5.10单因素协方差分析表

方差来源 回归 因素A 误差 总和 平方和 SSR SSA SSE 16.2222 自由度 1 m-1 n-m-1 n-1 平均平方和 SSR SSA/(m-1) SSE/(n-m-1) F 临界值 SSR(n-m-1)/SSE SSA(n-m-1)/SSE(m-1) F0.05 F0.05

以下只对只有1个因子1个协变量的协方差分析模型，介绍第二种理论：采用示性变量uj化协方差分析模型为一般线性模型

y(k)b0bx(k)1u1

(k)...mum(k)(k)k1,...n （5.9）

i0(k)~N(0,2)同方差分析的讨论一样，对此一般线性模型用最小二乘法求出b0,i的估计及3型平方和，因为设计矩阵不满秩，b0,a1,...am的估计不惟一；但是b的最小二乘估计和3型平方和是唯一的。由（4.9）易见：协变量只有一个，回归的自由度是1；示性变量有m个，相互的只有m-1个，因素A的自由度是m-1；误差的自由度是n-1-(m-1)=n-m。

协方差分析可用SAS解，仍用GLM过程计算，在model语句等号后跟因子和协变量，再加斜杠，后跟solution可求出b,b0,a1,...am的估计值。

例4.4的协方差分析解：用数学模型

yijbxijiiji1,2,3j1,2,...8

i0ij~N(0,2)具体计算采用SAS程序 data pig； input a$@ do i＝1 to 8； inputx@；

inputy@； output； end； cards；

A1 15 85 13 S3 11 65 12 76 12 80 16 91 14 84 17 90 A2 17 97 16 90 18 100 18 95 21 103 22 106 19 99 18 94 A3 22 24 91 20 83 23 95 25 100 27 102 30 105 32 110 Proc glm ； Class a；

Model y＝a x／solution； 1smeans a； means a；

得到主要输出是 General Linear Models Procedure Dependent Variable: Y Sum of Mean Source DF Squares Square F Value Pr > F Model 3 2328.343765 776.114588 68.20 0.0001 Error 20 227.614568 11.380728 Corrected Total 23 2555.958333 R-Square C.V. Root MSE Y Mean 0.910947 3.658599 3.373534 92.20833

上表为一般线性模型的协方差分析表，误差是227.614568，有20个自由度。F检验表明具有高度显著的线性关系。 General Linear Models Procedure Dependent Variable: Y Source DF Type I SS Mean Square F Value Pr > F A 2 1317.583333 658.791667 57. 0.0001 X 1 1010.760432 1010.760432 88.81 0.0001 Source DF Type III SS Mean Square F Value Pr > F A 2 707.218765 353.609382 31.07 0.0001 X 1 1010.760432 1010.760432 88.81 0.0001

上表与前表结合给出误差、A因素、回归的自由度分别为20，2和1；3型平方和为227．614568、707.218765和1010.760432。A因素、回归的F值分别为31.07和88.81，都是高度显著的。 协方差分析表就是表5.14。

表5.11

方差来源 回归 因素A 误差 总和 General Linear Models Procedure Dependent Variable: Y T for H0: Pr > |T| Std Error of Parameter Estimate Parameter=0 Estimate INTERCEPT 35.93518188 B 5.47 0.0001 6.577140 A A1 12.79324180 B 3.75 0.0013 3.4047 A2 17.33559201 B 7.20 0.0001 2.40915113 A3 0.00000000 B . . . X 2.40156919 9.42 0.0001 0.283321 NOTE: The X'X matrix has been found to be singular and a generalized inverse was used to solve the normal equations. Estimates followed by the letter 'B' are biased, and are not unique estimators of the parameters.

平方和 1010.76 707.218 227.615 2555.96 自由度 1 2 20 23 平均平方和 1010.76 353.609 11.381 F Prob>F 0.0001 0.0001 88.8* 31.07 上表给出,i,b的估计，英语说明指出，用一般线性模型算出的,i的估计是用广义逆解出的，是有偏的。 Least Squares Means A Y LSMEAN A1 94.9586305 A2 99.5009807 A3 82.1653887

上表给出i的最小二乘估计，正如协方差分析第二种理论所指出，它是唯一的。 General Linear Models Procedure Level of --------------Y-------------- --------------X-------------- A N Mean SD Mean SD A1 8 81.7500000 8.34522960 13.7500000 2.12132034 A2 8 98.0000000 5.12695956 18.6250000 1.99553072 A3 8 96.8750000 8.99900788 25.3750000 4.06860805

上表给出y在各个水平的均值：81.75，98，96.875。它们也称为未修正均值，因为是未考虑协变量的影响而算出的。i的最小二乘估计称为修正的均值，因为是考虑协变量的影响而算出的，修正公式为 iyib(xix)。

当给出水平和协变量的值后，可用公式yyib(xxi)进行预报。

从修正均值的表中看出扣除猪的初始重量后，第二种饲料的平均增重是99.5；仍是最大，第一种饲料的平均增重是94.96次之；第三种饲料的增重是82.17，效果最差。当考虑到饲料价格因素时，这种修正的增重更有用处。

从表5.17可见，对于单因素的协方差分析，与不考虑协变量的方差分析相比较．因素A的自由度不变，误差的自由度减少了k(k是协变量个数)；回归部分具有k个自由度，这一规律对于多因素协方差分析是普遍适用的。

对于多个因素，考虑交互作用和不考虑交互作用的，多个协变量的协方差分析问题，SAS-GLM过程采用第二种理论，建立一般线性模型。

对于例5.6可用数学模型

yijkbxijkijijki1,2,3,4j1,2,3,4k1,2

i0,j0ijk~N(0,2)并采用下列程序

data fab; do b=1 to 4; do a=1 to 4; do rep=1 to 2; input x @; input y @;

output; end; end; end; cards;

49.0 71 49.2 73 49.8 73 49.8 75 49.9 76 49.8 73 49.7 75 49.8 73 49.5 72 49.3 73 49.9 76 49.8 74 50.2 79 50.1 77 49.4 73 49.4 72 49.7 75 49.5 73 50.1 78 50.0 77 49.7 74 50.0 75 49.5 70 49.6 71 49.9 77 49.7 75 49.6 74 49.3 74 49.5 74 49.2 73 49.0 69 48.9 69 ; proc glm; class a b;

model y=a b a*b x/solution; lsmeans a b a*b;

run;

SAS计算结果主要如下 The GLM Procedure Dependent Variable: y Sum of Source DF Squares Mean Square F Value Pr > F Model 16 1.0087500 10.25069 9.49 <.0001 Error 15 16.2100000 1.0806667 Corrected Total 31 180.2187500 R-Square Coeff Var Root MSE y Mean 0.9100 1.407771 1.039551 73.84375

上表为一般线性模型的协方差分析表，误差是16.21，有15个自由度。F检验表明具有高度显著的线性关系。

Source DF Type I SS Mean Square F Value Pr > F a 3 70.59375000 23.53125000 21.77 <.0001 b 3 8.59375000 2.858333 2.65 0.0866 a*b 9 79.53125000 8.83680556 8.18 0.0002 x 1 5.29000000 5.29000000 4.90 0.0429 Source DF Type III SS Mean Square F Value Pr > F a 3 12.94841274 4.31613758 3.99 0.0283 b 3 2.32408298 0.77469433 0.72 0.5572 a*b 9 18.57792699 2.021411 1.91 0.1287 x 1 5.29000000 5.29000000 4.90 0.0429 上表与前表结合给出误差、A,B因素、A*B、回归的自由度分别为15、3、3、9和1；3型平方和为16.21、12.94841274、2.32408298、18.57792699、5.29。A,B因素、A*B、回归的F值分别为3.99、0.72、1.91、4.9。A因素、回归的影响是显著的,B因素、A*B的影响是不显著的。协方差分析表就是表5.11。

方差来源 回归 因素A 因素B 交互作 用A*B 误差 总和

Standard Parameter Estimate Error t Value Pr > |t| Intercept -156.1700000 B 101.7747156 -1.53 0.1457 a 1 3.0900000 B 2.05031 1.51 0.1526 a 2 2.7000000 B 1.4701474 1.84 0.0862 a 3 2.6600000 B 1.3312751 2.00 0.02 a 4 0.0000000 B . . . b 1 1.3200000 B 1.9614213 0.67 0.5112 b 2 1.4300000 B 1.3985731 1.02 0.3228 b 3 -1.2600000 B 1.6238309 -0.78 0.4498 b 4 0.0000000 B . . . a*b 1 1 -2.1000000 B 3.4478012 -0.61 0.5516 a*b 1 2 -3.0900000 B 2.2987953 -1.34 0.19 a*b 1 3 0.1800000 B 2.2198739 0.08 0.93 a*b 1 4 0.0000000 B . . . a*b 2 1 -2.9300000 B 1.7426053 -1.68 0.1134 a*b 2 2 -2.2700000 B 1.4738182 -1. 0.1443 a*b 2 3 2.0000000 B 1.4701474 1.36 0.1938 a*b 2 4 0.0000000 B . . . 平方和 5.29 12.94841274 2.32408298 18.57792699 16.21 180.2187500 自由度 1 3 3 9 15 31 平均平方和 5.29 4.31613758 0.77469433 2.021411 1.0806667 4.9 3.99 0.72 1.91 F Prob>F 0.0429 0.0283 0.5572 0.1287 a*b 3 1 -2.6200000 B 1.5969888 -1. 0.1217 a*b 3 2 -0.6100000 B 1.03841 -0.37 0.7152 a*b 3 3 -0.0400000 B 1.4847761 -0.03 0.97 a*b 3 4 0.0000000 B . . . a*b 4 1 0.0000000 B . . . a*b 4 2 0.0000000 B . . . a*b 4 3 0.0000000 B . . . a*b 4 4 0.0000000 B . . . x 4.6000000 2.0791024 2.21 0.0429 上表为参数的估计值，其中只有x的系数b=4.6是确定的。