第七章 方差分析和协方差分析
（医学统计之星）
上次更新日期：
方差分析和协方差分析在SAS系统中由SAS/STAT模块来完成，其中我们常用的有ANOVA过程和GLM过程。前者运算速度较快，但功能较为有限；后者运算速度较慢，但功能强大，我们做协方差分析时就要用到GLM过程。本章将首先介绍方差分析所用数据集的建立技巧，然后重点介绍这两个程序步。
其实，这里的速度快慢只是相对而言，SAS的处理速度是首屈一指的。举个例子，这个暑假我做了一个有6600条记录的，7因素的，交叉设计的方差分析（是不是已经有人喊头痛了？），我先是用SPSS FOR WIN95 来做，运行了大约10分钟才出结果。我又换用SAS FOR WIN95 来做，结果用了――秒！
§ 方差分析数据集的建立技巧
方差分析的数据集格式
统计分析所用的数据格式和我们在分析整理资料时所用的格式是不同的。一般来说，数据集中应至少有一个结果变量，用于记录不同处理因素水平下观察值的大小；至少有一个处理因素变量，用于记录处理因素的类型及其水平数。以单因素方差分析为例，就应有一个结果变量和一个处理因素变量；而两因素的方差分析应有一个结果变量和两个处理因素变量。
例 某职业病防治院对31名石棉矿工中的石棉肺患者、可疑患者及非患者进行了用力肺活量测定，请给出数据集的结构（卫统p44 例）。
解：数据集中应有两个变量，x和group。x记录肺活量的大小；group取值为1、2或3，分别代表石棉肺患者、可疑患者及非患者。
例 某厂医务室测定了10名氟作业工人工前、工中及工后4小时的尿氟浓度，请给出数据集的结构（卫统p46 例）。
解：数据集中应有三个变量，x、group和worker。x记录尿氟浓度；group取值为1、2或3，分别代表工前、工中及工后；worker取值为1到10，分别代表10名工人。
方差分析数据集的建立技巧
可见方差分析的数据集其变量取值有一定的规律，因此可以利用循环语句和判断语句来简化输入。
例 请建立例的数据集。
解：我们可以一个一个的输入变量group的值，但这里给出一种更酷的建立方法。
libname a 'c:\user';
data _1;
group=1;
if _n_>11 then
此处的_n_为SAS系统变量，表明读入的数据为第几条记录，由于石棉肺患者只有11人，因此第12条记录只能是可疑患者的记录，同理，从第20条记录开始的只能是非患者的记录。这里利用判断语句来完成从_n_到group的转换。
if _n_>20 then group=3;
else group=2;
;
input value @@;
cards;
... ...
;
proc print;
run;
例 请建立例的数据集。
解：此例中数据较有规律，各组的例数均相等，这可正是循环语句大显身手的时候。
libname a 'c:\user';
data _2;
do group=1 to 3;
这里的SAS语句非常类似于Foxbase语句，两个循环分别控制了group和worker的取值，这在输入数据时是常见的一种技巧。
do worker=1 to 10;
input value @@;
output;
两个变量只是循环变量，因此要用output语句写入数据集中。
end;
end;
cards;
... ...
请注意x变量和处理因素变量的对应关系，不要改变数据块的次序，否则建立的数据集是错的。
;
proc print;
run;
§ 引 例
例 请回答在例中三组石棉矿工的用力肺活量有无差别？
解：由于是三组做比较，因此要用方差分析法。在建立数据集后可用ANOVA过程或GLM过程进行分析，同时给出选项进行各组均数的两两比较（此处用SNK法），两过程的输出格式基本一致。对于单因素的方差分析，在比较之前一般要考虑方差齐性的问题。这里可以用GLM过程中的HOVTEST选项实现，实际上这也是在SAS中做方差齐性检验的唯一途径。
① 设定数据库环境：
LIBNAME A 'C:\USER';
② 数据步，建立数据集（略）：
③ ANOVA或GLM过程，进行方差分析并同时做两两比较（SNK法）：
PROC ANOVA DATA=_1;
PROC GLM DATA=_1;
CLASS GROUP;
CLASS GROUP;
MODEL VALUE=GROUP;
MODEL VALUE=GROUP /SS1;
MEAN GROUP /SNK;
MEAN GROUP /SNK HOVTEST;
RUN;
RUN;
这里的GLM过程中有两个选项，SS1要求只输出方差分析SS1模型的结果（默认输出SS1和SS3两种模型的结果），HOVTEST则要求进行方差齐性检验（该选项只能在SAS 以上的版本中实现）。
§ ANOVA过程
ANOVA过程是SAS系统中用于方差分析的几个程序步之一，它适用于平衡的数据。这意味着如果按所分析的变量将观测值分类，则每个格子内的观测值数目应相同。如果数据不满足这一要求，则ANOVA过程的分析结果并不可靠，分析者应改用GLM过程来处理数据。
语法格式
PROC ANOVA [DATA= <数据集名>
MANOVA
按多元分析的要求略去有任一缺失值的记录
OUTSTAT= <数据集名>] ;
指定统计结果输出的数据集名
CLASS <处理因素名列>;
必需，指定要分析的处理因素
MODEL <应变量名=处理因素名列> / [选项];
必需，给出分析用的方差分析模型
MEANS <变量名列> / [选项] ;
指定要两两比较的因素及比较方法
BY <变量名列>;
FREQ <变量名>;
MANOVA H= 效应 E= 效应 M= 公式...;
指定多元方差分析的选项
语法说明
Means语句的选项主要用来指定两两比较的方法和检验水准，主要有：
两两比较方法 BON（Bonferroni检验）、DUNCAN（新复极差法）、DUNNETT、GABRIEL、LSD（最小显著差法）、SCHEFFE、SIDAK 、SNK（q检验）、TUKEY 、WALLER
其它 ALPHA= p值 更改水准
结果解释
这里以例的数据为例，ANOVA过程的输出结果分两大部分，第一部分是分析变量的列表描述，第二部分是统计结果，先是所用方差分析模型的检验，然后是各个处理因素和交互作用的检验结果，最后则是两两比较的结果。
Analysis of Variance Procedure 方差分析过程
Class Level Information 处理因素取值情况
Class Levels Values
处理因素变量名 取值数 具体值
GROUP 3 1 2 3
WORKER 10 1 2 3 4 5 6 7 8 9 10
Number of observations in data set = 30 样本大小为30
--------------以下为方差分析的统计结果部分，首先是对方差分析所用模型的检验-------------
Analysis of Variance Procedure
Dependent Variable: VALUE 应变量名为VALUE
Source DF Sum of Squares Mean Square F Value Pr > F
变异来源 自由度 SS MS F p值
Model 11
Error 18
Corrected Total 29
R-Square . Root MSE 应变量名 Mean
所用模型的决定系数 变异系数 剩余标准差 应变量均数
-----------统计结果的第二部分，对各个具体处理因素及其交互作用的检验----------
Dependent Variable: VALUE
Source DF Anova SS Mean Square F Value Pr > F
变异来源 自由度 SS MS F p值
GROUP 2
WORKER 9
---------统计结果的第三部分，对GROUP变量做两两比较，只有MEAN语句中指定了才会出现------
Analysis of Variance Procedure
Student-Newman-Keuls test for variable: VALUE
注意这里我们用SNK法做例子，其余比较方法的输出结果与之类似
NOTE: This test controls the type I experimentwise error rate under the
complete null hypothesis but not under partial null hypotheses.
Alpha= df= 18 MSE=
α水准为，自由度为18，MS误差为
Number of Means 2 3
Critical Range
SNK法两两比较时的界值，两均数之差大于该界值时则两组有统计学差异
Means with the same letter are not significantly different.
SNK Grouping Mean N group(处理因素变量名)
A 10 2
A
B A 10 3
B
B 10 1
检验的分组结果 各组均数 例数 组别
注：为了便于理解，SAS将两两比较的结果直接用英文字母的形式标示出来。两两比较结果的最右侧是处理因素变量的取值，最左侧标以字母A、B、C等等，用以表示该处理组和其它组有无差异。如果两组有相同的字母（如2、3两组），则两者之间无差异；而如果两组间只有不同的字母，则表示两组间的差异有统计学意义。可以看出此次比较的结果为：1组和2组有区别，其余均两两无区别。
应用实例
例 现测定了某湖湖水不同季节氯化物的含量，问不同季节氯化物含量有无区别（卫统p226 题）？
解：
data _6;
do season=1 to 4;
input x @@;
output;
end;
cards;
...
注意输入次序
;
proc anova data=_6;
class season;
model x=season;
mean season/snk;
run;
在本例中我为了节省地方，将大部分数据块用省略号略去了，可在带习时我惊讶的发现：有的同学在上机时敲入的这个程序居然也带有省略号？！显然是我误导了大家，辜负了同志们的信任，在此仅致以深深的歉意。
你们到底要我怎么写才好 ．．．
§ GLM过程
GLM过程即广义线形模型(General Liner Model)过程，它使用最小二乘法对数据拟合广义线形模型。GLM过程中可以进行回归分析、方差分析、协方差分析、剂量――反应模型分析、多元方差分析和偏相关分析等等，其功能之强大可见一斑。
以上所提到的各种分析均可归入广义线形模型的范畴，因此均可用GLM过程来分析。
语法格式
非常走运，在我们所用到的范围里，GLM过程的语法结构和ANOVA过程完全相同，这可大大的方便了我们的学习。这里我们只解释协方差分析的做法：由于协变量并非我们研究的处理因素，因此在CLASS语句中可不能有它的位置；但是，另一方面协变量要影响结果变量的取值，因此必须在模型中引入，怎么办？只要在MODEL语句中将其写入即可（写在分析变量的首位）。瞧，就这么简单！
结果解释
这次我们仍然很走运，GLM过程的结果和ANOVA过程是雷同的――注意是雷同而不是相同。在处理因素、协变量及交互作用的检验结果部分，GLM过程会给出两种结果（Ⅰ类和Ⅲ类结果）。在我们用到的范围里，两种结果是完全相同的，所以用哪一种都可以。
应用实例
例 今用A、B两药治疗12名贫血病人，测得治疗一个月时的红细胞增加数（百万/mm3），问两药的治疗效果如何，两药同时使用的效果如何（医统第二版P92 例）？
解：该实验的设计为析因设计，在分析时要将交互作用考虑在内。程序如下：
data _5;
do b=1 to 0 by -1;
do a=1 to 0 by -1;
do tmpvar=1 to 3;
tmpvar为临时变量，使一个格子里有3个观测值
input x@@;
output;
end;
end;
end;
drop tmpvar;
删除临时变量tmpvar
cards;
;
proc glm data=_4;
class a b;
model x=a b a*b;
a*b表示a、b的交互作用
run;
例 请给出《医学统计学》第一版108页例的分析程序。
解：协方差分析的关键步骤是分析资料是否满足协方差分析的应用条件，程序如下：
data _1;
group=1;
if _n_>13 then group=2;
input y x@@;
cards;
...
...
;
proc gplot;
做散点图，比较两组的相关趋势
plot y*x=group ;
按group变量分组绘图
symbol1 i=rl v=star;
在散点图上加回归线
symbol2 i=rl v=plus;
proc reg;
做x、y的回归分析
model y=x;
by group;
proc glm data=_1;
class group;
model y= x group;
x为协变量
run;
