2.方差分析
对于n>20样本,当原H0为真时,统计量T=T+-T-接近于0,其方差为
建立检验统计量
近似于标准正态分布。
由于T=T+-T-=2T+-n(n+1)/2,故可将上式中的T改写为T+的形式:
标准正态分布使用显著水平α=0.05时,拒绝区域为z<-1.96和z>1.96,因为2.24>1.96,计算出z统计量的值,判断拒绝H0与否。
三、SAS实现(PROCUNIVARIATE过程步)
例1检验提高学生某种素质的训练是否有效。
随机地选取15名学生作为试验样本,在训练开始前做了一次测验,每个学生的素质按优、良、中、及、差打分,经过三个月训练后,再做一次测试对每个学生打分(素质提高用+表示,降低用-表示,无变化用0表示)。
表1训练前后的素质比较
学生编号
训练之前
训练之后
差异符号
1
中
优
+
2
及格
良
+
3
良
中
-
4
差
中
+
5
良
良
0
6
中
优
+
7
差
及格
+
8
良
优
+
9
中
差
-
10
差
中
+
11
中
优
+
12
及格
良
+
13
中
及格
-
14
中
优
+
15
差
中
+
为了处理,先对定性资料进行量化:
用1,2,3,4,5,分布表示差、及格、中、良、优。
代码:
datatraining;
inputbeforeafter@@;
d=after-before;
datalines;
3524431344
3512453113
3524323513
;
run;
procprintdata=training;
title'原始数据';
run;
procunivariatedata=training;
vard;
run;
运行结果及说明:
注意:
只能调用univariate过程,而不能调用means过程来进行符号检验。
分析变量为单样本数据集training中的d变量。
符号检验统计量M(Sign)=4,它是取正符号和负符号两者之间的小者作为检验统计量(?
)
Pr>=|M|计算的概率是二项分布的两尾概率之和,因此它是双侧检验,检验正符号和负符号是否相同,结果为0.0574。
在显著水平设定为0.1时,由于0.0574<α=0.1,拒绝原假设。
符号检验的缺点是丢失了差值d大小的信息,如果设定检验的显著水平为0.05,那么本例检验结果却由于0.0574>0.05,则变为不能拒绝原假设。
但是,如果用考虑差值d大小的信息的Wilcoxon符号秩检验,即SgnRank,由于0.0154<0.05,仍然得到拒绝原假设的检验结果。
例2某制造商想要比较两种不同的生产方法所花费的生产时间是否有差异。
随机选取了11个工人,每一个工人都分别随机地使用两种不同的生产方法来完成一项相同的任务。
任务完成时间的正差值表示生产方法1需要更多的时间,负差值表示生产方法2需要更多的时间。
表2两种不同生产方法完成任务的时间(分钟)
工人编号
n
生产方法M
差值D
绝对差值
秩次
R
符号秩次R
M1
M2
D=M1-M2
|D|
-
+
1
10.2
9.5
0.7
0.7
8
8
2
9.6
9.8
-0.2
0.2
2
2
3
9.2
8.8
0.4
0.4
3.5
3.5
4
10.6
10.1
0.5
0.5
5.5
5.5
5
9.9
10.3
-0.4
0.4
3.5
3.5
6
10.2
9.3
0.9
0.9
10
10
7
10.6
10.5
0.1
0.1
1
1
8
10.0
10.0
0
0
—
—
—
9
11.2
10.6
0.6
0.6
7
7
10
10.7
10.2
0.5
0.5
5.5
5.5
11
10.6
9.8
0.8
0.8
9
9
符号秩次总和T-=5.5,T+=49.5
5.5
49.5
代码:
datatime;
inputm1m2@@;
d=m1-m2;
datalines;
10.29.59.69.89.28.8
10.610.19.910.310.29.3
10.610.510.010.011.210.6
10.710.210.69.8
;
run;
procprintdata=time;
title'原始数据';
run;
procunivariatedata=timenormal;
vard;
run;
运行结果及说明:
“normal”选项,对差值作正态性检验。
差值D的正态性检验的结果为0.5339>0.05
配对资料如果其差值不是具体数字,只能用符号检验。
但如果差值有具体数字,而使用符号检验,相当于只利用了它的“+”、“-”,而对数字大小中所包含信息却未加利用。
此时,若符合正态分布则使用配对资料的t检验;若不符合正态分布则用Wilcoxon符号秩检验。
差值D的正态性检验的结果为0.5338>0.05,因此不能拒绝差值D具有正态性。
因为制造商拒绝相信差值D具有正态性,所以采用Wilcoxon符号秩检验。
Wilcoxon符号秩统计量S=22。
SAS建议在n≤20时,Pr>=|S|的概率由S的精确分布计算,而S的分布是尺度二项分布的卷积,所以精确结果为p值=0.0234<α=0.05,拒绝原假设H0,即两种不同的生产方法所花费的生产时间是有差异。
若
>20时,将符号秩统计量S标准化成自由度为
-1的t统计量来计算显著水平(注意跟前文的转换成标准正态分布略有不同),原因是当
较大时,t分布渐近标准正态分布。
另外,SAS系统在计算秩统计量S的方差时,用结值来修正方差。
p值=0.0194<α=0.05,也是拒绝原假设H0.
(三)Wilcoxon秩和检验
一、两样本的Wilcoxon秩和检验
Wilcoxon秩和检验,用来决定两个独立样本是否来自相同的或相等的总体。
如果这两个独立样本来自正态分布和具有方差齐性(相同方差),则可以采用t检验比较均值。
但当这两个条件都不能确定时,常用Wilcoxon秩和检验。
Wilcoxon秩和检验是基于样本数据秩和。
先将两样本看成是单一样本(混合样本)然后由小到大排列观察值统一编秩。
若“原假设H0:
两个独立样本来自相同的总体”为真,则小的、中等的、大的秩值大约均匀分布在两个样本中。
若“备择假设H1:
两个独立样本来自不相同的总体”为真,则其中一个样本有更多的小秩值,这样就会得到一个较小的秩和;另一个样本将会有更多的大秩值,会得到一个较大的秩和。
设两个独立样本为:
第一个x样本容量为n1,第二个y样本容量为n2,在容量为n=n1+n2的混合样本中,x样本的秩和为Wx,y样本的秩和为Wy,且有
定义
以x样本为例,若它们在混合样本中享有最小的n1个秩,则Wx取到最小值n1(n1+1)/2;同样Wy可能取的最小值为n2(n2+1)/2。
那么,Wx的最大取值等于混合样本的总秩和减去Wy的最小值,即
同样,Wy也同理。
所以W1和W2均为取值在0与
之间的变量。
当原假设H0为真,所有的xi和yi相当于从同一总体中抽得的独立随机样本,可看成一排n个球随机地指定n1个为x球另n2个为y球,共有
种可能(且是等可能的)。
基于这样分析,在原假设H0为真的条件下可求出W1和W2的概率分布为为样本大小为n1和n2的Mann-Whitney-Wilcoxon分布。
一个具有实际价值的方法是,对于每个样本中的观察数≥8的大样本来说,我们可以采用标准正态分布Z来近似检验。
由于W1的中心点为n1n2/2,故Wx中心点μ为
Wx的方差σ2为
若样本中存在结值,需要对方差做修正:
其中,τj为第j个结值的个数(结值的存在将使方差变小)。
标准化后Wx为
其中分子±0.5是为了对离散变量进行连续性修正,对于Wx-μ>0减0.5修正,对于Wx-μ<0加0.5修正。
二、PROCNPAR1WAY过程步(单因子非参数方差分析)
NPAR1WAY过程,是分析变量的秩,并计算几个基于经验分布的函数和通过一个单因子分类变量的响应变量确定的秩得分的统计量。
秩的得分计算有:
Wilcoxon得分、中位数得分、Savage得分和VanderWaerden得分等。
然后再由秩得分计算简单的线性秩统计量,由这个秩统计量可以检验一个变量的分布在不同组中是否具有相同的位置参数,或者在EDF检验下,检验这个变量分布在不同组中是否分布相同。
秩得分的统计量也可以先用procrank过程计算秩得分,然后用procanova过程分析这些秩得分而得到。
秩得分计算,用线性秩统计量:
其中Ri为第i个观察的秩,a(Ri)为秩得分,Ci是一个指示向量(由0和1组成),它表示了第i个观察所属的类,n是观察的总数。
下面介绍NPAR1WAY过程的四种不同的a(Ri)秩得分的计算:
(1)Wilcoxon得分
a(Ri)=Ri
它对Logistic分布的位置移动是局部最优的。
在计算两样本情况下的Wilcoxon秩和统计量时,过程对零假设下的渐近标准正态分布的Z统计量进行一个连续的±0.5校正。
(2)Median得分
又称为中位数得分。
当观察的秩大于中位点时,中位数得分为1,否则为0.对于双指数分布,中位数得分是局部最优。
(3)VanderWaerden得分
简称为VW的得分,是对正态分布的次序统计量的期望值的近似:
a(Ri)=F-1(Ri/(n+1))
其中F-1(x)是标准正态的累积分布函数的反函数,这个得分对正态分布是最优的。
(4)Savage得分
是指数分布的次序统计量的期望值,减去1使得得分以0为中心:
它在指数分布中比较尺度的不同性或在极值分布中的位置移动上是最优的。
基本语法:
PROCNPAR1WAYdata=数据集<可选项>;
BY变量;
CLASS变量;
EXACT统计量选项;
FREQ变量;
OUTPUT<统计量选项>;
VAR变量列表;
说明:
(1)可选项:
ANOVA——方差分析
CONOVER——协方差分析
D——运用Kolmogorov-Smirnov(D)统计量评分进行分析
KLOTZ——运用Klotz评分进行分析
MEDIAN——运用中位数评分进行分析
MOOD——运用Mood评分进行分析
SAVAGE——运用Savage评分进行分析(指数分布)
SCORES=DATA——以原始数据为评分值进行分析
ST——运用Siegel-Tukey评分进行分析
VW/NORMAL——运用VanderWaerden评分进行分析(通过应用反正态分布累积函数得到近似的正态得分)
WILCOXON——Kruskal-Wallis秩和检验
EDF——计算基于经验分布函数的统计量
(2)EXACT语句,对指定的统计量(选项)进行精确概率的计算。
例3某航空公司的CEO注意到飞离亚特兰大的飞机放弃预定座位的旅客人数在增加,他想知道,是否从亚特兰大起飞的飞机比从芝加哥起飞的飞机有更多的放弃预定座位的旅客。
获得一个从亚特兰大起飞的9次航班和从芝加哥起飞的8次航班上放弃预定座位的旅客人数样本。
表3放弃预定座位的旅客人数及统一秩值
航班
次数
亚特兰大(
组)
芝加哥(
组)
放弃人数
统一编秩
放弃人数
统一编秩
1
11
5.5
13
7
2
15
9
14
8
3
10
3.5
10
3.5
4
18
12
8
1
5
11
5.5
16
10
6
20
13
9
2
7
24
16
17
11
8
22
15
21
14
9
25
17
秩和
96.5
56.5
代码:
datanoshows;
dogroup=1to2;
inputn;
doi=1ton;
inputx@@;
output;
end;
end;
dropin;
datalines;
9
111510181120242225
8
13141081691721
;
run;
procprintdata=noshows;
title'原始数据';
run;
procnpar1waydata=noshowswilcoxon;
classgroup;
varx;
run;
运行结果及说明:
选项wilcoxon要求进行wilcoxon秩和检验。
要注意,若两组样本是配对样本,应该使用配对t检验或wilcoxon符号检验,因为使用wilcoxon秩和方法,将损失配对信息。
组1和组2的秩和分别为96.50和56.50。
原假设H0为真时(组1和组2的总体分布相同),期望秩值分别为
(96.50+56.50)×9/(9+8)=81.0
(96.50+56.50)×8/(9+8)=72.0
标准差为10.3795614,每组平均得分分别为
96.50/9=10.7222222
56.50/8=7.0625000
Wilcoxon两样本秩和统计量(较小的秩和)S=56.5000,正态近似检验统计量Z=-1.44515(连续性修正因子为0.5,加在分子上),正态分布的双尾p值之和为0.1484>α=0.05,不能拒绝原假设H0.
同时还给出了近似t检验和卡方检验的结果:
近似t检验的p值=0.1677,近似卡方检验统计量为2.2300,自由度为1,p值=0.1354。
结果都是相同的,不能拒绝原假设H0.
(四)完全随机设计的Kruskal-Wallis秩和检验
一、概述
方差分析,可以检验三个或更多总体的均值是否相等的问题,数据是被假设成具有正态分布和方差齐性(相等的方差),此时F检验才能奏效。
但有时数据不能完全满足这些条件,不妨将数据转换成秩统计量(秩统计量的分布与总体分布无关),可以摆脱总体分布的束缚。
在比较两个以上的总体时,广泛使用非参数的Kruskal-Wallis秩和检验,它是对两个以上的秩样本进行比较,本质上它是两样本时的Wilcoxon秩和检验方法在多于两个样本时的推广。
Kruskal-Wallis秩和检验,首先要求从总体中抽取的样本必须是对立的,然后将所有样本的值混合在一起看成是单一样本,再把这个单一的混合样本中值从小到大排序,序列值替换成秩值,最小的值给予秩值1,多个相同值时用平分秩值。
将数据样本转换成秩样本后,再对这个秩样本进行方差分析,但此时构造的统计量KW不是组间平均平方和除以组内平均平方和,而是组间平方和除以全体样本秩方差。
这个KW统计量是我们判定各组之间是否存在差异的有力依据。
二、基本原理
设有k组样本,ni是第i组样本中的观察数,n是所有样本中的观察总数,Ri·是第i组样本中的秩和,Rij是第i组样本中的第j个观察值的秩值。
需要检验的原假设H0为各组之间不存在差异,或者说各组的样本来自的总体具有相同的中心或均值或中位数。
在H0为真时,各组样本的秩平均应该与全体样本的秩平均比较接近:
所以组间平方和为
恰好是刻划该接近程度的一个统计量,除以全体样本秩方差消除量纲的影响。
样本方差的自由度为n-1,所以全体样本的秩方差为:
因此,Kruskal-Wallis秩和统计量KW为:
若样本中存在多个相同值(结值)则需要调整KW公式,校正系数C为:
其中,其中τj为第j个结值的个数。
调整后的KWc统计量为:
KWc=KW/C
如果每组样本中的观察数目至少有5个,那么样本统计量KWc非常接近自由度为k-1的卡方分布。
因此,我们将用卡方分布来决定KWc统计量的检验。
三、SAS实现(PROCNPAR1WAY过程步)
例4某制造商从来自3个大学的雇员中随机地抽取了3个独立样本,想知道来自这3个不同大学的雇员在管理岗位上的表现是否有所不同。
表4来自三个不同大学雇员得分及统一秩值
雇员
大学A
统一编秩
大学B
统一编秩
大学C
统一编秩
1
25
3
60
9
50
7
2
70
12
20
2
70
12
3
60
9
30
4
60
9
4
85
17
15
1
80
15.5
5
95
20
40
6
90
18.5
6
90
18.5
35
5
70
12
7
80
15.5
75
14
秩和
组A秩和
95
组B秩和
27
组C秩和
88
代码:
datacolleges;
dogroup=1to3;
inputn;
doi=1ton;
inputx@@;
output;
end;
end;
datalines;
7
25706085959080
6
602030154035
7
50706080907075
;
run;
procnpar1waydata=collegeswilcoxon;
classgroup;
varx;
run;
运行结果:
组1、组2和组3的秩和分别为95.0、27.0和88.0;原假设H0(组1、组2、组3的总体分布相同)为真时,期望秩值分别为
(95+27+88)×7/(7+6+7)=73.50
(95+27+88)×6/(7+6+7)=63.00
(95+27+88)×7/(7+6+7)=73.50
各组的标准差分别为12.5718985、12.0786894、12.5718985。
每组平均得分分别为95/7=13.5714286、27/6=4.50、88/7=12.5714286。
原假设H0:
各组的总体均值相等。
按修正公式修正后的多样本的Kruskal-Wallis秩和检验统计量为8.9839,用自由度为DF=3-1=2的卡方分布近似,得到大于近似卡方检验统计量8.9839的概率为p=0.0112<α=0.05,拒绝H0,表明各组的总体分布的差异是有统计学意义的。
根据平均秩和的结果,组1的最高,组2的最低,因此至少组1和组2的差异是显著的。
例5对于例4,也可以用freq过程,在tables语句中选项用scores=rank和cmh,查看第二项统计量即为Kruskal-Wallis检验。
代码:
procfreqdata=collegesformachar='|