第3章方差分析精品课程

第三章方差分析（初稿）3.1 单因素方差分析的数学模型首先让我们看两个例子：例 3.1 设甲、乙、丙三块麦田的基本苗数（按面积大小抽取样本点数）得表 3.1，问三块地的基本苗数是否有显著差别？表 3.1 三块麦田的基本苗数甲 乙 丙21 20 2429 25 2224 25 2822 23 2525 29 2130 31 2627 24 .26 26 . 20 . 21 .例 3.2 为了研究淬火温度和等温温度对铣刀硬度的影响三种不同淬火温度和三种等温温度淬火，测得铣刀平均硬度如表 3.2，检验淬火温度及等温温度是否对硬度 y 有显著影响。表 3.2 淬火温度与等温温度对硬度的影响淬火温度等温温度B1 B2 B3A1 64 66 68A2 66 68 67A3 65 67 68这两个例子和以前回归分析的问题不同首先，它们都只考察某种因素(地块、淬火温度、等温温度)在一系列试验中对产品某个指标 (寿命、得率 )的影响是否显著，而不要求建立回归方程；其次，这些因素可以不是定量的(如地块) ，即或这些因素是定量 (如淬火温度与等温温度)但其数值也不作为回归中变量的观察值，而只是代表一种处理（试验方案） ；最后，当因素确定后，可以作反复的试验。这两个例子和以前均值假设检验也不同，均值假设检验不考虑因素问题，而方差分析要考虑。在许多科学研究中都遇到和这两个例子类似的问题。尤其是科学研究中常涉及许多因素，例如研究作物栽培时，要考虑播种期、品种、土质、施肥方式、灌溉方式对产量的影响；在化学反应中要观察原料成分、剂量、催化剂、温度、压力，搅拌速度等对得率的影响。这些因素中要选出影响大的，以进一步安排更细致的试验，而判断一个因素的影响“是否大”的主要方法就是方差分析。我们所考察的。影响产品指标的因素(如产地，温度等) 也称为因子，用大写字母 A,B,C表示，例 1 有一个因子（地块） ，例 2 有 2 个因子（如淬火温度与等温温度） 。因素所能处的状况，如甲、乙、丙；60，65，70，75，称为因子的水平，简称为水平。例 1 的因子有三个水平（甲、乙、丙） ，例 2 每个因子恰也有 3 个水平，水平常以表示。,.,21BA因子也可以看成是一种变量，其取值不是数，而是水平。例如“产地”是一个变量，它取的值是“北京” 、 “上海” 、 “南京”等。这种变量称为属性变量，定性变量或分类变量本节只讨论一个因子，即一个分类变量的方差分析单因子方差分析。方差分析的目的在于找出自变量与因变量之间的线性关系，或自变量对因变量的实验效果。这种实验效果可分为：主效果、交互效果、镶嵌效果。Qualitative Variable(自变量，又称独立变量、定性变量)，Classification Variable(分类变量，其数值多半是不连续的。Response Variable（因变量，又称反应变量，其数值则是连续的）一般地，假设因素 A 有 k 个水平： 。第 j 个水平做实验 次，得指标 y，的kA,.1 jn个数据 。例 3.1 中 。通常作如下假设：jnjnjjy,.21 6;02;83n（1） 同一个水平 下得到的观测值 ，是由于实验过程中各种偶然因素j jnjy,.1的干扰及测量误差所致，每次实验中这些偶然因素的总和称为实验误差，它们是方差相同的零均值正态随机变量；（2） 所有误差相互独立；（3） 由于水平的不同，可能会给 一个定量的确定性的影响，其大小是未知的。ijy于是我们建立单因子方差分析数学模型（3.1）它 们 相 互 独 立),0(,.21,.Nkjniyij jjii其中 相互独立， 。(4.1)式称为单因素方差分析的数学模型。ij ),(2ij判断这个因素的影响是否显著就是要检验假设：, 不全相等kH.:10 k.:1令 ， ，jnjjijny/)(nyji/容易证明 是 的最小二乘估计jyj作方差分解 kjniijTjyS12)(kjni jjijj yy1 2)(kjnijijjy12kjnijj12)(kjnijijj12kjjyn12)(并令，kjnijijjySE12)(kjjynSA12)(即: -总的误差平方和kjniijTj12)(组内差，反映试验误差影响的大小。kjnijijEjyS12)(-组间差，反映因素 A 的各个水平不同引起kjjjAyn12)(的误差，若 A 的水平引起的误差显著时， 就比较大，反之就比较小。AS则有： EATS分别称为组内差和组间差。组内差又称为误差，用以估计实验误差影响的大小；组间差反映因素 A 的水平不同引起的系统差异。若 A 的水平不同引起的系统差异 (即组间的差距)显著时，SSA 就比较大；反之，当 A 引起的系统差异不显著时， SSA 就比较小而SSE 主要是由试验误差引起的。SSA 由 k 个平方之和形成，但有一个恒等式约束，只有 k-1 个自由度；同理 SSE 有 n-k 个自由度。下列定理给出kjjyn10)(SSA/SSE 的分布定理： 对于所给的模型，若 成立时，kH.:10则 ),()/(nFknSFEA如果 F 的值超过临界值（通常取为 ） ，就否定 ；当 F 超过 时，就称为高05. 0H01.度显著。当 成立时，F 的值不应太大，若 F 的值大于临界值时，就kH.:10应否定 ，即认为 间存在显著差异。,一般总用方差分析表表示计算结果，其形式为：方差来源 平方和 自由度 均方 F 值 临界值因素 A ASK-1 )1/(kSAF 05.误差 EN-K E总和 TN-13.2 方差分析的计算方差分析的计算一般都很复杂，可用 SAS 软件计算。SAS 中有 GLM，ANOVA 和NESTED 等过程可用方差分析。其中 GLM 过程和 ANOVA 最常用。3.2.1 PROC ANOVAPROC ANOVA 执行对均衡数据的方差分析，其应变量取连续值。所谓均衡数据是指对每个类变量的组合有同样多的响应变量观察值数。如果数据的非均衡的就要用 PROC GLM,做方差分析。 ANOVA 的语法为： PROC ANOVA DATA= SAS-data-setMANOVAMULTIPASSOUTSTAT= 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 ，给出应变量和自变量。其后有若干选项。322 PROC GLMPROC 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 TOLERANCEE E1 E2 E3 E4SS1 SS2 SS3 SS4 ALPHA= pCLM CLI P XPX INVERSESINGULAR= 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第五类选项 可用来调整统计的精度：ZETAGLM 过程主要有四个语句： PROC GLM，CLASS，MODEL 和 LSMEANS 语句PROC GLM 语句用以调用 GLM 过程，有许多选项，其中 DATA选项用以说明GLM 过程所加工的数据集。CLASS 语句说明哪些变量是分类变量。方差分析中的因素都是分类变量，响应变量和协方差分析中的协变量不是分类变量。例如 class x z;就指示计算机把因子 x ,z 作为分类变量。MODEL 语句中有等号，等号前是响应变量.LSMEANS 语句用以求待估参数的最小二乘估计。323 计算实例例 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: ySum ofSource DF Squares Mean Square F Value Pr FModel 2 6.7666667 3.3833333 0.32 0.7314Error 21 223.7333333 10.6539683Corrected Total 23 230.5000000R-Square Coeff Var Root MSE y Mean0.029356 13.18805 3.264042 24.75000Source DF Type I SS Mean Square F Value Pr Farea 2 6.76666667 3.38333333 0.32 0.7314Source DF Type III SS Mean Square F Value Pr Farea 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 的方差分析表方差来源 平方和 自由度 均方 F 值 F 分布随机变量大于 1.14 的概率因素 area 6.76666667 2 3.38333333 0.32 0.7314误差 223.7333333 21 10.6539683总和 230.5000000 23表 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 2429 25 2224 25 2822 23 2525 29 2130 31 2627 24 .26 26 . 20 . 21 .;proc glm;class area;model y=area;run;例 33 本文件（ 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.03dok5 17.7 3dok5 24.8 3dok5 27.9 3dok5 25.2 3dok5 24.33dok4 17.0 3dok4 19.4 3dok4 9.1 3dok4 11.9 3dok4 15.83dok7 20.7 3dok7 21.0 3dok7 20.5 3dok7 18.8 3dok7 18.63dok13 14.3 3dok13 14.4 3dok13 11.8 3dok13 11.6 3dok13 14.2compos 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 过程采用采用哑变量方法，以便作方差分析和协方差分析。该过程的数学原理是建立一般线性模型（类似回归模型）：（3.2）mxbby.10并用最小二乘法求解。在方差分析问题中，引进的变量 都是示性变量，即只取 0 或 1 的变量。mx.1GLM 过程对每一因子的每一水平，通过 class 语句产生 1 个示性变量。对于例 4.1 的计算，一般线性模型是3210xbxby其中对于地块甲 =1，否则为 0；对于地块乙 =1，否则为 0；对于地块丙 =1，否则为 0。从而设计2 3x矩阵是1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 对于一般线性模型，由于设计矩阵不满秩， 的解不是唯一的，但残差平方和是唯一确定mb,.的，从而F检验可行。这从以下例3.1的一般线性模型计算可见。令 。为使 S 最小，分别对2)()(3)(2)(101( kkknk yxbxS3210,b偏导得、0)()(3)(2)(101 kkknk yxxb)( (1()(3)(2)(101 kkkknk xb0)( (2()(3)(2)(101 kkkknk yxb)( (3()(3)(2)(101 kkkknk xb注意到 是是示性函数，= +, 所以 的解不是唯一的。由、分别可得ix 3210,b10by20by3代入可得最小值 。231)(ynSjjSAS-GLM 过程的假设检验总要用到 3 型平方和(Type III SS)。3 型平方和定义如下：设方差分析问题有 m 个因子，不包含某个因子一般线性模型的残差平方和减去包含全体因子一般线性模型的残差平方和之差称为该因子的 3 型平方和。例 3.1 中 Type III SS 计算步骤如下：先求不含地块因子的最小值=223.73333332312)(01 )(injiiinjkk yybi再求求不含地块因子的最小值= 230.50000002)()(3)(2)(10241(mi kkkk yxbx 231)(yniiType III SS= =230.5-223.7333333=6.76666667231)(jiiinjyi 21)(ii参见复旦大学编，概率论，第二册，第一分册 29 页，定理 4(Cochran)，37 页定理5。3.3 多因子方差分析实际问题中往往要考虑多因素的影响，例如，检验药效时，要用多种剂量的药对不同品种动物试验观察其效果；化工反应要考虑反应温度、浓度、反应时间等对得率的影响。和单因子方差分析一样，需要检验这些因子的影响是否显著，这就是多因子方差分析问题。例 3.2 中的因子就是淬火温度及等温温度，是双因子分析问题。由于因素越多，公式越复杂，我们着重介绍两因素的方差分析，更多因素的方差分析的原理是一样的，用 SAS 计算的程序也是相同的。331 一般双因素无重复试验无交互作用的模型设有两个因子 A、B，各有 p、q 个水平， 时 y 的观测值为 。假设：),(jiBAijy（1） A 的第 I 个水平使 y 增加 ；0,1pi（2） B 的第 j 个水平使 y 增加 ；,1qjj（3） 试验误差是方差相等均值为零的正态随机变量；（4） 试验误差相互独立。则得到模型为： 且 相 互 独 立),0(,.1;,.211Nqjpiyijjji ijjij要检验因素 A、B 的影响是否显著，就相当于分别检验以下假设：， 不全为零；0.10pHpH.:1， 不全为零。:2qq2令 ， ，pijjy1./)(jijiy1./)(piijjy1)/(可以证明： 是 的最小二乘估计； 是 的最小二乘估计。作方差分解.iij.jpiqjijTjyS12)( piqj jiijij yyy1 2 )()()(piqj ijij12)(pii12.)(qjjp12.)(并令 ， ，piqj ijijyySE12)( iiySA12.)(qjjypSB12.)(其中 SSA 称为 A 的组间差，它由 p 个平方和组成，有一等式 pii1.0)(约束，自由度是 p 一 1，反映 A 引起的系统差异；SSB 称 B 的组间差，反映 B 引起的系统差异，自由度是 q 一 1；SSE 称为误差，反映各种随机因素引起的差异，由 pq 个平方和组成，包含(p+q-1) 个等式 qj ijij piyy1 ,.10)(pi ijij pjy1 ,.0)(等式 可用上述等式中前 p 个之和减去后 q-1 个之和而得)所以自pi iqi1)由度为(p-1(q-1)即有： qpji jiijEjjBpiiAqjiijT yySyS,1, 22.1., 2)(),(,)则有方差分解： EBATS其中 反映了 A 引起的系统误差， 反映了 B 引起的系统误差， 反映了其它随S ES机误差。我们有定理：定理 3。2： 当 为真时，有：01H)1(,1()1(/ qpFqpSFEAA当当 为真时，有：02H)1(,1()(1/ qpqFpSqFEBB可由 是否超过临界值，判断响应的假设是否成立。BA,上述结果可以类似地用方差分析表的形式给出。例 5.2（续） ，可用下列 SAS 程序data quench;do eqt=1 to 3;do qut=1 to 3;input hard ;output;end;end;cards;64 66 6866 68 6765 67 68;proc glm;class eqt qut;model hard=eqt qut;run;执行此程序得到的主要输出有Dependent Variable: hardSum ofSource DF Squares Mean Square F Value Pr FModel 4 13.11111111 3.27777778 4.21 0.0962Error 4 3.11111111 0.77777778Corrected Total 8 16.22222222R-Square Coeff Var Root MSE hard Mean0.808219 1.325084 0.881917 66.55556Source DF Type I SS Mean Square F Value Pr Feqt 2 1.55555556 0.77777778 1.00 0.4444qut 2 11.55555556 5.77777778 7.43 0.0450Source DF Type III SS Mean Square F Value Pr Feqt 2 1.55555556 0.77777778 1.00 0.4444qut 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方差来源 平方和 自由度 平均平方和 F 临界值因素 A（等温温度）1.5556 2 0.7778 1.0 =6.9405.因素 B（淬火温度）11.5556 2 5.7778 7.43 =6.94.误差 3.1110 4 0.7778总和 16.2222 8由表 3.4 可见：不否定 ，否定 从而 B 因素的作用是显著的，A 因素的作用不显著，10H20即淬火温度对铣刀硬度有显著影响，而等温温度(即闷火温度) 没有显著影响模型式的优点是不需做重复试验，即对 A,B 不同水平，每一组合只需做一次试验，这大大节省人力、物力，但也有缺点，如例 4.3 所示。 对于多因子方差分析，GLM 过程对每一因子的每一水平，通过 class 语句产生 1 个示性变量。例 3.2 中有 6 个示性变量（3 个示性变量 表示等温温度，3 个示性变量21,A表示沾火温度） 。一般线性模型是321,B 3625143210 BbbAby由于设计矩阵不满秩， 解不唯一。但残差平方和是唯一确定的，从而610,.F 检验可行。 例 3.2 中 Type III SS 计算步骤如下：先求包含所有因子的最小值=3.11112)()(36)(25)(14)(3)(2)(1091(min kkkkkkk yBbbAbAb B再求不含等温温度因子的最小值=4.666672)()(36)(25)(14091(i kkkk yB等温温度因子的 3 型平方和就是Type III SS=4.66667-3.11111=1.55556然后求不含沾火温度因子的最小值=14.666672)()(3)(2)(1091(min kkkk yAbAb沾火温度因子的 3 型平方和就是Type III SS=14.66667-3.11111=11.55556332 双因素有重复试验有交互作用模型例 3.4 对三种推进器（ 、 、 ）和四种燃料（ 、 、 、 ）试验火1M231F234箭的射程，得表 5.5。推进器 M 与燃料 F 的作用是否显著?表 3.5 火箭的射程推进器燃料 123M1F58.2，52.6 56.2，41.2 65.3，60.8249.1，42.8 54.1，50.5 51.69，48.4360.1，58.3 70.9，73.2 39.2，40.74F75.8，71.5 58.2，51.0 48.7，41.4以 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.849.1 42.8 54.1 50.5 51.6 48.460.1 58.3 70.9 73.2 39.2 40.775.8 71.5 58.2 51.0 48.7 41.4;proc glm;class f m;model r=f m;run;从所得结果可见推进器和燃料引起的射程差异不显著，是否这些推进器和燃料对射程的影响不大呢? 并不是这样的，从表（4.5）中可见， 和 的组合射程很远； 和1M4F1M的组合射程很近， 不同水平的射程值作了平均后就不大不小了的缘故。因而模型3F1M(4.4)不能反映推进器和燃料的搭配对射程的影响。为了反映两种因素搭配后的影响，必须引入“交互作用”的概念。设有两个因子 A、B，各有 p、q 个水平，对 A、B 的每个组合 作 k 次重复试),(jiBA验（为了检验交互作用必须有重复实验） 。y 的观测值为。假设：ktjpiyjt ,.1;,.;,.1（1） 的第 I 个水平使 y 增加 ；0,1pi（2） B 的第 j 个水平使 y 增加 ；,1qjj（3） 的交互作用使 Y 增加 。),(jiA qjpiqjipijj ,.1;,.0,11 （4） 试验误差是方差相等均值为零的正态随机变量；（5） 试验误差相互独立。于是有数学模型： 且 相 互 独 立),0(,.1;,.1,0, ,;,.1211Nqjpiktqjpiyijtqjijpijjjii itjijjijt其中 分别称为因子 A、B 的主效应； 称为 A 与 B 的交互效应。当考虑,jiij交互作用时，方差分析问题就要检验， 不全为零；0.:10pHpH.:1， 不全为零；2qq2， 不全为零。.:130pp.:13qpjiijpiijjqjijiktijtij yyyy,1,1. 1.1,piqjijkktTjS12)( piqj jiijijijijkkt yyyy1 2 )()()()(piqj ijijk12)(piikq12.)(qjjy12.)(pijijijkkt y12)(并令 ， ，piqj ijijyykSAB12)( piiykqSA12.)(，qjj12.)(piqjijijkktSE12)(其中 SSA 称为 A 的组间差，反映 A 引起的发散程度，统计它们包含等式个数容易看出，自由度是 p 一 1；SSB 称 B 的组间差，反映 B 引起的发散，自由度是 q 一 1；SSAB 反映AB 交互作用效应引起的发散；自由度是(p-1)(q-1),SSE 称为误差，反映各种随机因素引起的发散，自由度是 pq(k-1)。qjjBpiiAkqptjiijtTykSyS12.12.,1, 2)(,)(kqptji ijijtEqpji jiijAByS y,1,12,1 2)( )则有方差分解： EABATSS其中 反映了 A 引起的系统误差， 反映了 B 引起的系统误差， 反映了 A 与 BS BS的交互作用， 反映了其它随机误差。我们有定理：E定理： 当 为真时，有：01H)1(,()1(/ kpqFkpqSFEAA当当 为真时，有：02H)1(,()1(/ kpqFkpqSFEBB当 为真时，有03H )1(),1()1(/ kpqpFkpqSFEAB可由 是否超过临界值，判断响应的假设是否成立。而且这 3 个检验是相互AB,独立的。所以当 之一大于临界值时，A,B 因素的主效应或交互作用效应显ABF,著。上述结果可以类似地用方差分析表的形式给出。表 5.8 方差分析方差来源平方和 自由度 平均平方和 F因子 A SSA p 一 1 SSA(p 一 1) )1(/kpqSEA因子 B SSB q 一 1 SSB(q 一 1) )(/B交互作用SSAB (p-1)(q-1) SSAB/(p-1(q-1) )1(/kpqSEA误差 SSE pq(k-1)总和 TSS Pq（k-1 ）用 SAS 软件计算交互作用时，只需在有关因子中加*号。例如 MODEL z=a b c a*b;指示计算机把 a,b,c 作为因子;考虑 a,b 的交互作用。例 3.4 的解：考虑交互作用的方差分析可用下述模型2,13,4,.1 tjiyijtijjiijt 其中 ，0i， ， ，j,.iij ,0jij且相互独立。),(2Nijt并用程序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.849.1 42.8 54.1 50.5 51.6 48.460.1 58.3 70.9 73.2 39.2 40.775.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 ProcedureDependent Variable: RSum of MeanSource DF Squares Square F Value Pr FModel 11 2401.348333 218.304394 11.06 0.0001Error 12 236.950000 19.745833Corrected Total 23 2638.298333R-Square C.V. Root MSE R Mean0.910188 8.080549 4.443628 54.99167从 Error行DF列可查得误差自由度为12，从Som of Squares列可查得误差（残差平方和）为236.95，从 Pr F的值为0.0001可见，一般线性模型的线性关系是高度显著的。Dependent Variable: RSource DF Type I SS Mean Square F Value Pr FF 3 261.675000 87.225000 4.42 0.0260M 2 370.980833 185.490417 9.39 0.0035F*M 6 1768.692500 294.782083 14.93 0.0001Source DF Type III SS Mean Square F Value Pr FF 3 261.675000 87.225000 4.42 0.0260M 2 370.980833 185.490417 9.39 0.0035F*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 值 PrFF(燃料) 3 261.657000 87.225000 4.42 0.0260M(发动机) 2 370.980933 185.49017 4.39 0.0035F*M(交互作用) 6 1768.692500 294.782086 14.93 0.0001SSE(误差) 12 236.9500 19.745833总和 23 2638.29833从方差分析表可见：F(燃料)对应的概率在 0.01 与 0.05 之间，说明燃料的作用显著；M( 发动机)及 F*M(交互作用)对应的概率小于 0.01，说明发动机与交互作用的影响高度显著。Least Squares MeansF RLSMEAN1 55.71666672 49.41666673 57.06666674 57.7666667上表给出不同燃料的射程均值（的最小二乘估计）。Least Squares MeansM RLSMEAN1 58.55000002 56.91250003 49.5125000上表给出不同发动机的射程均值（的最小二乘估计）。Least Squares MeansF M RLSMEAN1 1 55.40000001 2 48.70000001 3 63.05000002 1 45.95000002 2 52.30000002 3 50.00000003 1 59.20000003 2 72.05000003 3 39.95000004 1 73.65000004 2 54.60000004 3 45.0500000上表给出不同燃料与发动机组合的射程均值（的最小二乘估计）。由表可见，最远的射程是由第四种燃料与第一种发动机组合形成的，平均射程为 7365 公里。注意：考虑交互作用时，不要写赋值语句333 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。例 35title randomized complete block;data rcb;input block trtment $ yield worth;cards;1 a 32.6 1121 b 36.4 1301 c 29.5 1062 a 42.7 1392 b 47.1 1432 c 32.9 1123 a 35.3 1243 b 40.1 1343 c 33.6 116proc anova;class block trtment;model yield worth = block trtment;run;例 36*-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.0141 39.5112 37.9111 35.4121 36.7122 38.2132 36.4131 34.8221 42.7222 41.6212 40.3211 41.6241 44.5242 47.6231 43.6232 42.8proc anova;class block a b;model response