SAS学习系列29. 方差分析Ⅱ—ANOVA,GLM过程步.docx

28. 方差分析Ⅱ—ANOVA,GLM过程步

SAS提供了ANOVA和GLM过程步进行方差分析。ANOVA过程步主要处理均衡数据（分类变量的每个水平的观察数是相等），该过程考虑到均衡设计的特殊构造，处理起来速度更快更省内存，也可以处理拉丁方设计、若干不完全的均衡区组设计数据等。
若试验设计不均衡，也不是前面几种实验设计数据，则应该使用GLM过程。

（一）PROC ANOVA过程步

一、基本语法
PROC ANOVA data=数据集 <可选项> ;
CLASS 分类变量列表;
MODEL 因变量=效应变量列表 </可选项>;
<MEANS 效应变量列表 </可选项> ;>
<TEST <H=效应变量列表> E=效应变量列表;>
说明：
（1）CLASS语句是必不可少的，必须放在MODEL语句之前，用来指定分类、区组变量（单因素方差分析只有一个变量）；
（2）MODEL语句也是必不可少的，该语句用来规定因变量和自变量效应（单因素方差分析的自变量就是分类变量）。若没有规定自变量的效应，则只拟合截距，假设检验为因变量的均值是否为0. Model语句的主要形式有4种：
① 主效应模型
model y=a b c;
② 含有交叉因素的模型
model y=a b c a*b a*c b*c a*b*c;
③ 嵌套模型
model y=a b c(a b);
④ 包含嵌套、交叉和主效应的模型
model y=a b(a) c(a) b*c(a);
（3）MEANS语句必须出现在MODEL语句之后，用来计算在效应变量所对应的因变量均值，但这些均值没有针对模型中的效应进行修正。若要计算修正的均值需要用GLM过程步的LSMEANS语句；
（4）MEANS语句的可选项主要有两个内容，一是选择多重比较的检验方法，二是设定这些检验的参数（只能用于主效应）；
bon——对所有主效应均值之差进行Bonferroni的t检验；
duncan——对所有主效应均值进行Duncan的多重极差检验；
smm|gt2——当样本量不等时，基于学生化最大模和Sidak不相关t不等式，等到Hochberg的GT2方法，对主效应均值进行两两对比检验；
snk——对所有主效应均值进行Student-Newman-Keuls的多重极差检验；
t|lsd——对所有主效应均值进行两两t检验，它相当于在单元观察数相等时Fisher的最小显著差检验；
tukey——对所有主效应均值进行Tukey的学生化极差检验；
waller——对所有主效应均值进行Waller-Duncan的k比率检验；
……
alpha=p——设置显著水平；
clm——对变量的每个水平的均值按置信区间形式输出；
e=效应变量——指定在多重对比检验中所使用的误差均方。默认使用残差均方。指定的效应变量必须是在model语句中出现过的；
kratio=值——给出Waller-Duncan检验的类型1/类型2的误差限制比例。Kratio的合理值为50、100（默认）、500，大约相当于两水平时alpha值为0.1、0.05、0.01.
hovtest——要求输出组间方差齐性的Levene检验；
……
（5）TEST语句指定效应变量（H=）和误差变量（E=）做F检验，误差变量必须要指定且只能指定1个效应变量。默认是用残差均方作为误差项对所有平方和（SS）计算F值。

例1 来自四个地区少女篮球队队员身高的数据，变量包括地区、身高（C:\MyRawData\GirlHeights.dat）：

读入数据，做单因素（身高）方差分析，判断她们的身高是否存在显著性地区差异。
代码：
data heights;
infile 'c:\MyRawData\GirlHeights.dat';
input Region $ Height @@;
run;
* Use ANOVA to run one-way analysis of variance;
proc anova data = heights;
class Region;
model Height = Region;
means Region / SCHEFFE;
title "Girls' Heights from Four Regions";
run;
运行结果及说明：

CLASS语句中分类变量有4个不同的水平值，共64个观测值；

因变量Height的方差分析表，因变量的总平方和（1030.000）、属于模型部分的平方和（196.625）、属于误差部分的平方和（833.375），自由度为（3,60,63），模型的均方MS（65.541667），误差的均方MSE（13.889583），F值=MS/MSE=4.72, P值=0.0051<α=0.05, 否定原假设，即不同地区Height的均值不全相等（至少有两个不相等）。
R2=196.625/1030.000=0.90898,
变异系数CV=6.134771=100×根MSE / Height均值 (%),
因变量的标准差（根MSE）为3.726873
效应变量Region的方差分析表，同因变量的方差分析表中“模型”行。下面是默认输出的盒形图：

Levene的方差齐性检验结果（0.4514>0.05）表明：不能拒绝不同地区身高的方差是相等的原假设。

MEANS语句中的SCHEFFE选项，比较不同区域的平均身高，Scheffe分组A、B，在显著水平α=0.05下，认为同组内身高没有差异。

例2 接 例1 四个地区的Height均值不同，但可能存在某2个或某3个或地区的身高均值相同。除了用SCHEFFE选项，还可以对均值做多重比较和置信区间分析。
代码（部分）：
proc anova data = heights;
class Region;
model Height = Region;
means Region / DUNCAN;
means Region / LSD CLM CLDIFF;
运行结果及说明：

DUNCAN选项，输出组间均值比较的多重极差检验，各组均值按从小到大排列，3个均值间的比较，就看3个地区最大和最小均值之差是否大于临界值2.773，North、East、West均值之差60.750-58.688=2.062<2.773, 故这三个地区均值没有显著差异（α=0.05）

各地区Height均值t检验的置信区间：均值±1.863714.

LSD最小显著差检验，0.05显著水平下，两两比较的最小显著差为2.6357，若显著则被标上“***”，例如，South与North均值之差为2.750>2.6357，故有显著差异。

（二）PROC GLM过程步

GLM过程步分析符合一般线性模型（General Linear Models）的数据，因此取名GLM。可用在简单回归、多元回归、方差分析、协方差分析、加权回归、多项式回归、偏相关分析、多元方差分析等。
GLM过程步的语法与ANOVA过程步基本相同。区别是GLM过程多了些MODEL模型，并可以多三条语句：contrast、estimate和lsmeans.
1. MODEL模型（a、b、c表示分类变量；y1、y2、x1、x2代表连续变量）：
Model y=x1; ——线性回归
Model y=x1 x2; ——多元线性回归
Model y=x1 x1*x1; ——多项式回归
Model y1 y2=x1 x2; ——多元回归
Model y=a; ——单因素方差分析
Model y=a b c; ——主效应模型
Model y=a b a*b; ——交叉因素模型
Model y=a b(a) c(b a); ——嵌套模型
Model y1 y2=a b; ——多元方差分析模型
Model y=a x1——协方差分析模型
……
e1/e2/e3/e4——输出模型中每一效应的类型1/类型2/类型3/类型4的可估函数，并计算相应的平方和；
ss1/ss2/ss3/ss4——对每个效应，输出与类型1/类型2/类型3/类型4的可估函数相关的平方和；
cli/clm——打印每一观察的预测值/预测均值的置信限，两者不能同时使用；
p——打印自变量没有缺失值的每一观察值、预测值、残差值，以及Durbin-Waston统计量；

2. contrast语句
用来检验均值的线性组合关系的原假设。有三个基本参数，一是标签，二是分类变量名，三是效应均值线性组合的系数表（系数的次序是匹配分类变量按字母数字次序的水平值）。
示例：
contrast 'US vs NON-U.S.' brand 2 2 2 -3 -3;
检验H0：2μ1+2μ2+2μ3-3μ4-3μ5=0

3. estimate语句
用来估计效应均值的线性组合的值，格式同contrast语句。
示例：（分数系数的表示）
estimate '1/3(a+b)－2/3c' Man 1 1 -2 /divisor=3;

4. lsmeans语句
用来计算效应变量修正后的均值，最小二乘均值（LSM），这是针对非均衡数据设计的。可选参数：
stderr——输出LSM的标准差和H0：LSM=0的概率值；
tdiff——输出假设检验H0：LSM(i)=LSM(j) 的t值和相应的概率值；
slice=效应变量——通过规定的这个效应来分开交叉的LSM效应。例如，假定交叉项A*B是显著的，如果想对B的每个效应检验A的效应，使用下面语句：
lsmeans A*B /slice=B;

例3考虑在5种不同品牌的人工合成胶合板材料上进行磨损时间测试，每种品牌的材料做四次试验，且都是采用的同一种磨损措施，所有的试验都是在完全随机的顺序下在相同的机器上完成的。

品牌ACMX、AXAX和CHAMP来自美国制造商，而品牌TUFFY和XTRA来自非美国制造商。我们想要比较美国品牌的均值与非美国品牌的均值是否有差异。
代码：
data veneer;
input brand $ wear @@;
datalines;
ACME 2.3 ACME 2.1 ACME 2.4 ACME 2.5
CHAMP 2.2 CHAMP 2.3 CHAMP 2.4 CHAMP 2.6
AJAX 2.2 AJAX 2.0 AJAX 1.9 AJAX 2.1
TUFFY 2.4 TUFFY 2.7 TUFFY 2.6 TUFFY 2.7
XTRA 2.3 XTRA 2.5 XTRA 2.3 XTRA 2.4
;
run;
proc glm data = veneer;
class brand;
model wear=brand;
contrast 'US vs NON-U.S.' brand 2 2 2 -3 -3;
estimate 'US vs NON-U.S.' brand 2 2 2 -3 -3;
title 'Wear Tests for five brands';
run;
运行结果：

程序说明：
（1）根据题意，原假设
H0: (μACME+μAJAX+μCHAMP)/3=(μTUFFY+μXTRA)/2
等价于H0: 2(μACME+μAJAX+μCHAMP)-3(μTUFFY+μXTRA)=0, 故contrast语句的系数表为2,2,2,-3,-3. （注意到均值对应关系是按字母顺序排列）；
（2）美国品牌均值与非美国品牌均值比较的平方和为0.27075，F值为13=0.27075/0.020833，P值=0.0026<α=0.05，拒绝原假设H0，说明美国品牌均值与非美国品牌均值是不同的；
（3）效应线性组合的参数估计为
-1.425=3×(2.325+2.050+2.375)-2×(2.600+2.375)
对于原假设H0参数是否为0的t检验，t值=-3.60，P值=0.0026<α=0.05，拒绝原假设（注意到t检验的p值与F检验的p值相同，这是因为两种检验是相同的，F值等于t值的平方）。

例4 （随机单位组试验设计的方差分析）
某食品公司对一种食品设计了四种包装。为了考察哪种包装最受欢迎，选了10个有近似相同销售量的商店作试验，其中两种包装各指定两个商店，另两种包装各指定三个商店销售。在试验期间各商店的货架排放位置、空间都尽量一致，营业员的促销方法也基本相同。观察在一定时期的销售量（数据见下表）。试比较四种包装的销售量是否一致。
表 四种包装在10个商店中的销售量
包装类型
（treat）
商店（block）
商店数
n
1
2
3
A1
12
18

2
A2
14
12
13
3
A3
19
17
21
3
A4
24
30

2

注意，包装类型A1和A4在商店3里没有进行试验，所以这是有不平衡数据集的随机区组设计。
代码：
data pack;
input treat $ n;
do block=1 to n;
input y @@;
output;
end;
datalines;
A1 2
12 18
A2 3
14 12 13
A3 3
19 17 21
A4 2
24 30
;
run;
proc print data = pack;
title 'Sales for Four Different Pack';
run;
proc glm data = pack;
class block treat;
model y = block treat;
means block treat / SNK;
means block treat / DUNNETT('1');
means block treat;
run;
运行结果及说明：

读入数据，用n商店数控制每次读入数据数目（output不能缺），并输出原始数据集。

有两个分组变量，一是包装类型treat，包含四个水平A1、A2、A3、A4；二是商店名block，包含三个水平1、2、3. 共10个观测。CLASS语句，指定分组变量：包装类型treat，商店名block.

总模型方差分析结果：P值=0.0515，基本上有显著意义；R2=0.884868=269/304, 模型变异基本反映了总变异。
对于单因素不平衡数据的方差分析，类型Ⅰ和类型Ⅲ的平方和就不相同了，分组变量的变异计算应该采用类型Ⅲ的平方和。
分组变量block的方差分析结果p=0.5789>α=0.05，不具有显著意义，说明食品在3家不同商店进行销售时，销售量的均值没有显著差异；分组变量treat的方差分析结果p=0.0256<α=0.05，具有显著意义，说明4种不同包装食品的销售量的均值具有显著差异，但没有指出具体哪几种包装之间有显著差异。

MEANS语句的snk选项，指定采用多极差检验法对均值进行多级比较。3个组比较时，大均值与小均值之差的临界值为8.607705，而2个组比较时，临界值为6.7057385. “SNK”分组结果表明：3个商店（2, 1, 3）标有相同字母“A”，说明了3个商店的销售量均值没有显著差异。

对treat组进行snk多极差检验， “SNK”分组结果显示，包装A3，A1，A2出现了标有相同的字母 “B”，没有显著差异，它们与包装A4有显著差异。若看任意两种包装的差异，例如，A4与A2为27－13=14>10.992537，有显著差异。
结论： A4包装的销售量均值最高，其他三种包装销售量基本相同。另外，区组观察数的调和均数为2.4=4/(1/2+1/3+1/3+1/2)。

DUNNETT (‘1’) 选项，要求所有分组均值分别与对照组均值进行比较，采用dunnett的双尾t检验；也可用dunnetl（单尾t检验，分组的均值是否显著地小于对照组的均值）或dunnetu（单尾t检验，分组的均值是否显著地大于对照组的均值）。对照组在括号内规定为‘1’，即分组变量的第1个水平分组，第1家商店和A1包装。
用 Dunnett双侧检验的t临界值为3.33563，A2组与A1组均值之差为2<2×3.33563，无显著意义；A3组与A1组均值之差为-0.25<2×3.33563，无显著意义；另外也输出了均值之差的置信限。

第三个MEANS语句，用来输出各个分组的均值和标准差。

例5 （双因素实验设计的方差分析）
研究饮食和健美操对减肥的作用。饮食对减肥肯定有一定作用，适当的健美操对减肥也有效果。那么哪一种饮食配上哪一样健美操最为有效呢？因为饮食与饮食这两种减肥手段之间存在着交互作用，会加强减肥的效果。现有三套饮食方案称为a、b、c，五种不同的健美操标记为1、2、3、4、5。构成成了3×5=15种水平组合，选择了情况基本相同的90个肥胖人进行试验，将他们随机地指派到这15个组中且每组6人。经过一段时间后，体重的下降结果如下表所示：
表 3×5双因素设计的试验结果
饮食方案
food
健美操train
1
2
3
4
5

a

22.1
24.1
19.1
22.1
25.1
18.1
27.1
15.1
20.6
28.6
15.1
24.6
22.3
25.8
22.8
28.3
21.3
18.3
19.8
28.3
26.8
27.3
26.8
26.8
20.0
17.0
24.0
22.5
28.0
22.5

b

13.5
14.5
11.5
6.0
27.0
18.0
16.9
17.4
10.4
19.4
11.9
15.4
15.7
10.2
16.7
19.7
18.2
12.2
15.1
6.5
17.1
7.6
13.6
21.1
21.8
22.8
18.8
21.3
16.3
14.3

c

19.0
22.0
20.0
14.5
19.0
16.0
20.0
22.0
25.5
16.5
18.0
17.5
16.4
14.4
21.4
19.9
10.4
21.4
24.5
16.0
11.0
7.5
14.5
15.5
11.8
14.3
21.3
6.3
7.8
13.8

代码：
data fatness;
do i=1 to 3;
Input food $ ;
do train=1 to 5;
do j=1 to 6;
input y @@;
output;
end;
end;
end;
datalines;
a
22.1 24.1 19.1 22.1 25.1 18.1
27.1 15.1 20.6 28.6 15.1 24.6
22.3 25.8 22.8 28.3 21.3 18.3
19.8 28.3 26.8 27.3 26.8 26.8
20.0 17.0 24.0 22.5 28.0 22.5
b
13.5 14.5 11.5 6.0 27.0 18.0
16.9 17.4 10.4 19.4 11.9 15.4
15.7 10.2 16.7 19.7 18.2 12.2
15.1 6.5 17.1 7.6 13.6 21.1
21.8 22.8 18.8 21.3 16.3 14.3
c
19.0 22.0 20.0 14.5 19.0 16.0
20.0 22.0 25.5 16.5 18.0 17.5
16.4 14.4 21.4 19.9 10.4 21.4
24.5 16.0 11.0 7.5 14.5 15.5
11.8 14.3 21.3 6.3 7.8 13.8
;
run;
proc print data = fatness;
title 'Weight-loss Programs Based on Food and Train';
run;
proc glm data = fatness;
class food train;
model y = food train food*train;
lsmeans food train food*train;
lsmeans food*train / SLICE = food SLICE = train;
Contrast 't1 vs t4 in f1' train 1 0 0 -1 0 food*train 1 0 0 -1 0;
Contrast 't2 vs t4 in f1' train 0 1 0 -1 0 food*train 0 1 0 -1 0;
Contrast 't3 vs t4 in f1' train 0 0 1 -1 0 food*train 0 0 1 -1 0;
Contrast 't4 vs t5 in f1' train 0 0 0 1 -1 food*train 0 0 0 1 -1;
Contrast 't2 vs t5 in f3' train 0 1 0 0 -1 food*train 0 0 0 0 0 0 0 0 0 0 0 1 0 0 -1 ;
run;
运行结果及说明：
原始数据集（部分）如下

共有两个因素food和train，故CLASS语句有这两个分组变量名。除了要考察这两个因素的主效应外，还要考察这两个因素的交互效应，表示为food*train，所以需要在MODEL语句的后面加上这个交互效应。
用LSMEANS语句替代MEANS语句的主要原因是，对于非均衡的试验数据需要计算最小二乘均值，它是一种调整后的均值。
第二个LSMEAN语句的作用，考虑到交叉项food*train是显著情况时，通过SLICE选项规定的food效应和train效应来分开交叉的food*train效应。
CONTRAST语句是作更进一步的对比，前四个CONTRAST语句是把因素food固定在第一个水平a上，然后对food因素有显著交互作用的train因素的某两个水平之间进行比较；最后一个CONTRAST语句是把因素food固定在第三个水平c上，对train因素的第二个水平均值和第五个水平均值进行比较。
要注意food*train交叉效应的参数化形式的规则为：先变右下标，即f1*t1, f1*t2, ……, f1*t5, f2*t1、……, f3*t5.

总的模型方差分析结果：F值=4.87，P值=0.0001，模型效应是显著的。模型中有两个主效应food和train及一个交互效应food*train，其中主效应food和交互效应food*train是显著的，而主效应train，F值=0.14，P值=0.9648，是不显著的。
基本结论：饮食控制和健美操对减肥是有作用的，3种不同的饮食控制方案对减肥效果是有区别的，而5种不同的健美操对减肥效果是没有区别的，同时饮食方案和健美操的不同组合对减肥效果也是有区别的。

由于主效应food是显著的，说明三种饮食方案对减肥的效果是不同的，再通过查看三种饮食方案减肥体重的最小二乘均值均值，可以得出a方案最好，c方案最差，且a方案和c方案的差异应该是显著的，至于a与b的比较及b与c比较，可以采用多重比较的方法进一步分析。
要知道交互效应food*train显著的具体原因，需要比较两因素在各种组合时的均值差异，以便寻找最好的组合方案。先作切片（slice）分析，分别固定food因素在三个水平a、b、c上，再对交互效应food*train的五种不同组合均值进行分析，其中在切片a、b上无显著性（0.4546>0.05，0.2546>0.05），而在c上有显著性（0.0414<0.05），即只有当选择c饮食方案时，选择不同的健美操才会存在减肥效果区别。
对train的五种水平作切片分析，结果都是显著的（0.0383、0.0341、0.0100、0.0001、0.0010），说明无论采用哪种健美操，选择不同的饮食方案对减肥效果都存在区别。

进一步的分析，把因素food固定在减肥效果最好的第一个水平a上，然后把train因素的每个水平分别与第四个水平进行比较，结果显示都是无显著性（0.1050、0.1119、0.2718、0.1599），与前面的切片分析是一致的。
最后把因素food固定在第三个水平c上，对train因素的最大水平均值t2和最小水平t5均值进行比较，结果显示是显著的，同样证实了前面的切片分析。
综上得到结论：最佳效果的减肥措施是选择a饮食方案搭配5种健美操中的任何一种都可以。

例6（拉丁方设计的方差分析）
研究5种不同的防护服（A、B、C、D、E）对脉搏数的影响。采用5×5拉丁方试验设计，选用5个受试者，在5个不同日期进行试验，在行、列与字母上分别安排3个因素（日期、受试者、防护服），得到试验结果数据如下表所示：
表 5×5拉丁方试验设计的数据
日期
date
受试者person
1
2
3
4
5
1
A 129.8
B 116.2
C 114.8
D 104.0
E 100.6
2
B 144.4
C 119.2
D 113.2
E 132.8
A 115.2
3
C 143.0
D 118.0
E 115.8
A 123.0
B 103.8
4
D 133.4
E 110.8
A 114.0
B 98.0
C 110.6
5
E 142.8
A 110.6
B 105.8
C 120.0
D 109.8

代码：
data pulse;
do date=1 to 5;
do person=1 to 5;
input cloth $ y @@;
output;
end;
end;
datalines;
A 129.8 B 116.2 C 114.8 D 104.0 E 100.6
B 144.4 C 119.2 D 113.2 E 132.8 A 115.2
C 143.0 D 118.0 E 115.8 A 123.0 B 103.8
D 133.4 E 110.8 A 114.0 B 98.0 C 110.6
E 142.8 A 110.6 B 105.8 C 120.0 D 109.8
;
run;
proc print data = pulse;
title 'Pulse for Different Protective Suit';
run;
proc anova data = pulse;
class date person cloth;
model y = date person cloth;
run;
运行结果及说明：

拉丁方设计应该用ANOVA过程步。实际上拉丁方试验设计是一种特殊类型的3个因素试验设计，其水平数必须相同，因此在CLASS语句中有3个分类变量名（date、person、cloth）。在3个水平交叉的单元上只有一次试验，且不存在3个分类变量的交互效应，所以在MODEL语句等号的右边也只有这3个分类变量名。
所作的三个原假设为：①各种防护服的平均脉搏数相同；②各个受试者的平均脉搏数相同；③不同日期的平均脉搏数相同。
结果分析：5种防护服对脉搏数影响的差异在统计学上无显著意义（0.3445>0.05）。而受试者之间的差异是有统计学意义（0.0001<0.05）,说明试验的差异主要来自受试者个体的情况，如身体状况、心理状况等，而与5种不同的防护服无关，基本上也与试验在哪一天发生无关。
若要进一步比较某个因素的任两个水平的平均脉数是否相同，可增加MEANS或CONTRAST语句。