SAS系统和数据分析全随机设计KruskalWallis秩和检验.doc

第二十九课 完全随机设计Kruskal-Wallis秩和检验一、 完全随机设计的Kruskal-Wallis秩和检验方差分析过程关注三个或更多总体的均值是否相等的问题，数据是被假设成具有正态分布和相等的方差，此时F检验才能奏效。但有时采集的数据常常不能完全满足这些条件。在两两样本比较时，我们不妨尝试将数据转换成秩统计量，因为秩统计量的分布与总体分布无关，可以摆脱总体分布的束缚。在比较两个以上的总体时，广泛使用非参数的Kruskal-Wallis秩和检验，它是对两个以上的秩样本进行比较，本质上它是两样本时的Wilcoxon秩和检验方法在多于两个样本时的推广。Kruskal-Wallis秩和检验，首先要求从总体中抽取的样本必须是独立的，然后将所有样本的值混合在一起看成是单一样本，再把这个单一的混合样本中的值从小到大排序，序列值替换成秩值，最小的值给予秩值1，有结值时平分秩值。将数据样本转换成秩样本后，再对这个秩样本进行方差分布，但此时我们构造的统计量KW不是组间平均平方和除以组内平均平方和，而是组间平方和除以全体样本秩方差。这个KW统计量是我们判定各组之间是否存在差异的有力依据。设有组样本，是第组样本中的观察数，是所有样本中的观察总数，是第组样本中的秩和，是第组样本中的第个观察值的秩值。需要检验的原假设为各组之间不存在差异，或者说各组的样本来自的总体具有相同的中心或均值或中位数。在原假设为真时，各组样本的秩平均应该与全体样本的秩平均比较接近。所以组间平方和为组间平方和(29.1)恰好是刻画这种接近程度的一个统计量，除以全体样本秩方差，可以消除量纲的影响。样本方差的自由度为。所以，全体样本的秩方差为全体样本的秩方差=(29.2)因此，Kruskal-Wallis秩和统计量KW为：(29.3)如果样本中存在结值，需要调整公式(4.3.3)中的KW统计量，校正系数C为：(29.4)其中，第j个结值的个数。调整后的KWc统计量为：(29.5)如果每组样本中的观察数目至少有5个，那么样本统计量KWc非常接近自由度为的卡方分布。因此，我们将用卡方分布来决定KWc统计量的检验。例29.1某制造商雇用了来自三个本地大学的雇员作为管理人员。最近，公司的人事部门已经收集信息并考核了年度工作成绩。从三个大学来的雇员中随机地抽取了三个独立样本，见表29.1中的第2、4、6列所示。制造商想知道来自这三个不同的大学的雇员在管理岗位上的表现是否有所不同。表29.1 来自三个不同大学的雇员得分及统一秩值雇员大学A统一编秩大学B统一编秩大学C统一编秩12536095072701220270123609304609485171518015.5595204069018.569018.5355701278015.57514秩和组A秩和95组B秩和27组C秩和88为了计算KW统计量，我们首先必须将来自三个大学的20名雇员统一按考核成绩编排秩值，见表29.1中的第3、5、7列所示。本例中，7，6，7，20，95，27，88，3，三个总体的考核成绩分布是相同的。我们用(29.3)式计算KW统计量为：用(29.4)式计算校正系数C，从表29.1中我们可以发现，相等成绩值和相等的个数分别为60分3个，70分3个，80分2个，90分2个。所以：调整后的KWc为：查表可知道，自由度为的卡方分布，在显著水平下，分布的上尾临界值为5.99，由于8.985.99，因此拒绝原假设。所以，秩和最低的B组至少与秩和最高的A组是不同的。二、 freq频数过程Freq频数过程可以生成单向和n向的频率表和交叉表。对于双向表（二维表），该过程计算检验统计量和关联度。对于n向表，该过程进行分层分析，计算每一层和交叉层的统计量。这些频数也能够输出到SAS数据集里。1. freq过程说明proc freq过程一般由下列语句控制：proc freq data=数据集 ;by 变量列表 ;tables 交叉表的表达式 ;weight 变量;output ;run ;该过程proc freq语句是必需的。其余语句是供选择的。另外，该过程只能使用一个output语句。2. proc freq语句的选项l order=freq/data/internal/formatted规定变量水平的排列次序。freq表示按频数下降的次序，data表示按输入数据集中出现的次序，internal表示按非格式化值的次序（缺省值），formatted按格式化值的次序。l formachar(1,2,7)=三个字符规定用来构造列联表的轮廓线和分隔线的字符。缺省值为formachar(1,2,7)=|+，第一个字符用来表示垂直线，第二个字符用来表示水平线，第三个字符用来表示水平与垂直的交叉线。l page要求freq每页只输出一张表。否则，按每页行数允许的空间输出几张表。l noprint禁止freq过程产生所有输出。3. by语句一个by语句能够用来得到由by变量定义的分组观察，并分别进行分析。过程要求输入的数据集已按by变量排序。4. tables语句可以包括多个tables语句。如果没有tables语句，对数据集中的每个变量都生成一个单向频数表。如果tables语句没有选项，则计算tables语句中规定变量每个水平的频数、累计频数、占总频数的百分比及累计百分比。Tables语句中的交叉表的表达式，请参见第二章第二节proc tabulate过程中的table语句的用法。Tables语句中的主要选项如下：l all要求计算所有选项的检验和度量，包括chisq、measures和cmh。l chisq要求对每层是否齐性或独立性进行卡方检验，包括pearson卡方、似然比卡方和Mantel-Haenszel卡方。并计算依赖于卡方统计量的关联度，包括phi系数、列联系数和Cramer V。对于22联列表还自动计算Fisher的精确检验。l cmh要求Cochran-Mantel-Haenszel卡方统计量，用于2维以上表时，检验行变量和列变量是否有线性相关。exact要求对大于22表计算Fisher的精确检验。Fisher的精确检验是假设行与列的边缘频数固定，并且在零假设为真时，各种可能的表的超几何概率之和。l measures要求计算若干个有关相关的统计量及它们的渐近标准误差。alpha=p设定100（1p）%置信区间。缺省值为alpha=0.05。l scores=rank/table/ridit/modridit定义行/列得分的类型以便用于cmh统计量和pearson相关中。在非参数检验中，一般常用scores=rank，用于指定非参数分析的秩得分。cellchi2要求输出每个单元对总卡方统计量的贡献。cumcol要求在单元中输出累计列百分数。l expected在独立性（或齐性）假设下，要求输出单元频数的期望值。l deviation要求输出单元频数和期望值的偏差。missprint要求所有频数表输出缺失值的频数。l missing要求把缺失值当作非缺失值看待，在计算百分数及其他统计量时包括它们。out=输出数据集建立一个包括变量值和频数的输出数据集。l sparse要求输出在制表要求中变量水平的所有可能组合的信息。l list以表格形式打印二维表。l nocum/norow/nocol/nofreq/noprint分别不输出累计频率数、行百分率、列百分率、单元频数、频数表。l5. 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 75;proc 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 E Wilcoxon Scores (Rank Sums) for Variable X Classified by Variable GROUP Sum of Expected Std Dev Mean GROUP N Scores Under H0 Under H0 Score 1 7 95.0 73.5000000 12.5718985 13.5714286 2 6 27.0 63.0000000 12.0786894 4.5000000 3 7 88.0 73.5000000 12.5718985 12.5714286 Average Scores Were Used for Ties Kruskal-Wallis Test (Chi-Square Approximation) CHISQ = 8.9839 DF = 2 Prob CHISQ = 0.0112结果说明：组1、组2和组3的秩和（Sum of Scores）分别为95.0、27.0和88.0。原假设（组1、组2和组3的总体分布相同）为真时，期望秩值（Expected）分别为（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，各组的标准差（Std Dev）分别为12.5718985、12.0786894、12.5718985。每组平均得分（Mean Score）分别为95/7=13.5714286、27/6=4.50和88/7=12.5714286。按公式(29.5)调整后多样本的Kruskal-Wallis秩和检验统计量为8.9839 ，用自由度为DF=31=2的卡方分布近似，得到大于近似卡方检验统计量8.9839的概率为p =0.01120.05，因此不能拒绝group和x不线性相关。第二项cmh统计量为ANOVA统计量，首先要求列变量x是有序的。原假设为每一层的3个行的x平均得分是相等的，本例只有一层，且得分指定scores=rank选项，即用秩得分方法，因此就是Kruskal-Wallis秩和检验统计量，自由度为行数减1，即31=2，渐近自由度为2的卡方分布，KW= 8.984，p=0.0110.05，不能拒绝不相关。第三十课 Spearman等级相关分析四、 秩相关的Spearman等级相关分析前面介绍了使用非参数方法比较总体的位置或刻度参数，我们同样也可以用非参数方法比较两总体之间的相关问题。秩相关（rank correlation）又称等级相关，它是一种分析和等级间是否相关的方法。适用于某些不能准确地测量指标值而只能以严重程度、名次先后、反应大小等定出的等级资料，也适用于某些不呈正态分布或难于判断分布的资料。设和分别为和各自在变量X和变量Y中的秩，如果变量X与变量Y之间存在着正相关，那么X与Y应当是同时增加或减少，这种现象当然会反映在（，）相应的秩（，）上。反之，若（，）具有同步性，那么（，）的变化也具有同步性。因此：(30.1)具有较小的数值。如果变量X与变量Y之间存在着负相关，那么X与Y中一个增加时，另一个在减小，具有较大的数值。既然由（，）构成的样本相关系数反映了X与Y之间相关与否的信息，那么在参数相关系数的公式中以和分别代替和，不是同样地反映了这种信息吗？基于这种想法，Charles Spearman秩相关系数应运而生：(30.2)与形式上完全一致，但在中的秩，不管X与Y取值如何，总是只取1到之间的数值，因此它不涉及X与Y总体其他的内在性质，例如，秩相关不需要总体具有有限两阶矩的要求。由于：因此，公式(30.2)可以化简为：(30.3)显然在=时，秩相关系数达到最大值1。又因为：而在每对+=时达到最小值，最小值求法为：所以，最小的为：最大的为：故秩相关系数的最小值为12=1。在原假设和不相关的情况为真时，即秩相关系数为0时，的期望值为0，样本的方差为(30.4)自由度为且分布关于零点对称。当10时，的样本分布可以标准化为近似的t分布：(30.5)例30.1某公司想要知道是否职工期望成为好的销售员而实际上就能有好的销售记录。为了调查这个问题，公司的副总裁仔细地查看和评价了公司10个职工的初始面试摘要、学科成绩、推荐信等材料，最后副总裁根据他们成功的潜能给出了单独的等级评分。二年后获得了实际的销售记录，得到了第二份等级评分，见表30.1中的第1到4列所示。统计问题为是否职工的销售潜能与开始二年的实际销售成绩一致。表30.1 职工的销售潜能与销售成绩的秩相关分析职工编号潜能等级销售成绩成绩等级12400111243603113730052441295652556280711633504117102001000892608119822091110538523944Spearman秩相关系数的计算过程见表30.1中的第5到6列所示，最后计算结果为表明潜能与成绩之间是较强的正相关，高的潜能趋向于好的成绩。秩相关系数原假设为0的t检验统计量为：查表自由度为8，t=3.05的双侧p=0.0158。在0.05显著水平上，t分布的上临界点为2.30，由于3.052.30，因此，拒绝秩相关系数为0的原假设，接受潜能与成绩之间存在秩相关。五、 Corr相关过程Corr相关过程用于计算变量之间的相关系数，包括Pearson（皮尔逊）的乘积矩相关和加权乘积矩相关。还能产生三个非参数的关联测量：Spearman的秩相关，Kendall的tau-b和Hoeffding的相关性度量D。该过程也可以计算偏相关等一些单变量的描述性统计量。1. Corr过程说明proc corr过程一般由下列语句控制：proc corr data=数据集 ;var 变量列表;with 变量列表 ;partial 变量列表 ;weight 变量 ;freq 变量 ;By 变量列表 ;run ;proc corr语句调用corr过程，且是唯一必需的语句。如果只使用proc corr这一条语句，过程计算输入数据集中所有数值变量之间的相关系数。其余语句是供选择的。2. proc corr语句的选项l outp数据集名产生含有Pearson相关系数的一个新数据集。l outs数据集名产生含有Spearman等级相关系数的一个新数据集。l outk数据集名产生含有Kendallb相关系数的一个新数据集。l outh=数据集名产生含有Hoeffding D统计量的一个新数据集。l pearson要求计算通常的pearson乘积矩相关系数，是缺省值。l hoeffding要求计算并输出Hoeffding 的D统计量。l kendall要求计算并输出Kendallb相关系数。l spearman要求计算并输出Spearman等级相关系数。l vardefdf | weight | wgt | wdf指定计算方差时的除数：df（自由度n1），weight或wgt（权重之和），n（观察数），wdf（权重之和1）。缺省值为df。l cov计算协方差方差矩阵。l sscp 要求输出平方和与交叉积和。l csscp 要求输出偏差平方和与交叉积和。l best=数值对每个变量输出指定个数的绝对值最大的相关系数。l noprint禁止所有打印输出。l noprob禁止输出同这些相关有联系的显著性概率。l nosimple对原始数据执行标准方差分析。l rank要求按绝对值从高到低的次序对每个变量输出相关系数。l nocorr抑制Pearson相关的计算及输出。l nomiss 将带有某一变量缺失值的观测值从所有计算中除去。l nosimple不输出每个变量的简单描述性统计量。3. var语句该语句列出要计算相关系数的变量。例如，var a b c;则计算a和b，a和c，b和c三对变量之间的相关系数。4. with语句为了得到变量间的特殊组合的相关系数，该语句和var语句联合使用。用var语句列出的变量在输出相关阵的上方，而用with语句列出的变量竖在相关阵左边。例如，var a b;with x y z;则生成x和a，y和a，z和a，x和b，y和b，z和b。5. partial语句为了计算Pearson偏相关，Spearman偏秩相关，Kendall偏tau-b，用该语句给出偏出去（即固定）的变量名。6. weight语句为了计算加权的乘积矩相关系数，用该语句给出权数变量名。该语句仅用于Pearson相关。7. freq语句当规定freq语句时，输入数据集中的每个观察假定代表n个观察，其中n是该观察中freq变量中的值。观察的总数规定为freq变量值的和。8. by语句使用by语句能够获得用by变量定义的分组观察的独立分析结果。六、 实例分析例30.1的SAS程序如下：data study.persons ;input x y ;y=400-y; cards;2 400 4 360 7 300 1 295 6 280 3 350 10 200 9 260 8 220 5 385;proc corr data=study.persons spearman; var x;with y;run;程序说明：建立输入数据集persons