带有多终止事件的复发事件数据可加可乘比率模型

Additive and Multiplicative Ratio Modeling of Recurrent Event Data with Multiple Termination Events

摘要: 在已有的联合建模方法基础上，深入分析了具有两类终止事件的复发事件数据。研究中采用含脆弱变量的可加可乘比率模型来描述复发事件过程，同时运用含脆弱变量的Cox型模型刻画终止事件过程。两类事件通过共同的脆弱变量建立联系，从而揭示终止事件与复发事件之间的关联性。进一步提出了联合模型中边缘参数和关联参数的估计方程，并探讨了估计量的渐进性质。

Abstract: Recurrent event data with two types of termination events were analyzed in depth on the basis of existing joint modeling approaches. In the study, an additive and multiplicative ratio model with fragile variables is used to describe the recurrent event process, while a Cox-type model with fragile variables is applied to portray the termination event process. The two types of events are linked through a common vulnerability variable, thus revealing the correlation between termination and recurrence events. The estimation equations for the edge and association parameters in the joint model are further proposed and the asymptotic nature of the estimates is explored.

1. 引言

复发事件数据在多个学科领域，如医学、经济学和生物学等，是一种常见的特殊数据类型。这种数据类型下，研究中的个体可能在特定时期内多次经历同一事件，具体表现包括：肿瘤再次出现、多次入院治疗、重复感染、连续服药以及周期性经济下滑等现象[1]。基于复发事件的强度或速率函数，已有多种方法用于分析复发事件数据[2]-[4]，Cook和Lawless (2007) [1]等人对复发事件方法进行了出色的综述。

在许多应用中，复发事件个体可能会出现死亡、治愈、接受其他治疗等终止事件，从而停止随访。通常终止事件与研究复发事件往往存在着紧密联系。分析带有复发事件的终止事件的方法可以分为强度模型、边际模型和偏边际模型[9]。强度模型通过引入一个共同的脆弱性变量，描述复发事件与终止事件之间的相互依赖。这些模型假设事件强度由可观测的协变量和未观测的脆弱变量共同决定，从而提供了理解事件关系的新视角[6]。作为强度模型的替代方案，一些作者提出了边际模型[7]-[9]，边际模型通过存活者与死亡个体比例的平均值来定义复发事件的发生概率，具有假设灵活、对泊松分布偏差稳健和解释直观等优点，适用于复杂关系情境。然而，其简化计算也导致模型参数的解释力减弱[18]。

偏边际模型主要研究存活个体的事件复发率，并可通过脆弱变量来描述复发和终止事件间的关联系[2] [5] [6] [11]-[13]。例如，Cook和Lawless (1997) [11]提出了个体在特定时间点的复发事件平均值和复发率概念。Ye等(2007) [12]构建了一个基于脆弱变量的联合半参数建模框架，在保持边际模型特性的同时，有效刻画了终止事件与复发事件的相依关系，即仅以协变量和脆弱变量为条件，而不考虑过程的先前历史。

Kalbfleisch等(2013) [10]提出了一种基于估计方程的新方法，用于估计联合脆弱性模型中的边际参数和关联参数。这种方法具有三个主要优点：无需泊松过程假设即可进行参数估计，提高了鲁棒性；直接估计共同脆弱性分布；能够估计事件间的依赖程度。然而实际应用中，加法模型作为乘法模型的重要替代选择，也具有其独特价值。Pan和Schaubel (2009) [5]采用加性模型来估计条件复发事件率，并使用比例风险模型来估计终止事件风险。这种方法可以评估治疗对幸存者生存率和复发事件率的影响，有助于解释研究结果间的差异。Zeng和Cai (2010) [2]以及Sun和Kang (2013) [13]分别研究了加法率模型和加法、乘法率模型，其乘法率模型体现在使用比例风险率模型对终止事件进行建模，由于基线率函数依赖于非参数的脆弱变量，因此复发和终止事件之间关联的脆弱变量被视为一干扰项。Chen等(2016) [14]考虑了部分Aalen加性模型，该模型对复发事件条件比率和终止事件风险率均采用了乘性脆弱变量。但他们并没有建立所提估计值的大样本性质，Sun等(2017) [15]提出了一种新模型，利用共享脆弱变量同时分析复发和终止事件，阐释二者的关联性，并给出了大样本估计特性。Lin和Ying (1995) [16]则指出，实际应用中，协变量可能同时呈现加性和乘性效应。基于这一观点，Sun和Kang (2013) [13]开发了一个复发事件的可加可乘模型，引入非参数形式的脆弱变量。然而，这种处理方式使得模型在解释复发事件和终止事件相关性方面存在冗余。孙琴等(2019) [17]构建了一个联合模型：复发事件采用含脆弱变量的可加和可乘结构，而终止事件则使用带脆弱变量的乘法模型。但因其复发和终止事件是唯一的，不能应用于更复杂的情况，在实际情况中，一个复发事件往往伴随着多个终止事件，就如上文所述，而研究者也十分关心它们之间的联系。

在孙琴等(2019) [17]这篇文章的展望提出的新的复发事件与终止事件的联合分析模型的基础上，继续对带有两种终止事件相关的脆弱变量的复发事件加法乘法比率模型进行分析，进一步提出了联合模型中边缘参数和关联参数的估计方程，最后再给出估计量的参数估计及其证明。

2. 模型假设和参数估计

定义${\tilde{N}}^R\left(t\right)$ 为$\left(0,t\right]$ 时间窗内观察到的复发事件累积计数过程，$J\left(t\right)={\left(Z{\left(t\right)}^{\text{T}},W{\left(t\right)}^{\text{T}}\right)}^{\text{T}}$ 为$p+q$ 维外生协变量[18]，其中$Z\left(t\right)$ 为$p\times 1$ 维协变量，$W\left(t\right)$ 为$q\times 1$ 维协变量。定义$D_i$ 和$C_i$ 分别代表第i个观察对象的终止事件时间和删失时间。其中终止事件会停止未来复发事件的发生，即当$t>D$ 时，$\text{d}{\tilde{N}}^R\left(t\right)=0$。$Y_i$ 表示终止事件类型，个体历经$k$ 种终止事件，则$Y_i=k$ ，$k$ 取值为1,2。记$T=\left(C\wedge D\right)$ ，$\text{Δ}\left(t\right)=I\left(T\ge t\right)$ ，其中$I(\cdot )$ 是示性函数，由于数据删失，${\tilde{N}}^R\left(t\right)$ 未能完全观测，令$N^R\left(t\right)={\tilde{N}}^R\left(t\wedge T\right)$ 表示复发事件在区间$\left(0,t\right]$ 上实际数量，同理，令$N^{Dk}\left(t\right)={\tilde{N}}^{Dk}\left(t\wedge T\right)$ 表示终止事件在区间$\left(0,t\right]$ 上实际数量，其中${\tilde{N}}^{Dk}\left(t\right)=I\left(D\le t\right),$ $k=1,2$ 。考虑包含$n$ 个独立同分布个体的样本，其观测数据形式为 $\left\{N_i^R\left(t\right),N_i^{D1}\left(t\right),N_i^{D2}\left(t\right),T_i,{\text{Δ}}_i\left(t\right),J_i\left(t\right),0\le t\le T_i,i=1,\cdots ,n\right\}$ 。

令$v_1$ 和$v_2$ 是两个非负不可观测的脆弱变量且分别独立于$J\left(t\right)$ 。根据Ye [12]和Kalbfleisch [10]。我们考虑一个给定$J\left(t\right)$ ，$D=s$ 及脆弱变量$v$ 的复发事件的偏边际比率模型

${\text{dΛ}}_R\left(t|v\right)=P\left\{\text{d}{\tilde{N}}^R\left(t\right)=1|J\left(t\right),D=s,v\right\},s\ge t$ (1)

特别是，${\text{dΛ}}_R\left(t|v\right)$ 可能会受协变量$J\left(t\right)$ 和脆弱变量$v$ 影响，不过与终止事件$D=s\ge t$ 无关。说明当协变量给定时，$v$ 对复发事件与终止事件进行联合建模。

则表达式(1)可以表示成

${\text{Λ}}_R\left(t|v\right)=P\left\{{\tilde{N}}^R\left(t\right)=1|J\left(t\right),D\ge t,v\right\}$ (2)

这表明给定$J\left(t\right)$ 和$v$ ，${\text{dΛ}}_R\left(t|v\right)$ 指定了存活到时间$t$ 的受试者中复发事件的边际比率。为了方便分析，假设存在两种终止事件，因此考虑如下的可加可乘比率模型

${\text{dΛ}}_R\left(t|v\right)=v_1{v}_2\left[\text{exp}\left\{{\gamma }^{\text{T}}Z\left(t\right)\right\}{\text{dΛ}}_0^R\left(t\right)+{\beta }^{\text{T}}W\left(t\right)\text{d}t\right]$ (3)

其中回归参数$\gamma$ 和$\beta$ 分别是$p\times 1$ 维向量和$q\times 1$ 维向量，${\text{dΛ}}_0^R\left(t\right)$ 是一个未指定的基线比率函数。

令${\text{dΛ}}_{Dk}\left(t|v\right)=P\left\{\text{d}{\tilde{N}}^{Dk}\left(t\right)=1|J\left(t\right),D\ge t,v\right\}$为给定协变量$J\left(t\right)$ 和脆弱变量$v$ 下的风险函数。我们指定以下对于终止事件的乘法风险率模型：

${\text{dΛ}}_{Dk}\left(t|v\right)=v_k\text{exp}\left\{{\alpha }^{\text{T}}J\left(t\right)\right\}{\text{dΛ}}_{0k}^D\left(t\right),k=1,2$ (4)

其中回归参数$\alpha$ 是$\left(p+q\right)\times 1$ 维度向量，${\text{dΛ}}_{0k}^D\left(t\right)$ 是一个未指定的基线风险函数，因而两种终止事件的风险函数可以用(4)式表达，为了便于表示，模型(3)和(4)的$J\left(t\right)$ 假定为同一组协变量，所提出的估计方法可进一步推广，适用于处理两个模型中存在差异的协变量集合。

此外，类似Ye (2007) [12]，假设脆弱变量$v_1$ 和$v_2$ 分别服从伽马分布$G\left(1,{\theta }_1\right)$ 和$G\left(1,{\theta }_2\right)$ ，为了可识别这个模型，固定$E\left(v_1\right)=1$ ，$E\left(v_2\right)=1$ ，且脆弱变量分布假设的稳健性已在孙琴[17]和Qu [15]中得到验证。我们在下文中假设给定$J(\cdot )$ ，删失时间$C$ 与$\left\{{\tilde{N}}^R(\cdot ),{\tilde{N}}^{D1}(\cdot ),{\tilde{N}}^{D2}(\cdot ),D,v\right\}$ 独立。

3. 估计方法

令$\eta ={\left({\gamma }^{\text{T}},{\beta }^{\text{T}}\right)}^{\text{T}}$ 。注意到，如果脆弱变量$v_1$ 和$v_2$ 已知，可以利用Liu [19]以及Gill [3]的方法分别应用到模型(3)和(4)，从而得到$\eta$ 、${\alpha }_1$ 和${\alpha }_2$ 的估计方程。但在实际情况中，$v_1$ 和$v_2$ 通常是未知的，这使得直接应用上述方法变得不可行。对此，通过考虑一个包含$\eta$ ，${\alpha }_1$ 和${\alpha }_2$ 的边际模型，其通过给定$D\ge t$ 和$X\left(t\right)$ 的条件下，取模型(3)和(4)的条件期望。假设脆弱变量$v_1$ 和$v_2$ 服从伽马分布情况下有：

$\begin{array}{c}{\text{dΛ}}_R\left(t\right)=P\left\{\text{d}{\tilde{N}}^R\left(t\right)=1|J\left(t\right),D\ge t\right\}\\ ={\psi }_1\left(t\right){\psi }_2\left(t\right)\left[\text{exp}\left\{{\gamma }^{\text{T}}Z\left(t\right){\text{dΛ}}_0^R\left(t\right)+{\beta }^{\text{T}}W\left(t\right)\text{d}t\right\}\right]\end{array}$ (5)

$P\left\{\text{d}{\tilde{N}}^{Dk}\left(t\right)=1|D\ge t,J\left(t\right)\right\}={\psi }_k\left(t\right)\text{exp}\left\{{\alpha }_k^{\text{T}}J\left(t\right)\right\}{\text{dΛ}}_{0k}^D\left(t\right)$ (6)

其中

${\psi }_k\left(t\right)=E\left[v_k|D\ge t,Z\left(t\right)\right]={\left[1+{\theta }_k\text{exp}\left\{{\alpha }_k^{\text{T}}J\left(t\right)\right\}{\text{dΛ}}_{0k}^D\right]}^{-1},k=1,2$

${\psi }_3\left(t\right)=E\left[v_1{v}_2|D\ge t,X\left(t\right)\right]={\psi }_1\left(t\right){\psi }_2\left(t\right)$

给其定义

$\text{d}{M}_i^R\left(t\right)={\psi }_{3i}\left(t\right)\text{d}{N}_i^R\left(t\right)-{\text{Δ}}_i\left(t\right)\left[\text{exp}\left\{{\gamma }^{\text{T}}{Z}_i\left(t\right)\right\}+{\beta }^{\text{T}}{W}_i\left(t\right)\text{d}t\right]$

$\text{d}{M}_i^{Dk}\left(t\right)={\psi }_{ki}\left(t\right)\text{d}{N}_i^{Dk}\left(t\right)-{\text{Δ}}_i\left(t\right)\exp \left\{{\alpha }_k^{\text{T}}{J}_i\left(t\right)\right\}{\text{dΛ}}_{0k}^D\left(t\right),k=1,2.$

其中

${\psi }_{ki}\left(t\right)={\left[1+{\theta }_k{{\displaystyle \int }}_0^t\exp \left\{{\alpha }_k^{\text{T}}{J}_i\left(t\right)\right\}{\text{dΛ}}_{0k}^D\left(t\right)\right]}^{-1},k=1,2.$

在假设下，模型(3)和(4)可以得出$M_i^R\left(t\right)$ 和$M_i^{Dk}\left(t\right)$ 是零均值随机过程。

当${\psi }_1\left(t\right),{\psi }_2\left(t\right)$ 已知时，相当于知道${\psi }_3\left(t\right)$ ，可以用以下方程估计${\text{Λ}}_0^R\left(t\right),{\text{Λ}}_{0k}^D\left(t\right),\eta ,{\alpha }_k$ ，其中$k=1,2$

${\displaystyle \sum }_{i=1}^n\left[{\psi }_{3i}\left(t\right)\text{d}{N}_i^R\left(t\right)-{\text{Δ}}_i\left(t\right)\left\{\text{exp}\left\{{\gamma }^{\text{T}}{Z}_i\left(t\right)\right\}{\text{dΛ}}_0^R\left(t\right)+{\beta }^{\text{T}}{W}_i\left(t\right)\text{d}t\right\}\right]=0,0\le t\le \tau ,$

${\displaystyle \sum }_{i=1}^n\left[{\psi }_{ki}\left(t\right)\text{d}{N}_i^{Dk}\left(t\right)-{\text{Δ}}_i\left(t\right)\text{exp}\left\{{\alpha }_k^{\text{T}}{J}_i\left(t\right)\right\}{\text{dΛ}}_{0k}^D\left(t\right)\right]=0,0\le t\le \tau ,k=1,2,$

${\displaystyle \sum }_{i=1}^n{{\displaystyle \int }}_0^{\tau }\left[D_i\left(t;\eta \right)-\overline{D}\left(t;\eta \right)\right]\left[{\psi }_{3i}\left(t\right)\text{d}{N}_i^R\left(t\right)-{\text{Δ}}_i\left(t\right){\beta }^{\text{T}}{W}_i\left(t\right)\text{d}t\right]=0,k=1,2$

${\displaystyle \sum }_{i=1}^n{{\displaystyle \int }}_0^{\tau }\left\{J_i\left(t\right)-{\overline{J}}_k\left(t\right)\right\}{\psi }_{ki}\left(t\right)\text{d}{N}_i^{Dk}\left(t\right)=0,k=1,2$

其中$\tau$ 是常数，使得$P\left(T_i\ge \tau \right)>0$ ，并且

${\overline{J}}_k\left(t\right)={\displaystyle \sum }_{i=1}^n{\text{Δ}}_i\left(t\right)\text{exp}\left\{{\alpha }_k^{\text{T}}{J}_i\left(t\right)\right\}{J}_i\left(t\right)/{\displaystyle \sum }_{i=1}^n{\text{Δ}}_i\left(t\right)\text{exp}\left\{{\alpha }_k^{\text{T}}{J}_i\left(t\right)\right\}$

根据Lin和Ying [16]，我们选取

$D_i\left(t;\eta \right)={\left(Z_i^{\text{T}}\left(t\right),W_i^{\text{T}}\left(t\right)/\text{exp}\left\{{\left({\alpha }_1+{\alpha }_2\right)}^{\text{T}}{J}_i\left(t\right)\right\}\right)}^{\text{T}}$

定义$\overline{D}\left(t;\eta \right)={\left(\overline{Z^{\text{T}}}\left(t;\eta \right),\overline{W^{\text{T}}}\left(t;\eta \right)\right)}^{\text{T}}$，其中

$\overline{Z}\left(t;\eta \right)={\displaystyle \sum }_{i=1}^n{\text{Δ}}_i\left(t\right)\text{exp}\left\{{\gamma }^{\text{T}}{Z}_i\left(t\right)\right\}{Z}_i\left(t\right)/{\displaystyle \sum }_{i=1}^n{\text{Δ}}_i\left(t\right)\text{exp}\left\{{\gamma }^{\text{T}}{Z}_i\left(t\right)\right\}$

$\overline{W}\left(t;\eta \right)={\displaystyle \sum }_{i=1}^n{\text{Δ}}_i\left(t\right)W_i\left(t\right)/{\displaystyle \sum }_{i=1}^n{\text{Δ}}_i\left(t\right)\text{exp}\left\{{\gamma }^{\text{T}}{Z}_i\left(t\right)\right\}$

而${\psi }_k\left(t\right)$ 是未知的，为了估计${\theta }_k$ ，令

$w_{1i}^k\left(t\right)=E\left[N_i^R\left(t\right)|\text{d}{N}_i^{Dk}\left(t\right)=1,J_i\left(t\right)\right]$

和

$w_{2i}\left(t\right)=E\left[N_i^R\left(t\right)|{\text{Δ}}_i\left(t\right)=1,J_i\left(t\right)\right]$

在模型假设下，可得

$w_{1i}^k\left(t\right)=\left({\theta }_k+1\right){\psi }_{3i}{\left(t\right)}^{-1}{\displaystyle {\int }_0^t\left[\text{exp}\left\{{\gamma }^{\text{T}}{Z}_i\left(t\right)\right\}{\text{dΛ}}_0^R\left(u\right)+{\beta }^{\text{T}}{W}_i\left(u\right)\text{d}u\right]}$

和

$w_{2i}\left(t\right)={\psi }_{3i}{\left(t\right)}^{-1}{\displaystyle {\int }_0^t\left[\text{exp}\left\{{\gamma }^{\text{T}}{Z}_i\left(t\right)\right\}{\text{dΛ}}_0^R\left(u\right)+{\beta }^{\text{T}}{W}_i\left(u\right)\text{d}u\right]}$

因此有

$w_{1i}^k\left(t\right)/w_{2i}\left(t\right)={\theta }_k+1$ (7)

当$k=1,2$ ，${\theta }_1$ 和${\theta }_2$ 描述了复发事件和各个终止事件的关联性。

如Kalbfleisch [10]所讨论的，为${\theta }_k$ 构造了以下估计方程：

${\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^{\tau }\left\{N_i^R\left(t\right)-\left({\theta }_k+1\right)G_k\left(t\right)w_{2i}\left(t\right)\right\}\text{d}{N}_i^{Dk}\left(t\right)}=0,k=1,2$

其中

$G_k\left(t\right)={\displaystyle \sum }_{i=1}^n{w}_{2i}{\left(t\right)}^{-1}{\text{Δ}}_{ki}^{*}\left(t\right)N_i^R\left(t\right)/{\displaystyle \sum }_{i=1}^n{\text{Δ}}_{ki}^{*}\left(t\right)$

${\text{Δ}}_{ki}^{\text{*}}\left(t\right)={\text{Δ}}_i\left(t\right)\left\{1-N_i^{Dk}\left(t\right)\right\}$ 代表个体$i$ 在$t$ 时刻仍存活且未经历终止事件，并且在之后的某个时间点发生了第$k$ 类终止事件。

令$\gamma ={\left({\eta }^{\text{T}},{\alpha }_1^{\text{T}},{\alpha }_2^{\text{T}},{\theta }_1,{\theta }_2,{\text{Λ}}_0^R,{\text{Λ}}_{01}^D,{\text{Λ}}_{02}^D\right)}^{\text{T}}$ ，通过解估计方程$U\left(\gamma \right)={\left(U_1^{\text{T}},U_{21}^{\text{T}},U_{22}^{\text{T}},U_{31},U_{32},U_4,U_{51},U_{52}\right)}^{\text{T}}=0$ 来获得$\gamma$ 的估计值$\widehat{\gamma }={\left({\widehat{\eta }}^{\text{T}},{\widehat{\alpha }}_1^{\text{T}},{\widehat{\alpha }}_2^{\text{T}},{\widehat{\theta }}_1,{\widehat{\theta }}_2,{}_{\widehat{\text{Λ}}}^0,{}_{\widehat{\text{Λ}}}^{01},{}_{\widehat{\text{Λ}}}^{02}\right)}^{\text{T}}$，其中

$U_1={\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^{\tau }\left[D_i\left(t;\eta \right)-\overline{D}\left(t;\eta \right)\right]}\left[{\psi }_{1i}\left(t\right){\psi }_{2i}\left(t\right)\text{d}{N}_i^R\left(t\right)-{\text{Δ}}_i\left(t\right){\beta }^{\text{T}}{W}_i\left(t\right)\text{d}t\right],$

$U_{2k}={\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^{\tau }\left\{J_i\left(t\right)-{\overline{J}}_k\left(t\right)\right\}{\psi }_{ki}\left(t\right)\text{d}{N}_i^{Dk}}\left(t\right),k=1,2$

$U_{3k}={\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^{\tau }\left\{N_i^R\left(t\right)-\left({\theta }_k+1\right)G_k\left(t\right)w_{2i}\left(t\right)\right\}\text{d}{N}_i^{Dk}}\left(t\right),k=1,2$

$U_4={\displaystyle \sum }_{i=1}^n\left[{\psi }_{3i}\left(t\right)\text{d}{N}_i^R\left(t\right)-{\text{Δ}}_i\left(t\right)\left\{\text{exp}\left\{{\gamma }^{\text{T}}{Z}_i\left(t\right)\right\}{\text{dΛ}}_0^R\left(t\right)+{\beta }^{\text{T}}{W}_i\left(t\right)\text{d}t\right\}\right],0\le t\le \tau$

$U_{5k}={\displaystyle \sum }_{i=1}^n\left[{\psi }_{ki}\left(t\right)\text{d}{N}_i^{Dk}\left(t\right)-{\text{Δ}}_i\left(t\right)\text{exp}\left\{{\alpha }^{\text{T}}{J}_i\left(t\right)\right\}{\text{dΛ}}_{0k}^D\left(t\right)\right],0\le t\le \tau ,k=1,2$

由于每个参数的估计式取决于其他参数的估计式，因此可以通过递归过程的方获得估计方程的解。因此，提出以下的迭代算法来求解$U\left(\gamma \right)=0$ 。

步骤0 选取初值${\theta }_1^{\left(0\right)}$ ，${\theta }_2^{\left(0\right)}$ ，${\alpha }_1^{\left(0\right)}$ ，${\alpha }_2^{\left(0\right)}$ 和${\text{Λ}}_{01}^{D\left(0\right)}\left(t\right)$ ，${\text{Λ}}_{02}^{D\left(0\right)}\left(t\right)$ 。

步骤1 令${\psi }_{ki}^{\left(0\right)}\left(t\right)={\psi }_{ki}\left(t;{\theta }_k^{\left(0\right)},{\alpha }_k^{\left(0\right)},{\Lambda }_{0k}^{D\left(0\right)}\right)$ ，由上述方程得

${\psi }_{3i}^{\left(0\right)}\left(t\right)={\psi }_{3i}\left(t;{\theta }_1^{\left(0\right)},{\theta }_2^{\left(0\right)},{\alpha }_1^{\left(0\right)},{\alpha }_2^{\left(0\right)},{\Lambda }_{01}^{D\left(0\right)}\left(t\right),{\Lambda }_{01}^{D\left(0\right)}\left(t\right)\right)$ ，把${\psi }_{ki}^{\left(0\right)}\left(t\right)$ 、${\psi }_{3i}^{\left(0\right)}\left(t\right)$ 代入$U_1=0$ ，$U_{2k}=0$ ，$U_4=0$ 和$U_{5k}=0$ ，其中$k=1,2$ ，则求解得到更新估计${\eta }^{\left(1\right)},{\alpha }_1^{\left(1\right)},{\alpha }_2^{\left(1\right)},{\text{Λ}}_0^{R\left(1\right)}\left(t\right)$ 和${\text{Λ}}_{01}^{D\left(1\right)}\left(t\right)$ ，${\text{Λ}}_{02}^{D\left(1\right)}\left(t\right)$ 。

步骤2 给定${\eta }^{\left(1\right)},{\text{Λ}}_0^{R\left(1\right)}\left(t\right)$ ，${\psi }_{3i}^{\left(0\right)}\left(t\right)={\psi }_{3i}\left(t;{\theta }_1^{\left(0\right)},{\theta }_2^{\left(0\right)},{\alpha }_1^{\left(1\right)},{\alpha }_2^{\left(1\right)},{\text{Λ}}_{01}^{D\left(1\right)}\left(t\right){\text{,Λ}}_{02}^{D\left(1\right)}\left(t\right)\right)$

解出$U_{31}=0$ ，$U_{32}=0$ ，得到估计${\theta }_1^{\left(1\right)}$ 和${\theta }_2^{\left(1\right)}$ 。

步骤3 重复步骤1和2，更新估计值，直到收敛。

注意，初始的估计值${\theta }_k^{\left(0\right)}$ ，${\alpha }_k^{\left(0\right)}$ 和${\text{Λ}}_{0k}^{D\left(0\right)}\left(t\right)$ ，$k=1,2$ 可以有多种的选取方法，为了简单起见，因而选择${\theta }_1^{\left(0\right)}=1$ ，${\theta }_2^{\left(0\right)}=1$ ，${\alpha }_1^{\left(0\right)}=0$ ，${\alpha }_2^{\left(0\right)}=0$ ，${\text{Λ}}_{01}^{D\left(0\right)}\left(t\right)$ 和${\text{Λ}}_{02}^{D\left(0\right)}\left(t\right)$ 为累计基线风险函数的Nelson-Aalen类型估计。对于收敛性，该算法大多数时候会收敛，但偶尔会发生非收敛，具体取决于初始值设置。

令${\xi }_1={\left({\eta }^{\text{T}},{\alpha }_1^{\text{T}},{\theta }_1\right)}^{\text{T}}$ ，${\xi }_2={\left({\eta }^{\text{T}},{\alpha }_2^{\text{T}},{\theta }_2\right)}^{\text{T}}$ ，$\xi ={\xi }_1\cup {\xi }_2={\left({\eta }^{\text{T}},{\alpha }_1^{\text{T}},{\alpha }_2^{\text{T}},{\theta }_1,{\theta }_2\right)}^{\text{T}}$ ，则${\widehat{\xi }}_1={\left({\widehat{\eta }}^{\text{T}},{\widehat{\alpha }}_1^{\text{T}},{\widehat{\theta }}_1\right)}^{\text{T}}$，${\widehat{\xi }}_2={\left({\widehat{\eta }}^{\text{T}},{\widehat{\alpha }}_2^{\text{T}},{\widehat{\theta }}_2\right)}^{\text{T}}$，$\widehat{\xi }={\widehat{\xi }}_1\cup {\widehat{\xi }}_2={\left({\widehat{\eta }}^{\text{T}},{\widehat{\alpha }}_1^{\text{T}},{\widehat{\alpha }}_2^{\text{T}},{\widehat{\theta }}_1,{\widehat{\theta }}_2\right)}^{\text{T}}$，${\xi }_1$ ，${\xi }_2$ 和$\xi$ 的真值分别为${\xi }_{10}={\left({\eta }_0^{\text{T}},{\alpha }_{10}^{\text{T}},{\theta }_{10}\right)}^{\text{T}}$ ，${\xi }_{20}={\left({\eta }_0^{\text{T}},{\alpha }_{20}^{\text{T}},{\theta }_{20}\right)}^{\text{T}}$ ${\xi }_0={\xi }_{10}\cup {\xi }_{20}={\left({\eta }_0^{\text{T}},{\alpha }_{10}^{\text{T}},{\alpha }_{20}^{\text{T}},{\theta }_{10},{\theta }_{20}\right)}^{\text{T}}$ 。下面描述了所提出的估计器的渐近属性。首先考虑$\widehat{\xi }$ 、${}_{\widehat{\text{Λ}}}^0\left(t\right)$、${}_{\widehat{\text{Λ}}}^{01}\left(t\right)$和${}_{\widehat{\text{Λ}}}^{02}\left(t\right)$的存在性、唯一性和强相合性。

4. 估计量的渐进性质

为了研究所提出的估计的渐近性质，需要以下正则性条件：

(C1) $\left\{N_i^R(\cdot ),N_i^{D1}(\cdot ),N_i^{D2}(\cdot ),T_i,J_i(\cdot )\right\},i=1,\cdots ,n$ 是独立同分布的。

(C2) $E\left\{N_i^R\left(\tau \right)\right\}$ 是有界的，$\text{Pr}\left(T\ge \tau \right)>0$ 。

(C3) $Z_i\left(t\right)$ 在$\left[0,\tau \right]$ 上几乎处处为有界变差随机过程。

(C4) 存在${\eta }_0$ 的一个紧集$ℬ$ ，满足$\eta \in ℬ$ ，$\text{Γ}\left(\eta \right)$ 是非奇异的，其中$\text{Γ}\left(\eta \right)$ 是$-\partial \tilde{U}\left(\eta \right)/\partial {\eta }^{\text{T}}$的极限值，其中$\tilde{U}\left(\eta \right)$ 的定义在(A.4)中给出。

定理1 $\widehat{\xi }$ 、${}_{\widehat{\text{Λ}}}^0\left(t\right)$、${}_{\widehat{\text{Λ}}}^{01}\left(t\right)$和${}_{\widehat{\text{Λ}}}^{02}\left(t\right)$几乎处处存在且唯一。

证明 令${\psi }_{ki}\left(t;{\text{Λ}}_{0k}^D\right)={\psi }_{ki}\left(t;{\theta }_k,{\alpha }_k,{\text{Λ}}_{0k}^D\right)$ ，其中$k=1,2$ ，

${\text{Λ}}_0^R\left(t;{\text{Λ}}_0^D\right)={\text{Λ}}_0^R\left(t;{\alpha }_1,{\alpha }_2,{\theta }_1,{\theta }_2,{\text{Λ}}_{01}^D,{\text{Λ}}_{02}^D\right)$ 。对任意向量${\gamma }^{*}$ 和变量$Z^{*}$ ，定义

$S_{k{Z}^{*}}\left(t;{\gamma }^{*}\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\text{Δ}}_i\left(t\right)\text{exp}\left\{{\gamma }^{*\text{T}}{Z}_i^{*}\left(t\right)\right\}{Z}_i^{*\otimes k}\left(t\right),k=0,1,$

其中对于任意向量$a$ ，有$a^{\otimes 0}=1,a^{\otimes 1}=a$ 。令

${\overline{N}}^{Dl}\left(t\right)=n^{-1}{\displaystyle \sum }_{i=1}^n{\psi }_{li}\left(t;{}_{\widetilde{\text{Λ}}}^{0l}\right)N_i^{Dl}\left(t\right),{\overline{N}}^R\left(t\right)=n^{-1}{\displaystyle \sum }_{i=1}^n{\psi }_{1i}\left(t;{}_{\widetilde{\text{Λ}}}^{01}\right){\psi }_{2i}\left(t;{}_{\widetilde{\text{Λ}}}^{02}\right)N_i^R\left(t\right)$

及$s_{k{Z}^{*}}\left(t;{\gamma }^{*}\right)$ 为$S_{k{Z}^{*}}\left(t;{\gamma }^{*}\right)$ 的极限$\left(k=0,1\right),\left(l=1,2\right)$ 。

对给定的$\xi$ ，令${}_{\widetilde{\text{Λ}}}^0\left(t;\xi \right)$及${}_{\widetilde{\text{Λ}}}^{01}\left(t;{\xi }_1\right)$和${}_{\widetilde{\text{Λ}}}^{02}\left(t;{\xi }_2\right)$分别为$U_4=0$ 及$U_{51}=0$ 和$U_{52}=0$ 的解。则

${}_{\widetilde{\text{Λ}}}^0\left(t;\xi \right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^t\frac{{\psi }_{1i}\left(t;{}_{\widetilde{\text{Λ}}}^{01}\right){\psi }_{2i}\left(t;{}_{\widetilde{\text{Λ}}}^{02}\right)\text{d}{N}_i^R\left(u\right)-{\text{Δ}}_i\left(u\right){\beta }^{\text{T}}{W}_i\left(u\right)\text{d}u}{S_{0Z}\left(u;\gamma \right)}}$

且${}_{\widetilde{\text{Λ}}}^{01}\left(t;{\xi }_1\right)$和${}_{\widetilde{\text{Λ}}}^{02}\left(t;{\xi }_2\right)$均满足方程

$\begin{array}{c}{}_{\widetilde{\text{Λ}}}^{0k}\left(t;{\xi }_k\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^t\frac{{\psi }_{ki}\left(t;{}_{\widetilde{\text{Λ}}}^{0k}\right)\text{d}{N}_i^{Dk}\left(u\right)}{S_{0X}\left(u;{\alpha }_k\right)}}\end{array}$ (8)

其中$k=1,2$ ，(8)是Volterra积分方程，且存在唯一解。

固定$\xi$ ，令

${\text{Λ}}_k^D\left(t;{\xi }_k\right)={\displaystyle {\int }_0^t\frac{E\left\{\text{d}{\overline{N}}^{Dk}\left(u\right)\right\}}{s_{0X}\left(u;{\alpha }_k\right)}}$ (9)

的解为${\text{Λ}}_k^D\left(t;{\xi }_k\right)$ ，其中$k=1,2$ 。

这是Volterra积分方程且有唯一解，满足${\text{Λ}}_1^D\left(t;{\xi }_{10}\right)\equiv {\text{Λ}}_{01}^D\left(t\right)$ ，${\text{Λ}}_2^D\left(t;{\xi }_{20}\right)\equiv {\text{Λ}}_{02}^D\left(t\right)$ ，令

$K_l\left(t;{\xi }_l\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^t\exp \left\{{\alpha }_l^{\text{T}}{J}_i\left(u\right)\right\}}\frac{{\theta }_l\text{d}{N}_i^{Dl}\left(u\right)}{S_{0X}\left(u;{\alpha }_l\right)}$

利用(8)和(9)可得

${}_{\widetilde{\text{Λ}}}^{0l}\left(t;{\xi }_l\right)-{\text{Λ}}_l^D\left(t;{\xi }_l\right)={\displaystyle {\int }_0^t\left[{}_{\widetilde{\text{Λ}}}^{0l}\left(t;{\xi }_l\right)-{\text{Λ}}_l^D\left(t;{\xi }_l\right)\right]\text{d}{K}_l\left(t;{\xi }_l\right)}+{ℰ}_{nl}\left(t;{\xi }_l\right)$

其中

${ℰ}_{nl}\left(t;{\xi }_l\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^t\frac{\text{d}{M}_i^{Dl}\left(u\right)}{S_{0X}\left(u;{\alpha }_{l0}\right)}}$

$l=1,2$ ，由强大数定律可知，在$t\in \left[0,\tau \right]$ 和$\xi \in ℬ$ 条件下，几乎必然一致地有${ℰ}_{nl}\left(t;{\xi }_l\right)\to 0$ 。

因为上式也是一个Volterra积分方程，求解可得

${}_{\widetilde{\text{Λ}}}^{0l}\left(t;{\xi }_l\right)-{\text{Λ}}_l^D\left(t;{\xi }_l\right)=\frac{1}{P_{nl}\left(t\right)}{\displaystyle {\int }_0^L{P}_{nl}\left(u-\right)\text{d}{ℰ}_{nl}\left(t;{\xi }_l\right)}$ (10)

其中

$P_{nl}\left(t\right)={\displaystyle \prod }_{s\le \text{t}}\left\{1-\text{d}{K}_l\left(t;{\xi }_l\right)\right\}$

$l=1,2$ ，是${ℰ}_{nl}\left(t;\xi \right)$ 在$\left[0,t\right]$ 上的乘积积分[1]及$P_n\left(u-\right)$ 表示$P_n\left(u\right)$ 的左连续版本。利用乘积积分的渐进性质[20]，一致强大数定律[21]及Lin和Ying (2011) [9]得出在$t\in \left[0,\tau \right]$ 和$\xi \in ℬ$ 时，${}_{\widetilde{\text{Λ}}}^{01}\left(t;{\xi }_1\right)$和${}_{\widetilde{\text{Λ}}}^{02}\left(t;{\xi }_2\right)$分别几乎处处收敛到${\text{Λ}}_1^D\left(t;{\xi }_1\right)$ 和${\text{Λ}}_2^D\left(t;{\xi }_2\right)$ ，设

${\text{Λ}}^R\left(t;\xi \right)={\displaystyle {\int }_0^t\frac{E\left\{\text{d}{\overline{N}}^R\left(u\right)\right\}-{\beta }^{\text{T}}{s}_{1W}^{*}\left(u\right)\text{d}u}{s_{0Z}\left(u;\gamma \right)}}$

令${\text{Λ}}^R\left(t;{\eta }_0\right)\equiv {\text{Λ}}_0^R\left(t\right)$ ，其中，$s_{1W}^{*}\left(u\right)$ 为$S_{kW}\left(t;0\right)$ 的极限值，在$t\in \left[0,\tau \right]$ 和$\xi \in ℬ$ 上，${}_{\widetilde{\text{Λ}}}^0\left(t;\xi \right)$几乎处处一致收敛到${\text{Λ}}^R\left(t;\xi \right)$ 。因此证明$\widehat{\xi }$ 、${}_{\widehat{\text{Λ}}}^{01}\left(t\right)$、${}_{\widehat{\text{Λ}}}^{02}\left(t\right)$及${}_{\widehat{\text{Λ}}}^0\left(t\right)$的存在性和唯一性，只需证明：

$\tilde{U}\left(\xi \right)=\left({\tilde{U}}_1{\left(\xi \right)}^{\text{T}},{\tilde{U}}_{21}{\left({\xi }_1\right)}^{\text{T}},{\tilde{U}}_{22}{\left({\xi }_2\right)}^{\text{T}},{\tilde{U}}_{31}\left(\xi \right),{\tilde{U}}_{32}\left(\xi \right)\right)=0$ (11)

有唯一解，则

${\tilde{U}}_1\left(\xi \right)={\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^{\tau }\left\{D_i\left(t;\eta \right)-\overline{D}\left(t;\eta \right)\right\}\left\{{\psi }_{1i}\left(t;{}_{\widetilde{\text{Λ}}}^{01}\right){\psi }_{2i}\left(t;{}_{\widetilde{\text{Λ}}}^{02}\right)\text{d}{N}_i^R\left(t\right)-{\text{Δ}}_i\left(t\right){\beta }^{\text{T}}{W}_i\left(t\right)\text{d}t\right\}}$

${\tilde{U}}_{2k}\left({\xi }_k\right)={\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^{\tau }\left\{J_i\left(t\right)-\overline{J}\left(t\right)\right\}{\psi }_{ki}\left(t;{}_{\widetilde{\text{Λ}}}^{0k}\right)\text{d}{N}_i^{Dk}\left(t\right)}$

${\tilde{U}}_{3k}\left(\xi \right)={\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^{\tau }\left\{N_i^R\left(t\right)-\left({\theta }_k+1\right){\tilde{G}}_k\left(t\right){\tilde{\omega }}_{2i}\left(t\right)\right\}\text{d}{N}_i^{Dk}\left(t\right)}$

其中$k=1,2$ ，则

${\tilde{\omega }}_{2i}\left(t\right)={\psi }_{1i}^{-1}\left(t;{}_{\widetilde{\text{Λ}}}^{01}\right){\psi }_{2i}^{-1}\left(t;{}_{\widetilde{\text{Λ}}}^{02}\right){\displaystyle {\int }_0^t\left[\text{exp}\left\{{\gamma }^{\text{T}}{Z}_i\left(u\right)\right\}{\text{d}\widetilde{\text{Λ}}}_0^R\left(u\right)+{\beta }^{\text{T}}{W}_i\left(u\right)\text{d}u\right]}$

${\tilde{G}}_k\left(t\right)$ 为$G_k\left(t\right)$ 的${\omega }_{2i}\left(t\right)$ 替换为${\tilde{\omega }}_{2i}\left(t\right)$ 可得。

设

$\widehat{\text{Γ}}\left(\xi \right)=-n^{-1}\partial \tilde{U}\left(\xi \right)/\partial {\xi }^{\text{T}}$

根据强大数定律以及${}_{\widetilde{\text{Λ}}}^{01}\left(t;{\widehat{\xi }}_1\right)$，${}_{\widetilde{\text{Λ}}}^{02}\left(t;{\widehat{\xi }}_2\right)$和${}_{\widetilde{\text{Λ}}}^0\left(t;\widehat{\xi }\right)$的一致收敛性，则对于$\xi \in ℬ$ 的$\widehat{\text{Γ}}\left(\xi \right)$一致收敛到非随机过程$\text{Γ}\left(\xi \right)$ 。可证：

$n^{-1}\tilde{U}\left({\xi }_0\right)\to 0$

所以，$\widehat{\text{Γ}}\left(\xi \right)$的一致相合性及(C4)得：当$n$ 足够大，$\xi \in ℬ$ 时，$\widehat{\text{Γ}}\left(\xi \right)$非奇异。根据逆函数定理[18]，则$\tilde{U}\left(\xi \right)=0$ 存在唯一解$\widehat{\xi }$ ，即存在唯一估计

$\widehat{\xi }$ ，${}_{\widehat{\text{Λ}}}^{01}\left(t\right)\equiv {}_{\widetilde{\text{Λ}}}^{01}\left(t;{\widehat{\xi }}_1\right)$，${}_{\widehat{\text{Λ}}}^{02}\left(t\right)\equiv {}_{\widetilde{\text{Λ}}}^{02}\left(t;{\widehat{\xi }}_2\right)$和${}_{\widehat{\text{Λ}}}^0\left(t\right)\equiv {}_{\widetilde{\text{Λ}}}^0\left(t;\widehat{\xi }\right)\left(0\le t\le \tau \right)$

定理1证毕。

定理2 $\widehat{\xi }$ 强相合于${\xi }_0$ ，且在$t\in \left[0,\tau \right]$ 上，${}_{\widehat{\text{Λ}}}^0\left(t\right)$几乎处处收敛到${\text{Λ}}_0^R\left(t\right)$ ，${}_{\widehat{\text{Λ}}}^{01}\left(t\right)$和${}_{\widehat{\text{Λ}}}^{02}\left(t\right)$几乎处处分别收敛到${\text{Λ}}_{01}^D\left(t\right)$ 和${\text{Λ}}_{02}^D\left(t\right)$ 。

证明 对(11)进行一阶泰勒展开，可得

$n^{-1}\tilde{U}\left(\widehat{\xi }\right)-n^{-1}\tilde{U}\left({\xi }_0\right)=-\widehat{\text{Γ}}\left({\xi }_0\right)\left(\widehat{\xi }-{\xi }_0\right)+o\left(\Vert \widehat{\xi }-{\xi }_0\Vert \right)$

因此几乎处处有

$\widehat{\text{Γ}}\left({\xi }_0\right)\left(\widehat{\xi }-{\xi }_0\right)+o\left(\Vert \widehat{\xi }-{\xi }_0\Vert \right)=o\left(1\right)$

由于$\xi \in ℬ$ ，$\widehat{\text{Γ}}\left(\xi \right)$非奇异，上述等式可知$\widehat{\xi }$ 是强相合的。由${}_{\widetilde{\text{Λ}}}^{01}\left(t;\widehat{\xi }\right)$，${}_{\widetilde{\text{Λ}}}^{02}\left(t;\widehat{\xi }\right)$和${}_{\widetilde{\text{Λ}}}^0\left(t;\xi \right)$的一致收敛性表明：对于$t\in \left[0,\tau \right]$ ，${}_{\widehat{\text{Λ}}}^{01}\left(t\right)$，${}_{\widehat{\text{Λ}}}^{02}\left(t\right)$和${}_{\widehat{\text{Λ}}}^0\left(t\right)$几乎处处一致地分别收敛到

${\text{Λ}}_1^D\left(t;{\xi }_{10}\right)\equiv {\text{Λ}}_{01}^D\left(t\right)$ ，${\text{Λ}}_2^D\left(t;{\xi }_{20}\right)\equiv {\text{Λ}}_{02}^D\left(t\right)$ 和${\text{Λ}}^R\left(t;{\xi }_0\right)\equiv {\text{Λ}}_0^R\left(t\right)$ 。

定理2证毕。

定理3 $\sqrt{n}\left(\widehat{\xi }-{\xi }_0\right)$ 按分布收敛到具有均值为0，协方差矩阵为${\text{Γ}}^{-1}\text{Σ}{\left({\text{Γ}}^{\text{T}}\right)}^{-1}$ 的正态随机变量。

证明 令${}_{\widetilde{\text{Λ}}}^{01}\left(t\right)={}_{\widehat{\text{Λ}}}^{01}\left(t;{\xi }_{10}\right)$，${}_{\widetilde{\text{Λ}}}^{02}\left(t\right)={}_{\widehat{\text{Λ}}}^{02}\left(t;{\xi }_{20}\right)$，由定理1得

${}_{\widetilde{\text{Λ}}}^{0l}\left(t\right)-{\text{Λ}}_{0l}^D\left(t\right)={\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^t\left[{}_{\widetilde{\text{Λ}}}^{0l}\left(t\right)-{\text{Λ}}_{0l}^D\left(t\right)\right]\text{d}{K}_l\left(u;{\xi }_{l0}\right)}+\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^t\frac{\text{d}{M}_i^{Dl}\left(u\right)}{S_{0X}\left(u;{\alpha }_{l0}\right)}}$

解得

${}_{\widetilde{\text{Λ}}}^{0l}\left(t\right)-{\text{Λ}}_{0l}^D\left(t\right)=\frac{1}{{\widehat{P}}_{nl}\left(t;{\xi }_{l0}\right)}{\displaystyle {\int }_0^t{\widehat{P}}_{nl}\left(u-;{\xi }_{l0}\right)}\frac{{{\displaystyle \sum }}_{i=1}^n\text{d}{M}_i^{Dl}\left(u\right)}{n{S}_{0X}\left(u;{\alpha }_{l0}\right)},$

其中$l=1,2$ ，

${\widehat{P}}_{nl}\left(t;{\xi }_{l0}\right)={\displaystyle \prod }_{s\le t}\left\{1-\text{d}{K}_l\left(s;{\xi }_{l0}\right)\right\}$

表示关于$K_l\left(s;{\xi }_0\right)$ 在$\left[0,t\right]$ 的乘积积分。根据${}_{\widetilde{\text{Λ}}}^{01}\left(t\right)$和${}_{\widetilde{\text{Λ}}}^{02}\left(t\right)$的一致收敛性，结合一致强大数定律[21]及Lin和Ying (2011) [9]研究，则$t\in \left[0,\tau \right]$ 时，有

${}_{\widetilde{\text{Λ}}}^{0l}\left(t\right)-{\text{Λ}}_{0l}^D\left(t\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\varphi }_{1i}^l\left(t\right)+o_p\left(n^{-1/2}\right)$ (12)

其中

${\varphi }_{1i}^l\left(t\right)=\frac{1}{P_l\left(t;{\xi }_{l0}\right)}{\displaystyle {\int }_0^t{P}_l\left(u-;{\xi }_{l0}\right)\frac{\text{d}{M}_i^{Dl}\left(u\right)}{s_{0X}\left(u;{\alpha }_{l0}\right)}}$

$P_l\left(t;{\xi }_0\right)$ 是${\widehat{P}}_{nl}\left(t;{\xi }_0\right)$ 极限，$l=1,2$ 。

令

${\text{K}}^{*}\left(t;{\xi }_0\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^t\exp }\left\{{\left({\alpha }_{10}+{\alpha }_{20}\right)}^{\text{T}}{J}_i\left(u\right)\right\}\frac{{\theta }_{10}{\theta }_{20}\text{d}{N}_i^R\left(u\right)}{S_{0Z}\left(u;{\gamma }_0\right)}.$

由(12)得

${}_{\widetilde{\text{Λ}}}^0\left(t\right)-{\text{Λ}}_0^R\left(t\right)={\displaystyle {\int }_0^t\left\{{}_{\widetilde{\text{Λ}}}^{01}\left(u\right){}_{\widetilde{\text{Λ}}}^{02}\left(u\right)-{\text{Λ}}_{01}^D\left(u\right){\text{Λ}}_{02}^D\left(u\right)\right\}{\text{dK}}^{*}\left(u;{\xi }_0\right)}+\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^t\frac{\text{d}{M}_i^R\left(u\right)}{S_{0Z}\left(u;{\gamma }_0\right)}}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\varphi }_{2i}\left(t\right)+o_p\left(n^{-\frac{1}{2}}\right)$ (13)

其中

${\varphi }_{2i}\left(t\right)={\displaystyle {\int }_0^t{\varphi }_{1i}^1\left(u\right){\varphi }_{1i}^2\left(u\right){\text{dK}}^{*}\left(u;{\xi }_0\right)}+{\displaystyle {\int }_0^t\frac{\text{d}{M}_i^R\left(u\right)}{s_{0Z}\left(u\right)}}$

在(11)中，对${\tilde{U}}_1\left({\xi }_0\right)$ 可分解为

$\begin{array}{c}{\tilde{U}}_1\left({\xi }_0\right)={\displaystyle \sum }_0^t{\displaystyle {\int }_0^{\tau }\left\{D_i\left(t;{\eta }_0\right)-\overline{D}\left(t;{\eta }_0\right)\right\}}\left\{{\psi }_{1i}\left(t;{}_{\widetilde{\text{Λ}}}^{01}\right){\psi }_{2i}\left(t;{}_{\widetilde{\text{Λ}}}^{02}\right)\text{d}{N}_i^R\left(t\right)-{\text{Δ}}_i\left(t\right){\beta }^{\text{T}}{W}_i\left(t\right)\text{d}t\right\}\\ ={\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^{\tau }{\text{Δ}}_i\left(t\right)\text{exp}\left\{{\left({\alpha }_{10}+{\alpha }_{20}\right)}^{\text{T}}{J}_i\left(t\right)\right\}}\left\{D_i\left(t;{\eta }_0\right)-\overline{D}\left(t;{\eta }_0\right)\right\}\left\{{}_{\widetilde{\text{Λ}}}^{01}\left(u\right){}_{\widetilde{\text{Λ}}}^{02}\left(u\right)-{\text{Λ}}_{01}^D\left(u\right){\text{Λ}}_{02}^D\left(u\right)\right\}\text{d}{N}_i^R\left(t\right)\\ +{\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^{\tau }\left\{D_i\left(t;{\eta }_0\right)-\overline{D}\left(t;{\eta }_0\right)\right\}\text{d}{M}_i^R\left(t\right)}\end{array}$

令

${\tilde{H}}_1\left(t\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^t{\theta }_{10}{\theta }_{20}{\text{Δ}}_i\left(u\right)\text{exp}\left\{{\left({\alpha }_{10}+{\alpha }_{20}\right)}^{\text{T}}{J}_i\left(t\right)\right\}\left\{D_i\left(u;{\eta }_0\right)-\overline{D}\left(u;{\eta }_0\right)\right\}\text{d}{N}_i^R\left(u\right)}$

且$\tilde{d}\left(t;{\eta }_0\right)$ 和$H_1\left(t\right)$ 分别是$\overline{D}\left(t;{\eta }_0\right)$ 和${\tilde{H}}_1\left(t\right)$ 的极限。可以证

${\tilde{U}}_1\left({\xi }_0\right)={\displaystyle \sum }_{i=1}^n{\vartheta }_{1i}+o_p\left(n^{\frac{1}{2}}\right)$ (14)

其中

${\vartheta }_{1i}={\displaystyle {\int }_0^{\tau }\left\{D_i\left(t;{\eta }_0\right)-\overline{d}\left(t;{\eta }_0\right)\right\}\text{d}{M}_i^R\left(t\right)}+{\displaystyle {\int }_0^t{\varphi }_{1i}^1\left(t\right){\varphi }_{1i}^2\left(t\right)\text{d}{H}_1\left(t\right)}$

类似地

${\tilde{U}}_2\left({\xi }_0\right)={\displaystyle \sum }_{i=1}^n{\vartheta }_{2i}^l+o_p\left(n^{\frac{1}{2}}\right)$ (15)

其中

${\vartheta }_{2i}^l={\displaystyle {\int }_0^{\tau }\left\{J_i\left(t\right)-\overline{j}\left(t\right)\right\}\text{d}{M}_i^{Dl}\left(t\right)}+{\displaystyle {\int }_0^{\tau }{\varphi }_{1i}^l\left(t\right)\text{d}{H}_{2l}\left(t\right)}$

$\overline{j}\left(t\right)$ 和$H_{2l}\left(t\right)$ 分别是$\overline{J}\left(t\right)$ 和${\tilde{H}}_{2l}\left(t\right)$ 的极限，$l=1,2$ ，且

${\tilde{H}}_{2l}\left(t\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^t{\theta }_{l0}{\text{Δ}}_i\left(u\right)\text{exp}\left\{{\alpha }_{l0}^{\text{T}}{J}_i\left(u\right)\right\}\left\{J_i\left(u\right)-\overline{J}\left(u\right)\right\}\text{d}{N}_i^{Dl}\left(u\right)}$

设

${\omega }_{2i}^{*}\left(t\right)={\psi }_{1i}^{-1}\left(t;{\text{Λ}}_{01}^D\right){\psi }_{2i}^{-1}\left(t;{\text{Λ}}_{02}^D\right){\displaystyle {\int }_0^t\left[\text{exp}\left\{{\gamma }^{\text{T}}{Z}_i\left(u\right)\right\}{\text{dΛ}}_0^R\left(u\right)+{\beta }^{\text{T}}{W}_i\left(u\right)\text{d}u\right]}$

$\tilde{Q}\left(t\right)=\frac{{{\displaystyle \sum }}_{i=1}^n{\text{Δ}}_i^{*}\left(t\right){\tilde{\omega }}_{2i}^{-1}\left(t\right)N_i^R\left(t\right)}{{{\displaystyle \sum }}_{i=1}^n{\text{Δ}}_i^{*}\left(t\right)}$

且$q\left(t\right)$ 为$\widehat{Q}\left(t\right)$ 的极限。经计算可得

$\begin{array}{c}{\tilde{U}}_{3l}\left({\xi }_0\right)={\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^{\tau }\left\{N_i^R\left(t\right)-\left({\theta }_{l0}+1\right){\omega }_{2i}^{*}\left(t\right)q\left(t\right)\right\}\text{d}{N}_i^{Dl}\left(t\right)}\\ -\left({\theta }_{l0}+1\right){\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^{\tau }\tilde{Q}\left(t\right)\left\{{\tilde{\omega }}_{2i}\left(t\right)-{\omega }_{2i}^{\text{*}}\left(t\right)\right\}\text{d}{N}_i^{Dl}\left(t\right)}\\ -\left({\theta }_{l0}+1\right){\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^{\tau }\left\{\tilde{Q}\left(t\right)-q\left(t\right)\right\}{\omega }_{2i}^{*}\left(t\right)\text{d}{N}_i^{Dl}\left(t\right)}\end{array}$ (16)

其中$l=1,2$ ，和(A.7)的证明一样：(A.9)右边的第二项化简得

$\left({\theta }_{l0}+1\right){\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^{\tau }\left[{\varphi }_{1i}^l\left(t\right)\text{d}{H}_{3l}\left(t\right)-{\varphi }_{2i}\left(t\right)\text{d}{H}_{4l}\left(t\right)\right]}+o_p\left(n^{\frac{1}{2}}\right)$ (17)

其中$H_{3l}\left(t\right)$ 和$H_{4l}\left(t\right)$ 分别是${\tilde{H}}_{3l}\left(t\right)$ 和${\tilde{H}}_{4l}\left(t\right)$ 的极限，且

${\tilde{H}}_{3l}\left(t\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^t{\theta }_{l0}}\left[\text{exp}\left\{{\alpha }_{l0}^{\text{T}}{J}_i\left(s\right)\right\}\tilde{Q}\left(s\right){\psi }_{1i}^{-1}\left(t;{}_{\widetilde{\text{Λ}}}^{01}\right){\psi }_{2i}^{-1}\left(t;{}_{\widetilde{\text{Λ}}}^{02}\right){\omega }_{2i}^{*}\left(s\right)\right]\text{d}{N}_i^{Dl}\left(s\right)$

和

${\tilde{H}}_{4l}\left(t\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^t\exp \left\{{\gamma }_0^{\text{T}}{Z}_i\left(s\right)\right\}\tilde{Q}\left(s\right){\psi }_{1i}^{-1}\left(t;{}_{\widetilde{\text{Λ}}}^{01}\right){\psi }_{2i}^{-1}\left(t;{}_{\widetilde{\text{Λ}}}^{02}\right)\text{d}{N}_i^{Dl}\left(s\right)}$

设${\text{dΦ}}_l\left(t\right)=E\left\{{\omega }_{2i}^{*}\left(t\right)\right\}\text{d}{N}_i^{Dl}\left(t\right)$ ，同样的，(A.9)右边的第三项等于

$\begin{array}{l}\left({\theta }_{l0}+1\right){{\displaystyle \sum }}_{i=1}^n{\displaystyle {\int }_0^{\tau }\left[{\varphi }_{1i}^l\left(t\right)H_{5l}\left(t\right)+{\varphi }_{2i}\left(t\right)H_6\left(t\right)\right]\text{d}\Phi \left(t\right)}\\ -\left({\theta }_{l0}+1\right){{\displaystyle \sum }}_{i=1}^n{\displaystyle {\int }_0^{\tau }\left[\frac{{\omega }_{2i}^{*-1}\left(t\right){\Delta }_i^{*}\left(t\right)N_i^R\left(t\right)}{E\left\{{\Delta }_i^{*}\left(t\right)\right\}}-\frac{q\left(t\right)}{E\left\{{\Delta }_i^{*}\left(t\right)\right\}}{\Delta }_i^{*}\left(t\right)\right]\text{d}\Phi \left(t\right)}+o_p\left(n^{1/2}\right)\end{array}$ (18)

其中$H_{5l}\left(t\right)$ 和$H_6\left(t\right)$ 分别是${\tilde{H}}_{5l}\left(t\right)$ 和${\tilde{H}}_6\left(t\right)$ 的极限，且

${\tilde{H}}_{5l}\left(t\right)=\frac{{\theta }_{l0}{{\displaystyle \sum }}_{i=1}^n\exp \left\{{\alpha }_{l0}^{\text{T}}{J}_i\left(t\right)\right\}{\psi }_{1i}^{-1}\left(t;{}_{\widetilde{\text{Λ}}}^{01}\right){\psi }_{2i}^{-1}\left(t;{}_{\widetilde{\text{Λ}}}^{02}\right){\tilde{\omega }}_{2i}^{-1}\left(t\right){\text{Δ}}_i^{*}\left(t\right)N_i^R\left(t\right)}{{{\displaystyle \sum }}_{i=1}^n{\text{Δ}}_i^{*}\left(t\right)}$

${\tilde{H}}_6\left(t\right)=\frac{{{\displaystyle \sum }}_{i=1}^n\exp \left\{{\gamma }_0^{\text{T}}{Z}_i\left(t\right)\right\}{\psi }_{1i}^{-1}\left(t;{}_{\widetilde{\text{Λ}}}^{01}\right){\psi }_{2i}^{-1}\left(t;{}_{\widetilde{\text{Λ}}}^{02}\right){\tilde{\omega }}_{2i}^{-1}\left(t\right){\omega }_{2i}^{*-1}\left(t\right){\text{Δ}}_i^{*}\left(t\right)N_i^R\left(t\right)}{{{\displaystyle \sum }}_{i=1}^n{\text{Δ}}_i^{*}\left(t\right)}$

由上面计算，则

${\tilde{U}}_{3l}\left({\xi }_0\right)={\displaystyle \sum }_{i=1}^n{\vartheta }_{3i}^l+o_p\left(n^{1/2}\right)$ (19)

$l=1,2$ ，其中

$\begin{array}{c}{\vartheta }_{3i}^l={\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^{\tau }\left\{N_i^R\left(t\right)-\left({\theta }_{l0}+1\right){\omega }_{2i}^{*}\left(t\right)q\left(t\right)\right\}\text{d}{N}_i^{Dl}\left(t\right)}\\ -\left({\theta }_{l0}+1\right){\displaystyle \sum }_{i=1}^n{\displaystyle {\int }_0^{\tau }\left[{\varphi }_{1i}^l\left(t\right)\text{d}{H}_{3l}\left(t\right)-{\varphi }_{2i}\left(t\right)\text{d}{H}_{4l}\left(t\right)\right]}\\ +\left({\theta }_{l0}+1\right){{\displaystyle \sum }}_{i=1}^n{\displaystyle {\int }_0^{\tau }\left[{\varphi }_{1i}^l\left(t\right)H_{5l}\left(t\right)+{\varphi }_{2i}\left(t\right)H_6\left(t\right)\right]\text{d}\Phi \left(t\right)}\\ -\left({\theta }_{l0}+1\right){{\displaystyle \sum }}_{i=1}^n{\displaystyle {\int }_0^{\tau }\left[\frac{{\omega }_{2i}^{*-1}\left(t\right){\Delta }_i^{*}\left(t\right)N_i^R\left(t\right)}{E\left\{{\Delta }_i^{*}\left(t\right)\right\}}-\frac{q\left(t\right)}{E\left\{{\Delta }_i^{*}\left(t\right)\right\}}{\Delta }_i^{*}\left(t\right)\right]\text{d}\Phi }\end{array}$ (20)

令

${\vartheta }_i={\left({\vartheta }_{1i}^{\text{T}},{\vartheta }_{2i}^1{}^{\text{T}},{\vartheta }_{2i}^2{}^{\text{T}},{\vartheta }_{3i}^1,{\vartheta }_{3i}^2\right)}^{\text{T}}$

及$\Gamma =\Gamma \left({\xi }_0\right)$ ，根据泰勒一阶展开可得

$n^{1/2}\left(\widehat{\xi }-{\xi }_0\right)=n^{-1/2}{\Gamma }^{-1}{{\displaystyle \sum }}_{i=1}^n{\vartheta }_i+o_p\left(1\right)$

通过多元中心极限理论，$n^{1/2}\left(\widehat{\xi }-{\xi }_0\right)$ 是零均值和协方差矩阵为${\Gamma }^{-1}\Sigma {\left({\Gamma }^{\text{T}}\right)}^{-1}$ 的渐进正态分布，其中$\Sigma =E\left\{{\vartheta }_i{\vartheta }_i^{\text{T}}\right\}$ 。定理3证毕。

定理4 当$t\in \left[0,\tau \right]$ ，$\left\{\sqrt{n}\left({}_{\widehat{\text{Λ}}}^0\left(t\right)-{\text{Λ}}_0^R\left(t\right)\right)\right\}$，$\left\{\sqrt{n}\left({}_{\widehat{\text{Λ}}}^{01}\left(t\right)-{\text{Λ}}_{01}^D\left(t\right)\right)\right\}$和$\left\{\sqrt{n}\left({}_{\widehat{\text{Λ}}}^{02}\left(t\right)-{\text{Λ}}_{02}^D\left(t\right)\right)\right\}$均弱收敛到零均值随机过程。

证明 由于证明$n^{1/2}\left\{{\widehat{\Lambda }}_{01}^D\left(t\right)-{\Lambda }_{01}^D\left(t\right)\right\}$ 和$n^{1/2}\left\{{\widehat{\Lambda }}_{02}^D\left(t\right)-{\Lambda }_{02}^D\left(t\right)\right\}$ 的弱收敛性类似，为方便起见，只需证明其一即可，不妨证明$n^{1/2}\left\{{\widehat{\Lambda }}_{01}^D\left(t\right)-{\Lambda }_{01}^D\left(t\right)\right\}$ 和$n^{1/2}\left\{{}_{\widehat{\text{Λ}}}^0\left(t\right)-{\text{Λ}}_0^R\left(t\right)\right\}$的弱收敛性，注意到

${}_{\widehat{\text{Λ}}}^{01}\left(t\right)-{\text{Λ}}_{01}^D\left(t\right)=\left\{{}_{\widehat{\text{Λ}}}^{01}\left(t;{\widehat{\xi }}_1\right)-{}_{\widehat{\text{Λ}}}^{01}\left(t;{\xi }_{10}\right)\right\}+\left\{{}_{\widehat{\text{Λ}}}^0\left(t;{\xi }_{10}\right)-{\text{Λ}}_0^D\left(t\right)\right\}$ (21)

令

${}_{\widehat{\text{Υ}}}\left(t;{\xi }_1\right)=\partial {}_{\widehat{\text{Λ}}}^{01}\left(t,{\xi }_1\right)/\partial {\xi }_1$

根据一致强大数定律，当$t\in \left[0,\tau \right]$ 和$\xi \in ℬ$ 时，可证明${}_{\widehat{\text{Υ}}}\left(t;{\xi }_1\right)$几乎处处收敛到非随机函数${\text{Υ}}_{11}\left(t;{\xi }_1\right)$ ，由泰勒一阶展开和(12)，(20)，(21)得

$n^{1/2}\left\{{}_{\widehat{\text{Λ}}}^{01}\left(t\right)-{\text{Λ}}_{01}^D\left(t\right)\right\}=n^{-1/2}{\displaystyle \sum }_{i=1}^n{\text{Ψ}}_{1i}\left(t\right)+o_p\left(1\right)$

其中

${\text{Ψ}}_{1i}^1\left(t\right)=-{\text{Υ}}_{11}\left(t;{\xi }_{10}\right){\text{Γ}}^{-1}{\vartheta }_i+{\varphi }_{1i}^1\left(t\right)$

类似的，由(13)和(20)得

$n^{1/2}\left\{{\widehat{\Lambda }}_0^R\left(t\right)-{\Lambda }_0^R\left(t\right)\right\}=n^{-1/2}{{\displaystyle \sum }}_{i=1}^n{\Psi }_{2i}\left(t\right)+o_p\left(1\right)$

其中

${\text{Ψ}}_{2i}\left(t\right)=-{\text{Υ}}_2\left(t;{\xi }_0\right){\text{Γ}}^{-1}{\vartheta }_i+{\varphi }_{2i}\left(t\right)$

及${\text{Υ}}_2\left(t;{\xi }_0\right)$ 为$\partial {}_{\widehat{\text{Λ}}}^0\left(t;\xi \right)/\partial \xi$极限，给定

${\text{Ψ}}_i\left(t\right)={\left({\text{Ψ}}_{1i}^1\left(t\right),{\text{Ψ}}_{1i}^2\left(t\right){\text{,Ψ}}_{2i}^2\left(t\right)\right)}^{\text{T}}$

在时刻$t$ 的条件下，${\text{Ψ}}_i\left(t\right)\left(i=1,\cdots ,n\right)$ 独立同分布且具有零均值。由多元中心极限定理，$n^{1/2}\left\{{\widehat{\Lambda }}_{01}^D\left(t\right)-{\Lambda }_{01}^D\left(t\right)\right\}$ ，$n^{1/2}\left\{{\widehat{\Lambda }}_{02}^D\left(t\right)-{\Lambda }_{02}^D\left(t\right)\right\}$ 和$n^{1/2}\left\{{}_{\widehat{\text{Λ}}}^0\left(t\right)-{\text{Λ}}_0^R\left(t\right)\right\}$的有限维分别弱收敛于零均值过程。此外，${\text{Ψ}}_i\left(t\right)$ 可表示为$t$ 的单调函数的累加或乘积形式，故此过程是紧密的[22]。则$n^{1/2}\left\{{\widehat{\Lambda }}_{01}^D\left(t\right)-{\Lambda }_{01}^D\left(t\right)\right\}$ ，$n^{1/2}\left\{{\widehat{\Lambda }}_{02}^D\left(t\right)-{\Lambda }_{02}^D\left(t\right)\right\}$ 和$n^{1/2}\left\{{}_{\widehat{\text{Λ}}}^0\left(t\right)-{\text{Λ}}_0^R\left(t\right)\right\}$是紧的，并在$\left(s,t\right)$ 处的协方差函数为$E\left\{{\text{Ψ}}_i\left(s\right){\text{Ψ}}_i{\left(t\right)}^{\text{T}}\right\}$ 。定理4证毕。

5. 结束语

针对已存在双重终止事件的复发事件联合模型，继续对其进行研究，通过建立估计方程，推导出方程参数的相合性与渐近正态性。

因为模型(1)和(2)允许脆弱变量$v$ 与复发事件与终止事件之间相关性是正向，也可以是负向。对于这种情况，可以采用以下模型：

${\text{dΛ}}_R\left(t|v\right)=v^{-1}\exp \left\{{\gamma }^{\text{T}}Z\left(t\right){\text{dΛ}}_0^R\left(t\right)+{\beta }^{\text{T}}W\left(t\right)\text{d}t\right\},$

其中，模型(2)不变，$v$ 为单位均值，方差小于1的伽马分布。可得

${\psi }^{*}\left(t\right)=E\left[v^{-1}|X\left(t\right),D\ge t\right]=\psi \left(t\right)/\left(1-\theta \right)$

和

${\omega }_{1i}\left(t\right)/{\omega }_{2i}\left(t\right)=1-\theta .$

这表明复发事件的发生率与终止事件呈反向关联，因为在时间$t$ 有一个终止事件的个体预计，会比在时间$t$ 之后有一个终止事件的个体有更少的复发事件。此外，通过在$U_1$ ，$U_3$ 和$U_4$ 中分别用${\psi }^{*}\left(t\right)$ 和$1-\theta$ 代替$\psi \left(t\right)$ 和$1+\theta$ ，可以构造与前几节相同的估计方程。

然而，在分析终止事件和复发事件数据的联合模型时，协变量的选择标准与确定方法应当遵循以下原则：首先需要通过Wald统计检验和效应量评估，即标准化系数绝对值大于1.96来识别协变量的显著性；对于显著协变量，应进一步通过似然比检验比较其在加法模型和乘法模型中的拟合优度，以确定其更适合纳入可加部分$W\left(t\right)$ 还是可乘部分$Z\left(t\right)$ 。当协变量维度较高时，则采用LASSO惩罚回归方法进行变量筛选；对于存在高度删失的数据，应当使用逆概率加权法或多重插补技术来校正选择偏差。此外，需要通过Bootstrap重抽样来验证最终模型的稳定性，要求核心参数的变化率不超过15%。这一系统化的协变量选择流程既保证了统计严谨性，又能适应不同类型的数据特征。

采用广义估计方程法进行参数估计，在半参数模型框架下可能存在效率提升空间。未来研究可着眼于开发更简洁高效的推断方法。另外，为评估模型(1)和(2)对数据的拟合程度，可考虑基于残差$\text{d}{\widehat{M}}_i^R\left(t\right)$和$\text{d}{\widehat{M}}_i^D\left(t\right)$的图形和数值分析方法[13]。其中，$\text{d}{\widehat{M}}_i^R\left(t\right)$和$\text{d}{\widehat{M}}_i^D\left(t\right)$通过将$\text{d}{M}_i^R\left(t\right)$ 和$\text{d}{M}_i^D\left(t\right)$ 中的参数替换为相应估计值获得。这一问题的理论深化和数值模拟值得进一步探索。

基金项目

国家自然科学基金项目(11561010)，广西高校中青年教师科研基础能力提升项目(2022KY0707)资助。

NOTES

^*通讯作者。