聚类复发事件加性模型的分段常数化估计

Estimation of an Additive Rate Model for Clustered Recurrent Events Based on the Piecewise Constant Method

摘要: 研究大型的医疗数据具有重要意义，但大量数据会导致计算复杂度变高，尤其涉及到反复发生的疾病数据时。本文对大型数据的聚类复发事件提出了一个半参数边际加性比率函数模型。在进行参数估计时，随着数据量级的增加而导致未知参数增加，为了提高计算效率，对基准比率函数进行分段常数化。利用估计方程方法给出了模型参数的估计量的表达式。并证明估计量满足相合性和渐近正态性等大样本性质。通过数值模拟，验证了提出的估计方法，模拟结果表明对参数部分的估计结果较好。与对基准比率函数未做分段常数化的方法进行比较，发现分段常数化的方法偏差更小，而且计算时间显著缩短。最后将模型和方法应用到慢性肉芽肿病的数据中，将病人按照医院进行分类，找到影响各医院病人病情复发的显著因素和因素的影响方式。

Abstract: Studying large medical databases is important. However, the use of large databases may introduce computational difficulties, particularly when the event of interest is recurrent. A semiparametric marginal additive rate function model is proposed for large databases clustered recurrence events. In the parameter estimation, the number of baseline rate functions increases as the number of classes increases, leading to an increase in the parameters. The baseline rate functions were re-garded as piecewise constant during the estimation, which improved the computational efficiency. The expressions of the estimate of the model parameters were given using the estimating equation method. The estimators were proved consistent and asymptotic normality. The proposed estima-tion method is verified by numerical simulation, and the simulation results show that the estima-tion results of the parameter part are good, compared with the method without piecewise constant of the risk ratio function. It is found that the bias of the piecewise constant method is smaller, and the calculation time is significantly shortened. Finally, the model and method were applied to the data of chronic granulomatous disease, and the patients were classified by hospital to find out the significant factors affecting patients' disease recurrence in each hospital and the influencing ways of the factors.

1. 引言

在一段时间内，若某个感兴趣事件重复发生，我们称之为复发事件。复发事件在临床试验以及工业领域中很常见，例如癌症的反复发生，反复感染病毒，机器反复发生故障等等。在实际研究中，我们遇到的复发数据可能是来自不同医疗机构的病人，或者是家族遗传病，这时，来自于同一个机构或家庭的病人就具有相似性，我们需要根据这些属性(例如医疗机构，家族)，将研究的多个个体进聚类，得到聚类复发事件数据。

复发事件是近些年来的一个热门研究领域。根据复发事件建模对象的不同，主要有两种研究：一种是对复发事件过程的强度函数建模 [1] [2] ；另一种是对边际均值或比率函数进行建模 [3] 。均值模型和比率模型具有计算简单，真值易估计、稳健性强等优点，因此吸引了大批学者进行研究。有些学者假设协变量效应和基准比率之间是乘积关系的乘性模型，例如Lin [4] 或者Lawless [5] 。相应的，也有学者研究协变量效应和基准比率之间是相加关系的加性模型，例如Liu和Wu [6] 等人。Cai [7] 等人则在平均寿命函数中同时考虑了加性效应和乘性效应。死亡事件也受到了学者的关注，Sun [8] 等人提出了一种特定复发事件均值的半参数方法，该方法结合了死亡事件的加性模型和条件复发事件的加性模型，在这篇论文的基础上，马燕 [9] 所提出的方法则是基于乘性模型。Xu [10] 等人还研究了复发事件和死亡事件的联合尺度变化模型，并用脆弱变量来解释复发事件和死亡事件之间的相关性，Chan和Wang [11] 通过以死亡时间为源头，通过逆向计数的方法来研究复发事件。在估计间歇性观测的时间依赖协变量的比率模型时，Sun等人在2020 [12] 年提出了一种基于逆速率加权和核平滑的估计方法，Lyu Tianmeng等人2021 [13] 年则提出了使用加性–乘性比率的半参数回归模型来解决上诉问题。对被调查的个体可能会经历几种不同类型的多类型复发事件，杜彦斌 [14] 等人提出了一种半参数转移模型，作者利用广义估计方程的思想对参数进行了估计。杨青龙 [15] 用广义半参数风险模型研究了多类型复发事件，并用重采样方法计算估计参数的方差。

对于聚类复发事件的研究主要集中在边际均值或者边际比率函数上。Schaubel和Cai [16] 对聚类复发事件提出了一个半参数乘积比率模型。Liu [17] 推广了该模型，假设来自于一个中心的数据为一个类，不同的类具有不同的基准比率函数，并且考虑了死亡事件，当中心数较大时，作者提出对基准比率函数分段常数化，提高计算效率，但是Liu [17] 只研究了乘积比率模型，没有考虑加性的情况。Liu等人 [18] 在聚类复发事件的比率函数加入了一个起乘积作用的类固定效应项，但该方法的计算效率不高。He [19] 等人提出了一个加性比率函数模型，在模型中加入带加性作用的类固定效应项。此外，Sun [20] 等人还分别考虑了加性比率模型和加性风险回归模型的核加权估计过程。

本论文借鉴了Liu [17] 等人的基准比率函数分段常数化的方法，但他只研究了乘性比率函数模型的情况，我们在此基础上研究加性比率函数模型的估计，并且进一步研究了估计量大样本性质。

本文第一部分给出记算符号，并建立模型和估计方法，第二部分给出估计量的相合性和渐近正态性等大样本理论，第三部分进行数值模拟验证提出的方法，并与基准比率函数不做分段化处理时的结果进行了比较。第四部分分析了一个慢性肉芽肿病的数据，第五部分是论文的总结和展望，渐近性质的证明在附录里。

2. 模型和估计方法

2.1. 符号说明及模型建立

假设有n个观测的个体，将这些研究的个体分成K个类，用 $k\left(k=1,2,\cdots ,K\right)$ 表示类的指标值，用 $n_k$

代表第k个类所包含的样本数，满足 ${\displaystyle \underset{k=1}{\overset{K}{\sum }}{n}_k}=n$ ，同时用 $i\left(i=1,2,\cdots ,n\right)$ 表示第i个个体，对于个体i，我们

用 $G_i$ 表示该个体的类别， $G_i=k$ 表示个体i在第k个类。 $C_i$ 表示个体i的右删失时间，假设删失时间与复发事件过程是独立的。 $Y_i\left(t\right)=I\left(t\le C_i\right)$ 表示个体处于删失的风险，这里的 $I(\cdot )$ 是一个示性函数。用

$N_i^{\ast }\left(t\right)$ 表示在t时刻之前，个体i潜在的事件复发次数， $N_i\left(t\right)={\displaystyle {\int }_0^t{Y}_i\left(s\right)\text{d}{N}_i^{\ast }\left(s\right)}$ 表示在t时刻之前观测到

的复发次数。用 $\tau$ 表示观测截止时间。假设给定协变量的条件下，右删失时间与复发事件过程是独立的。

用 $\text{d}{N}_i^{\ast }\left(t\right)=N_i^{\ast }\left(t\right)-N_i^{\ast }\left(t-\right)$ 表示第i个个体在时刻t的瞬间复发事件发生次数，它只能取1 (有事件发生)或0 (没有事件发生)。在复发事件数据分析中，定义 $E\left\{\text{d}{N}_i^{\ast }\left(t\right)|N_i^{\ast }\left(u\right),Z_i\left(u\right),0\le u<t,G_i=k\right\}$ 为给定历史数据下，t时刻复发事件的风险或强度函数。我们研究复发事件的边际强度函数 $E\left\{\text{d}{N}_i^{\ast }\left(t\right)|Z_i\left(t\right),G_i=k\right\}$ ，也叫做比率函数。对比率函数定义如下的加性模型：

$E\left\{\text{d}{N}_i{}^{*}\left(t\right)|Z_i\left(t\right),G_i=k\right\}=\left({\rho }_{0k}\left(t\right)+Z_i{\left(t\right)}^{{\rm T}}\beta \right)\text{d}t,$ (1)

式中： $Z_i\left(t\right)$ 表示与时间有关的协变量， ${\rho }_{0k}\left(t\right)$ 表示与类别有关的基准比率函数， $\beta$ 是待估计的参数。

定义 $G_{ik}=I\left(G_i=k\right)$ ， $Y_{ik}\left(t\right)=G_{ik}{Y}_i\left(t\right)$ ， $\text{d}{N}_{ik}^{\ast }\left(t\right)=G_{ik}\text{d}{N}_i^{\ast }\left(t\right)$ ， $\text{d}{N}_{ik}\left(t\right)=G_{ik}\text{d}{N}_i\left(t\right)$ ，在右删失与复发事件独立的假设下，有 $E\left\{\text{d}{N}_i^{\ast }\left(t\right)|Z_i\left(t\right),G_i=k,Y_i\left(t\right)=1\right\}=E\left\{\text{d}{N}_i\left(t\right)|Z_i\left(t\right),G_i=k\right\}$ ，结合

$N_i\left(t\right)={\displaystyle {\int }_0^t{Y}_i\left(s\right)\text{d}{N}_i^{\ast }\left(s\right)}$ ，模型(1)变成

$E\left\{\text{d}{N}_{ik}\left(t\right)|Z_i\left(t\right),Y_{ik}\left(t\right)\right\}=Y_{ik}\left(t\right)\left({\rho }_{0k}\left(t\right)+Z_i{\left(t\right)}^{{\rm T}}\beta \right)\text{d}t.$ (2)

2.2. 估计方法

我们利用分段常数化的思想，将 $\left[0,\tau \right]$ 切割成L个区间， $a_0<a_1<\cdots <a_L$ 是L个区间的切割点。设 ${\rho }_{0k}\left(t\right)$ 在区间 $\left(a_{l-1},a_l\right]$ 上为一个常值函数，记成 ${\rho }_{kl}$ 。同样的，将 $z_i\left(t\right)$ 进行分段常值化，用 $z_{il}$ 表示 $z_i\left(t\right)$ 在

时间区间 $\left(a_{l-1},a_l\right]$ 上的值。记 $d_{ikl}={\displaystyle {\int }_{a_{l-1}}^{a_l}\text{d}{N}_{ik}\left(t\right)}$ ， $t_{ikl}={\displaystyle {\int }_{a_{l-1}}^{a_l}{Y}_{ik}\left(t\right)\text{d}t}$ ，根据以上的记法，模型(2)以及估计方

程的思想，得到如下的估计方程：

${\displaystyle \underset{i=1}{\overset{n}{\sum }}\left[d_{ikl}-t_{ikl}\left({\rho }_{kl}+{\beta }^{\text{T}}{Z}_{il}\right)\right]}=0,l=1,2,\cdots ,L$ (3)

${\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}{\displaystyle {\int }_{a_{l-1}}^{a_l}{Z}_{il}\left[\text{d}{N}_{ik}\left(t\right)-Y_{ik}\left(t\right)\cdot \left({\rho }_{kl}+Z_{il}^{\text{T}}\beta \right)\text{d}t\right]}}}}=0$ (4)

由方程(3)可得

${\rho }_{kl}=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{d}_{ikl}}-{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}{\beta }^{\text{T}}{Z}_{il}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}},l=1,2,\cdots ,L$ (5)

将(5)式代入方程(4)，可得

$\hat{\beta }={\left[{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}{t}_{ikl}{\left(Z_{il}-{\bar{Z}}_{kl}\right)}^{\otimes 2}}}}\right]}^{-1}\cdot {\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}\left(Z_{il}-{\bar{Z}}_{kl}\right)}}}\cdot d_{ikl},$ (6)

其中 ${\bar{Z}}_{kl}=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}{Z}_{il}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}}$ 。在上式中，对于每一个向量 $\nu$ ， $v^{\otimes 2}=v{v}^{\prime }$ 。将 $\beta$ 的估计量，代入到(5)式，得到 ${\rho }_{kl}$ 的估计结果如下：

${\hat{\rho }}_{kl}\left(\hat{\beta }\right)=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{d}_{ikl}}-{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}{\hat{\beta }}^{\text{T}}{Z}_{il}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}},$ (7)

根据分段常值化的假定，我们可以进一步求出累积基准比率函数 ${\mu }_{0k}\left(t\right)={\displaystyle {\int }_0^t{\rho }_{0k}\left(u\right)\text{d}u}$ 的估计：

${\hat{\mu }}_{0k}\left(t;\hat{\beta }\right)={\displaystyle \underset{l=1}{\overset{L}{\sum }}{\hat{\rho }}_{kl}\left(\hat{\beta }\right)\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)}$ (8)

式中： $a\wedge b=\min \left(a,b\right)$ ，表示两个数中的较小者。

3. 渐近性质

记 ${\bar{z}}_{kl}=\frac{E\left(t_{ikl}{Z}_{il}\right)}{E\left(t_{ikl}\right)}$ 是 ${\bar{Z}}_{kl}=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}{Z}_{il}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}}$ 的极限， ${\bar{t}}_{kl}=E\left(t_{ikl}\right)$ 是 $\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}$ 的极限。另外记 $\text{d}{M}_{ikl}=d_{ikl}-t_{ikl}\left({\rho }_{kl}+{Z^{\prime }}_{il}{\beta }_0\right)$ 。

为了研究估计量的渐近性质，我们需要以下正则性条件：

a) $P\left(C_{ij}\ge \tau \right)>0$ 对于一切 $\tau >0$ 都成立；

b) $E\left(t_{ikl}\right)>0$ ，对于所有的 $l=1,2,\cdots ,L$ 均成立。(为了保证 ${\bar{z}}_{kl}$ 的分母不等于0)；

c) $d_{ikl},Z_{il},l=1,\cdots ,L$ ，都是有界的；

d) A是正定矩阵，这里 $A=E\left[{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}{t}_{ikl}{\left(Z_{il}-{\bar{z}}_{kl}\right)}^{\otimes 2}}}\right]$ 。

3.1. 定理一

在上述正则化条件下， $\hat{\beta }$ 是一个强相合估计。收敛关系 $\hat{\beta }\to {\beta }_0$ 几乎处处成立。而且 $\sqrt{n}\left(\hat{\beta }-{\beta }_0\right)$ 收敛到一个正态分布，这个正态分布的均值是0，方差是 $A^{-1}\Omega A^{-1}$ ，可以由 ${\hat{A}}^{-1}\hat{\Omega }{\hat{A}}^{-1}$ 估计，这里

${\eta }_{\text{i}}={\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}\left(Z_{il}-{\bar{z}}_{kl}\right)\text{d}{M}_{ikl}}}$ 。

3.2. 定理二

在上述条件下，对于任意给定的k, $k=1,2,\cdots ,K$ ， ${\hat{\mu }}_{0k}\left(t,\hat{\beta }\right)$ 关于 $t\in \left(0,\tau \right]$ 几乎处处一致收敛于 ${\mu }_0\left(t\right)$ ， $\sqrt{n}\left[{\hat{\mu }}_{0k}\left(t,\hat{\beta }\right)-{\mu }_0\left(t\right)\right]$ 弱收敛于零均值的高斯过程，其在 $\left(s,t\right)$ 处的协方差函数为

${\Gamma }_k\left(s,t\right)={\displaystyle \underset{j=1}{\overset{n}{\sum }}E\left({\varphi }_{jk}\left(s\right){\varphi }_{jk}\left(t\right)\right)}$ ，这里

${\hat{\varphi }}_{jk}\left(t\right)=\left(-\frac{1}{\sqrt{n}}{\displaystyle \underset{l=1}{\overset{L}{\sum }}{{\bar{z}}^{\prime }}_{kl}\left(A^{-1}{\eta }_j\right)}\right)\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)+\sqrt{n}{\displaystyle \underset{l=1}{\overset{L}{\sum }}\frac{\text{d}{M}_{jkl}}{{\bar{t}}_{kl}}\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)}$

4. 数值模拟

4.1. 数据生成

在本节中，我们按照模型随机生成有限的样本，对参数进行模拟研究。首先，我们假设n个个体来自于K个类( $K=10,50$ )，每个类包含 $n_k$ 个样本( $n_k=2,5,10,20$ )，总个体 $n=K\times n_k$ ，对每一个个体 $i\left(i=1,2,\cdots ,n\right)$ 我们生成两个协变量 $Z_{1i},Z_{2i}$ ，对于第k类生成一个脆弱变量 $W_k\left(k=1,2,\cdots ,K\right)$ ，假设 $Z_{1i}$ 来自区间(0, 1)上的均匀分布， $Z_{2i}$ 来自参数为0.5的两点分布， $W_k$ 来自于均值和方差都为1的伽玛分布，并假设来自同一个类的个体共用一个相同的 $W_k$ 。假设删失时间C服从区间(u, v]上的均匀分布，u在[0, 0.5]之间变化，v在[2, 6]之间变化，试验终止时间 $\tau =3$ ，这些设置保证复发事件的平均发生次数在2~3次之间。我们从下列模型中生成复发事件：

$E\left\{\text{d}{N}_i^{\ast }\left(t\right)|Z_i\left(t\right),G_i=k,W_k\right\}=\left({\rho }_{0k}\left(t\right)+Z_i{\left(t\right)}^{{\rm T}}{\beta }_0+W_k\right)\text{d}t$ (9)

我们分别考虑 ${\rho }_{0k}\left(t\right)=1$ 、 ${\rho }_{0k}\left(t\right)=t$ 两种情况， ${\beta }_0={\left(0.5,0.5\right)}^{\prime }$ ，为了对基准比率函数实行分段常数化，将事件发生时间的整个区间 $\left(0,\tau \right]$ 分别分成 $L=10$ 或 $L=15$ 个小区间，在每一个小区间上 ${\rho }_{0k}\left(t\right)$ 为一个常

数。我们根据提出的方法来估计模型中的回归参数 $\beta$ 和累积基准比率函数 ${\mu }_{0k}\left(t\right)={\displaystyle {\int }_0^t{\rho }_{0k}\left(u\right)\text{d}u}$ 。

4.2. 模拟结果

$\beta$ 的估计结果见表1和见表2。表中BIAS表示 $\beta$ 的估计值 $\hat{\beta }$ 和真值 ${\beta }_0$ 之间的偏差在1000次中的平均值，ASE表示 $\hat{\beta }$ 的计算的标准差估计在1000次模拟中的平均值，ASD表示1000次模拟中所得到的 $\hat{\beta }$ 的观测值的样本标准差，CP表示在正态分布下，真值 ${\beta }_0$ 的95%置信区间内，覆盖到的 $\hat{\beta }$ 的百分比。模拟结果可以看出：估计值和真值的偏差都在0.05以内，样本标准差和估计量的标准差估计比较接近， ${\beta }_0$ 的95%置信区间内，覆盖真值的比例在95%附近，估计效果较好。

Table 1. The parameter estimation results of repeated simulations 1000 times when μ 0 k = t

表1. ${\mu }_{0k}=t$ 时，重复模拟1000次的参数估计结果

Table 2. The parameter estimation results of repeated simulations 1000 times when μ 0 k = 0.5 t 2

表2. ${\mu }_{0k}=0.5t^2$ 时，重复模拟1000次的参数估计结果

另外，我们也绘制了 $K=50$ ， $n_k=10$ ， $L=15$ 时， ${\mu }_{01}\left(t\right)$ 在不同时间点处的估计曲线以及置信度为95%的置信带的估计，见图1和图2。 ${\mu }_{01}\left(t\right)$ 的估计的95%置信带内基本上包含了真实的 ${\mu }_{01}\left(t\right)$ 。图1和图2如下所示。

Figure 1. Estimate when ${\mu }_{01}\left(t\right)=t$

图1. ${\mu }_{01}\left(t\right)=t$ 时的估计

Figure 2. Estimate when ${\mu }_{01}\left(t\right)=0.5t^2$

图2. ${\mu }_{01}\left(t\right)=0.5t^2$ 时的估计

在评价加入分段常数化的计算效率的时候，我们使用联想YOGA14S笔记本电脑，该电脑基于AMD R7-5800H处理器。在设置模型参数时，我们设置 $n=800$ ， $\beta =\left[0.5,0.5\right]$ ， ${\rho }_{0k}\left(t\right)=t$ ，用本文的分段常数化方法与不加入分段常数化的传统方法进行了比较。在不进行分段常数化时，根据计算， $\beta$ 的估计表达式为：

$\hat{\beta }={\left({\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle {\int }_0^{\tau }{Y}_{ik}\left(t\right){\left(Z_i\left(t\right)-{\bar{Z}}_k\left(t\right)\right)}^{\otimes 2}\text{d}t}}}\right)}^{-1}\cdot {\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle {\int }_0^{\tau }\left(Z_i\left(t\right)-{\bar{Z}}_k\left(t\right)\right)\text{d}{N}_{ik}\left(t\right)}}}$ (10)

其中 ${\bar{Z}}_k\left(t\right)=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{Z}_i\left(t\right)Y_{ik}\left(t\right)}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{Y}_{ik}\left(t\right)}}$ ，相应的 ${\bar{z}}_k\left(t\right)=\frac{E\left(Z_i\left(t\right)Y_{ik}\left(t\right)\right)}{E\left(Y_{ik}\left(t\right)\right)}$ 。

同理，此时也可以计算与标准差 $\hat{\beta }$ 有关的值，分别如下：

$A=E\left({\displaystyle {\int }_0^{\tau }{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{Y}_{ik}\left(t\right){\left(Z_i\left(t\right)-{\bar{z}}_k\left(t\right)\right)}^{\otimes 2}\text{d}t}}}\right)$ (11)

${\eta }_i={\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle {\int }_0^{\tau }\left(Z_i\left(t\right)-{\bar{z}}_k\left(t\right)\right)\text{d}{M}_{ik}\left(t\right)}}$ (12)

比较的结果在表3中，从结果中发现使用传统的估计方法，估计量的偏差大于分段常数化的偏差结果。无论使用哪种方法，在估计方差时，样本观测的标准差和模拟计算的结果均十分接近。

两种方法的差别主要体现在计算时间上，对于恒定的个体 $n=800$ ，分别设置 $K=20,40,80$ 计算500次模拟结果，结果记录在表3中。可以看到，用分段常数化的方法运行时间显著低于传统方法，传统方法所用的时间是分段常数的4倍以上；并且随着K的增大，这一比值也在不断增大。由于计算机性能的限制，包含更多个体的实验并没有展示。但分段常数化方法偏差更小，计算速度更快是毋庸置疑的，相较于传统方法优势很明显。

了解释分段常数化计算效率更高的原因，假设所有个体n的复发次数之和为 $N_r$ ，当我们使用分段常数化的方法计算文中的 $d_{ikl}$ ， $t_{ikl}$ 和 ${\rho }_{kl}$ 时，我们将时间区间一共分成了L段，这里的L由我们设置，此时我们设置的L一般远小于总个体n的复发次数之和 $N_r$ 。但当我们使用传统方法时，我们需要计算所有个体相邻两个发病时间点的 $d_{ikl}$ ， $t_{ikl}$ 和 ${\rho }_{kl}$ ，此时的 $L=N_r+1$ 。随着L的增大。计算所涉及的矩阵维度也随之增大，计算的数据量同样增大，因此计算时间大幅增加，每次计算都带有一定的误差，随着计算量级的增大，计算次数变多，最后所得到的偏差也随之增大。

Table 3. Computational time comparison results

表3. 运算时间比较结果

5. 实例分析

在这部分，我们分析了一个慢性肉芽肿病(CGD)数据 [16] ，这种疾病是免疫功能方面的遗传性疾病，大多发生在婴幼儿中。数据集包括来自不同医院的(13个医院)128个的病人，记录了他们的发病的间隔时间，是否接受Gamma干扰素治疗，基因信息，年龄，身高，体重，性别，是否删失，是否用过抗生素等信息。根据以前学者的研究 [4] ，我们主要关心患者年龄和接受Gamma干扰素治疗对疾病的重复感染比率函数的影响。用协变量 $Z_{i1}$ 表示是否使用Gamma干扰素， $Z_{i1}=1$ 表示使用Gamma干扰素， $Z_{i1}=2$ 表示使用安慰

剂。用 $Z_{i2}$ 表示年龄，这里我们对年龄做了规范化变换 $Z_{i2}^{\ast }=\frac{Z_{i2}-\underset{i=1,\cdots ,128}{\min }\left(Z_{i2}\right)}{\underset{i=1,\cdots ,128}{\max }\left(Z_{i2}\right)-\underset{i=1,\cdots ,128}{\min }\left(Z_{i2}\right)}$ 。假设重复感染的比

率函数满足模型：

$E\left\{\text{d}{N}_i^{\ast }\left(t\right)|Z_i,G_i=k\right\}=\left({\rho }_{0k}\left(t\right)+Z_i^{\text{T}}\beta \right)\text{d}t$ (13)

应用提出的估计方法对参数 $\beta$ 进行估计， $\beta$ 的估计结果见表4中。从表四的结果可以看出，协变量 $Z_{i1}$ 的估计系数 ${\beta }_1$ 为0.25%，且该系数的检验P值小于0.0001，说明系数是显著的，因此说明接受Gamma干扰素治疗能够显著降低重复感染的比率函数。且固定年龄时，使用Gamma干扰素治疗的患者比使用安慰剂的患者平均降低0.25%。协变量 $Z_{i2}$ 的估计系数 ${\beta }_2$ 为−0.14%，但是检验的P值却大于0.05，所以年龄的影响不显著。

Table 4. Parameter estimation of CGD data in model (1)

表4. CGD数据在模型(1)下的参数估计

为了对模型进行检验，我们绘制了残差 $r_{ik}={\displaystyle {\int }_0^{\tau }\text{d}{N}_{ik}\left(t\right)-Y_{ik}\left(t\right)\left({\rho }_{0k}\left(t\right)+Z_i{\left(t\right)}^{{\rm T}}\beta \right)\text{d}t}$ 对观测个体序号i

的残差图，见图3。在图3中，横轴表示个体序号，纵轴表示该个体在模型下的残差值。我们发现残差值在[−2, 2]之间随机变化。因此，没有理由拒绝提出的模型。

Figure 3. Scatterplot of residuals vs. individuals under additive ratio model

图3. 加性比率模型下残差对个体的散点图

6. 总结与展望

本文提出了一个聚类复发事件的加性比率模型，该模型将具有相似特征或来源的个体进行了分类，应用估计方程和非参部分分段常数化的方法对参数进行了估计，证明了估计量的相合性和渐近正态性。在数值模拟中获得了较好的估计效果，在与传统方法进行对比时，计算效率远远高于传统方法。并且对慢性肉芽肿病数据进行了分析，得到了与实际情况相吻合的结论。我们提出的模型和方法具有一定的应用前景，特别是可以应用到多中心医疗数据，家族遗传病医疗数据，以及社区医疗数据中，同时，我们的模型计算效率高，处理大规模医疗数据非常得心应手。

属于同一个类别的个体可能具有相关性，未来我们会考虑在模型中加入一个随机效应项，表示同一类个体之间的关系，那就需要对随机效应项的分布做一些假设，对分布中包含的参数进行估计。另外复发事件往往伴随着死亡事件的发生，我们会进一步研究这两类事件之间的联合建模，比如考虑加入随机效应项或者考虑对反向均值函数进行建模。加性模型可以推广到更加灵活的模型上，比如加性乘积模型，变系数模型，转移函数模型，部分线性模型等等，这些模型会用到一些非参数估计方法，统计推断会更复杂一些，更具有挑战性，适用范围也更广。

附录

定理一证明

因为 $\begin{array}{c}\hat{\beta }-{\beta }_0=\hat{\beta }-{\hat{A}}^{-1}\hat{A}{\beta }_0=\frac{1}{n}{\hat{A}}^{-1}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}\left(Z_{il}-{\bar{Z}}_{kl}\right)\left[d_{ikl}-t_{ikl}{\left(Z_{il}-{\bar{Z}}_{kl}\right)}^{\prime }{\beta }_0\right]}}}\\ =\frac{1}{n}{\hat{A}}^{-1}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}\left(Z_{il}-{\bar{Z}}_{kl}\right)\text{d}{M}_{ikl}}}}\end{array}$ 。

对于任意有限的 $k,l$ $E\left(\text{d}{M}_{ikl}|Z_{il}\right)=0$ ， $E\left(Z_{il}\text{d}{M}_{ikl}|Z_{il}\right)=Z_{il}E\left(\text{d}{M}_{ikl}|Z_{il}\right)=0$ ， $E\left(\left(Z_{il}-{\bar{z}}_{kl}\right)\text{d}{M}_{ikl}|Z_{il}\right)=0$ 。由大数定理可知， $\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}\left(Z_{il}-{\bar{Z}}_{kl}\right)\text{d}{M}_{ikl}}}}\stackrel{P}{\to }0$ ，所以，由 $\hat{A}$ 的收敛性以及矩阵A的正定性，可以得到 $\hat{\beta }-{\beta }_0=\frac{1}{n}{\hat{A}}^{-1}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}\left(Z_{il}-{\bar{Z}}_{kl}\right)\text{d}{M}_{ikl}}}}\stackrel{p}{\to }0$ ，相合性得证。

下面证明渐近正态性，

$\begin{array}{c}\sqrt{n}\left(\hat{\beta }-{\beta }_0\right)={\hat{A}}^{-1}\frac{1}{\sqrt{n}}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}\left(Z_{il}-{\bar{Z}}_{kl}\right)\text{d}{M}_{ikl}}}}+o_p\left(1\right)\\ ={\hat{A}}^{-1}\frac{1}{\sqrt{n}}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}\left(Z_{il}-{\bar{z}}_{kl}\right)\text{d}{M}_{ikl}}}}-{\hat{A}}^{-1}\frac{1}{\sqrt{n}}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}\left({\bar{Z}}_{kl}-{\bar{z}}_{kl}\right)\text{d}{M}_{ikl}}}}+o_p\left(1\right)\\ ={\hat{A}}^{-1}\cdot \frac{1}{\sqrt{n}}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\eta }_i}+o_p(\; 1\; )\end{array}$

其中 ${\eta }_i={\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}\left(Z_{il}-{\bar{z}}_{kl}\right)\text{d}{M}_{ikl}}}$ 。这里的 ${\eta }_i$ 是均值为0，且独立同分布的随机变量，所以

$Var\left(\sqrt{n}\left(\hat{\beta }-{\beta }_0\right)\right)\stackrel{p}{\to }{A}^{-1}Var\left(\frac{1}{\sqrt{n}}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\eta }_i}\right){A^{\prime }}^{-1}=A^{-1}Var\left({\eta }_i\right){A^{\prime }}^{-1}=A^{-1}E\left({\eta }_i^{\otimes 2}\right){A^{\prime }}^{-1}=A^{-1}\Omega {A^{\prime }}^{-1}$ 。

应用中心极限定理，可得 $n^{\frac{1}{2}}\left(\hat{\beta }-{\beta }_0\right)$ 的分布收敛到一个正态分布，这个正态分布的均值是0，方差是 $A^{-1}\Omega A^{-1}$ 。 $\Omega =E\left({\eta }_i^{\otimes 2}\right)=E\left[{\left({\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}\left(Z_{il}-{\bar{z}}_{kl}\right)\text{d}{M}_{ikl}}}\right)}^{\otimes 2}\right]$

定理一证明

因为 ${\mu }_{0k}\left(t,\hat{\beta }\right)-{\mu }_{0k}\left(t\right)={\hat{\mu }}_{0k}\left(t,\hat{\beta }\right)-{\hat{\mu }}_{0k}\left(t,{\beta }_0\right)+{\hat{\mu }}_{0k}\left(t,{\beta }_0\right)-{\mu }_{0k}\left(t\right)$ ，

$\begin{array}{c}{\hat{\mu }}_{0k}\left(t;\hat{\beta }\right)-{\hat{\mu }}_{0k}\left(t;{\beta }_0\right)={\displaystyle \underset{l=1}{\overset{L}{\sum }}\left({\hat{\rho }}_{kl}\left(\hat{\beta }\right)-{\hat{\rho }}_{kl}\left({\beta }_0\right)\right)\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)}\\ ={\displaystyle \underset{l=1}{\overset{L}{\sum }}\left[\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{d}_{ikl}}-{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}{\hat{\beta }}^{\text{T}}{Z}_{il}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}}-\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{d}_{ikl}}-{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}{\hat{\beta }}_0^{\text{T}}{Z}_{il}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}}\right]\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)}\\ =-{\displaystyle \underset{l=1}{\overset{L}{\sum }}\left[\left({\hat{\beta }}^{\text{T}}-{\hat{\beta }}_0^{\text{T}}\right)\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}{Z}_{il}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}}\right]\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)}\\ =-{\displaystyle \underset{l=1}{\overset{L}{\sum }}{\left(\hat{\beta }-{\beta }_0\right)}^{\prime }{\bar{Z}}_{kl}\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)}\end{array}$

$\begin{array}{c}{\hat{\mu }}_{0k}\left(t,{\beta }_0\right)-{\mu }_{0k}\left(t\right)={\displaystyle \underset{l=1}{\overset{L}{\sum }}\left[\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{d}_{ikl}}-{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}{\hat{\beta }}^{\text{T}}{Z}_{il}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}}-{\rho }_{0k}\left(t\right)\right]}\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)\\ ={\displaystyle \underset{l=1}{\overset{L}{\sum }}\left[\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{d}_{ikl}}-{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}{\hat{\beta }}^{\text{T}}{Z}_{il}}-{\rho }_{0k}\left(t\right){\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}}\right]}\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)\\ ={\displaystyle \underset{l=1}{\overset{L}{\sum }}\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{M}_{ikl}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}}}\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)\\ ={\displaystyle \underset{i=1}{\overset{n}{\sum }}\frac{{\displaystyle \underset{l=1}{\overset{L}{\sum }}{M}_{ikl}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}}}\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)\end{array}$

由上一个定理的证明我们得到 $\hat{\beta }-{\beta }_0\stackrel{p}{\to }0$ ，所以 ${\hat{\mu }}_{0k}\left(t;\hat{\beta }\right)-{\hat{\mu }}_{0k}\left(t;{\beta }_0\right)\stackrel{p}{\to }0$ ，同样的 ${\hat{\mu }}_{0k}\left(t,{\beta }_0\right)-{\mu }_{0k}\left(t\right)=\frac{\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{M}_{ikl}}}{\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}}\stackrel{p}{\to }0$ 所以得到 ${\hat{\mu }}_{0k}\left(t,\hat{\beta }\right)-{\mu }_{0k}\left(t\right)\stackrel{p}{\to }0$ ，所以相合性得证。

下面证明渐近正态性。

$\begin{array}{l}\sqrt{n}\left[{\mu }_k\left(t,\hat{\beta }\right)-{\mu }_{0k}\left(t\right)\right]=\sqrt{n}\left[{\hat{\mu }}_{0k}\left(t,\hat{\beta }\right)-{\hat{\mu }}_k\left(t,{\beta }_0\right)+{\hat{\mu }}_k\left(t,{\beta }_0\right)-{\mu }_{0k}\left(t\right)\right]\\ ={\displaystyle \underset{l=1}{\overset{L}{\sum }}\left[-\sqrt{n}{\left({\beta }_0-\hat{\beta }\right)}^{\prime }{\bar{Z}}_{kl}+\sqrt{n}\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}\text{d}{M}_{ikl}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}}\right]}\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)\\ ={\displaystyle \underset{l=1}{\overset{L}{\sum }}\left[-\frac{1}{\sqrt{n}}\cdot {\left({\hat{A}}^{-1}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\eta }_i}\right)}^{\prime }{\bar{Z}}_{kl}+o_P\left(1\right)+\sqrt{n}\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}\text{d}{M}_{ikl}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}}\right]}\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)+o_p(\; 1\; )\end{array}$

交换连加顺序，并应用大数定律，上式可以转化为如下关系：

$\begin{array}{l}\sqrt{n}\left({\mu }_{0k}\left(t,\hat{\beta }\right)-{\mu }_{0k}\left(t\right)\right)\\ ={\displaystyle \underset{i=1}{\overset{n}{\sum }}\left[\left(-\frac{1}{\sqrt{n}}{\left({\hat{A}}^{-1}{\displaystyle \underset{l=1}{\overset{L}{\sum }}{\eta }_i}\right)}^{\prime }{\bar{Z}}_{kl}\right)\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)+\frac{1}{\sqrt{n}}{\displaystyle \underset{l=1}{\overset{L}{\sum }}\frac{\text{d}{M}_{ikl}}{\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{t}_{ikl}}}}\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)\right]}+o_p\left(1\right)\\ ={\displaystyle \underset{i=1}{\overset{n}{\sum }}\left[\left(-\frac{1}{\sqrt{n}}{\displaystyle \underset{l=1}{\overset{L}{\sum }}{\left(A^{-1}{\eta }_i\right)}^{\prime }{\bar{z}}_{kl}}\right)\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)+\frac{1}{\sqrt{n}}{\displaystyle \underset{l=1}{\overset{L}{\sum }}\frac{\text{d}{M}_{ikl}}{{\bar{t}}_{kl}}}\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)\right]}+o_p\left(1\right)\\ ={\displaystyle \underset{i=1}{\overset{n}{\sum }}{\varphi }_{ik}}+o_p(\; 1\; )\end{array}$

这里， ${\varphi }_{ik}\left(t\right)=\left(-\frac{1}{\sqrt{n}}{\displaystyle \underset{l=1}{\overset{L}{\sum }}{\left(A^{-1}{\eta }_i\right)}^{\prime }{\bar{z}}_{kl}}\right)\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)+\frac{1}{\sqrt{n}}{\displaystyle \underset{l=1}{\overset{L}{\sum }}\frac{\text{d}{M}_{ikl}}{{\bar{t}}_{kl}}}\left({\alpha }_l\wedge t-{\alpha }_{l-1}\wedge t\right)$ 。

因为 $E\left({\eta }_i\right)=0$ ，所以 $E\left({\displaystyle \underset{l=1}{\overset{L}{\sum }}{\eta }_i}\right)=0$ ，则等式右边的第一项的期望为0，同样的 $E\left(\text{d}{M}_{ikl}\right)=0$ ，所以 $E\left({\displaystyle \underset{l=1}{\overset{L}{\sum }}\text{d}{M}_{ikl}}\right)=0$ ，则等式右边第二项的期望为0，所以有 $E\left({\varphi }_{ik}\right)=0$ 。又因为，对于固定的t，上式 $\sqrt{n}\left[{\hat{\mu }}_{0k}\left(t,\hat{\beta }\right)-{\mu }_{0k}\left(t\right)\right]$ 渐近于均值为零的独立同分布的随机变量之和，则由多元中心极限定理知， $\sqrt{n}\left[{\hat{\mu }}_{0k}\left(t,\hat{\beta }\right)-{\mu }_{0k}\left(t\right)\right]$ 以有限维分布弱收敛到均值为零的高斯过程，其在 $\left(s,t\right)$ 处的协方差函数为 ${\Gamma }_k\left(s,t\right)={\displaystyle \underset{j=1}{\overset{n}{\sum }}E\left({\varphi }_{jk}\left(s\right){\varphi }_{jk}\left(t\right)\right)}$ 。