如何用SAS软件正确分析生物医学科研资料XXIII.结果变量为多值有序变量的高维列联表资料的统计分析与SAS软件实现（一）

鲍晓蕾，王璐，胡良平

当观测结果是定性资料时，人们习惯将资料整理成列联表形式。比如“2×2 列联表资料”、“R×C 列联表资料”和“高维列联表资料”等。所谓高维列联表，也就是表中涉及到的定性变量的个数 k ≥ 3。对于高维列联表资料，根据结果变量的性质可将其分为以下三类：一是结果变量为二值变量的高维列联表；二是结果变量为多值有序变量的高维列联表；三是结果变量为多值名义变量的高维列联表。本文将详细介绍结果变量为多值有序变量的高维列联表及其用SAS 软件实现统计分析的内容。

对于结果变量为多值有序变量的高维列联表可使用CMH 校正的秩和检验或有序变量多重 logistic 回归分析等统计分析方法。对数线性模型无法考察资料的有序性，因此不宜选用。若采用 CMH 校正的秩和检验，应注意结合待考察的原因变量是否为多值有序变量而选择合适的统计量；若采用有序变量多重 logistic 回归分析，应注意结合原因变量中是否存在多值名义变量或多值有序变量而决定对原因变量的赋值方法。本文将详细介绍 CMH 校正的秩和检验。

CMH 检验作为扩展的 MH 检验，从一定意义上讲，是一种概括统计方法，它在考虑控制分层因素影响的前提下，根据 R×C 表格中行变量与列变量，即原因变量与结果变量属性的不同，给出三种检验统计量。当行变量与列变量均为多值有序变量时，可选择非零相关统计量；当行变量为名义变量而列变量为多值有序变量时，可选择行平均得分统计量，也称为方差分析统计量（注：有别于定量资料方差分析）；当行变量与列变量均为名义变量或行变量是多值有序变量而列变量为名义变量时，可选择一般关联统计量。

下面以简表形式表示第 h 层的 R×C 表（表 1），h =1、2 …… q。q 为层数（即分层因素的水平数），R 为行数，C 为列数。

上表中 nhij表示第 h 层第 i 行第 j 列所对应的频数；nhi+为第 h 层第 i 行的合计数，i = 1、2 …… R；nh+j为第 h 层第 j 列的合计数，j = 1、2 …… C；nh为第 h 层的合计数。根据上表可以得到第 i 行对应的概率为 Phi+=nhi+/nh，第 j 列对应的概率为 Ph+j= nh+j/nh。

为了方便起见，用向量的形式表示频数及各行与各列的概率：

其中 nh代表各个格子里的频数，是（R×C）× 1 的列向量是它的转置；Ph*+代表各行的概率，是 R×1 的列向量，是它的转置；Ph+*代表各列的概率，是 C×1的列向量，是它的转置。需要注意的是这里向量都是以粗体来表示的，特别是要注意区别用粗体表示的 nh与未加粗体的 nh之间的区别，nh表示所有格子的频数向量，而 nh为第 h 层的合计频数。

在各层之间相互独立，并且每一层中行合计与列合计固定的假定下，原假设 H0为：在任何一层中，原因变量与结果变量之间没有关联。在原假设成立的条件下，频数向量服从多重超几何分布，它的期望值与协方差矩阵分别为：

其中，DPh+*是以向量 Ph+*中元素作为主对角线元素的对角阵；DPh*+是以向量 Ph*+中元素作为主对角线元素的对角阵；⊗ 表示克罗内克积。

在原假设成立时，也可写出第 h 层第 i 行第 j 列所对应的频数 nhij的期望值为：

广义 CMH 统计量定义如下：

需要注意的是，当各层间效应方向不一致时，CMH 统计量的检验功效很低。

使用式 ④ 可以计算三种 CMH 统计量，当结果变量为多值有序变量时，根据原因变量为多值有序变量或名义变量，分别选择其中的非零相关统计量或行平均得分统计量。

在计算非零相关统计量时，列的评分阵 Ch是 1×C阵，行的评分阵 Rh是 1×R 阵，行与列的评分由 FREQ过程中的 SCORES 选项指定。非零相关统计量的自由度为1，它也被称为 Mantel-Haenszel 统计量。当行变量或列变量不是有序变量时，该统计量是没有意义的。非零相关统计量对应的备择假设为：至少在一层中，原因变量和结果变量之间存在线性相关。

在计算行平均得分统计量时，列的评分阵 Ch是 1×C阵，由 SCORES 选项指定；行的评分阵 Rh是 (R － 1)×R矩阵，由 FREQ 过程内部产生：

其中 IR－1是秩为 R － 1 的单位阵，JR－1是元素均为1 的（R － 1）× 1 的列向量。

行平均得分统计量的自由度为 R － 1，它所对应的备择假设为：至少在一层中，R 行之间的平均得分是不同的，也就是按原因变量分为 R 个组之后，不同组别之间关于结果变量的平均得分存在差异。

下面将通过实例向读者介绍如何通过 SAS 软件使用CMH 校正的秩和检验处理结果变量为多值有序变量的高维列联表资料。

【例 1】观察宫腔镜下分离宫腔粘连后两种抗粘连方法的治疗效果，在研究中也考虑了粘连程度对于疗效的影响，收集到的试验数据见表 2，试对该资料进行分析。

表2 不同粘连分度及治疗方法与治疗效果的关系

SAS 程序如下，程序名为 example 1。

D A T A e x a m p l e 1;d o a=1 t o 2;d o b=1 t o 2;d o c=1 t o 3;i n p u t f@@;o u t p u t;e n d; e n d; e n d;c a r d s;10 2 07 3 05 5 50 4 8;R U N;o d s h t m l;P R O C F R E Q;t a b l e s a*b*c/C M H;w e i g h t f;R U N;o d s h t m l c l o s e;

程序说明：首先建立数据集，程序中的 a 表示粘连分度，a = 1 表示 I、II 度，a = 2 表示 III、IV 度；b 表示治疗方法，b = 1 表示球囊组，b = 2 表示置环组；c 表示疗效，c = 1 表示治愈，c = 2 表示有效，c = 3 表示无效；变量 f 表示频数。数据的分析采用 FREQ 过程，在 tables语句中依次列出粘连分度、治疗方法和疗效，列在第一位的变量是需要控制的原因变量，列在第二位的变量是想要考察的原因变量，列在第三位的变量是结果变量。本程序中，a*b*c 表示控制 a 因素，考察 b 与 c 因素之间的关系。若换成 b*a*c，则表示控制 b 因素，考察 a 与 c 之间的关系，得到的结果是不一样的。Tables 语句中的 CMH 选项指定输出 CMH 统计量。ods html 语句则要求将结果以网页格式输出。

SAS 程序运行结果：

“b * c”的汇总统计量“a”的控制

程序运行结果的第一部分包括按粘连分度分级以后，治疗方法和疗效所形成的两个 2×3 列联表，其中包括频数、百分比、行百分比和列百分比。因与结果关系不大，鉴于篇幅关系，此处省略。结果的第二部分如上所示，输出了3 个 CMH 统计量，依次为非零相关统计量、行均值得分统计量和一般关联统计量。本例中结果变量是多值有序的，而原因变量是二值的，所以使用行平均得分统计量，此时的 CMH χ2检验也称为 CMH 校正的秩和检验。此处自由度 v = 1，= 5.3374，P = 0.0209 ＜ 0.05。总的样本含量为 49 例。

【例 2】在一项临床试验中，研究病程与依沙酰胺疗效的关系。试验在三所医院中同时进行，具体疗效数据见表 3，试对病程与治疗效果之间的关系进行分析。

表3 病程与依沙酰胺疗效的数据

SAS 程序如下，程序名为 example 2。

DATA example2;do a=1 to 3;do b=1 to 4;do c=1 to 3;input f@@;output;end; end; end; cards;24 10 3 10 4 235 26 10 10 8 425 12 3 11 4 230 22 8 10 10 520 8 2 10 4 120 17 5 12 10 6;RUN;ods html;PROC FREQ;tables a*b*c/CMH;weight f;RUN;ods html close;

程序说明：本程序与程序 example1 相似。首先建立数据集，程序中的 a 表示试验中心，a = 1 表示中心 1，a = 2表示中心 2，a = 3 表示中心 3；b 表示病程，b = 1 表示 ＜1 个月组，b = 2 表示 1～3 个月组，b = 3 表示 3 个月～ 5年组，b = 4 表示 ＞ 5年组；c 表示疗效，c = 1 表示治愈，c = 2 表示好转，c = 3 表示无效；变量 f 表示频数。调用FREQ 过程，a*b*c表示控制 a 因素，考察 b 与 c 因素之间的关系。Tables 语句中的 CMH 选项指定输出 CMH 统计量。

SAS 程序运行结果：

“b * c”的汇总统计量“a”的控制

程序运行结果的第一部分，即按实验中心分层后形成的二维列联表此处从略，只给出 CMH 检验的结果。本资料考察的病程和疗效均为多值有序变量，在分析时可以选择非零相关统计量考察原因变量与结果变量之间是否存在线性相关关系；也可以选用行均值得分差值统计量考察原因变量不同组之间关于结果变量的平均得分是否存在差异。非零相关统计量的自由度 v = 1，= 12.6274，P = 0.0004 ＜0.05，说明病程与疗效之间存在线性相关关系。行均值得分差值的自由度 v = 3，= 13.1908，P = 0.0004 ＜ 0.05，说明不同病程之间的疗效不同。总样本含量为 403 例。

统计与专业结论：C2MH1 χ = 12.6274，P = 0.0004 ＜ 0.05，说明至少在一层中，病程与疗效之间存在线性相关关系，根据数据可以看出随着病程增加，疗效有下降的趋势。= 13.1908，P = 0.0042 ＜ 0.05，说明不同病程之间的疗效不同。

[1] Hu LP.Medical statistics-analysis of quantitative and qualitative data applying the triple-type theory.Beijing: People’s Military Medical Press, 2009:354-363.(in Chinese)胡良平.医学统计学-运用三型理论分析定量与定性资料.北京:人民军医出版社, 2009:354-363.

[2] Hu LP.Statistics facing scientific problems -- (2) multi-factor designs and linear model analysis.Beijing: People’s Medical Publishing House, 2012:500-507.(in Chinese)胡良平.面向问题的统计学——(2)多因素设计与线性模型分析.北京: 人民卫生出版社, 2012:500-507.