双边数据风险差的一致性检验

Homogeneity Test of Risk Difference for Bilateral Data

摘要: 本文探讨了多组双边数据风险差的一致性假设检验问题及其检验过程，在一致性检验中，当样本量较小时，Wald统计量和Score统计量的第一类错误率接近于预设的显著性水平0.05，然而似然比统计量显示出了比较膨胀的第一类错误率。当样本量较大时，Score统计量检验效果更佳。随着样本量的增加，所有统计量检验效果趋于稳健。因此，在评估第一类错误率性能时，对于多组双边数据，建议采用统计量Score统计量进行风险差的一致性检验。

Abstract: This article explores the homogeneity testing of risk difference for multiple sets of bilateral data and its testing process. In the homogeneity test, when the sample size is small, the Type I error rates of Wald statistic and Score statistic are close to the preset significance level of 0.05, while Likelihood ratio statistic exhibits a relatively inflated Type I error rate. However, when the sample size is larger, the testing effect of Score statistic is better. As the sample size increases, the testing effects of all statistical measures tend to become more robust. Therefore, when evaluating the performance of Type I error rates for multiple sets of bilateral data, it is recommended to adopt Score statistic for the homogeneity test of risk difference.

1. 引言

在医学临床试验中，若研究对象是患者身体成对器官，如耳朵、眼睛等，便会产生相关双边数据，研究者通常将器官患病称之为出现响应。相比较于肝脏、心脏这些单一器官，成对器官产生的双边数据更加复杂。Morris [1] 发现对双边配对数据的统计检验如果忽略了内相关性的存在，可能会导致检验具有夸大的显著性水平。针对双边数据的研究，Ronser [2] 提出了一个假设当一侧器官出现响应时，另外一侧器官也出现响应的条件概率与无条件概率成正比的模型。在此模型基础上，研究者提出了Donner模型 [3] 和Dallal模型 [4] 等适用于研究双边数据的模型。一般来说，风险差经常被用于衡量分组双边数据组间响应率的差异性 [5] 。Zhang等人 [6] 认为，检测风险差一致性是临床试验中一个至关重要的问题。Lui等 [7] 研究了在缺失样本数据中两组双侧数据的风险差的一致性检验问题。Shen等 [8] 研究并推导了两组双边数据风险差一致性检验的三种检验方法。通常而言，鉴于对照组变量的差异性，多个观察组的设置较为常见。因此，在涉及相关配对的双边数据研究中，考虑包含多个观察组和一个对照组的情境是极具意义的 [9] 。综上所述，双边数据风险差的一致性假设检验问题具有重要的研究价值。本文不仅关注理论层面的探讨，更致力于解决实际应用中的问题。通过对双边数据风险差的一致性进行有效检验，更准确地评估风险差异，为决策提供科学依据。因此，本研究的开展具有重要的理论意义和实际应用价值。

2. Dallal模型

设 $m_i$ 为第i ( $i=1,2,\cdots ,g$ )组中的患者数量， $m_{hi}$ 为第i组中有h ( $h=0,1,2$ )个响应的患者数量， $p_{hi}$ 为第i组中无、单边及双边响应的概率，其中 $m_i=m_{0i}+m_{1i}+m_{2i}$ ， $p_{0i}+p_{1i}+p_{2i}=1$ ，具体数据结构如表1所示。

Table 1. Data structure

表1. 数据结构

记 $X_{ijk}$ 是第i组中第j( $j=1,\cdots ,m_i$ )个患者的第k ( $k=1,2$ )只器官响应情况的指标，若无响应，则记 $X_{ijk}=0$ ，否则 $X_{ijk}=1$ 。Dallal模型中包含两个假设：1) 第i组患者一侧器官有响应的概率为 $P\left(X_{ijk}=1\right)={\lambda }_i$ $\left(0\le {\lambda }_i\le 1\right)$ ；2) 患者一侧器官有响应，另一侧器官也有响应的概率为 $P\left(X_{ijk}=1|X_{ij\left(3-k\right)}=1\right)={\gamma }_i$ 。基于假设，可计算出第i组中无、单边及双边响应的概率分别为：

$p_{0i}=1-\left(2-{\gamma }_i\right){\lambda }_i,p_{1i}=2\left(1-{\gamma }_i\right){\lambda }_i,p_{2i}={\gamma }_i{\lambda }_i.$

3. 风险差的一致性检验

假设第一组为对照组，其余组实验组，那么风险差 ${\Delta }_i={\lambda }_i-{\lambda }_1\left(i=2,3,\cdots ,g\right)$ 。一致性检验考虑的问题是各实验组与对照组之间的风险差是否相等，即

如果不能否定原假设H₀，则认为各实验组与对照组之间的响应率无显著性差异。

3.1. H_a和H₀下的极大似然估计

记 $M=\left(m_1,m_2,\cdots ,m_g\right)$ ， $\lambda =\left({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_g\right)$ ， $\gamma =\left({\gamma }_1,{\gamma }_2,\cdots ,{\gamma }_g\right)$ ，根据数据结构，H_a下的对数似然函数可以写为：

$l_{11}\left(\lambda ,\gamma |M\right)={\displaystyle \underset{i=1}{\overset{g}{\sum }}\left(m_{0i}\log \left[1+{\lambda }_i\left({\gamma }_i-2\right)\right]+m_{1i}\log \left[2{\lambda }_i\left(1-{\gamma }_i\right)\right]+m_{2i}\log {\lambda }_i{\gamma }_i\right)}+\log C$ (1)

其中C为一个常数，设 ${\lambda }_i,{\gamma }_i\left(i=1,2,\cdots ,g\right)$ 和 ${\Delta }_i\left(i=2,3,\cdots ,g\right)$ 在H_a下的极大似然估计分别为 ${\hat{\lambda }}_i,{\hat{\gamma }}_i$ 和 ${\hat{\Delta }}_i$ ，则 ${\hat{\lambda }}_i$ 和 ${\hat{\gamma }}_i$ 的值是偏导方程组 $\partial l_{11}/\partial {\lambda }_i=0,\partial l_{11}/\partial {\gamma }_i=0$ 的解，求解方程组可得出 ${\hat{\lambda }}_i=\left(m_{1i}+2m_{2i}\right)/2m_i$ ， ${\hat{\gamma }}_i=2m_{2i}/\left(m_{1i}+2m_{2i}\right)$ 。又因为风险差 ${\Delta }_i={\lambda }_i-{\lambda }_1$ ，故风险差的估计值 ${\hat{\Delta }}_i=\left(m_{1i}+m_{2i}\right)/2m_i-\left(m_{11}+m_{21}\right)/2m_1$ 。在原假设H₀条件下，有 ${\lambda }_i={\lambda }_1+\Delta \left(i=2,3,\cdots ,g\right)$ ，则对数似然函数l₁₁等价于：

$\begin{array}{c}{l}_{10}\left({\lambda }_1,\Delta ,\gamma |M\right)={\displaystyle \underset{i=2}{\overset{g}{\sum }}\left(m_{0i}\log \left[1+\left({\lambda }_1+\Delta \right)\left({\gamma }_i-2\right)\right]+m_{1i}\log \left[2\left({\lambda }_1+\Delta \right)\left(1-{\gamma }_i\right)\right]+m_{2i}\log \left({\lambda }_1+\Delta \right){\gamma }_i\right)}\\ +m_{01}\log \left[1+{\lambda }_1\left({\gamma }_1-2\right)\right]+m_{11}\log \left[2{\lambda }_1\left(1-{\gamma }_1\right)\right]+m_{21}\log {\lambda }_1{\gamma }_1+\log C,\end{array}$ (2)

设 ${\lambda }_1,{\gamma }_i$ 和 $\Delta$ 在H₀下极大似然估计分别为 ${\tilde{\lambda }}_1,{\tilde{\gamma }}_i$ 和 $\tilde{\Delta }$ 。令l₁₀关于 ${\lambda }_1,{\gamma }_i$ 和 $\Delta$ 的偏导均为0，并求解方程组。然而上述方程组没有精确解，故选用费舍尔评分迭代算法计算 ${\tilde{\lambda }}_1,{\tilde{\gamma }}_i$ 和 $\Delta$ 的近似值，算法过程可简单描述为以下4步：

1) 定义各参数的初始值为：

${\lambda }_1^{\left(0\right)}={\hat{\lambda }}_1,{\gamma }_i^{\left(0\right)}={\hat{\gamma }}_i,{\Delta }^{\left(0\right)}=\frac{1}{g-1}{\displaystyle \underset{i=2}{\overset{g}{\sum }}{\Delta }_i},$

2) 第(t + 1)次迭代， ${\tilde{\Delta }}^{\left(t+1\right)}$ 的估计值更新为：

${\Delta }^{\left(t+1\right)}={\Delta }^{\left(t\right)}-{\left(\frac{{\partial }^2{l}_{10}}{\partial {\Delta }^2}\right)}^{-1}\frac{\partial l_{10}}{\partial \Delta },$

3) 第(t + 1)次迭代， ${\lambda }_1^{\left(t+1\right)}$ 和 ${\gamma }_i^{\left(t+1\right)}$ 的估计值更新为：

$\left[\begin{array}{c}{\lambda }_1^{\left(t+1\right)}\\ {\gamma }_1^{\left(t+1\right)}\\ ⋮\\ {\gamma }_g^{\left(t+1\right)}\end{array}\right]=\left[\begin{array}{c}{\lambda }_1^{\left(t\right)}\\ {\gamma }_1^{\left(t\right)}\\ ⋮\\ {\gamma }_g^{\left(t\right)}\end{array}\right]+I{\left({\lambda }_1^{\left(t\right)},{\gamma }_1^{\left(t\right)},\cdots ,{\gamma }_g^{\left(t\right)}\right)}^{-1}\left[\begin{array}{c}\frac{\partial l_{10}}{\partial {\lambda }_1}\\ \frac{\partial l_{10}}{\partial {\gamma }_1}\\ ⋮\\ \frac{\partial l_{10}}{\partial {\gamma }_g}\end{array}\right]{|}_{\Delta ={\Delta }^{\left(t+1\right)}}$

其中I是费舍尔信息矩阵。

4) 重复步骤1)~3)，直到所有参数的估计值趋于收敛。

3.2. 检验统计量

构造以下三个常见的检验统计量：似然比检验统计量、Wald检验统计量和Score检验统计量。

3.2.1. 似然比统计量

记 $\hat{\lambda }=\left({\hat{\lambda }}_1,{\hat{\lambda }}_2,\cdots ,{\hat{\lambda }}_g\right)$ ， $\hat{\gamma }=\left({\hat{\gamma }}_1,{\hat{\gamma }}_2,\cdots ,{\hat{\gamma }}_g\right)$ ， $\tilde{\lambda }=\left({\tilde{\lambda }}_1,{\tilde{\lambda }}_2,\cdots ,{\tilde{\lambda }}_g\right)$ ， $\tilde{\gamma }=\left({\tilde{\gamma }}_1,{\tilde{\gamma }}_2,\cdots ,{\tilde{\gamma }}_g\right)$ 。为了检验假设H₀，构造似然比统计量 $T_L^H$ 为：

$T_L=2\left[l_{11}\left(\hat{\lambda },\hat{\gamma }|M\right)-l_{10}\left({\tilde{\lambda }}_1,\tilde{\Delta },\tilde{\gamma }|M\right)\right]$

原假设H₀下， $T_L^H$ 渐近服从自由度为g-2的 ${\mathcal{X}}^2$ 分布。

3.2.2. Wald统计量

记 ${\beta }_1={\left({\lambda }_1,{\Delta }_2,{\Delta }_3,\cdots ,{\Delta }_g,{\gamma }_1,{\gamma }_2,\cdots ,{\gamma }_g\right)}_{1\times 2g}$ ，设 ${\hat{\beta }}_1$ 是 ${\beta }_1$ 在H_a下的全局极大似然估计。假设H₀等价于 ${\Delta }_2-{\Delta }_3={\Delta }_3-{\Delta }_4=\cdots ={\Delta }_{g-1}-{\Delta }_g$ ，即 $C_1{\beta }_1^{\text{T}}=0$ ，其中

$C_1={\left[\begin{array}{cccccccc}0& 1& -1& 0& 0& \cdots & \cdots & 0\\ 0& 0& 1& -1& 0& \cdots & \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& 0& 1& -1& \cdots & 0\end{array}\right]}_{\left(g-2\right)\times 2g}$

构造Wald检验统计量 $T_W^H$ 为：

$T_W^H=\left({\beta }_1{C}_1^{\text{T}}\right){\left(C_1{P}^{-1}{C}_1^{\text{T}}\right)}^{-1}\left(C_1^{\text{T}}{\beta }_1^{\text{T}}\right)|{}_{{\beta }_1={\hat{\beta }}_1},$

其中P是费舍尔信息矩阵，原假设H₀下， $T_W^H$ 渐近服从自由度为g-2的 ${\mathcal{X}}^2$ 分布。

3.2.3. Score统计量

记 ${\beta }_2={\left({\lambda }_1,\Delta ,{\gamma }_1,\cdots ,{\gamma }_g\right)}_{1\times \left(g+2\right)}$ ，设 ${\tilde{\beta }}_2$ 是 ${\beta }_2$ 在H₀下的极大似然估计。定义Score检验统计量 $T_S$ 为：

$T_S^H=U_1{P}^{-1}{U}_1^{\text{T}}|{}_{{\beta }_2={\tilde{\beta }}_2},$

其中 $U_1={\left(\frac{\partial l_{10}}{\partial {\lambda }_1},\frac{\partial l_{10}}{\partial {\Delta }_2},\frac{\partial l_{10}}{\partial {\Delta }_3},\cdots ,\frac{\partial l_{10}}{\partial {\Delta }_g},\frac{\partial l_{10}}{\partial {\gamma }_1},\frac{\partial l_{10}}{\partial {\gamma }_2},\cdots ,\frac{\partial l_{10}}{\partial {\gamma }_g}\right)}_{1\times 2g}$ ，P是费舍尔信息矩阵，原假设H₀下， $T_S^H$ 渐近服从自由度为g-2的 ${\mathcal{X}}^2$ 分布。

4. 数值模拟研究

在Monte-Carlo模拟中，设置每组的样本数量为 $m=50,100,150$ ，组数 $g=\text{3},\text{5}$ 。记 $\lambda =\left({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_g\right)$ ， $\gamma =\left({\gamma }_1,{\gamma }_2,\cdots ,{\gamma }_g\right)$ ，在 ${\lambda }_i$ 和 ${\gamma }_i\left(i=1,2,\cdots ,g\right)$ 的选择上，必须确保响应率 $p_{hi}$ 的取值在0.1到0.9之间，否则可能会导致出现数据为0的情况，从而产生不准确的结果，具体参数设置如表2所示。

为了评估三个检验统计量在检验中的性能，计算并比较了其在不同参数设置下的第一类错误率。在表2的每组参数配置下，随机模拟生成10,000个样本，并统计其中检验p值小于显著性水平的次数，通过拒绝次数除以10,000，计算出经验第一类错误率。如果检验的第一类错误率小于0.04或大于0.06，则意味检验表现过于保守或膨胀，否则是稳健的 [10] 。

Table 2. Parameter configuration

表2. 参数配置

在原假设H₀条件下，计算上述统计量在一致性检验中犯第一类错误的概率，取 $\Delta =0.05,0.075,0.1$ 且 ${\Delta }_2={\Delta }_3=\cdots ={\Delta }_g=\Delta$ 。在表2中的每个参数配置下随机生成10,000个样本，通过计算在显著性水平 $\alpha =0.05$ 下的拒绝H₀的比例求得第一类错误率，结果如表3和表4所示。结果表明， $T_W^H$ 和 $T_S^H$ 的第一类错误率接近于显著性水平，而 $T_L^H$ 的第一类错误率在m = 50时表现非常膨胀。同时可以发现，一致性检验统计量 $T_K^H\left(K=L,W,S\right)$ 的第一类错误率均随着样本量的增大而趋于稳健。

此外，在H₀假设下，取 $m=50,100,150$ ，在随机生成1000组参数 $\left(\lambda ,\gamma \right)$ 。对于每种参数设置，每个检验重复10,000次，然后计算第一类错误率。通过图1中的一组箱线图，比较了 $m=50,100,150$ 情况下，上述统计量在第一类错误率方面的表现。结果表明：在一致性检验中，当m = 50时， $T_W^H$ 和 $T_S^H$ 的第一类错误率接近于显著性水平0.05，而 $T_L^H$ 则产生了较为膨胀的第一类错误率。当样本量数较大时， $T_S^H$ 检验效果更好。所有统计量 $T_K^H\left(K=L,W,S\right)$ 随着样本量的增加也更加稳健。因此，基于三个统计量在第一类错误率的表现，对于多组相关配对数据，推荐构建统计量 $T_S^H$ 进行风险差的一致性检验。

Table 3. Type I Error Rate of Each Statistic when g = 3

表3. g = 3时各统计量第一类错误率

Table 4. Type I error rate of each statistic when g = 5

表4. g = 5时各统计量的第一类错误率

Figure 1. Box Plot of the Type I Error Rate of Each Statistic under 1000 Parameters

图1. 1000个参数下各统计量第一错误率箱线图

5. 结论

本文提出了双边数据风险差的一致性假设检验问题及其检验过程，模拟研究发现，当样本量较小时，Wald统计量和Score统计量优于似然比统计量。当样本量较大时，Score统计量检验效果更好。因此，针对第一类错误率性能的考量，对于多组相关配对数据，推荐构建Score统计量进行风险差的一致性检验。本研究还具有广阔的创新空间，未来可以深入研究其他统计量在风险差一致性检验中的表现，以寻找更优的检验方法。