部分区间删失数据下时间相依威布尔比例风险模型的参数估计

Parameter estimation using time-dependent Weibull proportional hazards model for survival analysis with partly interval censored data

摘要

目的 针对临床研究中常见的部分区间删失数据，提出构建时间相依威布尔比例风险模型的参数估计问题，同时探讨不同协变量对生存时间的影响。 方法 以威布尔分布作为比例风险模型的基准风险函数，同时在模型中引入时变协变量，建立时间相依威布尔比例风险模型。为了估计模型的参数，采用极大似然方法，并通过优化函数得到参数的估计值。 结果 数值模拟结果表明，在不同样本量及不同参数设置下，精确观测的比例越高，参数估计效果更好，其覆盖率均近似达到理论预期的95%。此外，随着样本量增大，各参数偏差均呈现减小趋势。 结论 将该方法运用到实例数据中进一步验证模型的有效性，相较于仅有精准观测个体的失效时间数据，具有额外的区间删失数据有助于给出有效的回归参数估计。此外，与含时变协变量的Cox模型进行对比，进一步表明采用时间相依威布尔比例风险模型可给出有效的估计结果。

Abstract

Objective To assess the validity and effectiveness of parameter estimation using a time-dependent Weibull proportional hazards model for survival analysis containing partly interval censored data and explore the impact of different covariates on the results of analysis. Methods We established a time-dependent Weibull proportional hazards model using the Weibull distribution as the baseline hazard function of the model which incorporated time-varying covariates. Maximum likelihood estimation was employed to estimate the model parameters, which were obtained by optimization of the likelihood function. Results and Conclusion Numerical simulation results showed that with higher proportions of precise observations across different sample sizes and parameter settings, the proposed model resulted in improved accuracy of parameter estimation with coverage probabilities approximating the theoretical expectation of 95%. As the sample sizes increased, the parameter biases of the model tended to decrease. Experiments with empirical data further validated the effectiveness of the model. Compared with the failure time data for each precisely observed individual, additional interval-censored data helped to obtain more effective estimates of the regression parameters. Comparison with the Cox model that included time-varying covariates further demonstrated the effectiveness of the time-dependent Weibull proportional hazards model for parameter estimation in survival analysis with partly interval censored data.

关键词

部分区间删失数据 / 时间相依协变量 / 威布尔比例风险模型 / 似然函数 / 极大似然估计

Key words

partly interval-censored data / time-dependent covariates / Weibull proportional hazards model / likelihood function / maximum likelihood estimation

在流行病学和生物医学研究中，某些研究对象感兴趣的事件发生时间能够被精准地捕捉到，而由于实验条件等某些客观因素，某些个体的失效时间则无法确切观测到的，只知道其发生在某一时间区间内，此时产生的数据类型为部分区间删失数据，该数据由失效时间的精确数据和区间删失数据组成。在艾滋病治疗方案的临床试验中，调查涉及患者对齐多夫定的耐药性，由于实验成本较高，对患者的评估较少，因此部分个体的耐药性发展是确切观测的，而其余患者的数据则是区间删失的，此时产生的数据为部分区间删失数据^［1］。关于部分区间删失数据，大多学者基于威布尔参数模型^{［2， 3］}、Cox模型、变换模型等经典半参数模型^［4］展开研究，采用极大似然方法求解参数^［5］。在此研究中，心脏病数据集和糖尿病病数据集是较为经典的实例数据集^{［6， 7］}，与完整的、无缺失的数据集相比，删失数据所提供的信息量相对较少，这一特性无疑增加了生存分析研究中相关方法的复杂性，且部分区间删失数据在各领域大量存在，因此围绕该类型数据展开研究具有深远意义。

在临床研究中通常会产生部分区间删失数据，这时需要建立一个合理的模型，以探讨感兴趣事件发生时间与风险因子之间的关系，Cox比例风险模型作为经典的半参数模型被广泛应用到医药统计研究中^［8-11］，该模型假设各风险因素对风险函数的影响相互独立，且风险因素对风险函数的影响是乘性的。例如，Pan^［12］应用Cox模型分析转移性结直肠癌患者在接受帕尼单抗联合FOLFIRI治疗过程中的疗效。上述学者基于部分区间删失数据构建模型时，通常假定协变量是固定的，然而在实际医学疾病数据分析中，往往会发现在随访时间内某些协变量的值是非恒定的，会随时间发生变化，这种协变量可以看作时间相依的协变量^{［13， 14］}。例如，研究发现吸烟对癌症风险存在一定的影响，然而在随访期间，吸烟状况一直在不断变化^［15］，这种变化在分析具体问题时不可忽视，需要在模型构建中加以考虑。目前，对于带有时变协变量的比例风险模型是对经典Cox模型的较为普遍的扩展^［16］。Crowley ^［17］扩展了传统模型，提出带有时变协变量的比例风险模型，时变协变量在一定程度上放宽具有时间固定协变量的Cox模型的比例风险假设，以便该模型更好的拟合数据。Cox模型的一个优势是在拟合右删失数据时可不必处理基准风险函数，但在部分区间删失数据中，不可避免的需要对基准风险函数进行估计。因此，在分析部分区间删失数据时如何处理模型中未知的基准风险函数是具有研究意义的。

在模型中估计非参数函数是较为困难的，威布尔分布形状参数与尺度参数的灵活性在生存分析中发挥极其重要的作用，考虑威布尔分布模型可以在一定程度上简化估计过程^［18］。在删失数据的分析问题中提出合适的模型来处理数据，对分析影响患者生存时间的风险因素是至关重要的。比例风险模型主要关注风险因素对生存时间的影响，而威布尔分布具有较好的灵活性，在拟合删失数据时是生存分析中使用最为广泛的寿命分布^{［19， 20］}。因此本文将它们结合起来构建时间相依威布尔比例风险模型，以更好的拟合部分区间删失数据。

1 方法

1.1 数据及模型

生存分析中得到的失效时间数据往往是删失的，设n_{\mathrm{1}}为能够精确观测到的样本个数，n_{\mathrm{2}}为一般区间删失的样本个数，总样本个数为n=n_{\mathrm{1}}+n_{\mathrm{2}}，其中前n_{\mathrm{1}}个样本的失效时间T_i(i=\mathrm{1,2},\dots ,n_{\mathrm{1}})是能够被精确观测到的，但对剩下的n_{\mathrm{2}}个样本，只知道感兴趣事件发生在时间区间\left.\left(L_i,R_i\right.\right]内，i=n_{\mathrm{1}}+\mathrm{1},n_{\mathrm{1}}+\mathrm{2},\dots ,n。这样就得到了本文的观测数据，即\left\{\left(T_i,Z_i\right),i=\mathrm{1},\dots ,n_{\mathrm{1}}\right\}和\left\{\left(L_i,R_i,Z_i\right),\right.i=\left.n_{\mathrm{1}}+\mathrm{1},\dots ,n\right\}，其中Z_i表示协变量向量。

在生存分析研究中，尽管大多数模型假设协变量是恒定的，但在实际案例应用中，经常会发现部分协变量呈现出随时间变化的特性。以病人血糖水平为例，这一指标与病人的生存时间存在显著的相关性，但其值并非固定，而是随着治疗进展或病情演变而发生变化。针对此类情况，将血糖水平视作时变协变量进行分析，将更为准确地捕捉其与生存时间之间的动态关系。因此，本文构建时间相依威布尔比例风险模型。

h\left(t\left|\mathit{z}\right.\right)=\lambda \gamma {\left(\lambda t\right)}^{\gamma -\mathrm{1}}exp\left\{{\mathit{\alpha }}^{\mathrm{\top }}\mathit{x}\left(t\right)+{\mathit{\beta }}^{\mathrm{\top }}\mathit{v}\right\},（1）

其中，\lambda >\mathrm{0}为尺度参数，\gamma >\mathrm{0}为形状参数，{\mathit{\alpha }}^{\mathrm{\top }}=\left({\alpha }_{\mathrm{1}},\dots ,{\alpha }_p\right)和{\mathit{\beta }}^{\mathrm{\top }}=\left({\beta }_{\mathrm{1}},\dots ,{\beta }_q\right)分别为时间相依协变量和时间独立协变量的回归系数向量，\mathit{z}={\left(\mathit{x}{\left(t\right)}^{\mathrm{\top }},{\mathit{v}}^{\mathrm{\top }}\right)}^{\mathrm{\top }}表示协变量向量，其中\mathit{x}\left(t\right)={\left(x_{\mathrm{1}}\left(t\right),x_{\mathrm{2}}\left(t\right),\dots ,x_p\left(t\right)\right)}^{\mathrm{\top }}为p维时间相依协变量向量，\mathit{v}={\left(v_{\mathrm{1}},v_{\mathrm{2}},\dots ,v_q\right)}^{\mathrm{\top }}为q维时间独立协变量向量。模型（1）对应的生存函数为

S\left(t\left|\mathit{z}\right.\right)=exp\left\{-H\left(t\left|\mathit{z}\right.\right)\right\},（2）

其中，累积风险函数H\left(t\left|\mathit{z}\right.\right)表示为

H\left(t\left|\mathit{z}\right.\right)={\int }_{\mathrm{0}}^t\lambda \gamma {\left(\lambda s\right)}^{\gamma -\mathrm{1}}exp\left\{{\mathit{\alpha }}^{\mathrm{\top }}\mathit{x}\left(s\right)+{\mathit{\beta }}^{\mathrm{\top }}\mathit{v}\right\}ds.

1.2 极大似然估计

极大似然估计作为一种统计方法，在于利用已有的数据，通过最大化似然函数来确定参数的估计值。本文定义参数向量\mathit{\theta }=(\lambda ,\gamma ,{\mathit{\alpha }}^{\mathrm{\top }},{\mathit{\beta }}^{\mathrm{\top }}{)}^{\mathrm{\top }}。对于部分区间删失数据，前n_{\mathrm{1}}个个体的失效时间是可精准观测到的，此时其似然函数可表示为

L_{\mathrm{1}}={\displaystyle \prod _{i=\mathrm{1}}^{n_{\mathrm{1}}}}f(t_i\left|{\mathit{z}}_i\right.)={\displaystyle \prod _{i=\mathrm{1}}^{n_{\mathrm{1}}}}h\left(t_i\left|{\mathit{z}}_i\right.\right)\cdot S\left(t_i\left|{\mathit{z}}_i\right.\right),（3）

其中，h\left(t_i{\left|\mathit{z}\right.}_i\right)与S\left(t_i\left|{\mathit{z}}_i\right.\right)分别为给定协变量时第i个个体的风险函数与生存函数。

对于其余n_{\mathrm{2}}个个体，无法观测到其失效时间，只知道感兴趣事件发生在\left.\left(L_i,R_i\right.\right]时间区间内，此时似然函数如下

L_{\mathrm{2}}={\displaystyle \prod _{i=n_{\mathrm{1}}+\mathrm{1}}^n}\left[S\left(l_i\left|{\mathit{z}}_i\right.\right)-S\left(r_i\left|{\mathit{z}}_i\right.\right)\right],（4）

从而

L_{\mathrm{2}}={\displaystyle \prod _{i=n_{\mathrm{1}}+\mathrm{1}}^n}exp\left\{-{\int }_{\mathrm{0}}^{l_i}\lambda \gamma {\left(\lambda s\right)}^{\gamma -\mathrm{1}}\mathrm{e}\mathrm{x}\mathrm{p}\left\{{\mathit{\alpha }}^{\mathrm{\top }}\mathit{x}\left(s\right)+{\mathit{\beta }}^{\mathrm{\top }}\mathit{v}\right\}\mathrm{d}s\right\}\cdot ｛\mathrm{1}-exp｛-{\int }_{l_i}^{r_i}\lambda \gamma {\left(\lambda s\right)}^{\gamma -\mathrm{1}}exp\left\{{\mathit{\alpha }}^{\mathrm{\top }}\mathit{x}\left(s\right)+{\mathit{\beta }}^{\mathrm{\top }}\mathit{v}\right\}\mathrm{d}s｝｝.（5）

因此，总观测数据的全似然函数为L=L_{\mathrm{1}}\cdot L_{\mathrm{2}}，故而

L={\displaystyle \prod _{i=\mathrm{1}}^n}{\left[f\left(t_i\left|{\mathit{z}}_i\right.\right)\right]}^{{\delta }_i}\cdot {\left[S\left(l_i\left|{\mathit{z}}_i\right.\right)-S\left(r_i\left|{\mathit{z}}_i\right.\right)\right]}^{\mathrm{1}-{\delta }_i},（6）

其中，当{\delta }_i=\mathrm{1}时，此时可观测到的数据结构为\left\{\left(t_i,{\mathit{z}}_i\right),i=\mathrm{1},\dots ,n_{\mathrm{1}}\right\}；当{\delta }_i=\mathrm{0}时，得到的数据结构为\left\{\left(l_i,r_i,{\mathit{z}}_i\right),\right.i=n_{\mathrm{1}}+\mathrm{1},\left.\dots ,n\right\}。

为了求解模型的参数，选择极大化对数似然函数，对数似然函数为：

l\left(\lambda ,\gamma ,\mathit{\alpha },\mathit{\beta }\right)={\displaystyle \sum _{i=\mathrm{1}}^n}{\delta }_i\cdot logf\left(t_i\left|{\mathit{z}}_i\right.\right)+(\mathrm{1}-{\delta }_i)\cdot log\left[S\left(l_i\left|{\mathit{z}}_i\right.\right)-S\left(r_i\left|{\mathit{z}}_i\right.\right)\right].（7）

接下来，将对数似然函数l\left(\lambda ,\gamma ,\mathit{\alpha },\mathit{\beta }\right)关于参数\mathit{\theta }=(\lambda ,\gamma ,{\alpha }_{\mathrm{1}},\dots ,{\alpha }_p,{\beta }_{\mathrm{1}},\dots ,{\beta }_q{)}^{\mathrm{\top }}求偏导，令偏导为0再联立，进而可得参数估计值\widehat{\mathit{\theta }}=(\widehat{\lambda },\widehat{\gamma },{\widehat{\alpha }}_{\mathrm{1}},\dots ,{\widehat{\alpha }}_p,{\widehat{\beta }}_{\mathrm{1}},\dots ,{\widehat{\beta }}_q{)}^{\mathrm{\top }}，即

\left\{\begin{array}{c}{U}_{\lambda }=\frac{\partial l\left(\lambda ,\gamma ,\mathit{\alpha },\mathit{\beta }\right)}{\partial \lambda }=\mathrm{0}\\ U_{\gamma }=\frac{\partial l\left(\lambda ,\gamma ,\mathit{\alpha },\mathit{\beta }\right)}{\partial \gamma }=\mathrm{0}\\ U_{{\alpha }_d}=\frac{\partial l\left(\lambda ,\gamma ,\mathit{\alpha },\mathit{\beta }\right)}{\partial {\alpha }_d}=\mathrm{0}\\ U_{{\beta }_j}=\frac{\partial l\left(\lambda ,\gamma ,\mathit{\alpha },\mathit{\beta }\right)}{\partial {\beta }_j}=\mathrm{0}\end{array}\right.（8）

对数似然函数的Fisher信息矩阵如下

\mathit{I}\left(\widehat{\mathit{\theta }}\right)=-{\left[\begin{array}{cccc}\frac{\partial U\left(\lambda \right)}{\partial \lambda }& \frac{\partial U\left(\lambda \right)}{\partial \gamma }& \frac{\partial U\left(\lambda \right)}{\partial {\alpha }_d}& \frac{\partial U\left(\lambda \right)}{\partial {\beta }_j}\\ \frac{\partial U\left(\gamma \right)}{\partial \lambda }& \frac{\partial U\left(\gamma \right)}{\partial \gamma }& \frac{\partial U\left(\gamma \right)}{\partial {\alpha }_d}& \frac{\partial U\left(\gamma \right)}{\partial {\beta }_j}\\ \frac{\partial U\left({\alpha }_d\right)}{\partial \lambda }& \frac{\partial U\left({\alpha }_d\right)}{\partial \gamma }& \frac{\partial U\left({\alpha }_d\right)}{\partial {\alpha }_d}& \frac{\partial U\left({\alpha }_d\right)}{\partial {\beta }_j}\\ \frac{\partial U\left({\beta }_j\right)}{\partial \lambda }& \frac{\partial U\left({\beta }_j\right)}{\partial \gamma }& \frac{\partial U\left({\beta }_j\right)}{\partial {\alpha }_d}& \frac{\partial U\left({\beta }_j\right)}{\partial {\beta }_j}\end{array}\right]}_{\mathit{\theta }=\widehat{\mathit{\theta }}}

由此，可根据Fisher信息矩阵的逆矩阵来计算参数\widehat{\mathit{\theta }}的渐近方差-协方差矩阵\widehat{\mathit{\Sigma }}，具体公式如下，方差-协方差矩阵\widehat{\mathit{\Sigma }}的对角线元素为各个参数的方差，而非对角线元素为各参数之间的协方差。

\widehat{\mathit{\Sigma }}=\mathit{I}{\left(\widehat{\mathit{\theta }}\right)}^{-\mathrm{1}}.（9）

2 结果

2.1 模拟研究

本文根据部分区间删失数据，分别基于协变量均独立和部分协变量为时变协变量的两种情况下进行模拟研究。对于每个个体的观测时间，设置观测次数服从K~Ber\left(\mathrm{0.75}\right)+\mathrm{1}，当K=\mathrm{1}时，生成观测时间U_{\mathrm{1}}~Unif\left(\mathrm{0,3}\tau /\mathrm{4}\right)，若T\le U_{\mathrm{1}}，则定义失效时间所在区间为\left.\left(L,R\right.\right]=\left.\left(\mathrm{0},U_{\mathrm{1}}\right.\right]，若T>U_{\mathrm{1}}，则定义\left.\left(L,R\right.\right]=\left(U_{\mathrm{1}},\mathrm{\infty }\right)。当K=\mathrm{2}时，两次观测时间分别为U_{\mathrm{1}}~Unif\left(\mathrm{0,3}\tau /\mathrm{4}\right)\mathrm{、}{U}_{\mathrm{2}}~min\left\{\mathrm{0.1}+U_{\mathrm{1}}+Exp\left(\mathrm{1}\right)\tau /\mathrm{2},\tau \right\},若T\le U_{\mathrm{1}}，则失效时间所在时间区间为\left.\left(L,R\right.\right]=\left.\left(\mathrm{0},U_{\mathrm{1}}\right.\right]；若，此时\left.\left(L,R\right.\right]=\left.\left(U_{\mathrm{1}},U_{\mathrm{2}}\right.\right]；若T>U_{\mathrm{2}}，\left.\left(L,R\right.\right]=\left(U_{\mathrm{2}},\mathrm{\infty }\right)。然后生成精准观测结果，当R=\mathrm{\infty }时，设\mathrm{\Delta }=\mathrm{0}；当时，生成\mathrm{\Delta }~Ber\left(pt\right)，此时\mathrm{\Delta }=\mathrm{1}，则假定失效时间T是准确观察到的，pt给出了除右删失观测值外，精确观测到失效时间的比例，本研究考虑pt=\mathrm{0.4,0.75}两种情况下的估计效果，令研究的时间长度\tau =\mathrm{6}。

2.1.1 协变量均与时间独立情况

此部分假设给定协变量x，且失效时间与观测时间相互独立时，采用威布尔比例风险模型拟合部分区间删失数据，此时设协变量均与时间无关，具体如下：

h\left(t\left|x\right.\right)=\lambda \gamma {\left(\lambda t\right)}^{\gamma -\mathrm{1}}exp\left\{{\alpha }_{\mathrm{1}}{x}_{\mathrm{1}}+{\alpha }_{\mathrm{2}}{x}_{\mathrm{2}}+{\alpha }_{\mathrm{3}}{x}_{\mathrm{3}}\right\}（10）

其中，\lambda为尺度参数，\gamma为形状参数，\alpha =({\alpha }_{\mathrm{1}},{\alpha }_{\mathrm{2}},{\alpha }_{\mathrm{3}}{)}^{\mathrm{\top }}为未知的回归参数向量，x={\left(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}}\right)}^{\mathrm{\top }}是与时间相互独立的协变量向量，这里令x_{\mathrm{1}}~Ber\left(\mathrm{0.3}\right)，x_{\mathrm{2}}\sim U\left(-\mathrm{1,1}\right)，x_{\mathrm{3}}\sim Ber\left(\mathrm{0.5}\right)。

此部分设置两组不同参数真值{\mathit{\theta }}_{\mathrm{1}}=(\lambda ,\gamma ,{\alpha }_{\mathrm{1}},{\alpha }_{\mathrm{2}},{\alpha }_{\mathrm{3}}{)}^{\mathrm{\top }}分别为：（1）{(\mathrm{1.1,1.2,0.6,1},\mathrm{1})}^{\mathrm{\top }}；（2）{(\mathrm{0.5,0.8},-\mathrm{0.7,0.4,0.7})}^{\mathrm{\top }},同时考虑pt为0.4、0.75，样本量n为200、400和600，各自重复1000次进行模拟。为考虑参数估计效果，选取偏差（BIAS）、样本标准差（SSE）、理论标准差（ESE）、95%的经验覆盖概率（CP）这四个感兴趣的评价指标（表1，2）。

从上表中可以看出，协变量均独立的情况时，在上述参数真值设置下，随着样本量的增长使得参数估计的偏差趋势性减小。此外，利用正态近似构建的95%置信区间，上述结果显示其覆盖率近似达到理论预期的95%。进一步观察发现，随着精确观测值比例的提升，参数估计的标准差显著降低，与此同时，样本标准差与理论标准差的差异也呈逐步缩小的态势，这与预期一致。

2.1.2 部分协变量与时间相依情况

本部分考虑构建时间相依威布尔比例风险模型，具体如下：

h\left(t\left|\mathit{z}\right.\right)=\lambda \gamma {\left(\lambda t\right)}^{\gamma -\mathrm{1}}exp\left\{\alpha x(t)+{\beta }_{\mathrm{1}}{v}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{v}_{\mathrm{2}}\right\}（11）

其中，v={\left(v_{\mathrm{1}},v_{\mathrm{2}}\right)}^{\mathrm{\top }}是与时间独立的协变量向量，v_{\mathrm{1}}\sim U\left(-\mathrm{1,1}\right)，v_{\mathrm{2}}\sim Ber\left(\mathrm{0.5}\right)，而令x(t)为与时间相依协变量，这里使x(t)=x\cdot logt，其中x\sim Ber(\mathrm{0.4})。

在部分协变量非独立时，即部分协变量为时变协变量的情况下，此部分设置3组不同参数真值{\mathit{\theta }}_{\mathrm{2}}=(\lambda ,\gamma ,\alpha ,{\beta }_{\mathrm{1}},{\beta }_{\mathrm{2}}{)}^{\mathrm{\top }}，分别为：（1）{(\mathrm{1.1,1.2,0.6},-\mathrm{0.2,0.2})}^{\mathrm{\top }}；（2）{\left(\mathrm{0.6},\right.\mathrm{0.8,1},-\mathrm{2},-\mathrm{2})}^{\mathrm{\top }}；（3）{(\mathrm{0.4,0.5},-\mathrm{0.3,0.5,0.7})}^{\mathrm{\top }}，同时考虑pt为0.4、0.75，样本量n为200、400和600，并且各自重复观测1000次进行模拟设置（表3~5）。

在3组不同参数真值的设置下，参数估计的结果均呈现出高度的无偏性，同时随着精确观测值在数据集中占比的增加，参数估计的标准差降低，这一发现与预期相符。进一步地，通过正态近似方法构建的95%置信区间，其实际覆盖率与理论预期的95%极为接近。另外，当样本量逐步增大时，参数估计的偏差呈现出明显的下降趋势。

验证时间相依威布尔比例风险模型的有效性，将其与带有时变协变量的Cox模型的估计结果进行了对比分析。设置感兴趣参数\alpha =\mathrm{0.2},{\beta }_{\mathrm{1}}=\mathrm{1.2},{\beta }_{\mathrm{2}}=-\mathrm{1}，同时考虑pt为0.4、0.75，样本量n为200和400，循环1000次。模拟结果表明，随着样本量的增加，两种模型感兴趣参数估计的偏差和标准差均有减小趋势。此外，相较于带有时变协变量的Cox模型，时间相依威布尔比例风险模型中时变协变量对应的系数的估计值的偏差与标准差均略小，进一步验证所提出的模型在估计时变协变量的回归系数是有效的（表6）。

2.2 实例分析

2.2.1 丹麦糖尿病数据集

在I型糖尿病的研究中发现，高死亡率主要集中在发展为肾衰竭的患者，因此，糖尿病肾病的发展可视为I型糖尿病患者生存的预后标志物。研究人员对1933~1972年转到Steno纪念医院并被确诊的I型糖尿病患者进行了随访，直到其死亡或移民。在进入本研究前或研究结束时，共有731例患者已诊断为糖尿病，这意味着数据集不包含右删失数据。其中有595名患者的糖尿病发病时间有确切记录，这些构成了精确观测数据。对于剩下的136例患者，仅掌握他们最后无病症状时的持续时间和首次出现肾病的持续时间，其构成了区间删失数据。在这731例患者中，有454例男性（X=1），其余的277例为女性（X=0）。针对该糖尿病实例数据，采用所提出的方法，重点关注患者性别对糖尿病的影响。表7展示了可精准观测到失效时间数据的个体和部分区间删失数据两种情况下的相关参数估计值（Estimate）、标准差（ESE）和P值。

当性别作为协变量时，在仅有精准观测数据下性别的回归参数估计值为-0.1439，标准差为0.0838，部分区间删失数据对应的性别回归参数估计值为-0.1293，标准差为0.0777。可以看出相较于仅有可精准观测到感兴趣事件失效时间的个体（595例），具有额外的区间删失数据（136例）有助于给出有效的回归参数估计。除此之外，其P值均小于0.05，因此认为性别协变量具有显著影响，性别对应的参数为负值，则表示相较于女性，男性患有糖尿病肾病的概率更低（表7）。

2.2.2 ACTG175数据集

本实例数据选自R语言speff2trial包的ACTG175数据集，来自艾滋病临床试验组175研究，在CD4细胞计数为200~500/mm³的2139名hiv感染者中，随访AIDS或死亡的复合临床终点的发生情况，以最先发生者为准。对受试者随机指定四种治疗方案，包括单独使用齐多夫定、联合使用齐多夫定和二腺苷、联合使用齐多夫定和扎西他滨以及单独使用二腺苷。在受试者中有1266例参与者在研究前接受了抗逆转录病毒治疗（ART），而873例参与者未接受抗逆转录病毒治疗。由于ART的使用经常导致艾滋病毒耐药性会影响抗逆转录病毒药物的功效，目前尚不清楚耐药性如何影响艾滋病的发生或死亡，因此仅对873例未接受抗逆转录病毒治疗的参与者进行了分析。在这873受试者中有157例是精准观测的，其余受试者仅知道感兴趣事件发生的时间区间。为探讨治疗方案对AIDS或死亡发生的影响，设3种治疗方式（联合使用齐多夫定和二腺苷、联合使用齐多夫定和扎西他滨以及单独使用二腺苷）的治疗指标值为1，仅齐多夫定的治疗指标值为0。CD4细胞计数反映了免疫状态，受试者在首次就诊和第20周、第96周时检查CD4细胞计数。感兴趣的是评估治疗和时间相依协变量lo{g}_{\mathrm{10}}(CD\mathrm{4})与AIDS发生或死亡时间的关系。

采用时间相依威布尔比例风险模型对ACTG175数据集进行拟合时，将治疗方式和CD4细胞计数作为协变量分析。结果显示，威布尔分布的尺度参数估计值为0.2082，形状参数估计值为2.0588，治疗方式和CD4细胞计数对应的回归参数估计值分别为-0.5697和-0.6914，标准差分别为0.1253和0.0466，且P值均小于0.05。因此，根据结果可以发现治疗方式和CD4细胞计数对患者的生存时间具有显著影响。这一结果表明，对于实验开始前未接受抗逆转录病毒治疗的受试者，相较于仅使用齐多夫定，联合治疗方案能显著降低艾滋病相关死亡风险，进一步验证了联合治疗方案的优势。

此外，使用含有时变协变量的Cox模型，将治疗方式视为不随时间变化的协变量，而CD4细胞计数作为时变协变量进行分析。治疗方式和CD4细胞计数对应的回归参数的估计值分别为-0.5751和-0.6043，标准差分别为0.1728和0.1000，治疗方式和CD4细胞计数同样显示具有显著影响，该模型的估计结果与本文模型一致。上述两种不同模型的结果比较进一步表明采用时间相依威布尔比例风险模型拟合ACTG175数据集可给出更有效的估计结果。

3 讨论

在生物学和医学领域，研究者通常聚焦于各风险因子对患者生存期的影响，这有助于为患者制定更有效的治疗方案，其中较为常见的是Cox比例风险回归模型。为了解搭桥手术对患者生存的影响，Guan^［21］ 基于Cox比例风险模型对斯坦福心脏移植数据进行研究。Cox比例风险回归模型的两个优点是它能够纳入时变协变量效应和时变协变量^［22-25］。Thackham ^［26］ 同样基于斯坦福心脏移植数据探讨年龄、搭桥手术以及心脏移植等因素对心脏病患者的生存期的影响，其中是否做过心脏移植手术被视为时变协变量，表示患者是否在特定时间进行了心脏移植。但在上述的研究中不可避免需要对Cox模型中的基准风险函数进行估计。因此，本文考虑将传统Cox比例风险模型中的基准风险函数用威布尔模型替代，进而将对非参数的估计转变为对参数的估计，采用极大似然方法得到模型的参数估计值，在一定程度上降低了计算难度。

本文以部分区间删失数据的统计分析为出发点，通过引入时间相依协变量，将Cox模型扩展至带有时变协变量的范畴，结合时间固定协变量和时间相依协变量充分考虑各风险因子的影响，将Cox比例风险模型中的基准风险函数替换为威布尔模型，建立时间相依威布尔比例风险模型。根据Zhou ^［27］的模拟设置生成部分区间删失数据，模拟研究结果显示，在不同样本量、不同参数设置、以及不同删失比率下，运用时间相依威布尔比例风险模型拟合部分区间删失数据具有较好的效果。此外，将模型应用于丹麦糖尿病数据集^［28］与ACTG175数据集^{［29， 30］}，实例分析的结果进一步表明所提方法不仅切实可行，且在实际应用中取得了显著的效果。

本文所提出的模型中充分探索了各个风险因子对生存函数与风险函数的线性影响，然而在实际问题中有时无法保证所有风险因子均呈现线性效应，该模型忽略了非线性影响情况，从而在实际问题的解释中存在一定的局限性。除此之外，对于部分区间删失数据并不局限于构建Cox比例风险模型，实际上还存在其他生存分析经典模型可供选择，关于部分区间删失数据的分析仍有众多潜在的研究方向值得探讨。

参考文献

[1] Li JQ, Ma J. On hazard-based penalized likelihood estimation of accelerated failure time model with partly interval censoring[J]. Stat Methods Med Res, 2020, 29(12): 3804-17.

[2] Fauzi NA, Elfaki F, Ali MY. Some method on survival analysis via weibull model in the present of partly interval censored: a short review[J]. Int J Comput Sci Net, 2015, 15(4): 48-51.

[3] Mohsin Saeed N, Elfaki FAM. Parametric weibull model based on imputations techniques for partly interval censored data[J]. Austrian J Stat, 2020, 49(3): 30-7.

[4] Jiang JJ, Wang CJ, Pan D, et al. Transformation models with informative partly interval-censored data[J]. Stat Comput, 2023, 34(1): 8.

[5] 董小刚, 彭小草, 蒋京京, 等.部分区间删失数据下广义指数分布的参数估计及应用[J]. 吉林大学学报(理学版), 2022, 60(03): 557-567.

[6] Odell PM, Anderson KM, D'Agostino RB. Maximum likelihood estimation for interval-censored data using a Weibull-based accelerated failure time model[J]. Biometrics, 1992, 48(3): 951-9.

[7] Ramlau-Hansen H, Chr Bang Jespersen N, Kragh Andersen P, et al. Life insurance for insulin-dependent diabetics[J]. Scand Actuar J, 1987, 1987(1/2): 19-36.

[8] Kim JS. Maximum likelihood estimation for the proportional hazards model with partly interval-censored data[J]. J R Stat Soc Ser B Stat Methodol, 2003, 65(2): 489-502.

[9] Elfaki FAM, Azram M, Usman M. Parametric cox's model for partly interval-censored data with application to AIDS studies[J]. Int J Appl Phys Math, 2012: 352-4.

[10] 李纯净, 赵昱榕, 张 淼. 部分区间删失数据下比例风险模型的贝叶斯变量选择[J]. 东北师大学报: 自然科学版, 2023, 55(3): 37-44.

[11] Harrell FE. Regression modeling strategies: with applications to linear models, logistic and ordinal regression, and survival analysis [M]. 2nd ed. Cham, Switzerland: Springer International AG, 2015.

[12] Pan C, Cai B. A Bayesian semiparametric mixed effects proportional hazards model for clustered partly interval-censored data[J]. Stat Model, 2024, 24(5): 459-79.

[13] Fisher LD, Lin DY. Time-dependent covariates in the Cox proportional-hazards regression model[J]. Annu Rev Public Health, 1999, 20: 145-57.

[14] Austin PC. Generating survival times to simulate Cox proportional hazards models with time-varying covariates[J]. Stat Med, 2012, 31(29): 3946-58.

[15] Pinsky PF, Zhu CS, Kramer BS. Lung cancer risk by years since quitting in 30+ pack year smokers[J]. J Med Screen, 2015, 22(3): 151-7.

[16] Webb A, Ma J. Cox models with time-varying covariates and partly-interval censoring-a maximum penalised likelihood approach[J]. Stat Med, 2023, 42(6): 815-33.

[17] Crowley J, Hu M. Covariance analysis of heart transplant survival data[J]. J Am Stat Assoc, 1977, 72(357): 27.

[18] 王淑影, 郭祥道, 李红伟, 等.区间删失数据下Weibull 比例优势模型的参数估计[J]. 山东科技大学学报(自然科学版), 2024, 43(01): 73-81.

[19] Love CE, Guo R. Application of weibull proportional hazards modelling to bad-as-old failure data[J]. Quality & Reliability Eng, 1991, 7(3): 149-57.

[20] 王淑影, 汪 童, 黄 鹤. 融合面板计数数据下相依区间删失数据的参数估计[J]. 系统科学与数学, 2024, 1-18.（网络首发）

[21] Guan Z. Maximum approximate Bernstein likelihood estimation in proportional hazard model for interval-censored data[J]. Stat Med, 2021, 40(3): 758-78.

[22] Austin PC, Fang J, Lee DS. Using fractional polynomials and restricted cubic splines to model non‐proportional hazards or time‐varying covariate effects in the Cox regression model[J]. Stat Med, 2022, 41(3): 612-24.

[23] Xu R, Adak S. Survival analysis with time-varying relative risks: a tree-based approach[J]. Methods Inf Med, 2001, 40(2): 141-7.

[24] Adeleke KA, Abiodun AA, Ipinyomi RA. Semi-parametric non-proportional hazard model with time varying covariate[J]. J Mod App Stat Meth, 2015, 14(2): 68-87.

[25] Therneau TM, Grambsch PM. Modeling Survival Data: Extending the Cox Model[M]. New York, NY: Springer New York, 2000.

[26] Thackham M, Ma J. On maximum likelihood estimation of the semi-parametric Cox model with time-varying covariates[J]. J Appl Stat, 2020, 47(9): 1511-28.

[27] Zhou QN, Sun YQ, Gilbert PB. Semiparametric regression analysis of partly interval-censored failure time data with application to an AIDS clinical trial[J]. Stat Med, 2021, 40(20): 4376-94.

[28] Green A, Borch-Johnsen K, Andersen PK, et al. Relative mortality of type 1 (insulin-dependent) diabetes in Denmark: 1933-1981[J]. Diabetologia, 1985, 28(6): 339-42.

[29] Hammer SM, Katzenstein DA, Hughes MD, et al. A trial comparing nucleoside monotherapy with combination therapy in HIV-infected adults with CD4 cell counts from 200 to 500 per cubic millimeter. AIDS Clinical Trials Group Study 175 Study Team[J]. N Engl J Med, 1996, 335(15): 1081-90.

[30] 刘颖欣, 康 佩, 许 军, 等. 条件推断森林在生存分析中的应用[J]. 南方医科大学学报, 2020, 40(4): 475-82.

基金资助

吉林省自然科学基金优秀青年项目(20230101371JC)

RIGHTS & PERMISSIONS

版权所有©《南方医科大学学报》编辑部2021