打开微信扫一扫，关注圣才：
sc100xuexi
认证官方微博
认证官方微博
现代统计学与SAS应用：R×C表资料的统计分析1
发布人：&&发布日期： 10:08&&共1833人浏览
　　当列联表的行数R与列数C中有１个为２，另１个大于２时，称为２&ｋ或ｋ&２表；当R，C都大于２时，称为R&C表。多数场合下，这２种２维列联表的统计分析方法是相同的，故一般不细区分。特殊情况将另作介绍。
　　第１节 R&C表资料的分类
　　根据２维列联表中２个分组变量的类型以及分析的目的，可对R&C表资料进行分类。因为不同类型的资料或不同的分析目的，有不同的分析方法。为直观起见，首先，请看几组实际资料。
　　１．双向无序列联表
　　表3.2.1 某地6094人按2种血型划分的结果　表3.2.2 心律失常种类与心肌梗塞部位关系
　　━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
　　缓慢心
　　梗塞例数
　　────────────────律失常─────────────────
　　血型 MN血型：M
　　MN　合计
　　种类部位：下壁前壁真后壁心内膜下合计
　　───────────────────────────────────────
　　431　490　902　1823　窦性过缓
　　388　410　800　1598　被动心律
　　495　587　950　2032　房室阻滞
　　137　179　325　641　束支阻滞
　　───────────────
　　───────────────
　　77　6094
　　合计　16　27
　　━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
　　注∶被动心律的全称为被动性交界性心律。
　　２．单向有序列联表
　　表3.2.3
　　3种药物疗效的观察结果
　　表3.2.4 患者痰液中SSXBXB含量与病型的关系
　　━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
　　嗜酸性
　　患者例数
　　────────────────白细胞─────────────────
　　效　药物∶ A
　　C　合计
　　含量病型∶A型　B型　C型 D型合计
　　───────────────────────────────────────
　　5　3　11
　　7　5　19
　　31　50　45　126
　　3　3　20
　　5　22　24
　　2　0　10
　　───────────────
　　───────────────
　　100　85　85　270
　　合计　17　15　17　11　60
　　━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
　　注∶SSXBXB代表嗜酸性白细胞。
　　３．双向有序且属性不同的列联表
　　表3.2.5　眼晶状体混浊度与年龄之关系　表3.2.6
　　肺门密度与矽肺期次之关系
　　━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
　　晶状体
　　患者例数
　　混浊────────────────密度─────────────────
　　程度年龄∶20～ 30～ 40～合计　级别矽肺期次∶Ⅰ期　Ⅱ期Ⅲ期合计
　　───────────────────────────────────────
　　215　131　148 494
　　1　6　50
　　67　101　128 296
　　96　17 301
　　44　63　132 239
　　72　55 141
　　───────────────
　　───────────────
　　326　295　408 1029
　　169　78 492
　　━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
　　４．双向有序且属性相同的列联表
　　表3.2.7 哩检查室壁收缩运动的符合情况表3.2.8 学生文化课与体育课成绩之关系
　　━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
　　对比法
　　冠心病人数
　　文化课
　　学生人数
　　测定────────────────
　　─────────────────
　　结果核素法∶正常减弱异常合计　成绩体育∶不及格及格良好优秀合计
　　───────────────────────────────────────
　　不及格
　　5　56　74　15　150
　　9　17　34
　　10　32　128　17　187
　　9　13　8　36
　　───────────────
　　───────────────
　　67　53　27 147
　　25　104　223　43　395
　　━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
　　[分类说明] 上面给出了８张有代表性的２维列联表，通过观察喝较，不难看出它们确有差别。设表中横向、纵向分组变量分别为X、Y，则根据X与Y的不同性质和不同分析目的，可得如下分类，见表3.2.9（提示∶定性指标的统计性质分为有序变量、名义变量）。
　　表3.2.9
　　列联表的分类及其适用的统计分析方法
　　━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
　　变量的统计性质及其专业属性列联表分类列联表的编号　适用的统计分析方法
　　━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
　　X、Y皆为名义变量且属性不同双向无序表表3.2.（1；2）　卡方检验，Fisher精确检验
　　X、Y之一名义变量且属性不同单向有序表表3.2.（3；4）　秩和检验，CPD、Ridit检验
　　X、Y皆为有序变量且属性不同双向有序表表3.2.（5；6）　相关分析，线性趋势检验
　　X、Y皆为有序变量且属性相同双向有序表表3.2.（7；8）　一致性检验，特殊模型分析
　　━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
　　注∶各种列联表资料几乎都可运用对数线性模型和logistic回归模型来分析。
　　第２节双向无序R&C表资料的统计分析
　　对于双向无序表而言，所选用的统计方法应当与分组变量各水平的先后顺序或取值大小无关，仅仅与表中总频数、各行的合计、各列的合计有关。符合这些要求的方法有Pearson的拟合优度&2检验、基于似然函数（或熵的分解）导出的似然比&2检验、对数线性模型等方法。
　　表3.2.1资料中各网格上的观察频数都较大，用拟合优度&2统计量进行独立性检验既适合有方便，其公式见式（3.2.1）；而表3.2.2资料中有许多很小的观察频数，最好用SAS软件中FREQ过程并取选择项/EXACT进行Fisher的精确检验。如果没有现成的程序，也可用似然比卡方检验，公式参见第３篇第１章式（3.1.5），但此式中的ｉ代表R&C表中每一网格上的观察频数。
　　１．双向无序R&C表资料的拟合优度卡方检验
　　当A，B２变量都是无序变量时，常假定A与B之间互相独立，在此假设成立的前提下可导出式（3.2.1）的卡方检验统计量，它由统计学家K.Pearson于1900年从检验的拟合优度中发息提出的，故常称为Pearson拟合优度卡方检验。
　　式中 df=（R-1）（C-1），i=1，2，&，R；j=1，2，&，C；n为总频数；ni为第i行
　　的合计频数；nj为第j列的合计频数；nij为第i行第j列的频数。&2服从自由度为df的卡，拒绝域∶&2&&2&（df），则P&&。
　　２．双向无序R&C表资料的对数线性模型分析
　　关于对数线性模型的理论和方法，留到第３篇第４章中专门介绍。这里顺便用SAS软件中的CATMOD过程给出对数线性模型分析的结果，以便使读者知道对同１个资料往往有多种可供选用的方法。
　　[例3.2.1] 根据表3.2.1资料，试分析ABO与MN２种血型系统之间是否独立。
　　[分析与解答] H0∶２种血型系统之间互相独立，H1∶２种血型系统之间不独立，&=0.05。
　　（1）用计算器实现统计计算
　　&2=6094（0...）=8.595， df=（4-1）（3-1）=6
　　查&2临界值表，得∶&20.05（6）=12.59，因&2=8.595&&20.05（6）=12.59，故P&0.05，接受H0。
　　（2）用SAS软件实现统计计算
　　[SAS程序]
　　[CHISQN5.PRG]
　　DATA a；
　　PROC FREQ；
　　DO a=1 TO 4；
　　WEIGHT f；
　　DO b=1 TO 3；
　　TABLES a*b / CHISQ；
　　INPUT f @@；
　　OUTPUT；
　　END； END；
　　PROC CATMOD；
　　CARDS；
　　WEIGHT f；
　　431　490　902
　　MODEL a*b=_RESPONSE_
　　388　410　800
　　/ NOGLS NOPARM
　　495　587　950
　　NORESPONSE ML；
　　137　179　325
　　LOGLIN a b a*b； RUN；
　　LOGLIN a b；
　　（程序的数据步）
　　（程序的2个过程步）
　　[CHISQN5.PRG修改] 第１个过程步是用FREQ过程实现一般的卡方检验，第２个过程步是用CATMOD过程实现对数线性模型分析。 _RESPONSE_是进行此种分析的关键词，必须写。４个选择项的含义分别为∶ NOGLS（不用加权最小平方法计算）、ML（用最大似然估计法计算）、NOPARM（不输出模型中各参数，因为这些参数值对一般的用户作用不大，只要有基于最大似然估计的方差分析表就已达到分析目的）、NORESPONSE（不输出建立对数线性模型的反应矩阵）。
　　[输出结果及其解释]
　　STATISTICS FOR TABLE OF A BY B
　　Statistic
　　Chi-Square
　　Likelihood Ratio Chi-Square　6
　　Mantel-Haenszel Chi-Square
　　Phi Coefficient
　　Contingency Coefficient
　　Cramer's V
　　Sample Size = 6094
　　这是FREQ过程计算的结果，&
　　0.05，接受H0。
　　MAXIMUM LIKELIHOOD ANALYSIS OF VARIANCE TABLE
　　这是模型一的计算结果这是模型二的计算结果
　　Source
　　DF Chi-Square　Prob DF Chi-Square　Prob
　　629.82 0.0000　3
　　680.58 0.0000
　　553.15 0.0000　2
　　647.92 0.0000
　　8.59 0.1981
　　LIKELIHOOD RATIO　0
　　8.67 0.1931
　　这是CATMOD过程计算的结果∶模型一是饱和模型，无法计算误差。其结果已表明∶A（即ABO血型）、B（即MN血型）２分组变量对各网格上理论频数的对数值均有非常显著的影响，而二者之间的交互作用不显著，说明A与B互相独立，即接受H0。从模型二更清楚地看到这一点，这里的似然比检验是反应由模型推算的各网格上的理论频数与观察到的实际频数的吻合程度的，卡方值越小，P值越大，说明吻合程度越高。因&
　　0.05，说明用模型二拟合此资料是满意的。
　　[专业结论] ２种血型系统之间互相独立，即具有ABO血型中某种血型的人用MN血型系统来划分时，属于M、N、MN血型的可能性几乎相等，没有确定的倾向性。
　　３．双向无序R&C表资料的Fisher精确检验
　　通常的统计学教科书上只介绍四格表资料的Fisher精确检验法，我国统计工作者曾在《中国卫生统计》杂志上发表了２&ｋ表资料的精确概率计算方法；SAS软件中FREQ过程可实现R&C表资料的Fisher精确概率计算，下面将给出用程序实现计算的方法。
　　[例3.2.2] 根据表3.2.2资料，试分析缓慢心律失常种类与发生心肌梗塞部位之间是否独立。
　　[分析与解答] H0∶２属性之间互相独立，H1∶２属性之间不独立，&=0.05。前面已说明了本例适合于用Fisher的精确检验或用似然比卡方检验。当然，也可试用对数线性模型。
　　[SAS程序] 只需用表3.2.2中的４行４列原始频数取代[CHISQN5.PRG]中的数据，并将第１个过程步中的选择项CHISQ改为EXACT。若想用对数线性模型计算，因零频数不能取对数，故需修改如下∶在INPUT与OUTPUT语句之间加一句&f=f+1E-05；&这个语句的作用是使所有的频数都增加１个小正数0.00001，以便计算能进行下去，第２个过程步可以不变。
　　[主要输出结果及其解释]
　　Fisher's Exact Test （2-Tail）
　　2.62E-05
　　这是FREQ过程计算的主要结果，精确概率P=0..001，拒绝H0，接受H1。
　　MAXIMUM LIKELIHOOD ANALYSIS OF VARIANCE TABLE
　　Source
　　DF　Chi-Square
　　40.06　0.0000
　　6.70　0.0351
　　LIKELIHOOD RATIO
　　26.42　0.0002
　　这是CATMOD过程计算的主要结果，A（缓慢心律失常种类）、B（发生心肌梗塞部位）２分组变量对各网格上的理论频数之对数的影响为非常显著和显著。因似然比&2=26.42，P=0.，说明用不含交互作用的模型拟合该资料不合适（因本例资料中小格频数太多，用含交互作用项的饱和模型计算时，有些项无法计算出来），从另一个角度也就说明了A与B之间不独立。
　　[专业结论] 缓慢心律失常种类与发生心肌梗塞部位之间不独立，即窦性过缓多发生在下壁和前壁２缚位；房室阻滞多发生在下壁；而束支阻滞多发生在前壁。
　　第３节单向有序R&C表资料的统计分析
　　对于单向有序表而言，所选用的统计方法应当与有序的那个分组变量各水平的先后顺序或取值大小有关。显然，通常的卡方检验已无能为力了。此时，把有序变量当作半定量指标在计算中加以考虑的统计方法，如∶秩和检验、Ridit分析、CPD分析、有序变量的logistic回归模型和有序变量的对数线性模型就显得十分有用了。
　　表3.2.3与表3.2.4资料在形式和应选用的统计分析方法上都是一样的，故只讲解表3.2.3资料，而表 3.2.4资料可作为练习之用。首先，用秩和检验、Ridit分析、CPD分析这３种常用统计方法处理3.2.3资料，并在比较的基础上作出评价。然后，用有序变量的logistic回归模型分析该组资料。
　　[例3.2.3] 根据表3.2.3资料，试分析A，B，C３种药物的疗效之间是否有显著差别。　[分析与解答] H0∶３种药物的疗效总体相同，H1∶３种药物的疗效总体不同或不完全相同，&=0.05。
　　１．单向有序R&C表资料的秩和检验
　　（1）方法的概述
　　秩和检验的基本思想∶把像疗效这样的有序指标的各分组标志看作一个特定的取值区间，接受不同药物治疗的患者之合计值便形成了该区间的宽度。第１个区间的起点从１开始，前１个区间的终点加１便是后１个区间的起点，依次类推。把在这些区间内所取的值称为秩（本质上是顺序号），落于同１个区间内的接受各药物治疗的患者的秩无法准确地区分开，只能都用该区间的平均秩代替，最后求出各药物组的秩和。显然，当疗效按&治愈&无效&顺序排列时，秩和越小越好；反之亦然。当然，秩和只是间接的统计量，需要相应的检验统计量，方可作出统计推断。
　　先计算检验统计量H见式（3.2.2），再计算校正值Hc见式（3.2.3）。Hc近似服从df=C-1的卡，这里C为实验分组的组数。
　　（3.2.2）
　　（3.2.3）
　　式中ni，Ri分别为第i个实验组的样本含量及秩和（i=1，2，&，C），n为总样本含量，tj为有序指标的第j个组内的合计频数（即重复的秩数，j=1，2，&，R）。
　　（2）用计算器实现统计计算
　　表3.2.10
　　用秩和检验法分析３种药物疗效需要的数据
　　━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
　　───────────────秩次范围平均秩次─────────────
　　效　药物∶ A
　　C　合计
　　────────────────────────────────────────
　　1- 20　10.5　157.5
　　42.0　10.5
　　73　21- 93　57.0 2793.0　513.0　855.0
　　31　50　45　126　94-219　156.5 4851.5　7825.0　7042.5
　　5　22　24
　　51　220-270　245.0 0.0　5880.0
　　────────────────────────────────────
　　100　85　85　270
　　- 70.0 13788.0
　　━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
　　注∶表中&药物&为实验分组因素，&疗效&为有序的指标分组变量。
　　H=534.287
　　Hc的分母=1-[（203-20）+（733-73）+（）+（513-51）]/（）=0.
　　Hc=534.287/0..146
　　因C=3，df=3-1=2，查卡方值表，得∶ =9.21。又因
　　9.21，故P＜0.005，拒绝H0，接受H1。
　　２．单向有序R&C表资料的Ridit分析
　　（1）方法的概述
　　Ridit一词的前３个字母Rid系Relative to an Identified Distribution的缩写，它的后２个字母it系Unit的字尾，意即&与特定相对应的单位&。它是一种关于指标分组变量为有序变量的资料进行对比组与标准组比较的假设检验法。
　　Ridit分析的基本思想是∶先确定１个标准组作为特定总体，求得各等级（即有序变量的各水平）的参照值R，用各实验分组的各等级的频数与各R值加权平均，进而求得各实验分组的平均参照值。可以证明∶标准组R的均数 R=0.5，若H0（对比组系来自标准总体的随机样本）成立，则对比组总体R值的100（1-&）％置信区间包括0.5的概率为（1-&）；否则，将根据检验水准&拒绝H0，认为对比总体与标准总体有显著差别。
　　确定标准组的方法通常有以下２种∶①用实验分组中样本含量最大的组作为标准组；
　　②当实验分组的各组样本含量接近时，按指标分组的各水平组分别对各实验组频数进行合并，用合并后的结果（即各等级的频数）作为标准组。
　　计算标准组的R和对比组的R值的方法分别见式（3.2.4）和（3.2.5）；计算标准组R值的方差见式（3.2.6），计算对比组R-值的标准误差见式（3.2.7），对比组总体R值的95％与99％置信区间按式（3.2.8）与（3.2.9）计算。关于标准组各等级的R值的算法将结合实例讲解。
　　（3.2.4）
　　（3.2.5）
　　（3.2.6）
　　（3.2.7A）
　　（3.2.7B）
　　式中F为标准组各等级的频数，N为标准组总样本含量，f为某对比组各等级的频数，n为某对比组的总样本含量。
　　（3.2.8）
　　（3.2.9）
　　（2）用计算器实现统计计算
　　表3.2.11
　　用Ridit分析方法分析3种药物疗效时计算各等级R值需要的数据
　　━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
　　疗效　合计频数（2）&0.5 （2）累计下移1格（3）+（4） R=（5）/N
　　────────────────────────────────────────
　　10.0　0.=（270-20）
　　56.5　0.209259　157=（250-20-73）
　　156.0　0.577778　-42=（157-73-126）
　　244.5　0.9=（-42-126-51）
　　━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
　　注∶N=270。第（7）列是用CPD分析法分析单向有序资料所需的y值及其计算方法。
　　=（20&0.&0.&0.&0.905556）/270=0.5（核对之用）
　　=（15&0.&0.&0.&0.905556）/100=0.332481
　　=（4&0.&0.&0.&0.905556）/85=0.598148
　　=（1&0.&0.&0.&0.905556）/85=0.598933
　　=[2703-（203+733+）]/[12&2702&（270-1）]=0.
　　A药组∶ =（0.0）1/2=0.，
　　B药组∶ =（0.）1/2=0.，
　　C药组∶其标准误差与B药组相同。
　　R-的95%置信区间为∶
　　R-的99%置信区间为∶
　　A药组∶（0.2796，0.3854）；
　　A药组∶（0.2629，0.4020）
　　B药组∶（0.5408，0.6555）；
　　B药组∶（0.5227，0.6736）
　　C药组∶（0.5415，0.6563）。
　　C药组∶（0.5235，0.6744）。
　　显然，R-的95％和99％置信区间都不包含0.5，说明３个药物组的疗效与标准组都有非常显著的差别，但A药的置信区间在0.5之下，B，C２药的置信区间在0.5之上，计算时是按&治愈&无效&顺序排列的，即R-值越小，疗效越高，故可认为A药疗效高，B，C药疗效差，它们之间无显著差别。
　　３．单向有序R&C表资料的CPD分析
　　（1）方法的概述
　　CPD是交差积差（Cross Product Difference）的英文缩写，用它作为１个统计量的符号，代表有序分组数据某实验组的交差积差和，它相当于秩和检验中的秩或Ridit分析中的R-值。
　　CPD分析的基本思想与Ridit分析十分相似，它是利用各实验组的合计频数求出各等级的y值，此值类似于Ridit分析中的R值，然后，用各实验组的各等级下的频数与各y值对应相乘，并求和，就得到各实验组的CPD值。当疗效按&治愈&无效&顺序排列时，CPD值越大越好；反之亦然。当然，CPD只是间接的统计量，需要相应的检验统计量，方可作出统计推断。
　　H0∶CPD/n=0， H1∶CPD/n&0， &=0.05（注∶对于任１个实验组而言）。其检验统计量见式（3.2.10），式中的U服从标准正态。
　　（3.2.10）
　　式中n.j为第j个实验组频数之和，ni.为指标分组的第i个水平组频数之和，N为R&C表中的总频数，Scpd为CPD的标准误差。U的临界值为∶｜U｜&U（&/2c），P&&。这里C为实验组的组数，（&/c）/2所代表的概率是指标准正态曲线下的单侧概率，即按此值直接查标准正态曲线下的面积表得U的临界值。为使用方便，这里给出几组∶
　　表3.2.12
　　用CPD分析法分析单向有序资料所需的U临界值
　　────────────────────────────────────────
　　U（&/2c）界值
　　&　─────────────────────────────────────
　　C∶　2
　　────────────────────────────────────────
　　1.96　2.39　2.49　2.57　2.63　2.68　2.73　2.77　2.80
　　2.58　2.93　3.02　3.09　3.14　3.19　3.23　3.26　3.29
　　────────────────────────────────────────
　　根据式（3.2.10）算出的U值的大小，可用&显著性水平将其划分成３部分∶
　　U&-U（&/2c）　| -U（&/2c）&U&U（&/2c）
　　|　U&U（&/2c）
　　──────────|───────────|───────────
　　（2）用计算器实现统计计算
　　利用表3.2.10喉3.2.11，介绍CPD分析的计算方法∶先看表3.2.11最后１列y是怎样算得的。y1=270-20=250；y2=250-（20+73）=157；y3=157-（73+126）=- 42；y4=-42-（126+51）=-219。有了各yi值，求各药物组的CPD值方法如下∶
　　A药组∶ CPD=15&250+49&157+31&（-42）+5&（-219）=9046
　　B药组∶ CPD=4&250+9&157+50&（-42）+22&（-219）=-4505
　　C药组∶ CPD=1&250+15&157+45&（-42）+24&（-219）=-4541
　　为计算式（3.2.10）中的U值，先计算W值∶
　　/[3&270（270-1）]=78.
　　UA=&（270-100）W]1/2=7.8195
　　UB=-&（270-85）W]1/2=-4.0490
　　UC=-&（270-85）W]1/2=-4.0814
　　查表3.2.12，当&=0.01、C=3时，得∶U（&/2c）=2.93，因
　　2.93，故P&0.01；因UB= -4.，故P&0.01；因UC=-4.，故P&0.01。在区间图上表示出来便一目了然∶
　　UC、UB
　　──────────|────────────|──────────&治愈
　　[评论] 上面用３种非参数分析方法对同１个单向有序列联表资料进行了分析，显然，秩和检验法只能给出３种药物的疗效之间差别显著这样一个笼统的结论来；而Ridit分析和CPD分析却能把它们之间的差别反映得更加具体。后２种方法比较，CPD分析在计算上显得更简捷一些。
　　４．用SAS软件实现上述３种方法的统计计算
　　SAS程序参见本书第１篇第４章第５节[SAS程序─KRCTEST.SAS]。这里只给出输出结果。
　　Kruskal-Wallis's Test
　　Adjusted Chi-Square（Hc）
　　Prob&Hc
　　61.1456
　　0.0000
　　Result of CPD and RIDIT Analysis
　　CPD Values CPD's U-Values RIDIT Values RIDIT's U-Values 药物组
　　7.8195
　　0.3325
　　6.2047
　　-4.0490
　　0.5981
　　-3.3516
　　-4.0814
　　0.5989
　　-3.3784
　　ALPHA= 0.05　CPD's Critical U-Values　-2.3940
　　2.3940
　　RIDIT's Critical U-Values -1.9600
　　1.9600
　　ALPHA= 0.01　CPD's Critical U-values　-2.9352
　　2.9352
　　RIDIT's Critical U-Values -2.5758
　　2.5758
　　[输出结果的解释] 首先给出的是秩和检验的结果∶Hc=61.1456，P&0.001。接着给出的是分别用CPD和Ridit法分析的结果∶前３行分别是根据A、B、C药资料算得的有关统计量的值，即①CPD值、②检验CPD所对应的U统计量的值、③Ridit值、④检验Ridit所对应的U 统计量的值。因U0.01=2.576，由第②、④列U值可知∶它们的绝对值均大于2.576，故P＜0.01。说明各CPD值与0之间、各Ridit值与0.5之间的差别都非常显著。再结合U值的符号和疗效的排列顺序（治愈&无效）可知∶A药疗效最好，B、C药的疗效接近且很差。
　　最后给出检验CPD和Ridit所对应的U统计量值的95％和99％置信区间，当U值低于（或高于）区间的下限（或上限）时，则在给定的&（即ALPHA）水准上认为差别显著。
　　５．单向有序R&C表资料的logistic回归分析
　　（1）方法的概述
　　通常，logistic回归模型适用于因变量Y为二值的回归分析的场合，即Y=0或Y=1。若Y=0代表所关心的某事件或现象发生，其概率记为P（Y=0）=p；则Y=1就代表该事件或现象不发生，其概率记为 P（Y=1），并有关系式P（Y=1）=1-P（Y=0）=1-p。如果有很多因素影响Y的取值，就称它们是自变量，记为X1，X2，&，Xn，它们中既有定性的变量，也有定量的变量。经过研究发现∶刻划自变量与p之间变化规律的一个十分有效的表达式为式（3.2.11），此式就是著名的logistic线性回归模型。从此式中解出p，便得到 logistic非线性回归模型，见式（3.2.12）。
　　（3.2.11）
　　（3.2.12）
　　式（3.2.11）等号左边常用logit（p）表示，叫做对概率p作了logit变换；式（3.2.12）中的exp（*）代表e（自然对数的底，e&2.）的*次方。
　　当因变量Y可取１组有序的值（通常取秩或等级编号）时，如何使用logistic回归分析呢？仅问题的途径是每次把Y的取值按某种规定划分成２级，如Y=1与Y&2，这样，就又回到通常的logistic回归分析中去了。
　　假定Y取c富同的值1，2，&，c，并且序与数字的大小相同，则Y取各值的概率为pi=P（Y=i），i=1，2，&，c。于是，可引入不同的概率之比，取对数后相应的意义也不同。其中较常用的是累积logit变换，见式（3.2.13）。
　　（3.2.13）
　　这就依次将c类合并为２类，对２类进行logistic回归分析。
　　模型中各参数需用最大似然法进行估计，具体计算方法请看有关专著。
　　（2）用SAS软件实现统计计算─应用举例
　　[例3.2.4] 在前例中，已用３种非参数法分析了表3.2.3资料，下面拟用&有序变量的logistic回归分析法&分析A，B，C３种药物的疗效之间是否有显著差别。
　　[分析与解答] 首先，给出２段SAS程序，第１段程序是调用logistic过程分析单向有序列联表资料；第２段程序是利用已求出的回归参数代入logistic回归方程计算各药物的疗效所对应的概率。
　　[SAS程序─DXYXLOG1.PRG]
　　[SAS程序─DXYXGL.PRG]
　　DATA dxyx；
　　DATA abc；
　　DROP z1-z4；
　　%MACRO prob（x1，x2）；
　　x1=0； x2=0；
　　w1=EXP（-3.*&x1+0.00275*&x2）；
　　IF _N_=1 THEN x1=1；
　　w2=EXP（-1.*&x1+0.00275*&x2）；
　　IF _N_=2 THEN x2=1；
　　w3=EXP（0.*&x1+0.00275*&x2）；
　　ARRAY z z1-z4；
　　P1=ROUND（w1/（1+w1），0.0001）；
　　INPUT z1-z4；
　　P2=ROUND（w2/（1+w2），0.0001）；
　　DO OVER z；
　　P3=ROUND（w3/（1+w3），0.0001）；
　　y=_I_；
　　%MEND prob；
　　wt=z；
　　%prob（1，0）
　　OUTPUT；
　　q1=p1； q2=p2； q3=p3；
　　%prob（0，1）
　　CARDS；
　　q4=p1； q5=p2； q6=p3；
　　15 49 31 5
　　%prob（0，0）
　　4 9 50 22
　　q7=p1； q8=p2； q9=p3；
　　1 15 45 24
　　FILE PRINT；
　　PUT #2 @10 'P（=1）' @25 'P（&=2）' @40 'P（&=3）'；
　　PROC LOGISTIC；
　　PUT #4 @10 q1 @25 q2 @40 q3；
　　WEIGHT wt；
　　PUT #6 @10 q4 @25 q5 @40 q6；
　　MODEL y=x1-x2；
　　PUT #8 @10 q7 @25 q8 @40 q9；
　　[DXYXLOG1.PRG修改指导] z1-z4代表疗效的4个等级（治愈&无效），也是数组z的4个元素，属中间变量，用drop语句阻止它们被写入数据集dxyx中去。X1， X2都是哑变量（其个数比实验分组的组数少１），X1=1，X2=0时，代表A药；X1=0，X2=1时，代表B药，X1=X2=0时，代表C药。do- end构成了１个do组，下标变量为隐含的变量名，它从１变到４；&y=_i_&产生反应变量，其取值依次为1，2，3，4。 wt是与每１个Y的取值相对应的频数。注意数据输入方式∶横向为实验分组变量、纵向为有序指标分组变量。
　　[DXYXGL.PRG修改指导]　%MACRO-%MEND是１个宏（相当于高级计算机语言中的子程序），用它来实现概率值的多次计算。宏变量名为prob（X1，X2），其中X1，X2 为参数，共有３组取值，故需３次调用宏进行计算。W1～W3的表达式中系数是左侧程序计算的结果，其表达式是logistic方程式。每次所算得的P1， P2， P3分别代表Y=1（即疗效为治愈）、Y&2（即疗效为显效和治愈）、Y&3（即疗效在好转以上）的概率。
　　[左侧程序的输出结果及其解释]
　　Score Test for the Proportional Odds Assumption
　　Chi-Square = 5.8492 with 4 DF （p=0.2107）
　　这是在拟合有序变量logistic模型前对资料是否满足其前提条件所作的检验，结果表明前提条件已满足（因&
　　0.05）。
　　Criteria for Assessing Model Fit
　　Intercept
　　Intercept
　　Criterion
　　Covariates
　　Chi-Square for Covariates
　　663.122
　　600.143
　　664.577
　　602.567
　　-2 LOG L
　　657.122
　　590.143
　　66.980 with 2 DF （p=0.0001）
　　61.373 with 2 DF （p=0.0001）
　　这部分是关于模型中所考虑的协变量有无意义的检验，从左开始的第①列是４种统计量的名称（如AIC，SC等）；第②列是模型中仅有截距项时统计量的取值；第③列是模型中既有截距项又有协变量时统计量的取值；第④列是检验协变量贡献大小的卡方统计量的取值及其概率。因P＜0.0001，说明在模型中引入协变量（本例的协变量为X1、X2）是有意义的。
　　Analysis of Maximum Likelihood Estimates
　　Parameter Standard
　　Pr &　Standardized
　　Variable　Estimate　Error　Chi-Square Chi-Square　Estimate
　　INTERCP1
　　-3.7483　0.3377
　　123.1817
　　0.0001
　　INTERCP2
　　-1.5521　0.2429
　　40.8235
　　0.0001
　　INTERCP3
　　0.9782　0.2236
　　19.1370
　　0.0001
　　2.0850　0.3088
　　45.5979
　　0.0001
　　2.750162
　　0.00275　0.2956
　　0.0001
　　0.9926
　　0.003493
　　这是按累积logit变换法拟合有序变量logistic回归模型所得的参数，设∶
　　W1=EXP（-3.X1+0.00275X2）；
　　W2=EXP（-1.X1+0.00275X2）；
　　W3=EXP（0.X1+0.00275X2）；
　　则与A，B，C药对应的logistic回归方程手别为∶
　　A药∶令X1=1，X2=0，代入下列３式∶
　　P（Y=1）=P（治愈）=W1/（1+W1）=0.1593； P（Y&2）=P（显效以上）=W2/（1+W2）=0.6302；
　　P（Y&3）=P（好转以上）=W3/（1+W3）=0.9553，显然，P（Y=4）=P（无效）=1-P（Y&3）=0.0447。
　　同理，可算出与B，C药对应的概率。用右侧的程序计算的结果如下∶
　　各种疗效所对应的概率
　　P（治愈）
　　P（显效以上）　P（好转以上）　P（无效）
　　0.1593
　　0.6302
　　0.9553
　　0.0447
　　0.0231
　　0.1752
　　0.7273
　　0.2727
　　0.0230
　　0.1748
　　0.7268
　　0.2732
　　从这个结果可以清楚地看出∶ A药疗效明显优于B，C药，而B，C药的疗效几乎相同。这个结论不仅与前面用３种非参数法分析所得的结论是吻合的，而且以数学模型的形式将结果表达出来，使得对问题的研究更深入了一步。
　　第４节双向有序且属性不同的R&C表资料的统计分析
　　像表3.2.5（表3.2.6也一样）那样的R&C表资料，横向变量&晶状体混浊程度&与纵向变量&年龄&都是有序变量，其属性在专业上是不同的。当把这２个变量看成是地位平等的相互关系时，常需考察它们之间是否存在线性相关关系，即需对资料进行相关分析；当把&年龄&看成是自变量，把&晶状体混浊程度&看成是因变量时，常需考察它们之间是否存在直线变化趋势，即需对资料进行线性趋势检验。
　　１．双向有序且属性不同的R&C表资料的简单相关分析
　　（1）方法的概述
　　变量虽然是有序的，但毕竟还不是定量的，故需给有序变量的各等级赋值方可进行相关分析。最简单的赋值法是按顺序赋给秩次（即得分），即给行变量的等级赋值1，2，&，R和给列变量的等级赋值1，2，&，C。这样（X，Y）的不同取值就有R&C对，表中的R&C个频数就是这R&C对取值所对应的频数，用式（3.2.14）计算Spearman秩相关系数，并作显著性检验，这是比较粗糙的分析方法。
　　（3.2.14）
　　式中T为纵向任何２列合计频数与CPD值（概念与计算方法参见第３篇第２章第３节）交叉乘积相减后求和，Tx为总频数的立方与横向各行合计频数的立方和之差，Ty为总频数的立方与纵向各列合计频数的立方和之差。
　　（2）应用举例
　　[例3.2.5] 对表3.2.5资料进行Spearman秩相关分析。
　　[分析与解答] H0：&s=0， H1：&s&0，&=0.05。
　　①用计算器实现统计计算
　　20～　30～
　　535=（）
　　296　-255=（535-494-296）
　　239　-790=（-255-296-239）
　　合计　326
　　408　1029
　　CPD　63180　-5440　-57740
　　各列CPD值=该列频数与y值对应相乘后求和
　　T=[326&（-5440）-295&63180]+[326&（-57740）-408&63180]
　　+[295&（-57740）-408&（-5440）]=-
　　Tx=10293-（93）=-9407350
　　Ty=10293-（83）=-1311726
　　将T，Tx，Ty的值代入式（3.2.14），得：ｒs=-0.25336，df=1027，查等级相关系数ｒs临界值表，得：ｒs0.02（100）=0.233，因｜ｒs｜
　　0.233，故P&0.02，接受H1。
　　②用SAS软件实现统计计算（程序见[SXYXFX1A.PRG]），下面是此程序第１个过程步输出的结果：
　　|R| under Ho： Rho=0
　　/ N = 1029 / FREQ Var = F
　　1.00000
　　0.25336
　　0.0001
　　0.25336
　　1.00000
　　0.0001
　　[输出结果的解释] SAS输出的是ｒs=0.253，P＜0.0001。这与用式（3.2.14）算出来的符号相反，原因是该式运用了CPD统计量，其值与横向、纵向有序变量各等级排列的顺序有关。若纵向不变，把横向中的第１与第３两行频数互换，则可得出ｒs=0.253。若把表中２变量之间的变化关系通过图形来表达，则写在横向的变量&晶&度&，相当于纵轴y，从下而上数值应由小到大；而写在纵向的变量&年龄&，相当于横轴x，从左至右数值应由小到大，才可正确的反映２变量之间变化关系的真实情况。
e书题库免费下载
推荐的课程