竞争风险数据下关于时变暴露的逆概率处理加权研究

Inverse Probability of Treatment Weighting for Time-Varying Exposure in Competing Risks Data

摘要: 在生存分析中，当研究个体的暴露因素随时间动态变化，且结局事件存在竞争风险时，传统Cox比例风险模型的估计会产生偏倚。本文旨在将逆概率处理加权方法与Fine-Gray竞争风险模型相结合，提出一种双重加权模型，以精准估计竞争风险数据下时变暴露对感兴趣事件的因果效应。首先，构建逆概率权重以调整时间依赖性混杂。然后，将逆概率权重与Fine-Gray模型的删失权重结合，形成双重权重并嵌入加权得分方程。最后，利用恶性黑色素瘤患者数据，分别建立时变分段暴露变量和静态二元暴露变量模型进行对比分析。实例分析表明，在双重加权竞争风险模型中，相较于静态暴露变量，时变暴露变量能更精确地捕捉暴露的动态效应，估计结果更为稳健。

Abstract: In survival analysis, when the exposure factors of study subjects change dynamically over time and the outcome events involve competing risks, traditional Cox proportional hazards models may yield biased estimates. This paper aims to integrate the inverse probability of treatment weighting method with the Fine-Gray competing risks model to propose a doubly-weighted model, in order to accurately estimate the causal effect of time-varying exposure on the event of interest in the presence of competing risks. First, inverse probability weights are constructed to adjust for time-dependent confounding. Then, these weights are combined with the censoring weights from the Fine-Gray competing risks model to form double weights, which are embedded into a weighted score equation. Finally, using data from malignant melanoma patients, models with time-varying segmented exposure variables and static binary exposure variables are established for comparative analysis. The empirical analysis demonstrates that, in the doubly-weighted competing risks model, time-varying exposure variables more accurately capture the dynamic effects of exposure and yield more robust estimates compared to static exposure variables.

1. 引言

生存分析作为统计学的重要分支，在医学、生物学、工程学等众多领域有着广泛应用。随访研究中个体的治疗或暴露水平常随时间动态变化，且感兴趣事件往往与一个或多个竞争事件并存，使传统以固定协变量为前提的Cox模型在推断上存在一定的局限性。1999年，Fine和Gray [1]提出了子分布比例风险模型(简称FG模型)，为竞争风险分析奠定了基础。该模型针对竞争风险数据中感兴趣事件的累积发生率函数进行建模，为分析时变暴露在竞争风险环境下的作用提供了有力工具。在FG模型的基础上，众多学者围绕时变暴露的处理展开了深入研究。2009年，Cortese和Andersen [2]对竞争风险模型中的时变暴露进行了系统且深入的讨论，严格区分了外部时变暴露与内部时变暴露两类关键类型，详细剖析了不同类型时变暴露在传统竞争风险模型中的适配局限——外部时变暴露因与个体生存状态独立可直接纳入模型，内部时变暴露则因“存活依赖观测”特性导致传统模型无法实现效应估计与累积发生率预测的统一，并提出引入中间瞬态状态、扩展吸收状态、里程碑分析三类多状态解决方案以优化时变暴露的处理流程。他们的研究强调，明确时变暴露类型并选择适配的建模策略，是提升竞争风险模型对时变暴露效应推断准确性与临床解释性的核心前提。

2000年，Robins等[3]首次将IPW法应用于时变暴露的生存分析中，他们提出了边际结构模型，通过构建逆概率权重来消除时间依存性混杂因素的影响。在边际结构模型框架下，IPW法能够有效地调整时变暴露与生存结局之间的关联，使得模型能够更准确地估计时变暴露的因果效应。2008年，Cole和Hernan [4]详细阐述了时变暴露IPW权重的构造方法，区分了不稳定权重和稳定权重。不稳定权重定义为治疗分配概率的倒数，而稳定权重在分子中加入了边际治疗概率，从而减少了权重的方差。2017年，Grafféo等[5]针对传统R包ipw在处理时变暴露时仅支持单次暴露切换的局限，提出了扩展ipw包功能的方法，实现了对多阶段、间歇性时变暴露的逆概率加权估计优化。本文为准确评估时变暴露对感兴趣事件发生的影响，将逆概率处理加权嵌入FG模型，构造了双重权重体系，在子分布风险层面实现稳健的因果效应估计，并对比时变分段暴露变量与传统二元暴露变量二者结果进行分析。

本文探讨下竞争风险数据下关于时变暴露的逆概率处理加权进行研究，结构安排如下：第二部分介绍逆概率加权方法，基于逆概率权重和删失权重的竞争风险模型的核心理论与方法；第三部分为模拟分析；第四部分实例分析，对恶性黑色素瘤患者数据，应用上述方法进行实证研究，并对比分析不同暴露定义下的结果，以验证所提方法的有效性与稳健性，第五部分总结全文。

2. 方法介绍

2.1. 逆概率加权方法介绍

当数据存在时间依赖性暴露(如个体开始或停止治疗的时间不同、治疗剂量随时间变化)时，传统的固定暴露分析方法(如基线倾向评分)无法处理随时间动态变化的混杂因素(如治疗期间病情恶化影响后续治疗选择)，而逆概率加权可通过构建“时间依赖的权重”解决这一问题，核心是使用IPW可以通过创建一个伪种群来调整混杂因素的暴露，在该伪种群中，暴露与测得的混杂因素无关，从而无偏估计暴露的因果效应。

对于给定的患者$i$ ，直到有序失效时间$k$ 的权重定义为：

$ip{w}_{ik}={\displaystyle \prod _{t=1}^k\frac{1}{P\left(A_t=a_{t,i}|{\overline{A}}_{t-1}={\overline{a}}_{t-1,i},{\overline{C}}_t={\overline{c}}_{t,i},V=v_i\right)}}$ (2.1)

其中，$a_t$ 表示患者在时间$t$ 接受的治疗，${\overline{a}}_{t-1}$ 表示截至时间$t-1$ 的观察治疗暴露史，${\overline{c}}_t$ 表示观察到的时变混杂因素史，$v_i$ 表示观察到的时间固定协变量。当上述权重公式分母的条件部分定义的某些层中发现很少的治疗暴露值时，可能会出现不稳定的权重。因此，Robins等人[6]提出了稳定权重：

$s{w}_{ik}={\displaystyle \prod _{t=1}^k\frac{P\left(A_t=a_{t,i}|{\overline{A}}_{t-1}={\overline{a}}_{t-1,i},V=v_i\right)}{P\left(A_t=a_{t,i}|{\overline{A}}_{t-1}={\overline{a}}_{t-1,i},{\overline{C}}_t={\overline{c}}_{t,i},V=v_i\right)}}$ (2.2)

通过稳定化得到的权重分布更窄，这会产生更窄的置信区间，从而提高统计效率。

2.2. 基于逆概率权重的竞争风险模型建模

本节先介绍Fine与Gray关于累积发生率函数风险建模的研究。该方法与Cox比例风险模型[7]存在诸多共同点。虽然理论推导采用渐进式展开(从无删失数据情形入手)，但本文将直接呈现适用于不完整数据的一般情况——这正是我们通常面对的数据类型：既包含发生目标事件的个体，也存在发生竞争风险事件的个体，同时还存在删失观测案例。该模型基于

$\gamma \left(t,x\right)={\gamma }_0\left(t\right){\text{e}}^{\beta x}$ (2.3)

其中，$\gamma$ 表示子分布风险，${\gamma }_0$ 为子分布基准风险，$x$ 是协变量向量，$\beta$ 为系数向量。

该偏似然函数的形式与Cox比例风险模型中的形式相似。若仅考虑单一协变量$x$ ，其偏似然函数的数学表达式为：

$L\left(\beta \right)={\displaystyle \prod _{j=1}^r\frac{\exp \left(\beta x_j\right)}{{\displaystyle \sum _{i\in R_j}{w}_{ji}\exp \left(\beta x_i\right)}}}$ (2.4)

该乘积运算取自所有观察到目标事件发生的时间点($t_1<t_2<\cdots <t_r$ )。Fine-Gray模型与Cox比例风险模型相比，Fine-Gray模型主要有两个差异，风险集$R_j$ 的定义方式不同，并且引入了删失权重$w_{ji}$ 。风险集由在时间$t$ 之前未经历任何事件的个体，以及在时间$t$ 之前经历了竞争风险事件的个体组成。

$R_j=\{i:T_i\ge t\text{ or }\left(T_i\le t且经历过一个竞争事件的个体\right)\}$

因此，经历了其他类型事件的个体将始终保留在风险集中。在存在竞争风险模型的前提下，一种类型$i$ 事件在时间$t$ 之前不发生的概率称为子生存函数，定义为$S_i\left(t\right)=P\left(T>t,C=i\right)$ 。

删失权重$w_{ji}$ 的定义如下：

$w_{ji}=\frac{\widehat{S}\left(t_j\right)}{\widehat{S}\left(\min \left(t_j,t_i\right)\right)}$ (2.5)

其中，$\widehat{S}$ 是删失分布生存函数的Kaplan-Meier估计量。删失分布由$\left(T_i,C_i\right)$ 定义，其中$T_i$ 是到首次观察到事件发生的时间，$C_i$ 是一个指示变量：如果未观察到任何事件，则取值为1；如果观察到了任何类型的事件，则取值为0。在每个观测到感兴趣事件的时间点发生的时间点(下标为$j$ )，处于风险中的个体集合包括那些在时间$t_j$ 之前未发生任何类型事件的个体，以及在$t_j$ 之前发生过竞争风险事件的个体。前者的权重为1，后者的权重则小于或等于1。

通过这种加权机制，那些经历了竞争风险事件的个体并不会完全参与部分似然函数的计算：竞争风险事件发生时间点距离当前时间点$t_j$ 越远，相应个体所分配的权重就越小。

为方便计算，对给出的部分似然函数取自然对数后得到：

$\log L\left(\beta \right)={\displaystyle \sum _{j=1}^r\{\beta x_j-\log \left({\displaystyle \sum _{i\in R_j}{w}_{ji}\exp \left(x_i\beta \right)}\right)\}}$ (2.6)

对上述部分似然函数的自然对数应用求导运算，即可得到得分统计量，其形式为：

$U\left(\beta \right)={\displaystyle \sum _{j=1}^r\{x_j-\frac{{\displaystyle \sum _{i\in R_j}{w}_{ji}{x}_i\exp \left(x_i\beta \right)}}{{\displaystyle \sum _{i\in R_j}{w}_{ji}\exp \left(x_i\beta \right)}}\}}$ (2.7)

在Fine-Gray模型的基础上，对于存在时间依赖性暴露的数据可能会出现时变混杂，为消除这一影响引入了逆概率权重$ipw$ ，这是为了调整风险集中个体的相对重要性和调整每个事件的贡献强度，它确保在伪总体中，事件的发生与暴露的因果路径保持一致且在每个时间点，暴露分配与协变量历史无关。包含双重权重的加权部分似然函数表达式为：

$L^{IPW}\left(\beta \right)={{\displaystyle \prod _{j=1}^r(\frac{\exp \left(\beta x_j\right)}{{\displaystyle \sum _{i\in R_j}{w}_{ji}ip{w}_i\exp \left(\beta x_i\right)}})}}^{ip{w}_j}$ (2.8)

取自然对数后得到：

$\log L^{IPW}\left(\beta \right)={\displaystyle \sum _{j=1}^{r}ip{w}_j\{\beta x_j-\log \left({\displaystyle \sum _{i\in R_j}{w}_{ji}ip{w}_i\exp \left(x_i\beta \right)}\right)\}}$ (2.9)

对每个事件点$j$ 分别求导并将所有事件点导数相加，得到加权得分统计量表达式：

$U^{IPW}\left(\beta \right)={\displaystyle \sum _{j=1}^{r}ip{w}_j\{x_j-\frac{{\displaystyle \sum _{i\in R_j}{w}_{ji}ip{w}_i{x}_i\exp \left(x_i\beta \right)}}{{\displaystyle \sum _{i\in R_j}{w}_{ji}ip{w}_i\exp \left(x_i\beta \right)}}\}}$ (2.10)

3. 数据模拟

为评估在竞争风险生存数据中，基于逆概率加权的时变暴露模型与静态暴露模型对暴露效应估计的准确性，本文设计并实施了一项模拟研究。

3.1. 模拟参数设定

设第$i$ 个个体的基线协变量为$X=\left(X_{i1},X_{i2}\right)$ ，其中$X_{i1}~U\left(0,3\right)$ ；$X_{i2}~N\left(5,1\right)$ 。为构造暴露变量，令起始治疗时间$t_{0i}~U\left(20,500\right)$ 。并定义时变暴露变量$A_{it}$ 为二值变量：

$A_{it}=\{\begin{array}{c}0t<t_{0i}\\ 1t\ge t_{0i}\end{array}$

同时，定义静态暴露变量$T$ ，若个体在随访期间任何时间接受治疗则$T=1$ ，否则$T=0$ 。

生存时间$t$ 的生成基于竞争风险模型。设$\gamma \left(t\right)$ 为总风险函数，${\lambda }_1\left(t\right)$ 和${\lambda }_2\left(t\right)$ 分别为主要事件与竞争事件的特定风险函数。基线风险函数设定如下：

$\begin{array}{l}{\gamma }_0\left(t\right)=0.001\exp \left(-0.001t/\log \left(1.3\right)\right)\\ {\lambda }_{1,0}\left(t\right)=0.001\\ {\lambda }_{2,0}\left(t\right)={\gamma }_0\left(t\right)-{\lambda }_{1,0}\left(t\right)-\frac{{{\gamma }^{\prime }}_0\left(t\right)}{{\gamma }_0\left(t\right)}+\frac{{{\lambda }^{\prime }}_{1,0}\left(t\right)}{{\lambda }_{1,0}\left(t\right)}\end{array}$

其中，${{\gamma }^{\prime }}_0\left(t\right)$ 和${{\lambda }^{\prime }}_{1,0}\left(t\right)$ 分别为对应函数的一阶导数，通过数值微分计算。

协变量对风险的影响通过比例风险模型引入。设真实参数向量为$\xi =\left({\xi }_1,{\xi }_2,{\xi }_3\right)=\left(\log \left(1.15\right),\log \left(0.9\right),\log \left(2\right)\right)$ ，分别对应$X_1$ 、$X_2$ 和时变治疗$A\left(t\right)$ 的对数风险比。在时间$t$ ，给定协变量向量$X=\left(X_1,X_2,A\left(t\right)\right)$ 时，风险函数定义为：

$\begin{array}{l}\gamma \left(t,\xi ,X\right)={\gamma }_0\left(t\right)\exp \left({\displaystyle \sum _{k=1}^3{\xi }_k{X}_k}\right)\\ {\lambda }_1\left(t,\xi ,X\right)={\lambda }_{1,0}\left(t\right)\exp \left({\displaystyle \sum _{k=1}^3{\xi }_k{X}_k\exp \left(-0.0005t\right)}\right)\\ {\lambda }_2\left(t,\xi ,X\right)=\gamma \left(t,\xi ,X\right)-{\lambda }_1\left(t,\xi ,X\right)-\frac{\partial \gamma \left(t,\xi ,X\right)/\partial t}{\gamma \left(t,\xi ,X\right)}+\frac{\partial {\lambda }_1\left(t,\xi ,X\right)/\partial t}{{\lambda }_1\left(t,\xi ,X\right)}\end{array}$

在时间$t$ 发生事件的概率为${\lambda }_1\left(t,\xi ,X\right)+{\lambda }_2\left(t,\xi ,X\right)$ ，事件类型(1或2)的发生与${\lambda }_1$ 和${\lambda }_2$ 成比例。个体生存时间通过离散时间近似方法模拟，最大迭代次数设为1000。

为适应时变协变量分析，将生存数据转换为长格式的长型数据集。若个体在治疗开始前发生事件($t\le t_0$ )，则生成一条记录；若在开始治疗后发生事件($t>t_0$ )，则生成两条记录，分别对应治疗前时期$[0,t_0)$ 和治疗后时期$\left[t_0,t\right]$ ，其中治疗后时期的$A\left(t\right)=1$ 。

在模拟中构建时变暴露模型($Z_i\left(t\right)=\left(X_1,X_2,A\left(t\right)\right)$ )和静态暴露模型($Z_i\left(t\right)=\left(X_1,X_2,T\right)$ )，两种模型均基于加权部分似然得分方程估计参数$\beta =\left({\beta }_1,{\beta }_2,{\beta }_3\right)$ ，分别对应$X_1$ 、$X_2$ 和暴露变量($A\left(t\right)$ 或$T$ )的效应。加权指数项为$W_i\left(t_j\right)=w_i\left(t_j\right)s{w}_i\exp \left({\beta }^{\prime }{Z}_i\left(t_j\right)\right)$ ，其中$w_i\left(t_j\right)$ 为对于右删失数据，采用逆概率删失加权进行调整，删失分布通过Kaplan-Meier方法估计得到的权重；$s{w}_i$ 为逆概率权重。以此构造得分函数$U\left(\beta \right)$ ，参数估计通过求解$U\left(\beta \right)=0$ 获得。

模拟研究设置样本量$n$ 为1000，模拟重复次数为100次。

3.2. 模拟结果

两种模型的估计结果如下表1和表2。

Table 1. Estimated value under the time-varying exposure model

表1. 时变暴露模型下的估计值

	均值	真实值	偏差	标准误	p值
$X_1$	0.1688	0.1398	0.0270	0.0083	0.0012
$X_2$	−0.1174	−0.1054	−0.0120	0.0071	0.0910
$A\left(t\right)$	1.0787	0.6931	0.3855	0.0144	0.0000

Table 2. Estimated value under the static exposure model

表2. 静态暴露模型下的估计值

	均值	真实值	偏差	标准误	p值
$X_1$	0.1574	0.1398	0.0176	0.0278	0.5270
$X_2$	−0.1596	−0.1054	−0.0542	0.0232	0.0198
$T$	−0.1437	0.6931	−0.8369	0.3403	0.0139

从模拟中可以看出，时变暴露模型估计更稳定，能更好捕捉暴露效应。

4. 实例分析

4.1. 数据描述

本小节实例分析了恶性黑色素瘤患者生存时间影响因素的案例。数据来源于R软件包survival中的Melanoma数据集，该数据集收集了丹麦哥本哈根大学医院皮肤科1962年至1977年间确诊的恶性黑色素瘤患者的临床随访资料。研究旨在评估不同临床病理特征对患者生存预后的影响，为黑色素瘤的临床诊疗和预后评估提供循证医学依据。数据集包含205个观测值以及相关信息，包括患者的生存时间、生存状态、肿瘤浸润深度、细胞类型、表皮溃疡情况等指标，具体研究指标变量见表3。

Table 3. Variable definitions and meanings

表3. 变量定义与意义

变量	定义	意义
time	从诊断到事件发生或末次随访的时间	响应变量(以天为单位)
status	事件状态	结局变量(0 = 删失，1 = 黑色素瘤死亡， 2 = 其他原因死亡)
event	事件类型详细描述	分类结局变量(具体死亡原因或删失)
invasion_num	肿瘤浸润深度	协变量(描述肿瘤侵袭的分级变量)
ici	细胞类型指数	协变量(表示干扰素治疗的不同水平)
epicel_num	表皮细胞化生情况	协变量(1 = 存在，0 = 不存在)
ulcer	表皮溃疡情况	协变量(1 = 存在，0 = 不存在)
thick	肿瘤厚度	协变量(描述肿瘤厚度的连续变量)
sex	性别	协变量(Female = 女性，Male = 男性)
age	年龄	协变量(患者年龄)
logthick	肿瘤厚度的对数变换	协变量(肿瘤厚度的对数转换值)
traettime	治疗时间	动态变量(表示每个个体的治疗时间)
treat	治疗状态	静态暴露变量(0 = 未治疗，1 = 治疗)
tvc	治疗状态	时间依赖性暴露变量

4.2. 时变暴露变量的建立

数据预处理：将分类变量(侵袭程度、表皮情况、溃疡、性别)转换为数值型变量，便于建模分析。基于肿瘤厚度、浸润深度、溃疡、年龄4个临床指标加权评分(0~10分)，按评分将患者分为极高危(≥8分)、高危(6~7分)、中危(4~5分)、低危(2~3分)、极低危(0~1分) 5个风险层。根据“风险越高、干预越早”原则，为不同风险层设定个体化治疗启动时间(treattime)：极高/高/中/低危组分别基于生存时间的 2%/8%/20%/60%分位数，并限定时间区间；极低危组设为无穷大(视为无需治疗)，缺失/无穷值则设为超生存时间值。

构建时间依赖性数据结构：对于治疗时间晚于或等于生存时间的患者，生成单段数据：[0, ytime]区间，tvc = 0；对于治疗时间早于生存时间的患者，生成两段数据：治疗前段：[0, treattime]，tvc = 0；治疗后段：[treattime, ytime]，携带最终事件状态，tvc = 1；定义二分类暴露变量(treat)和时间依赖性暴露变量(tvc)进行对比研究：

静态暴露(treat)：二值变量，反映“是否曾接受治疗”

${\text{treat}}_i=I\left({\text{treattime}}_i<{\text{tyears}}_i\right)$

时变暴露(tvc)：分段常数函数

$Z_i\left(t\right)=\{\begin{array}{l}0t<{\text{treattime}}_i\\ 1t\ge {\text{treattime}}_i\end{array}$

4.3. 结果分析

为了评估逆概率加权以及暴露类型的对比研究效果，对此数据集，建立如下回归模型：

$\lambda \left(t|x\right)={\lambda }_0\left(t\right)\exp \left({\beta }_1\text{tvc}+{\beta }_2\text{epicel\_num}+{\beta }_3\text{invasion\_num}\right)$

在含有时间依赖性暴露tvc的情况下，使用公式(2.1)进行逆概率加权，代入公式(2.10)，方程$U^{IPW}\left(\beta \right)=0$ 得到的估计结果如下表4。

Table 4. Model estimation parameter results

表4. 模型估计参数结果

变量	估计值	标准误	p值
tvc	2.9999	0.4300	0.0000
epicel_num	−0.6967	0.2816	0.0133
invasion_num	−0.1634	0.2084	0.4329

基于表4分析表明，在双重加权竞争风险模型中，存在表皮细胞化生情况(epicel_num)和接受治疗时间依赖性暴露(tvc)在0.05的显著性水平上是显著的，接受治疗的生存时间是不接受治疗情况下的2.9999倍。

在含有静态暴露treat的情况下，同样使用公式(2.1)进行逆概率加权，代入公式(2.10)，方程$U^{IPW}\left(\beta \right)=0$ 得到的估计结果如下表5。

Table 5. Model estimation parameter results

表5. 模型估计参数结果

变量	估计值	标准误	p值
treat	1.8153	0.3850	0.0000
epicel_num	−0.6654	0.2819	0.0182
invasion_num	0.1159	0.2066	0.5746

基于表5分析表明，在双重加权竞争风险模型中，存在表皮细胞化生情况(epicel_num )和接受治疗静态暴露(treat)在0.05的显著性水平上是显著的，接受治疗的生存时间是不接受治疗情况下的1.8153倍。对比表2，可以得出在存在时变暴露的生存数据中，在其他变量保持不变的情况下，时变分段暴露变量相较于传统二元治疗变量能更精准捕捉暴露的时间动态效应。

5. 结语

本文研究了在存在时间依赖性暴露与竞争风险的生存数据中，结合逆概率加权与Fine-Gray模型的分析方法。本文提出的双重加权模型通过逆概率权重平衡混杂因素，并通过删失权重准确处理竞争风险事件，使得对目标事件风险的估计结果更为精准稳健。实例分析表明，在存在时变暴露的竞争风险数据中，具有逆概率加权的模型，时变分段暴露比静态暴露更精准地估计暴露对目标事件的因果效应，有效控制偏倚，保证结果的稳健性。

基金项目

本项目由2022年度辽宁省研究生教育教学改革研究项目(项目编号：LNYJG2022395)资助。

NOTES

^*通讯作者。