系统和数据分析Wilcoxon秩和检验

第二十八课Wilcoxon秩和检验两样本的Wilcoxon秩和检验两样本的 Wilcox on秩和检验是由 Mann , Whit ney和Wilcox on三人共同设计的一种检验， 有时也称为 Wilcoxon秩和检验，用来决定两个独立样本是否来自相同的或相等的总体。如果 这两个独立样本来自正态分布和具有相同方差时，我们可以采用t检验比较均值。但当这两个条件都不能确定时，我们常替换t检验法为 Wilcox on秩和检验。Wilcox on秩和检验是基于样本数据秩和。先将两样本看成是单一样本(混合样本)然后 由小到大排列观察值统一编秩。如果原假设两个独立样本来自相同的总体为真，那么秩将大 约均匀分布在两个样本中，即小的、中等的、大的秩值应该大约被均匀分在两个样本中。如 果备选假设两个独立样本来自不相同的总体为真，那么其中一个样本将会有更多的小秩值， 这样就会得到一个较小的秩和；另一个样本将会有更多的大秩值，因此就会得到一个较大的 秩和。设两个独立样本为：第一个的样本容量为，第二个样本容量为，在容量为n =nj n2的混合样本(第一个和第二个)中，样本的秩和为，样本的秩和为，且有:Wx Wy =12 n =n(n 1)2(28.1)我们定义:(28.2)(28.3)n 1 (n1 +1)Wi =Wx _2n 2( n2 +1)W2 =Wy2 2以样本为例，若它们在混合样本中享有最小的个秩，于是Wx二n1(n11，也是可能取去的最小值，即罟n2(n2 ;同样，的最大取值等于2n(n 1)n1(n11)2 2所以,式(28.2)和式(28.3)中的和均为取值在0与n(n 1)25( n1 1)2nzg1)2=n1 n2的变量。的最小值；同样可能取的最小值为(“21)。那么，的最大取值等于混合样本的总秩和减2当原假设为真时，所有的和相当于从同一总体中抽得的独立随机样本，和构成可分辨的排列情况，可看成一排个球随机地指定个为球，另个为球，共有C：2种可能，而且它们是等可能的。基于这样的分析，在原假设为真的条件下不难求出和的概率分布，显然它们的分布还是相同的，这个分布称为样本大小为和的Mann-Whit ney-Wilcoxo n 分布。一个具有实际价值的方法是，对于每个样本中的观察数大于等于8的大样本来说，我们可以采用标准正态分布来近似检验。由于的中心点为nm，根据式(28.2)，中心点为:2nin2nn 1)njn 1 + n2+1)R =2 2 2的方差从数学上可推导出：2nm 2 (n1 n2 1)a =12(28.4)(28.5)如果样本中存在结值，将影响到公式(28.5)中的方差，按结值调整方差的公式为：2 n1n2(n1 n2 1)厲门八(j ”j)CT =1212(n_j + n2)(n + n2 -1)其中，为第j个结值的个数。结值的存在将使原方差变小，这是一个显然正确的事实。(28.6)Wx -_0.5 z =CTn (m + n2 +1)Wx _ X 121 -0.52qn2(ni n2 1) n1n2 ( j -)1212( n1 + n2)( n1 + n2 1) N(0,1)(28.7)标准化后为：其中，分子加0.5或减0.5是为了对离散变量进行连续性修正，对于W I大于0减0.5修正，对于 Wx-小于0加0.5修正。例28.1某航空公司的CEO注意到飞离亚特兰大的飞机放弃预订座位的旅客人数在增加， 他特别有兴趣想知道，是否从亚特兰大起飞的飞机比从芝加哥起飞的飞机有更多的放弃预订 座位的旅客。获得一个从亚特兰大起飞的9次航班和从芝加哥起飞的8次航班上放弃预订座位的旅客人数样本，见表28.1中的第2列和第4列。表28.1 放弃预订座位的旅客人数及统一秩值航班次数亚特兰大(组)芝加哥(组)放弃人数统编秩放弃人数统编秩1115.513721591483103.5103.541812815115.51610620139272416171182215211492517秩和96.556.5如果假定放弃预订座位旅客人数的总体是正态分布且有相等的方差，我们可以采用两样 本比较的t检验。但航空公司的CEO认为这两个假设条件不能满足，因此采用非参数的Wilcox on秩和检验。将组与组看成是单一样本进行编秩，见表28.1中的第3列和第5列。最小值是8，秩值为1，最大值是25,秩值为17,有两个结值10和11,两个10平均分享秩值 3和4为3.5,两个11平均分享秩值5和6为5.5。如果两组放弃预订座位的旅客人数是相同 的，那么我们期望的两组秩和和大约是相同的；如果两组放弃预订座位的旅客人数是不相同 的，那么我们期望的两组秩和和也是非常不相同的。注意到n1 =9，n2 =8, =96.5，=56.5，H0 :两组放弃预订座位旅客人数的分布是相同的。标准正态分布值的计算结果为96.5，(9 80.5z： 2=1.445159(8)(9 8 1) _ 9(8)(8 -2 8 - 2) 彳 1212(9+8)(9+8 _1)如果设定显著水平=0.05，我们知道标准正态分布在0.05显著水平时，上临界值为1.645,下临界值为1.645，由于1.445 1.645，因此得到的是相同的结果，不能拒绝原假设。另外，要特别注意的是由于在连续型分布中随机地抽出个样本，几乎极少可能存在有些 值相等的情况，但在社会经济中有很多离散变量，很可能存在数值相同的情况，即样本中存 在着“结”。我们处理“结”的方法采用分享平均秩，但当大量“结”存在时，将可能直接影 响的方差，因此需要把式(28.5)中的方差修正为式(28.6)。但在手工计算和结值不多的情况下, 常使用未修正方差来简化计算，因为与修正方差的计算结果比较只存在一些小差异，大多数 情况下不影响最终的推断结果。二、单因子非参数方差分析的 nparlway过程单因子非参数方差分析的 npar1way过程是分析变量的秩， 并计算几个基于经验分布的函 数(EDF)和通过一个单因子分类变量的响应变量确定的秩得分的统计量。秩的得分计算分 成四种：Wilcox on得分、中位数得分、Savage得分和Van der Waerde n得分。然后，再由秩得 分计算简单的线性秩统计量，由这个秩统计量可以检验一个变量的分布在不同组中是否具有 相同的位置参数，或者在 EDF检验下，检验这个变量分布在不同组中是否分布相同。秩得分 的统计量也可以先用 proc rank过程计算秩得分，然后用proc anova过程分析这些秩得分而得到。1. 四种不同的秩得分计算用以下公式定义的统计量:nS 二 Ga(Ri)(28.8)i 4称为线性秩统计量。其中，是第个观察的秩，a(R)是秩得分，是一个指示向量(由 0和1组成)，它表示了第个观察所属的类，是观察的总数。nparlway过程的四种不同的 a(Ri)秩得分计算为：(1) Wilcoxon 得分在Wilcox on得分中：a(R)=(28.9)它对Logistic分布的位置移动是局部最优的。在计算两样本情况下的 Wilcox on秩和统计量时，过程对零假设下的渐进标准正态分布的z统计量进行一个连续的 +0.5和0.5校正。(2) Median 得分Median得分又称为中位数得分。当观察的秩大于中位点时，中位数得分为1,否则为0,即：a(Ri) =1当 Ri (n 1)/2(28.10)a(Rj =0当R 兰(n+1)/2对于双指数分布，中位数得分是局部最优。(3) Van der Waerden 得分Van der Waerden得分简称为 VW的得分。它是对正态分布的次序统计量的期望值的近似，即：a(Ri) = F lR/(n+1)(28.11)其中，F (x)函数是标准正态的累积分布函数的反函数，这个得分对正态分布是最优的。(4) Savage 得分Savage得分是指数分布的次序统计量的期望值。减去1使得得分以0为中心，即：Ra(R) =L 1/(n i +1) T(28.12)i=1Savage得分在指数分布中比较尺度的不同性或在极值分布中的位置移动上是最优的。2. nparlway过程说明proc npprocnpar1waydata=数据集 ;class分类变量；var 变量列表；by变量列表；为了使用proc nparlway过程，必须调用 proc和class语句。其余语句是供选择的。（1）proc nparlway 语句的选项anova对原始数据执行标准方差分析。edf计算基于经验分布函数（EDF ）的统计量，如Kolmogorov-Smirnov、Cramer-Von Meses、Kuiper 统计量。missing把class变量的缺失值看作一个有效的分类水平。median执行一个中位数得分分析。对于两样本产生一个中位数检验，对于更多 样本产生一个 Brown-Mood检验。savage执行一个Savage得分分析。该检验适用于数据服从指数分布的组间比较。vw执行一个Van der Waerden得分分析。这是一个通过应用反正态分布累积函数得到近似的正态得分。对于两个水平情况，这是一个标准Van der Waerden检验。wilcoxo n对数据或 Wilcoxo n得分进行秩分布。对于两个水平，它与Wilcox on秩和检验一样；对于任何数量的水平，这是一个Kruskal-Wallis检验。对于两样本情况，该过程使用一个连续的校正。（2） class语句class语句是必需的，它指定一个且只能一个分类变量。该变量用来标识数据中的各个类。Class语句变量可以是字符型或数值型。（3）var语句var语句命名要分析的响应变量或自变量。如果省略var语句，过程分析数据集中除 class语句指定的数据变量外的所有数值型变量。（4）by语句一个by语句能够用来得到由by变量定义的几个观察组，并用proc npar1way过程分别进行分析。当一个by语句出现时，过程希望输入的数据集已按by变量排序。三、实例分析例28.1的SAS程序如下：data study.noshows ;do group=1 to 2;input n;do i=1 to n;input x ;output;end;end;cards;11 15 10 18 11 20 24 22 25813 14 10 8 16 9 17 21proc npar1way data=study.noshows wilcoxon;class group;var x;run;程序说明：建立输入数据集 noshows,数据的输入和成组t检验相同，先输入本组数据的总数，然后输入组中每个数据。分组变量为group，共有两组取值为1和2。输入变量为x，存放每组中的数据。过程步调用npar1way过程，后面用选择项wilcoxon要求进行 wilcoxon秩和检验。要注意，如果两组样本是配对样本，应该使用配对t检验或wilcoxon符号检验，因为使用wilcoxon秩和方法，将损失配对信息。class语句后给出分组变量名group，var语句后给出要分析的变量X。主要结果如表 28.2所示。表28.2用nparlway过程进行 Wilcoxon秩和检验的输出结果N P A R 1 W A Y P R O C E D U R EWilcoxon Scores (Rank Sums) for Variable XClassified by Variable GROUPSum of Expected Std DevMeanGROUP N Scores Under HO Under HO Score1 996.500000081.010.379561410.72222222 856.500000072.010.37956147.0625000Average Scores Were Used for TiesWilcoxon 2-Sample Test (Normal Approximation)(with Continuity Correction of .5)S = 56.5000 Z = -1.44515 Prob |Z| = 0.1484第二十九课 结果说明：组1和组2的秩和 (Sum of Scores)分别为 96.50 和 56.50。原假设(组1和组2的总体分布相同)为真时，期望秩值(Expected)分别为(96.50+56.50) X 9/ (9+8) =81.0 和(96.50+56.50) X 8/ (9+8) =72.0,标准差(Std Dev)按公式(28.6)计算为10.3795614。每组平均得分（MeanScore） 分别为 96.50/9=10.7222222 和 56.50/ 8=7.0625000。Wilcoxon两样本秩和统计 量（较小的秩和）S =56.5000，正态近似 检验统计量Z = 1.44515 （连续性修正因 子为0.5，加在分子上），正态分布的双尾 p值之和为0.1484，不能拒绝原假设。同 时，还给出了近似t检验和卡方检验的结 果：近似t检验的p=0.1677，近似卡方检 验统计量为2.2300,自由度为1,p=0.1354。 结果都是相同的，不能拒绝原假设。完全随机设计Kruskal-Wallis秩和检验四、完全随机设计的Kruskal-Wallis秩和检验方差分析过程关注三个或更多总体的均值是否相等的问题，数据是被假设成具有正态分 布和相等的方差，此时 F检验才能奏效。但有时采集的数据常常不能完全满足这些条件。在 两两样本比较时，我们不妨尝试将数据转换成秩统计量，因为秩统计量的分布与总体分布无 关，可以摆脱总体分布的束缚。在比较两个以上的总体时，广泛使用非参数的Kruskal-Wallis秩和检验，它是对两个以上的秩样本进行比较，本质上它是两样本时的Wilcoxon秩和检验方法在多于两个样本时的推广。Kruskal-Wallis秩和检验，首先要求从总体中抽取的样本必须是独立的，然后将所有样本 的值混合在一起看成是单一样本，再把这个单一的混合样本中的值从小到大排序，序列值替 换成秩值，最小的值给予秩值1,有结值时平分秩值。将数据样本转换成秩样本后，再对这个秩样本进行方差分布，但此时我们构造的统计量KW不是组间平均平方和除以组内平均平方和，而是组间平方和除以全体样本秩方差。这个KW统计量是我们判定各组之间是否存在差异的有力依据。设有组样本，是第组样本中的观察数，是所有样本中的观察总数，是第组样本中的秩和， 是第组样本中的第个观察值的秩值。需要检验的原假设为各组之间不存在差异，或者说各组 的样本来自的总体具有相同的中心或均值或中位数。在原假设为真时，各组样本的秩平均应1+2+ nn+1该与全体样本的秩平均 二-一-比较接近。所以组间平方和为组间平方和八nii 4Rinin+1 22(29.1)恰好是刻画这种接近程度的一个统计量，除以全体样本秩方差，可以消除量纲的影响。样本 方差的自由度为n 一1。所以，全体样本的秩方差为全体样本的秩方差=+JRn -1 i42n 1I21n -1n fZ i-i 二 I厂nl 2Z i2n(n +1)2 4(29.2)1n -1n(n +1)(2n + 1)n(n +1)2、n(n 1)12因此，Kruskal-Wallis秩和统计量 KW 为:KW=组间平方和全体样本的秩方差kZ ni n(n 1) y1212 k n(n 1) i mni2Ri 2 -3( n 1)(29.3)如果样本中存在结值，需要调整公式(433)中的KW统计量，校正系数 C为:C=1(3，)(29.4)3n - n其中，第j个结值的个数。调整后的 KWc统计量为:KWc = KW /C(29.5)如果每组样本中的观察数目至少有5个，那么样本统计量 KWc非常接近自由度为 k - 1的卡方分布。因此，我们将用卡方分布来决定KWc统计量的检验。例29.1某制造商雇用了来自三个本地大学的雇员作为管理人员。最近，公司的人事部门 已经收集信息并考核了年度工作成绩。从三个大学来的雇员中随机地抽取了三个独立样本， 见表29.1中的第2、4、6列所示。制造商想知道来自这三个不同的大学的雇员在管理岗位上 的表现是否有所不同。表29.1来自三个不同大学的雇员得分及统一秩值雇员大学A统编秩大学B统编秩大学C 1统编秩12536095072701220270123609304609485171518015.5595204069018.569018.5355701278015.57514秩和组A秩和95组B秩和27组C秩和88为了计算KW统计量，我们首先必须将来自三个大学的20名雇员统一按考核成绩编排秩值，见表 29.1 中的第 3、5、7 列所示。本例中，n!=7,n2=6,n3=7，n = 20，Rj.=95,只2厂27，R3.=88，k =3，Ho :三个总体的考核成绩分布是相同的。我们用(29.3)式计算KW统计量为:KW 1295- -3(20 1) =8.916320(21) 767 _用(29.4)式计算校正系数 C，从表29.1中我们可以发现，相等成绩值和相等的个数分别为60分3个，70分3个，80分2个，90分2个。所以：3333203 -20C -(3 -3 3 -3 2 -22 -2) =0.9925调整后的KWc为:KWc 二 KW/C =8.9163/0.9925 = 8.9839查表可知道，自由度为k-1=2的卡方分布，在=0.05显著水平下，分布的上尾临界值为5.99，由于8.985.99，因此拒绝原假设。所以，秩和最低的B组至少与秩和最高的A组是不同的。五、freq频数过程Freq频数过程可以生成单向和 n向的频率表和交叉表。对于双向表(二维表)，该过程计 算检验统计量和关联度。对于n向表，该过程进行分层分析，计算每一层和交叉层的统计量。这些频数也能够输出到 SAS数据集里。1. freq过程说明proc freq 过程一般由下列语句控制：procfreqdata=数据集 ;by变量列表；tables交叉表的表达式 ;weight 变量;output ;run ;该过程proc freq语句是必需的。其余语句是供选择的。另外，该过程只能使用一个output语句。2. proc freq语句的选项order=freq/data/internal/formatted 规定变量水平的排列次序。 freq表示按频数下 降的次序，data表示按输入数据集中出现的次序，internal表示按非格式化值的次序(缺省值)，formatted按格式化值的次序。formachar (1,2,7)=三个字符规定用来构造列联表的轮廓线和分隔线的字符。缺省值为formachar(1,2,7)= +,第一个字符用来表示垂直线，第二个字符用来表示水 平线，第三个字符用来表示水平与垂直的交叉线。page要求freq每页只输出一张表。否则，按每页行数允许的空间输出几张表。 noprint禁止freq过程产生所有输出。3. by语句一个by语句能够用来得到由 by变量定义的分组观察，并分别进行分析。过程要求输入 的数据集已按by变量排序。4. tables 语句可以包括多个tables语句。如果没有tables语句，对数据集中的每个变量都生成一个单向 频数表。如果tables语句没有选项，则计算tables语句中规定变量每个水平的频数、累计频数、占总频数的百分比及累计百分比。Tables语句中的交叉表的表达式，请参见第二章第二节proc tabulate过程中的table语句的用法。Tables语句中的主要选项如下：all要求计算所有选项的检验和度量，包括chisq、measures和cmh。chisq 要求对每层是否齐性或独立性进行卡方检验，包括pearson卡方、似然比卡方和Mantel-Haenszel卡方。并计算依赖于卡方统计量的关联度，包括phi系数、列联系数和Cramer V。对于2X 2联列表还自动计算 Fisher的精确检验。cmh要求Cochran-Mantel-Haenszel卡方统计量，用于 2维以上表时，检验行变量和列变量是否有线性相关。exact要求对大于2X 2表计算Fisher的精确检验。Fisher的精确检验是假设行与列的边缘频数固定，并且在零假设为真时，各种可能的表的超几何概率之和。measures要求计算若干个有关相关的统计量及它们的渐近标准误差。alpha=p设定100 (1 p) %置信区间。缺省值为alpha= 0.05。scores=rank/table/ridit/modridit定义行/列得分的类型以便用于cmh统计量和pearson相关中。在非参数检验中，一般常用scores=rank,用于指定非参数分析的秩得分。cellchi2要求输出每个单元对总卡方统计量的贡献。cumcol要求在单元中输出累计列百分数。expected在独立性(或齐性)假设下，要求输出单元频数的期望值。deviation要求输出单元频数和期望值的偏差。misspri nt要求所有频数表输出缺失值的频数。missing要求把缺失值当作非缺失值看待，在计算百分数及其他统计量时包括它们。sparse要求输出在制表要求中变量水平的所有可能组合的信息。list以表格形式打印二维表。nocum/norow/nocol/nofreq/noprint分别不输出累计频率数、行百分率、列百分率、单元频数、频数表。5. weight 语句通常每个观察对频数计数的贡献都是1。然而当使用 weight语句时，每个观察对频数计数的贡献为这个观察对应的权数变量的值。6. output 语句该语句用于创建一个包含由 proc freq过程计算的统计量的 SAS数据集。由output语句创 建的数据集可以包括在 tables语句中要求的任意统计量。当有多个tables语句时，output语句 创建的数据集的内容相应于最后要求的那个表。六、实例分析例29.1的SAS程序如下：data study.colleges ;do group=1 to 3;input n;do i=1 to n;input x ;output;end;end;cards;725 70 60 85 95 90 80660 20 30 15 40 35750 70 60 80 90 70 75proc npar1way data=study.colleges wilcoxon;class group;var x;run;程序说明：建立输入数据集colleges，数据的输入和完全随机化方差分析的数据输入完全相同，先输入本组数据的总数，然后输入组中每个数据。分组变量为group，共有三组取值为1、2和3。输入变量为x,存放每组中的数据。过程步调用npar1way过程，后面用选择项wilcoxon，当样本数大于两个时，自动进行多样本的Kruskal-Wallis秩和检验。class语句后给出分组变量名group，var语句后给出要分析的变量x。主要结果如表29.2所示。表29.2 用npar1way过程进行多样本比较的Kruskal-Wallis秩和检验输出结果N P A R 1 W A Y P R O C E D U R EWilcoxon Scores (Rank Sums) for Variable XClassified by Variable GROUPSum of Expected Std DevMeanGROUP N Scores Under HO Under HO Score1795.073.500000012.571898513.57142862627.063.000000012.07868944.50000003788.073.500000012.571898512.5714286Average Scores Were Used for Ties结果说明：组1、组2和组3的秩和(Sum of Scores)分别为95.0、27.0和88.0。原假设 (组1、组2和组3的总体分布相同)为真时，期望秩值(Expected)分别为(95+27+88 )x7/( 7+6+7)=73.50、( 95+27+88)X 6/( 7+6+7) =63.00 和(95+27+88 ) X 7/ ( 7+6+7) =73.50, 各组的标准差(Std Dev)分别为 12.5718985、12.0786894、12.5718985。每组平均得分(MeanScore)分别为 95/7=13.5714286、27/6=4.50 和 88/7=12.5714286。按公式(29.5)调整后 多样本的Kruskal-Wallis秩和检验统计量为8.9839，用自由度为DF=3 仁2的卡方分布近似，得到大于近似卡方检验统计量8.9839的概率为p=0.01120.05，因此不能拒绝 group和x不线性相关。第二项cmh统计量为 ANOVA 统计量， 首先要求列变量 x是有序的。原假设为每一层的3个行的x平均得分是相等的，本例只有一层，且得分指定 scores=rank选项，即用秩得分方法，因此就是Kruskal-Wallis秩和检验统计量，自由度为行数减 1,即3-仁2,渐近自由度为2的卡方分布，KW= 8.984, p=0.0110.05，不能拒绝不相关。