分层单双边配对数据优势比的同质性检验及样本量的确定

Homogeneity Test and Sample Size Determination of Odds Ratios for Stratified Unilateral and Bilateral Paired Data

摘要: 在临床试验中，会对成对身体器官或部位进行试验，获得单边或双边的配对数据；同时，由于患者年龄、性别等混杂因素的影响，相同的治疗方法可能会产生不同的治疗效果。因此，本文基于优势比这一衡量治疗效果的指标对分层单双边混合配对数据进行研究。在Donner模型下提出对分层两组的单双边混合数据中各层优势比的同质性假设。利用迭代算法得到最大似然估计量；构造似然比检验统计量(T_LR)、Wald检验统计量(T_W)和拉格朗日乘数检验统计量(T_SC)进行同质性检验；基于第一类错误率和功效比较三种检验统计量的效果。为临床试验设计和试验数据的检验方法选择提供建议，以提高检验效果。本文还通过对比各检验方法，预测在固定的显著性水平下达到目标功效所需的最小样本量，既减少临床试验成本，又使后续进行的试验有更稳健的第一类错误率和更好的功效。蒙特卡洛模拟结果显示：在同质性假设下，T_W在三个检验统计量中表现最为冒进；T_LR相对冒进，总体表现较为稳健；而T_SC在各种情况下均表现稳健，且功效较高，因此最为推荐使用。对于试验样本量的确定，为了达到相同的功效，T_W所需样本量较小；随着层数的增加，三个检验统计量所需样本量差异逐渐减小；随着优势比增加，越来越接近于1，三个检验统计量所需样本量也会逐渐增加。

Abstract: In clinical trials, pairs of body organs or parts are tested to obtain unilateral or bilateral paired data. At the same time, due to the influence of confounding factors such as patient age and gender, the same treatment method may produce different therapeutic effects. Therefore, this paper studied stratified unilateral and bilateral paired data based on the odds ratio which is an index to measure the treatment effect. In the Donner model, the homogeneity hypothesis of the odds ratios of each stratum in the stratified unilateral and bilateral paired data in two groups is proposed. The maximum likelihood estimator is obtained by iterative algorithm. The likelihood ratio test statistic (T_LR), Wald test statistic (T_W) and Lagrange multiplier test statistic (T_SC) are constructed for homogeneity test. The effects of the three test statistics are compared based on the type I error rate and power, to provide suggestions for clinical trial design and test method selection of test data, so as to improve test results. This paper also predicts the minimum sample size required to achieve the target power at a fixed significance level by comparing these testing methods, which not only reduces the cost of clinical trials, but also enables the subsequent trials to have a more robust type I error rate and better power. Monte Carlo simulation results show that T_W is the most aggressive among the three test statistics under the homogeneity hypothesis. T_LR is relatively aggressive, the overall performance is relatively stable; T_SC is robust in all situations and has high power, so it is most recommended. For the determination of the sample size of the test, in order to achieve the same power, the sample size required by T_W is small. With the increase of the number of stratums, the difference of sample size required for the three test statistics gradually decreases. As the odds ratio increases and gets closer to 1, the sample size required for the three test statistics will gradually increase.

1. 引言

在临床试验中，研究人员想要了解疾病的治疗情况，可以从患者的特定组织或者器官收集数据，并将患者根据治疗方案分组。对于眼睛、耳朵等成对的身体器官，会带来单边或双边的配对数据，在这种情况下，只研究单边或双边的检验不太合适，会使得从数据中获得分析结果不完整、不准确 [1] 。Rosner (1982) [2] 指出成对身体器官的治疗结果通常是高度相关的，它们之间的相关性是不可忽视的，任何忽略这一依赖特征的统计过程都将导致不正确的推论。Donner (1989) [3] 提出了一个配对数据的相关系数模型，Thompson (1993) [4] 证明了Donner的模型可以使充分利用单个或成对器官的数据来提高检测的有效性。Ma和Wang (2021) [5] 构造似然比检验、Wald检验和score检验统计量，检验了单双边混合配对数据比例的同质性，并比较了这些方法下的经验性第一类错误率和功效。治疗效果不仅取决于治疗方案，也可能与病人的性别、年龄、不同治疗环境等影响有关。Kahan和Morris (2013) [6] 指出，对这类配对数据的分析如果受到分层的影响可能会产生混杂因素，不能够忽略这些因素。因此，需要对数据进行分层，以判断不同层间是否存在分层因素的影响。Sun等(2022) [7] 基于单边和双边的分层配对数据进行了风险差的同质性检验，她进一步研究了使检验满足所需功效的样本量。样本量的确定也是试验设计的重要组成部分 [8] ，因此将其融入配对数据的研究。此外，优势比是一个比较两组治疗效果的指标，是一种首选的前瞻性联想测量方法，主要用于回顾性或横断面抽样设计 [9] ，并可对此进行同质性检验。

本文基于分层单双边混合配对数据得到的优势比对同质性假设的检验方法进行研究，可以为生物医学的临床试验数据在检验过程中检验方法的选择上提供建议；对于样本量的研究，可以为将来临床试验的方案设计提供建议，以减少临床实验中人力、物力等成本的消耗。

2. 数据结构和最大似然估计

2.1. 数据结构

假设 $m_{lij}^{\left(q\right)}\left(l=0,1,2,i=1,2,q=1,2,j=1,\cdots ,J\right)$ 是指第j层个体在接受组i类别的治疗后，治疗部位的响应个数。相应的， $p_{lij}^{\left(q\right)}$ 分别表示在第j层中接受组i类别的治疗后，没有响应，单边响应，双边响应的概率。

Table 1. Data structure of the jth stratum

表1. 第j层的单双边混合数据结构

当 $q=1$ 时为单边数据， $m_{ij}^{\left(1\right)}=\left(m_{01j}^{\left(1\right)},m_{11j}^{\left(1\right)}\right)$ 服从二项分布 $B~\left(1,p_{lij}^{\left(1\right)}\right)$ ；当 $q=2$ 时为双边数据， $m_{ij}^{\left(2\right)}=\left(m_{0ij}^{\left(2\right)},m_{1ij}^{\left(2\right)},m_{2ij}^{\left(2\right)}\right)$ 服从三项分布 $m_{ij}^{\left(2\right)}~\left(m_{+ij}^{\left(2\right)},p_{0ij}^{\left(2\right)},p_{1ij}^{\left(2\right)},p_{2ij}^{\left(2\right)}\right)$ 。令 $Z_{hij}^{\left(1\right)}$ 表示治疗单边部位的个体在第j层的治疗部位的响应的指示变量， $Z_{hij}^{\left(1\right)}=1$ 表示存在响应， $Z_{hij}^{\left(1\right)}=0$ 表示未响应。 $Z_{hijk}^{\left(2\right)}\left(k=1,2\right)$ 表示治疗双边部位的个体在第j层第k个部位的响应的指示变量， $Z_{hijk}^{\left(2\right)}=1$ 表示存在响应， $Z_{hijk}^{\left(2\right)}=0$ 表示未响应。在Donner模型中，假设 $P\left(Z_{hij}^{\left(1\right)}=1\right)=P\left(Z_{hijk}^{\left(2\right)}=1\right)={\pi }_{ij}\left(0\le {\pi }_{ij}\le 1\right)$ ， $Corr\left(Z_{hijk}^{\left(2\right)},Z_{hij\left(3-k\right)}^{\left(2\right)}\right)={\rho }_j\left(0\le {\rho }_j\le 1\right)$ [5] 。

基于表1中的数据结构和Donner模型的假设，得到 $p_{lij}^{\left(q\right)}$ 的概率：

$p_{01j}^{\left(1\right)}=1-{\pi }_{ij}$ ， $p_{11j}^{\left(1\right)}={\pi }_{ij}$ ；

$p_{01j}^{\left(2\right)}={\rho }_j\left(1-{\pi }_{ij}\right)+\left(1-{\rho }_j\right){\left(1-{\pi }_{ij}\right)}^2$ ，

$p_{11j}^{\left(2\right)}=2{\pi }_{ij}\left(1-{\rho }_j\right)\left(1-{\pi }_{ij}\right)$ ，

$p_{21j}^{\left(2\right)}={\rho }_j{\pi }_{ij}+\left(1-{\rho }_j\right){\pi }_{ij}^2$ ； $i=1,2$ ， $j=1,\cdots ,J$ 。

令 ${\pi }_i={\left({\pi }_{i1},{\pi }_{i2},\cdots ,{\pi }_{iJ}\right)}^{\text{T}}$ ， $\rho ={\left({\rho }_1,{\rho }_2,\cdots ,{\rho }_J\right)}^{\text{T}}$ ，又 $m_{lij}^{\left(q\right)}=\left(m_{0ij}^{\left(1\right)},m_{1ij}^{\left(1\right)},m_{0ij}^{\left(2\right)},m_{1ij}^{\left(2\right)},m_{2ij}^{\left(2\right)}\right)\left(i=1,2\right)$ 相互独立，因此可以得到它的似然函数为：

$L\left({\pi }_1,{\pi }_2,\rho |m_{lij}^{\left(q\right)}\right)=f\left(m_{0ij}^{\left(1\right)},m_{1ij}^{\left(1\right)},m_{0ij}^{\left(2\right)},m_{1ij}^{\left(2\right)},m_{2ij}^{\left(2\right)}|{\pi }_i,\rho \right)=\underset{l=1}{\overset{2}{{\displaystyle \prod }}}f\left(m_{lij}^{\left(1\right)}|{\pi }_i,\rho \right)\underset{l=1}{\overset{3}{{\displaystyle \prod }}}f\left(m_{lij}^{\left(2\right)}|{\pi }_i,\rho \right)$ ，

则可以得到如下对数似然函数：

$l\left({\pi }_1,{\pi }_2,\rho \right)=\underset{j=1}{\overset{J}{{\displaystyle \sum }}}{l}_j\left({\pi }_{1j},{\pi }_{2j},{\rho }_j\right)+\log C$ ，

其中 $\begin{array}{c}{l}_j\left({\pi }_{1j},{\pi }_{2j},{\rho }_j\right)={\displaystyle \underset{i=1}{\overset{2}{\sum }}\{m_{0ij}^{\left(1\right)}\log \left(1-{\pi }_{ij}\right)+m_{1ij}^{\left(1\right)}\log {\pi }_{ij}+m_{0ij}^{\left(2\right)}\log \left[{\rho }_j\left(1-{\rho }_j\right){\left(1-{\pi }_{ij}\right)}^2\right]}\\ +m_{1ij}^{\left(2\right)}\log \left[2{\pi }_{ij}\left(1-{\rho }_j\right)\left(1-{\pi }_{ij}\right)\right]+m_{2ij}^{\left(2\right)}\log \left[{\rho }_j{\pi }_{ij}+\left(1-{\rho }_j\right){\pi }_{ij}^2\right]\}\end{array}$ ， $C=\underset{j=1}{\overset{J}{{\displaystyle \prod }}}\underset{i=1}{\overset{2}{{\displaystyle \prod }}}\frac{\left(m_{+ij}^{\left(1\right)}\right)!}{m_{0ij}^{\left(1\right)}!m_{1ij}^{\left(1\right)}!}\frac{\left(m_{+ij}^{\left(2\right)}\right)!}{m_{0ij}^{\left(2\right)}!m_{1ij}^{\left(2\right)}!m_{2ij}^{\left(2\right)}!}$ 为一个与参数 ${\pi }_{ij}$ 和 ${\rho }_j$ 无关的常数。

选取优势比作为衡量两组间治疗效果的指标，在分层模型中，令 ${\theta }_j=\frac{{\pi }_{1j}/\left(1-{\pi }_{1j}\right)}{{\pi }_{2j}/\left(1-{\pi }_{2j}\right)}$ 为第j层中两组之间的优势比，反映了各层的治疗方案在治疗效果上的差异，从而建立治疗方案与治疗效果间的联系。本文进行如下假设： $H_0:{\theta }_1={\theta }_2=\cdots ={\theta }_J=\theta \triangleq {\theta }_0$ ， $H_1:\theta \ne {\theta }_0$ ，其中 ${\theta }_0$ 为一个大于零的常数。

2.2. 全局最大似然估计

令 ${\hat{\pi }}_{1j}$ ， ${\hat{\pi }}_{2j}$ ， ${\hat{\rho }}_j$ 为 ${\pi }_{1j}$ ， ${\pi }_{2j}$ ， ${\rho }_j$ 的估计值，对于 ${\pi }_{ij}$ 和 ${\rho }_j\left(i=1,2;j=1,\cdots ,J\right)$ ，通过如下等式求得 ${\pi }_{1j}$ ， ${\pi }_{2j}$ ， ${\rho }_j$ 的全局最大似然估计值：

$\frac{\partial l_j}{\partial {\pi }_{ij}}=0$ ， $\frac{\partial l_j}{\partial {\rho }_j}=0$ 。

计算得到全局最大似然估计值 ${\hat{\pi }}_{ij}$ ， ${\hat{\rho }}_j$ 。由于上述式子不存在显式解，因此需要使用迭代算法对方程进行求解。使用Fisher得分法求解得到最大似然估计值 ${\hat{\pi }}_{ij}$ ， ${\hat{\rho }}_j$ 。

首先设置初值：

${\pi }_{1j}^{\left(0\right)}=\frac{m_{11j}^{\left(1\right)}+m_{11j}^{\left(2\right)}+2m_{21j}^{\left(2\right)}}{m_{+1j}^{\left(1\right)}+2m_{+1j}^{\left(2\right)}}$ ， ${\rho }_j^{\left(0\right)}=\frac{4\left(m_{01j}^{\left(2\right)}+m_{02j}^{\left(2\right)}\right)\left(m_{21j}^{\left(2\right)}+m_{22j}^{\left(2\right)}\right)-{\left(m_{11j}^{\left(2\right)}+m_{12j}^{\left(2\right)}\right)}^2}{\left[m_{11j}^{\left(2\right)}+m_{12j}^{\left(2\right)}+2\left(m_{21j}^{\left(2\right)}+m_{22j}^{\left(2\right)}\right)\right]\left[m_{11j}^{\left(2\right)}+m_{12j}^{\left(2\right)}+2\left(m_{01j}^{\left(2\right)}+m_{02j}^{\left(2\right)}\right)\right]}$ ，

参数的第( $t+1$ )次迭代值可以通过Fisher得分法得到：

$\left[\begin{array}{c}{\pi }_{1j}^{\left(t+1\right)}\\ {\pi }_{2j}^{\left(t+1\right)}\\ {\rho }_j^{\left(t+1\right)}\end{array}\right]=\left[\begin{array}{c}{\pi }_{1j}^{\left(t\right)}\\ {\pi }_{2j}^{\left(t\right)}\\ {\rho }_j^{\left(t\right)}\end{array}\right]+{\left[I_{0j}\right]}^{-1}\left({\pi }_{1j},{\pi }_{2j},{\rho }_j\right)\times \left[\begin{array}{c}\frac{\partial l_j}{\partial {\pi }_{1j}}\\ \frac{\partial l_j}{\partial {\pi }_{2j}}\\ \frac{\partial l_j}{\partial {\rho }_j}\end{array}\right]$ ，

其中 $I_{0j}$ 为关于 ${\pi }_{1j}$ ， ${\pi }_{2j}$ ， ${\rho }_j$ 的信息矩阵， $I_{0j}=\left[\begin{array}{ccc}{I}_{011j}& I_{012j}& I_{013j}\\ I_{021j}& I_{022j}& I_{023j}\\ I_{031j}& I_{032j}& I_{033j}\end{array}\right]$ ， $I_{011j}=E\left(-\frac{{\partial }^2{l}_j}{\partial {\pi }_{1j}^2}\right)$ ， $I_{012j}=E\left(-\frac{{\partial }^2{l}_j}{\partial {\pi }_{1j}{\pi }_{2j}}\right)$ ， $I_{013j}=E\left(-\frac{{\partial }^2{l}_j}{\partial {\pi }_{1j}{\rho }_j}\right)$ ， $I_{022j}=E\left(-\frac{{\partial }^2{l}_j}{\partial {\pi }_{2j}^2}\right)$ ， $I_{023j}=E\left(-\frac{{\partial }^2{l}_j}{\partial {\pi }_{2j}{\rho }_j}\right)$ ， $I_{033j}=E\left(-\frac{{\partial }^2{l}_j}{\partial {\rho }_j^2}\right)$ 。

2.3. 同质性假设下的最大似然估计

在研究中，对每层的优势比做一致性检验后，如果优势比的假设成立，那么这个优势比的值是多少就需要关注。因此在等优势比假设成立的基础上，对优势比进行同质性检验。

在 $H_0:{\theta }_1={\theta }_2=\cdots ={\theta }_J=\theta \triangleq {\theta }_0$ 假设下， ${\pi }_{2j}=\frac{{\pi }_{1j}}{{\pi }_{1j}-{\pi }_{1j}{\theta }_0+{\theta }_0}$ $\left(j=1,\cdots ,J\right)$。此时的未知参数只有 ${\pi }_{ij}$ 和 ${\rho }_j$ 。则相应的对数似然函数可以写作： $l^c\left({\pi }_1,\rho ,{\theta }_0\right)=\underset{j=1}{\overset{J}{{\displaystyle \sum }}}{l}_j\left({\pi }_{1j},{\rho }_j,{\theta }_0\right)+\log C$ 。

令 ${\tilde{\pi }}_{1j}$ ， ${\tilde{\rho }}_j$ 为 ${\pi }_{1j}$ ， ${\rho }_j$ 在 $H_0$ 假设下的限制性最大似然估计值，即求下列式子的解：

$\frac{\partial l_j^c}{\partial {\pi }_{1j}}=0$ ， $\frac{\partial l_j^c}{\partial {\rho }_j}=0$ ， $j=1,\cdots ,J$ 。

使用Fisher得分法计算 $H_0$ 假设下的最大似然估计值 ${\tilde{\pi }}_{1j}$ ， ${\tilde{\rho }}_j$ 。

首先，对于第j层，设置 ${\pi }_{1j}$ 和 ${\rho }_j$ 的初值为：

${\pi }_{1j}^{\left(0\right)}=\frac{m_{11j}^{\left(1\right)}+m_{11j}^{\left(2\right)}+2m_{21j}^{\left(2\right)}}{m_{+1j}^{\left(1\right)}+2m_{+1j}^{\left(2\right)}}$ ， ${\rho }_j^{\left(0\right)}=\frac{4\left(m_{01j}^{\left(2\right)}+m_{02j}^{\left(2\right)}\right)\left(m_{21j}^{\left(2\right)}+m_{22j}^{\left(2\right)}\right)-{\left(m_{11j}^{\left(2\right)}+m_{12j}^{\left(2\right)}\right)}^2}{\left[m_{11j}^{\left(2\right)}+m_{12j}^{\left(2\right)}+2\left(m_{21j}^{\left(2\right)}+m_{22j}^{\left(2\right)}\right)\right]\left[m_{11j}^{\left(2\right)}+m_{12j}^{\left(2\right)}+2\left(m_{01j}^{\left(2\right)}+m_{02j}^{\left(2\right)}\right)\right]}$ ，

参数 ${\pi }_{1j}$ 和 ${\rho }_j\left(j=1,\cdots ,J\right)$ 的第 $\left(t+1\right)$ 次迭代值可以通过Fisher得分法得到：

$\left[\begin{array}{c}{\pi }_{1j}^{\left(t+1\right)}\\ {\rho }_j^{\left(t+1\right)}\end{array}\right]=\left[\begin{array}{c}{\pi }_{1j}^{\left(t\right)}\\ {\rho }_j^{\left(t\right)}\end{array}\right]+I_{1j}^{-1}\left({\pi }_{1j},{\rho }_j\right)\times \left[\begin{array}{c}\frac{\partial l_j^c}{\partial {\pi }_{1j}}\\ \frac{\partial l_j^c}{\partial {\rho }_j}\end{array}\right]$ ，

其中 $I_{1j}=\left[\begin{array}{cc}{I}_{111j}& I_{113j}\\ I_{113j}& I_{133j}\end{array}\right]$ ， $I_{111j}=E\left(-\frac{{\partial }^2{l}_j^c}{\partial {\pi }_{1j}^2}\right)$ ， $I_{113j}=E\left(-\frac{{\partial }^2{l}_j^c}{\partial {\pi }_{1j}{\rho }_j}\right)$ ， $I_{133j}=E\left(-\frac{{\partial }^2{l}_j^c}{\partial {\rho }_j^2}\right)$ 。

由于当样本量足够大时，即 $m\to +\infty$ 时，似然比检验统计量、Wald检验统计量、score检验统计量等均渐进服从自由度为 $J-1$ 的 ${\chi }^2$ 分布，选择这三种方法对假设 $H_0$ 进行检验，比较它们的第一类错误率和功效。

3. 同质性检验

3.1. 检验统计量的构造

似然比检验分别计算在约束和无约束条件下的参数估计值，然后计算二者的对数似然函数是否足够接近 [10] ，所构造的似然比检验统计量的具体公式如下：

$T_{LR}=2\ln \frac{{\hat{L}}^c}{\hat{L}}=-2\left[l^c\left({\tilde{\pi }}_1,\tilde{\rho },\tilde{\theta }\right)-l\left({\hat{\pi }}_1,{\hat{\pi }}_2,\hat{\rho }\right)\right]=-2\underset{j=1}{\overset{J}{{\displaystyle \sum }}}\left[l_j^c\left({\tilde{\pi }}_{1j},{\tilde{\rho }}_j,{\tilde{\theta }}_j\right)-l_j\left({\hat{\pi }}_{1j},{\hat{\pi }}_{2j},{\hat{\rho }}_j\right)\right]$ ，

其中 ${\hat{\pi }}_{1j}$ ， ${\hat{\pi }}_{2j}$ ， ${\hat{\rho }}_j$ 是全局最大似然估计， ${\tilde{\pi }}_{1j}$ ， ${\tilde{\rho }}_j$ ， ${\tilde{\theta }}_j$ 是 $H_0$ 下的最大似然估计。在给定显著性水平 $\alpha$ 的情况下， $T_{LR}>{\chi }_{1-\alpha }^2\left(J-1\right)$ 时拒绝原假设， $T_{LR}\le {\chi }_{1-\alpha }^2\left(J-1\right)$ 时接受原假设。

Wald检验是检验无约束的估计值满足约束条件的程度。令 $\beta ={\left(\theta ,{\pi }_{11},{\rho }_1,\cdots ,{\pi }_{1J},{\rho }_J\right)}_{2J+1}^{\text{T}}$ ，相应的 $H_0$ 下的最大似然估计值为 ${\hat{\beta }}_{}={\left(\hat{\theta },{\hat{\pi }}_{11},{\hat{\rho }}_1,\cdots ,{\hat{\pi }}_{1J},{\hat{\rho }}_J\right)}_{2J+1}^{\text{T}}$ ， $K={\left(1,0,0,\cdots ,0\right)}_{2J+1}$ ，那么 $K\times \hat{\beta }$ 为 $H_0$ 下的最大似然估计值，则构造在假设 $H_0:{\theta }_1={\theta }_2=\cdots ={\theta }_J=\theta \triangleq {\theta }_0$ 下的Wald检验统计量为：

$T_W=\frac{{\left(\hat{\theta }-{\theta }_0\right)}^2}{Var\left(\beta \right)}=\frac{{\left(K\hat{\beta }-{\theta }_0\right)}^2}{I_2^{-1}}=\frac{{\left(\hat{\theta }-{\theta }_0\right)}^2}{K{I}_2^{-1}{K}^{\text{T}}}=\frac{{\left(\hat{\theta }-{\theta }_0\right)}^2}{I_2^{-1}\left(1,1\right)}$ ，

$I_2=\left(\begin{array}{cc}{I}_{244}& B\\ B^{\text{T}}& I_{NJ}\end{array}\right)$ 为关于 $\beta$ 的Fisher信息矩阵， $I_2^{-1}\left(1,1\right)$ 为 $I_2$ 的逆矩阵对角线上的第一个元素，其中 $I_{244}=E\left(-\underset{j=1}{\overset{J}{{\displaystyle \sum }}}\frac{{\partial }^2{l}_j^c}{\partial {\theta }^2}\right)$ ， $B=\left(E\left(-\frac{{\partial }^2{l}_1^c}{\partial \theta {\pi }_{11}}\right),E\left(-\frac{{\partial }^2{l}_1^c}{\partial \theta {\rho }_1}\right),\cdots ,E\left(-\frac{{\partial }^2{l}_J^c}{\partial \theta {\pi }_{1J}}\right),E\left(-\frac{{\partial }^2{l}_J^c}{\partial \theta {\rho }_J}\right)\right)$ ， $I_{NJ}=diag\left(I_{n1},I_{n2},\cdots ,I_{nJ}\right)$ ， $I_{nj}=\left(\begin{array}{cc}\frac{{\partial }^2{l}_j^c}{\partial {\pi }_{1j}^2}& \frac{{\partial }^2{l}_j^c}{\partial {\pi }_{1j}{\rho }_j}\\ \frac{{\partial }^2{l}_j^c}{\partial {\pi }_{1j}{\rho }_j}& \frac{{\partial }^2{l}_j^c}{\partial {\rho }_j^2}\end{array}\right)$ 。

score检验所感兴趣的参数为 $\theta$ ，因此令计分矩阵 $U=\left(\frac{\partial l^c}{\partial \theta },0,0,\cdots ,0,0\right)$ ，score检验统计量就转化为：

$T_{SC}^C={\left(\underset{j=1}{\overset{J}{{\displaystyle \sum }}}\frac{\partial l_j^c}{\partial \theta }\right)}^2{I}_2^{-1}\left(1,1\right){|}_{\theta =\tilde{\theta }}$ ，

其中 $I_2^{-1}\left(1,1\right)$ 为 $I_2$ 的逆矩阵对角线上的第一个元素。

3.2. 仿真模拟

本节使用蒙特卡罗模拟来比较似然比统计量、Wald检验统计量和score检验统计量的经验性的一类错误率和功效，用于衡量三个统计量的性能。第一类错误率是原假设成立时拒绝原假设的概率。当第一类错误率介于0.04和0.06之间时，检验被认为是稳健的；大于0.06时，检验被认为是冒进的；小于0.04时，被认为是保守的。功效是当原假设为假时接受它的概率 [11] 。一个好的检验应该同时具有稳健的第一类错误率和尽可能高的功效。

设样本量 $m=\left(m_{ij}^{\left(1\right)},m_{ij}^{\left(1\right)}\right)=\left(25,25\right),\left(50,50\right),\left(100,100\right)$ ， $J=2,4,6$ 。首先，计算不同参数配置下检验统计量的经验性第一类错误率，具体参数配置如表2所示。

Table 2. Parameter configuration of Monte Carlo simulation

表2. 蒙特卡罗模拟的参数配置

模拟同质性检验的经验性第一类错误率时，在每组的样本量相同的情况下进行构造。对于表2中的每种参数组合进行10000次模拟，经验性第一类错误率即为在原假设为真的情况下，拒绝原假设的次数除以10000。表3和附表1、附表2展示了层数 $J=2$ 时， $\theta =0.2,0.4,0.6$ 的经验性第一类错误率大小。

从表中可以看出，当层数相同时，三种检验方法的结果都可以稳定在0.04至0.06的范围内，其中 $T_W$ 更为冒进，而 $T_{LR}$ 则更为保守。

为了更直观反应三种不同方法在经验性第一类错误率上的差异以及它们的渐近性质，下面分别在不同层与不同相关系数的情形下随机生成1000组参数并做蒙特卡洛模拟，每一层分别绘制出三种方法在样本量为25,50,100情况下的经验性第一类错误率的箱线图。各层情况下的箱线图如图所示。

Table 3. The empirical TIEs (%) for J = 2

表3. 经验性第一类错误率( $J=2$ )

图1可以看出，当J和m增加时，各检验统计量的经验性的一类错误率显著趋向于0.05。其中， $T_{SC}$ 在各种情况下均表现稳健； $T_W$ 在层数和每组人数较少的情况下表现冒进，当层数和每组人数都较高时，才表现稳

Figure 1. The box plots of empirical TIEs

图1. 不同样本量下第一类错误率的箱线图

Table 4.The empirical powers (%) for J = 2

表4.经验性功效( $J=2$ )

健； $T_{LR}$ 与 $T_{SC}$ 相比较为冒进，总体表现较为稳健。因此，基于第一类错误率的表现，最推荐的检验统计量为 $T_{SC}$ 。

在参数配置相同的备择假设 $H_1$ 下，计算三个统计量的经验性功效。具体结果见表4和附表3、附表4。可以看出 $T_{LR}$ ， $T_W$ 和 $T_{SC}$ 的功效随着J和m的增大而增大，且在样本量较大的情况下表现更为良好。

3.3. 样本量的确定

由检验结果可知，样本量对检验统计量的第一类错误率和功效有显著影响，如果能提前预知在固定的显著性水平下达到期望功效所需的样本总量，就能使后续进行的试验有更稳健的第一类错误率和更好的功效。于是，提出了检验统计量 $T_{LR}$ ， $T_W$ ， $T_{SC}$ 的渐近样本量计算方法。

在 $H_1$ 假设下，检验统计量 $T_{LR}$ ， $T_{LW}$ ， $T_{SC}$ 达到特定功效所需的渐近样本量可由下式计算得出：

$\mathrm{Pr}\left(T_{LR}\ge {\chi }_{1-\alpha }^2\left(J-1\right)|H_{11}\right)=1-\beta$ ，

$\mathrm{Pr}\left(T_W\ge {\chi }_{1-\alpha }^2\left(J-1\right)|H_{11}\right)=1-\beta$ ，

$\mathrm{Pr}\left(T_{SC}\ge {\chi }_{1-\alpha }^2\left(J-1\right)|H_{11}\right)=1-\beta$ 。

由于上述式子不存在显式解，因此可通过迭代的方法对渐进样本量大小进行求解：

样本总量 $N={\displaystyle {\sum }_{j=1}^J{N}_j}$ ， $N_j=m_{++j}^{\left(1\right)}+m_{++j}^{\left(2\right)}$ ， $N_j=r_{j}N$ ， $m_{++j}^{\left(1\right)}=t_{j}N$ ， $m_{++j}^{\left(1\right)}=t_j{N}_j$ ， $j=1,2,\cdots ,J$ 。在 $H_1$ 假设条件下，通过设置 $r_j,t_j,{\pi }_{ij},{\rho }_j$ 的参数值，随机生成样本量N；

基于(1)中生成的样本量N，代入上式中求解得到渐近功效；

若步骤(2)中所得到的渐近功效小于所期望的功效 $1-\beta$ ，则增加样本量N；若所得到的渐近功效大于 $1-\beta$ ，则减小样本量N。最终使得渐近功效接近 $1-\beta$ 。

令 $r_j=\frac{1}{J}$ ， $t_j=0.5$ ， ${\pi }_{1j}=\left(0.3,0.5\right)$ ， $\theta =\left({\theta }_0,{\theta }_j\right)$ ，由于当 $\theta =1$ 时表示一层间的两组治疗效果相同，因此选择 ${\theta }_0=1$ 。利用上述迭代方法，求得 ${\theta }_j$ 或 ${\rho }_j$ 变化时，达到显著性水平 $\alpha =0.05$ 或功效 $1-\beta =0.8$ 所需的渐近样本量(具体结果如表5所示)。

Table 5. Sample size with 80% power

表5. 80%功效下的渐进样本量

4. 实证研究

Mandel等(1982) [12] 进行了一项研究双盲随机对照试验，该试验对214名患有中耳炎的儿童随机分配服用头孢克洛或阿莫西林进行治疗，以比较这两种抗生素的效果。经过14天的治疗试验，共有203例有效的治疗数据，由表6所示。

基于表6中的数据进行实证分析来展示上面所提到的所有检验以及确定样本量的方法，进而得到一定的应用价值。假设 $H_0:{\theta }_1={\theta }_2={\theta }_3\triangleq \theta$ ， $H_1:\theta \ne {\theta }_0$ ，来检验不同年龄下哪种抗生素的治疗效果较好。通过样本量估计，Fisher得分法和牛顿迭代算法，我们可以得到样本量估计值 ${\pi }_{1j}^{\left(0\right)},{\pi }_{2j}^{\left(0\right)},{\rho }_j^{\left(0\right)}$ 全局最大似然估计量 ${\hat{\pi }}_{1j},{\hat{\pi }}_{2j},{\hat{\rho }}_j$ ， $H_0$ 假设下的最大似然估计量 ${\tilde{\pi }}_{1j},{\tilde{\pi }}_{2j},{\rho }_j$ 结果汇总在表7中。

Table 6. The number of children after receiving treatment

表6. 各年龄层患病儿童接受治疗后改善情况

Table 7. Maximum likelihood estimator of each parameter

表7. 各参数的最大似然估计量

基于三类最大似然估计量构造检验统计量，三种检验统计量和对应 值如表8所示。

Table 8. Statistic values and p-value of homogeneity test

表8. 同质性检验下的统计值和p值

结果表明，在原假设 $H_0$ 下，p值均小于0.05，认为拒绝原假设，各层的优势比并不等于1。优势比近似等于0.2389，表示两种抗生素对于同一年龄段的患者的治疗效果不同，头孢克洛对儿童中耳炎有更好的治疗效果。然后根据提出的方法确定样本量。

Tabl 9. Estimated sample size under the desired power

表9. 目标功效下的样本量

再根据所提出的方法进行样本量的确定。考虑 $r_1=r_2=r_3=\frac{1}{3}$ ， $t_1=t_2=0.5$ ，参数估计值与上述相应

的最大似然估计值相同。在期望功效为80%和90%的情况下，表9给出了显著性水平为0.05条件下的近似样本量。结果表明，在同质性假设下， $T_{LR}$ 达到期望功效所需的样本量最少， $T_W$ 达到期望功效所需的样本量最大，因此最推荐使用似然比检验。

根据以上结果，表明头孢克洛对儿童中耳炎有更好的治疗效果，因此更推荐使用头孢克洛治疗儿童中耳炎。而分层效应不明显，说明年龄对儿童中耳炎的治疗影响不大。此外，根据估计，选取312名儿童急性中耳炎患者能够在研究不同治疗方案对治疗效果的影响时，使用可能少的人力成本下达到期望的功效。

5. 总结

本文对于分层的单双边混合配对数据各层的优势比构造检验统计量进行同质性检验并估计渐近样本量的大小。首先，构造基于Donner模型的同质性假设下参数的三种最大似然估计量。构造同质性假设下的最大似然估计量。利用Fisher得分算法迭代得到所需要的最大似然估计量。再通过估计量构造似然比检验统计量( $T_{LR}$ )、Wald检验统计量( $T_W$ )和score检验统计量( $T_{SC}$ )，选择第一类错误率和功效作为判断指标，对三种统计量进行比较。蒙特卡洛模拟结果显示，在同质性假设下， $T_W$ 在三个检验统计量中表现最为冒进； $T_{LR}$ 相对冒进，总体表现较为稳健；而 $T_{SC}$ 在各种情况下均表现稳健，且功效较高，因此最推荐的检验统计量为 $T_{SC}$ 。对于临床试验中样本量的确定，为了达到相同的功效， $T_W$ 所需样本量较小；随着层数J的增加，三个检验统计量所需样本量差异逐渐减小；随着优势比的增加，越来越接近于1，三个检验统计量所需样本量也会逐渐增加。

附录

Table A1. The empirical TIEs (%) for J = 4

附表1. 经验性第一类错误率(J = 4)

Table A2. The empirical TIEs (%) for J = 6

附表2. 经验性第一类错误率(J = 6)

Table A3. The empirical powers (%) for J = 4

附表3. 经验性功效(J = 4)

Table A4. The empirical powers (%) for J = 6

附表4. 经验性功效(J = 6)