如何正确运用χ2检验——高维表资料齐性检验与SAS实现

胡纯严 ，胡良平 ，2*

（1.军事科学院研究生院，北京 100850；2.世界中医药学会联合会临床科研统计学专业委员会，北京 100029*通信作者：胡良平，E-mail：lphu927@163.com）

对高维表资料进行差异性分析的基本思路是将高维表降为二维表，降维的重要举措就是按一个因素的全部水平或多个因素的全部水平组合对资料进行分层，从而使每层中的资料都是一个二维表资料。一种特殊的高维表就是分层后的二维表为2×2表（即含一个二值的原因变量和一个二值的结果变量），简记为“g×2×2表”。此时，就可依据“队列研究设计（含横断面研究设计）”或“病例对照研究设计”，推导出两套对此种高维表资料进行差异性分析的方法。此法涉及到若干个具体的分析内容，概括地表述为“定性资料的Meta分析［1］”，因篇幅所限，本文只讨论g×2×2表资料的齐性检验问题。

1 高维表资料齐性检验的必要性

1.1 高维表（g×2×2表）资料的表达模式

高维表（g×2×2表）资料的表达模式见表1。

表1 高维表（g×2×2表）的第h层2×2表资料的表达模式

1.2 高维表资料齐性的含义

在高维表资料中，至少有2个因素（或自变量），1个定性的结果变量。除了采用回归分析可以同时考察多个因素对定性结果变量的影响之外，差异性分析的思路是将其他因素当作分层变量，只研究剩余的一个因素对二值定性结果变量的影响，这被称为将高维表降维后使其成为二维表。显然，在分层变量（它可以是1个因素，也可以是多个因素的水平组合）的每个水平下，都有一张二维表。假定分层变量有g（g≥2）个水平，则有g张2×2表（注：本文不考虑g张R×C表）。只有在它们的内部构成基本一致时，才能将它们按某种规则合并或压缩成一张2×2表。于是，人们将“内部构成基本一致”的g张2×2表资料称为满足“齐性的高维表资料”。

1.3 高维表资料齐性检验方法概述

如何度量高维表资料是否满足齐性？若资料取自队列研究设计，则分层后的各张2×2表的相对危险度（RR）或危险率差（RD）之间接近相等，就可以认定此种高维表资料满足齐性；若资料取自病例对照研究设计，则分层后的各张2×2表的优势比（OR，简称为优比）之间接近相等，就可以认定此种高维表资料满足齐性。由于对同一个资料而言，通常OR值与RR值是比较接近的（但必须注意：选择2×2表资料的第1列还是第2列来计算RR的数值是至关重要的，其中有一个与OR值接近，而另一个与OR值的倒数接近），故在SAS统计软件［2］中，仅给出基于OR为效应指标的“齐性检验”的计算公式和SAS软件实现方法。然而，在文献［3-7］中，先给出一个统一的用于齐性检验的“Q检验统计量”；然后再按照3种效应指标（分别为“相对危险度RR或lnRR”“危险率差RD”和“优势比OR或lnOR”）分别展开“Q检验统计量”。也就是说，实际存在3个不同的用于齐性检验的“Q检验统计量”。

在SAS/STAT的FREQ 过程［2］中，高维表资料（特指g×2×2表）优势比齐性检验的方法有以下5种 ，其 中 ，“Breslow-Day检 验 ”“Breslow-Day-Tarone检验”和“Q检验”在本质上都是χ2检验；“I2度量统计量及其不确定性限值”在本质上是Z检验（也可被视为自由度为1的χ2检验）；而“Zelen's精确检验”在本质上类似于分析二维表资料的Fisher's精确检验，即基于超几何分布而推导出的计算方法。

2 高维表资料齐性检验及SAS实现

2.1 高维表资料优势比齐性检验的具体算法

2.1.1 Breslow-Day检验统计量与Breslow-Day-Tarone检验统计量

Breslow-Day检验统计量QBD公式如下，见式（1）：

在式（1）中，QBD为服从自由度为q-1的χ2分布的检验统计量。其中，“ORMH”“E(.)”和“Var(.)”分别代表“基于Mantel-Haenszel方法计算的优势比”“第h层2×2表中第（1，1）网格上的期望频数”和“第h层2×2表中第（1，1）网格上的方差”，其计算公式分别见下式：

应用式（1）所需要满足的前提条件：各层2×2表中理论频数＞5的格子数应大于80.0%。若资料不满足前述的前提条件，SAS软件建议改用Tarone's校正的Breslow-Day检验统计量，简称为Breslow-Day-Tarone检验统计量QBDT，见式（5）：

在上式中，QBDT为服从自由度为g-1的χ2分布的检验统计量。

2.1.2 Q检验统计量

Q检验统计量Q的公式如下，见式（6）：

在式（6）中，Q渐近地服从自由度为g-1的χ2分布。分别见下式：

在上述各式的计算过程中，若某一层的表格中有0频数，就将该层的全部频数分别加上0.5。

2.1.3 I2度量统计量及其不确定性限值

Higgins和Thompson于2002年提出了I2统计量［2］，它被用于度量在分层2×2表中层间非齐性的一种测度。I2以百分数的形式来表达，它可以被解释为层间变异占总变异的比例。其计算公式见式（11）：

在式（11）中，“g”代表表的层数，Q的计算公式见前文式（6）。

Higgins和Thompson于2002年提出，通过构建H统计量的置信区间，可以转换成构建I2的不确定限值（类似于“置信区间”）［2］。首先，定义H统计量如下式：

其次，需要给出log（H）［即ln（H）］的标准误的计算公式，见式（13）和式（14）：

最后，构建I2的不确定性限值（类似于“置信区间”）。

基于式（16），通过转换式（15）间接获得I2的不确定性限值：

当I2为0时，SAS/STAT中FREQ过程将置信限的下限设置为0，并以显著性水平α取代α/2来决定置信限的上限（注：在本质上，此时所求的就是单侧置信区间）。

2.1.4 Zelen's精确检验

Zelen's精确检验是Breslow-Day近似检验的精确版本。用于Zelen's精确检验的“参考集”包括所有可能的g×2×2表，它们的行、列和层合计频数以及各个2×2表中第（1，1）网格上的合计频数与观测到的高维表是相同的。在固定边缘合计的条件下，这个检验统计量是观测到的g×2×2表出现的概率（它可被称为“理论概率”），它是超几何分布概率的乘积。

Zelen's精确检验的“P值”是前述提及的“理论概率”中小于等于观测到的“表概率”的“理论概率之和”。这个检验与二维表时基于Fisher's精确检验［8-10］是相同的。

2.2 高维表资料相对危险度或危险率差齐性检验的具体算法

2.2.1 高维表资料相对危险度RR齐性检验的具体算法

基于“g×2×2表资料”估计共同相对危险度的前提条件是高维表资料应满足齐性。针对“相对危险度”的齐性检验（也称为异质性检验）的具体方法如下［3-5，11］：

在式（20）中，RRMH为共同相对危险度，其计算公式见式（21）；Q为服从自由度为df=g-1的χ2分布的检验统计量。

基于Mantel-Haenszel估计量（简称MH估计量）估计高维表资料共同相对危险度［2］，见式（21）：

2.2.2 高维表资料危险率差RD齐性检验的具体算法

基于“g×2×2表资料”估计共同危险率差的前提条件是高维表资料应满足齐性。针对“危险率差”的齐性检验（也称为异质性检验）的具体方法如下［3-5，11］：

在式（25）中，共同危险率差RDMH是基于MH法计算的，见下式：

在式（26）中，权重w′h由下式定义：

2.3 高维表资料优势比齐性检验的SAS实现

2.3.1 问题与数据

【例1】文献［3］提供了如下资料，试分析5项研究的OR值之间是否满足齐性。见表2。

表2 吸烟与肝细胞癌关系的5项病例对照研究结果

【例2】文献［3］提供了如下资料，试分析6项研究的OR值之间是否满足齐性。见表3。

表3 鼻咽癌与EB病毒感染关系的6项病例对照研究结果

2.3.2 多项研究的OR值之间齐性检验的SAS实现

【例3】沿用例1中的“问题与数据”，试检验5项研究的OR值之间是否满足齐性。

【分析与解答】设所需要的SAS程序如下：

【程序说明】“tables语句”中的选项“cmh”要求对高维表资料进行CMHχ2检验；其后括号内的两个选项的含义分别为：“I2”代表要求计算“I2测度统计量及其不确定性限值”的数值；“QOR”代表要求对各层OR是否满足齐性进行Q检验（注意：系统默认进行“Breslow-Day检验”）。“exact语句”中的选项“eqor”，要求对各层OR是否满足齐性进行Zelen's精确检验。

【SAS输出结果及解释】

以上是“优比齐性的Breslow-Day检验”的输出结果，QBD=0.6277，df=4，P=0.9599，说明5项研究的OR值之间满足齐性。

以上是“优比齐性的Q检验”的输出结果，Q=0.6273，df=4，P=0.9600，说明5项研究的OR值之间满足齐性。

以上输出的是“优比齐性检验I2测度统计量及其不确性限值”的计算结果，I2=0.00%；其不确定性限值的下限与上限都是0.00%。说明5项研究的OR值之间满足齐性。

【说明】本例实际上没有进行Zelen's精确检验，因为SAS的“日志窗口”中显示了一条警告信息，告知用户未进行Zelen's精确检验的理由（就本例而言，显示“频数太大”）。

【例4】沿用例2中的“问题与数据”，试检验6项研究的OR值之间是否满足齐性。

【分析与解答】设所需要的SAS程序如下：

【程序说明】在SAS程序的“tables语句”的CMH选项后的括号内，加上选项“BDT”，要求输出Tarone's校正的Breslow-Day检验统计量QBDT的数值及其P值。

【SAS输出结果及解释】

以上为两种检验的输出结果：其一，“优比齐性的Tarone's校正的Breslow-Day检验”的结果，QBDT=21.9424，df=5，P=0.0005，说明6项研究的OR值之间不满足齐性；其二，“优比齐性的Zelen's精确检验”的结果中，第1行上的“P＜0.0001”是指观测到的高维表发生的概率，而第2行上的“P=0.0004”是指各种符合特定条件的分层表的概率之和。此值与前面的校正结果“P=0.0005”是比较接近的。

以上为“优比齐性的Q检验”的输出结果，Q=19.2304，df=5，P=0.0017，说明6项研究的OR值之间不满足齐性。

以上输出的是“优比齐性检验I2测度统计量及其不确性限值”的计算结果，I2=74.00%＞50.00%；其不确定性限值的下限与上限分别是40.69%与88.60%，说明6项研究的OR值之间不满足齐性。

【说明】因篇幅所限，“高维表资料相对危险度或危险率差齐性检验的SAS实现”从略。

3 讨论与小结

3.1 讨论

本文介绍的统计分析方法主要适用于g×2×2表，而不适用于g×R×C表（R与C中至少有一个大于2）；当分层后的2×2表资料均来自病例对照研究设计时，才适合以“优势比OR”为效应指标；SAS中的FREQ过程只交待“优势比OR值的齐性检验”，而未提及“相对危险度RR的齐性检验”。本文补充介绍了“相对危险度RR”和“危险率差RD”的齐性检验方法。本文例3的计算结果（QBD=0.6277，df=4，P=0.9599）与文献［3］的结果（Q=0.60，df=4，P＞0.05）非常接近；本文例4的计算结果（QBDT=21.9424，df=5，P=0.0005）与文献［3］的相应结果（Q=12.64，df=5，P＜0.05）相差较大，但最终的结论相同。检验统计量Q值上存在的差距，主要原因在于所采取的计算公式不同所致。

【说明】在对g×2×2表资料进行齐性检验时，SAS软件［2］并没有严格区分资料所对应的研究设计类型是队列研究设计、横断面研究设计，还是病例对照研究设计，只是基于不同的数学原理推导出稍有区别的齐性检验公式或算法而已，但在输出结果中，一律显示“优比齐性检验”的字样；而在统计学教科书（例如：文献［3-5］）和RevMan软件［6-7］中，区分得比较清楚。

3.2 小结

本文围绕g×2×2表资料齐性检验问题，解释了进行齐性检验的必要性，介绍了多种齐性检验的计算公式，并结合两个实例，基于SAS软件给出了优势比齐性检验的结果，对输出结果给出了解释，并据此做出了统计结论和专业结论。