查看原文
其他

Stata:回归标准化估计因果效应-standsurv

连享会 连享会 2023-10-24

👇 连享会 · 推文导航 | www.lianxh.cn

理论 + 实证:从「读懂模型」到「折腾模型」
🎦 理论模型构建专题
📅 2022 年 4 月 23-24 日 (周六-周日)
🔑 郭凯明副教授 (中山大学)
🍓 课程主页https://gitee.com/lianxh/emodel

作者:陈文静 (中山大学)
邮箱:chenwj266@163.com

编者按:本文主要摘译自下文,特此致谢!
Sourse:Syriopoulou E, Mozumder S I, Rutherford M J, et al. Estimating causal effects in the presence of competing events using regression standardisation with the Stata command standsurv[J]. arXiv preprint arXiv:2109.03628, 2021. -PDF-


目录

  • 1. 背景介绍

  • 2. 命令介绍

  • 3. 实战演示

    • 3.1 估计总效应

    • 3.2 直接影响

    • 3.3 可分离的影响

  • 4. 参考资料

  • 5. 相关推文



温馨提示: 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码,直达原文:

原文是一篇研究因果效应的优秀文章,作者在 Github 上提供了所有的图表和 Stata 程序,使得读者能够完美重现文中的结果。本推文侧重点在于结合原文了解后估计命令 standsurv,以及对采用回归标准化估计竞争风险模型 (Comepeting Risk Model) 的因果效应有更深的认识。

1. 背景介绍

当我们对一个时间—事件 (time to event,TTE) 的结果感兴趣时,可能会出现阻止感兴趣事件发生的竞争性事件,即研究中的结局事件可能有多个。例如在治疗对前列腺癌死亡 (感兴趣事件) 的因果效应研究中,死于其他原因的患者便不可能再死于前列腺癌。因此,其他原因导致的死亡被认为是竞争事件。

在竞争事件环境下定义因果效应可能很复杂,需要特别考虑如何处理竞争事件。当竞争事件被容纳时,因果效应定义为治疗对感兴趣事件的总影响。当竞争事件被剔除时,因果效应定义为治疗对感兴趣事件的直接影响 (Stensrud 等,2020)。在治疗通过不同因果途径对感兴趣事件和竞争事件产生影响的情况下,可分离效应 (Separable Effects) 被定义 (Young 等,2020)。

具体地,我们来看一下相关概念:

  • 总效应 (Total Effects):在不需要消除竞争事件的环境下,治疗和感兴趣事件之间所有因果途径产生的总效应;
  • 直接效应 (Direct Effects):在消除竞争事件干预下,治疗对感兴趣事件的影响;
  • 可分离的直接效应 (Separable Direct Effects):治疗对感兴趣事件的影响,不以其对竞争事件的影响为中介;
  • 可分离的间接效应 (Separable Indirect Effects):治疗对感兴趣事件的影响只通过其对竞争事件的影响来实现。

在竞争风险模型下,研究中缺少使用正式的因果框架来描述估计结果,使得对其解释很麻烦。standsurv 命令应用了回归标准化,将估计值计算为所有个人特定预测的平均值。为了估计回归标准化的平均因果效应,需要事先拟合一个生存模型,对研究中每个治疗组下的每个个体进行预测,计算出针对针对个人估计值的平均值,并计算出治疗组之间的相关对比 (如:治疗组与治疗组之间的差异)。

2. 命令介绍

standsurv 是一个后估计命令,用于标准化估计参数生存模型的各种 (边际) 估计值,包括标准化生存函数、标准化生存函数的各种函数、以及这些函数的对比。这些函数包括失效函数 (failure function) 和危险函数 (hazard function)。

在标准化时,standsurv 命令可以保持特定的协变量不变,并在不同群体之间进行对比,例如差异和比率。同时, standsurv 还允许在竞争风险设置中指定多个生存模型,以计算标准化的特定原因累积发病率函数 (cause-specific cumulative incidence functions)。

// 命令安装
ssc install standsurv, replace
// 语法介绍
standsurv [, options]

主要 options 如下:

  • crmodels(modellist):设定竞争性风险模型;
  • at1()...atn():为每个原因设置具体的协变量值;
  • atvars():每个 atn() 选项对应的新变量;
  • cif:计算竞争风险模型的累积发病率函数;
  • contrast():对比 atn() 选项定义的协变量模式;
  • contrastvars():每个模式对比对应的新变量;
  • lincom(numlist)options 之间的线性组合;
  • timevar(varname):用于预测的时间变量 (默认 _t);
  • rmst:计算限制平均生存时间;
  • rmft:计算限制平均失败时间。

3. 实战演示

我们使用一个关于前列腺的公开试验数据 prostate.dta,来演示如何使用 standsurv 获得标准化的因果效应。相关数据和代码可通过以下方式获得:

cnssc install lxhget, replace
lxhget standsurv_example.zip, replace
unzipfile standsurv_example.zip, replace

示例数据包括 502 个被随机分配到雌激素治疗的个体,根据雌激素剂量被分为四个治疗组。为了简单起见,我们分析仅限于高剂量雌激素治疗组 (diethylstilbestrol,DES) 和安慰剂 (placebo)。我们感兴趣的是治疗对前列腺死亡的因果效应,但其它原因导致的死亡会阻止我们观察因前列腺死亡的发生。以下的变量将在我们的分析中使用:

  • rx:治疗组 (1:DES,0:安慰剂)
  • hyBinary:血红蛋白水平 (1:< 12 (g/100ml),0:≥ 12)
  • ageCat:年龄 (0:0-59 岁,1:60-74 岁,2:75-100 岁)
  • hx:心血管疾病史 (0:无,1:有)
  • normalAct:日常活动功能 (0:正常活动,1:其它)
  • dtime:随访月份
  • eventType:死因 (0:活着,1:因前列腺癌死亡,2:因其他原因死亡)

下图显示了各治疗组全因死亡的 Kaplan-Meier 失败曲线:

与安慰剂组相比,随机化后的前几个月,DES 组因任何原因死亡的概率更高。然而,在随机化后约 20 个月,曲线首次出现交叉,并在 60 个月内保持接近,这表明治疗对全因死亡概率的影响几乎可忽略不计。如果想了解更多关于生存曲线 Kaplan-Meier curve,可以参考连享会推文「Stata:生存分析一文读懂」。

3.1 估计总效应

在存在竞争事件下,特定病因的危险函数被定义为:

其中, 表示死因,即如果是由于前列腺癌导致的死亡 (感兴趣事件),则 ;如果是由于其他原因导致的死亡,则

3.1.1 特定病因的累积发病率函数 (CIF)

治疗对感兴趣事件的总影响,我们用特定病因的累计发生函数 (Cause-specific Cumulative Incidence Functions,CIFs) 来定义。 表示前列腺癌导致死亡的反事实累计发生率。在有其他原因导致的死亡作为竞争事件的情况下,前列腺死亡的边际反事实累积发生率定义为:

其中,

是全因生存率, 是前列腺癌的危险率。平均因果差异被定义为:

其含义是治疗 (通过所有因果途径) 对前列腺癌死亡的影响,包括那些可能由竞争事件介导的影响。同样,在有前列腺癌导致的死亡作为竞争事件的情况下,其他原因导致的死亡的边际反事实累积发生率为:

其中, 是其他原因的危险率 (Young 等,2020)。

通过 standsurv 命令的回归标准化,我们可以得到 DES 和安慰剂下的标准化特定病因的累计发生率函数 (CIFs)。为此,我们使用选项 cif 要求估计特定病因的 CIF (默认为总生存率),选项 crmodel(cancer other) 给出特定病因模型估计的名称。每个模型估计值都事先用 estimate store 存储在内存中。

/* obtain predictions of standardised CIFs as well as their difference
using standsurv with the option cif and crmodels() */
standsurv, crmodels(prostate other) cif at1(rx 0) at2(rx 1) timevar(timevar) ///
contrast(difference) ci atvars(CIF0 CIF1) contrastvars(CIF_diff)

上述命令语句的含义如下:

atn():根据其指定的固定协变量值创建一个标准化的 CIF。通过 at1(),我们强制将所有受试者的协变量 rx 设置为 0 (安慰剂);通过 at2(),我们强制将所有受试者的协变量 rx 设置为 1 (DES)。关键的一点是,在 DES 和安慰剂下,其余的混杂因素的分布被强制为相同,任何没有在 atn() 指定的协变量都保持其观察值。

例如,如果我们不使用观察到的年龄分布,而是想获得每个患者都属于最老的年龄组 (ageCat3) 的预测结果,那么可以通过在 at1() 选项中设置年龄组 ageCat2 取值为 0 ,年龄组 ageCat3 取值为 1 来获得。

contrast():要求对两个 CIFs (DES 下和安慰剂下) 进行比较,其差异参数要求取 CIFs 之间的差异。默认情况下,at1() 是参考,即对比是 ,这也可以通过 atref() 选项改变。atvar():给出了每个 atn() 创建的新变量的名称。contrastvar():给出了 contrast() 创建的新变量的名称。如在上面的例子中,创建了以下变量:

  • CIF0_prostate:在安慰剂下,由于前列腺癌导致死亡的 CIF;
  • CIF1_prostate:在 DES 治疗组下,由于前列腺癌导致死亡的 CIF;
  • CIF_diff_prostate:在 DES 治疗组和安慰剂下,由于前列腺癌导致死亡的差异;
  • CIF0_other:在安慰剂下,由于其他原因导致死亡的 CIF;
  • CIF1_other:在 DES 治疗组下,由于其他原因导致死亡的 CIF;
  • CIF_diff_other:在 DES 治疗组和安慰剂下,由于其他原因导致死亡的差异。

下图显示了 DES 和安慰剂下,前列腺癌死亡的标准化 CIF 和其他死因的 CIF,以及它们在随机化后的差异。随机化后 60 个月 (5 年),DES 下前列腺癌死亡的 CIF 为 21.3% (95% CI: 15.3%-29.5%);而安慰剂下则更高,为 27.7% (95% CI: 21.2%-36.2%)。

对于其他原因导致的死亡的 CIFs,情况正好相反,而且要高得多。在 DES 下,为 53.5% (95% CI: 45.9%-62.2%),而在安慰剂下,为 43.1% (95% CI: 35.9%-51.7%)。上述估计值考虑了竞争性事件也存在的事实。

3.1.2 因某一死因导致的预期寿命损失

治疗总效果的另一个估计值是时间 之前的预期寿命损失 (Mozumder 等,2021),当治疗设定为 时,在时间 前的边际反事实的预期寿命损失,也称为限制性平均失败时间 (Restricted Mean Failure Time,RMFT),被定义为:

其中, 表示死亡原因,RMFT 可进一步划分为在时间 之前由于每个特定原因 所造成的寿命损失:

当设定 时,在时间 由于原因 造成的预期寿命损失的边际反事实差异可以定义为:

预期寿命损失将研究人口与一个不死队列进行比较,即在时间点 上,随访期结束,所有人仍然活着。尽管预期寿命损失是探索不同原因影响的有用措施,但这种比较会使措施的解释具有挑战性,,因为比较涉及一个假设的结构。对于示例中预期寿命损失的估计,我们需要选择一个时间点 。估计结果将因 的选择而不同。这里我们选择 60 个月:

// generate time variable at 60 months (5 years) as t*
gen t_rmft60 = 60 in 1

拟合特定原因的模型后,通过使用选项 rmftcrmodels()cif 得到每个原因导致的标准化预期寿命损失:

/* obtain the standardised expected loss in life due to a cause of death
before 60 months under placebo and under DES as well as their difference */
standsurv, crmodels(prostate other) cif rmft ///
at1(rx 0) at2(rx 1) timevar(t_rmft60) contrast(difference) ci ///
atvars(RMFT0 RMFT1) contrastvars(RMFT_diff)

列出在 60 个月内因为前列腺癌而损失的寿命估计值:

// list the estimates of the life years lost due to prostate cancer before 60 months
// under placebo
list t_rmft60 RMFT0_prostate* in 1, noobs abb(22)
+---------------------------------------------------------------------+
| t_rmft60 RMFT0_prostate RMFT0_prostate_lci RMFT0_prostate_uci |
|---------------------------------------------------------------------|
| 60 10.113965 7.5143802 13.612872 |
+---------------------------------------------------------------------+

//under DES
list t_rmft60 RMFT1_prostate* in 1, noobs abb(22)
+---------------------------------------------------------------------+
| t_rmft60 RMFT1_prostate RMFT1_prostate_lci RMFT1_prostate_uci |
|---------------------------------------------------------------------|
| 60 6.9161976 4.7224429 10.129035 |
+---------------------------------------------------------------------+

// their difference
list t_rmft60 RMFT_diff_prostate* in 1, noobs abb(22)
+---------------------------------------------------------------------------------+
| t_rmft60 RMFT_diff_prostate RMFT_diff_prostate_lci RMFT_diff_prostate_uci |
|---------------------------------------------------------------------------------|
| 60 -3.1977678 -7.2035117 .80797615 |
+---------------------------------------------------------------------------------+

同理,列出在 60 个月内因其他原因而损失的寿命估计值:

// To list the estimates of the life years lost due to other causes before 60 months
// under placebo
list t_rmft60 RMFT0_other* in 1, noobs abb(22)
+------------------------------------------------------------+
| t_rmft60 RMFT0_other RMFT0_other_lci RMFT0_other_uci |
|------------------------------------------------------------|
| 60 15.63767 12.644827 19.338875 |
+------------------------------------------------------------+

// under DES
list t_rmft60 RMFT1_other* in 1, noobs abb(22)
+------------------------------------------------------------+
| t_rmft60 RMFT1_other RMFT1_other_lci RMFT1_other_uci |
|------------------------------------------------------------|
| 60 19.812743 16.498442 23.79284 |
+------------------------------------------------------------+

// their difference
list t_rmft60 RMFT_diff_other* in 1, noobs abb(22)
+------------------------------------------------------------------------+
| t_rmft60 RMFT_diff_other RMFT_diff_other_lci RMFT_diff_other_uci |
|------------------------------------------------------------------------|
| 60 4.1750728 -.54819852 8.8983441 |
+------------------------------------------------------------------------+

当选择 60 个月时,其它原因导致死因的损失月数比前列腺癌多。在 DES 下,由于其他原因导致死亡的损失月数也更多:安慰剂下为 19.8 个月,DES 下为 15.6 个月,差异为 4.2 个月。然而,因前列腺癌而损失的月数在安慰剂下更多:在 DES 下损失 6.9 个月,在安慰剂下损失 10.1 个月,它们之间的差异为 -3.2 个月。

我们也可以计算总的预期寿命损失,即每个原因造成的月数损失的总和,它量化了病人从时间点 0 到预先定义的时间点 所损失的平均寿命月数。尽管这也可以在拟合全因模型后获得,但在此我们展示如何在拟合特定病因模型后获得估计值:

/* total months lost (from both prostate and other causes) under placebo,
by setting the first two # in lincom() to 1 and the following two in 0 */
// with standsurv use options rmft and cif and crmodels()
standsurv, crmodels(prostate other) cif rmft ///
at1(rx 0) at2(rx 1) timevar(t_rmft60) lincom(1 1 0 0) ci ///
atvar(RMLT0b RMLT1b) lincomvar(RMLT_total0)

lincom(#…#):计算 atn() 选项的线性组合,而不用上述 contrast() 选项,来获得所有原因造成的总损失月数。对于安慰剂下的总损失寿命,lincom() 中前两个 # 对应 at1(),应设置为 1 (指由于前列腺癌导致的损失月数和由于其他原因导致的损失月数之和) 。

在安慰剂下,自随机化后 60 个月内,总的寿命损失为 25.8 个月:

// list
li RMLT_total0* in 1, noobs abb(22)
+-------------------------------------------------+
| RMLT_total0 RMLT_total0_lci RMLT_total0_uci |
|-------------------------------------------------|
| 25.751636 22.256293 29.246978 |
+-------------------------------------------------+

同样,通过将 lincom() 中最后两个 # 设置为1,对应 at2() ,可以得到在 DES 下总损失寿命 (指由于前列腺癌导致的损失月数和由于其他原因导致的损失月数之和):

/* total months lost (from both prostate and other causes) under DES,
by setting the first two # in lincom() to 0 and the following two in 1 */
standsurv, crmodels(prostate other) cif rmft ///
at1(rx 0) at2(rx 1) timevar(t_rmft60) lincom(0 0 1 1) ci ///
atvar(RMFT0c RMFTc) lincomvar(RMFT_total1)

在 DES 下,自随机化后 60 个月内,总的寿命损失为 26.7 个月:

// list
li RMFT_total1* in 1, noobs abb(22)
+-------------------------------------------------+
| RMFT_total1 RMFT_total1_lci RMFT_total1_uci |
|-------------------------------------------------|
| 26.728941 23.225487 30.232395 |
+-------------------------------------------------+

这导致安慰剂和 DEs 之间的总损失相差约 1 个月,与前面为每个具体原因计算的差异之和相同 (-3.2 和4.2)。

3.2 直接影响

考虑一种假设的干预措施,设定 ,即一种消除其他原因造成的竞争性死亡的干预措施。在这种干预下,不同水平的治疗之间的反事实对比是受控的直接效应,它量化了治疗对感兴趣的事件的影响,而不是由竞争事件来介导。

在剔除竞争事件的干预措施下,当治疗为 时,边际反事实的死亡概率为:

其中, 表示死于前列腺癌的净概率。如果消除了竞争事件,前列腺癌死亡净概率的平均因果差异就被定义为:

如果排除了竞争事件,在 DES 和安慰剂的作用下,前列腺癌死亡的净概率以及它们之间的差异可以通过回归标准化得到。只考虑前列腺癌死亡的特定病因模型 (该模型的估计值已储存在内存的 prostate 中) 如下:

estimates restore prostate

/* obtain the standardised net probability of prostate cancer death under
placebo and under DES as well as their difference */
// with standsurv use option failure
standsurv, failure at1(rx 0) at2(rx 1) timevar(timevar) contrast(difference) ci ///
atvars(F_net_prostate0 F_net_prostate1) contrastvars(F_net_prostate_diff)

failure: 用于计算因前列腺癌死亡的标准化失效函数 ,即前列腺癌的标准化死亡概率。

下图显示了 DES 和安慰剂的标准化净死亡概率,以及它们在随机化后的时间差异。随机化 60 个月后,DES 的标准化净死亡概率为 34%,安慰剂为 38%,差异为 4%。与前面相比,这些估计假设前列腺癌是唯一可能的死因,而且不可能死于其他原因。这样的解释可能具有挑战性,然而它可以捕捉到治疗对前列腺癌死亡率的直接影响,这种影响不受竞争性事件的中介。

3.3 可分离的影响

在治疗效果可以分解为不同成分的情况下,可以估计可分离的效果。假设治疗 可以被概念化为有两个二元成分,通过不同的因果途径发挥作用:一个是影响感兴趣事件的治疗成分 ,另一个是影响竞争事件的治疗成分 。则治疗对癌症死亡概率的可分离的直接影响定义为:

也就是说,当影响竞争事件的治疗成分 被设定为常数 时,影响感兴趣事件的治疗成分的效果。类似地,我们可以将治疗对癌症死亡概率的可分离间接效应定义为:

也就是说,当影响感兴趣事件的治疗成分被设定为一个常数时,影响竞争事件的治疗成分的效果。可分离效应虽然不涉及消除竞争事件的假设性干预,但也假设了一种干预措施,即在治疗的每个组成部分中分配不同的治疗值。为了估计前列腺癌例子的可分离效应,我们需要一个治疗变量的副本,使得在 standsurv 中分别操作这些变量。

// make copy of treatment variable, so can manipulate separately in standsurv
gen rx_c = rx
gen rx_o = rx

拟合包含变量 rx_crx_o 的特定原因生存模型:

// fit cause specific models (exit at 60)
// for prostate
stset dtime, failure(eventType==1) exit(time 60)
stpm2 normalAct ageCat2 ageCat3 hx hgBinary rx_c, scale(hazard) df(4) tvc(rx_c) dftvc(2)
estimates store prostate

// for other causes
stset dtime, failure(eventType==2) exit(time 60)
stpm2 normalAct ageCat2 ageCat3 hx hgBinary rx_o, scale(hazard) df(3)
estimates store other

命令 stpm2 的介绍如下:

// 命令安装
ssc install stpm2 //拟合灵活的参数生存模型

// 命令语法
stpm2 [varlist] [if] [in] [, options]

其中,

  • scale():指定生存模型拟合的比例;
  • df(#):基线危险函数的自由度;
  • eform:取幂系数;
  • tvc(varlist):定义与时间有关的变量;
  • dftvc(df_list):每个时间依赖效应的自由度。

灵活的参数生存模型 (Flexible Parametric Survival Models,FPMs) 不支持因子变量,必须在拟合模型前生成虚拟变量。参数估计与之前的模型相同,故不再展示。支持的生存模型估计命令包括 stpm2strcsstreg,如果想了解更多关于生存模型的介绍,可参考连享会相关推文「Stata: 生存分析一文读懂」。

使用 standsurv 命令,并加入更多的 atn() 选项,得到可分离的直接和间接效应:

/* obtain the standardised CIFs under DES, under placebo and if the effect of
other causes of death is removed as well as the differences (separable effects) */
// use option cif and crmodels()
standsurv, crmodels(prostate other) cif timevar(timevar) contrast(difference) ci ///
at1(rx_c 1 rx_o 1) at2(rx_c 1 rx_o 0) at3(rx_c 0 rx_o 0) ///
atvars(F_rx11 F_rx10 F_rx00) contrastvars(F_diff_indirect F_diff_total)

下图显示了在 DES 和安慰剂治疗下,标准化前列腺癌死亡累积发生率分别为 14.5% 和 21.7%,在去除其他死因影响的假设治疗下为 15.6% (实线)。当去除其他死因的影响后,前列腺癌死亡的累积发生率 (蓝线) 与 DES 下前列腺癌死亡的累积发生率 (黑线) 非常接近,表明治疗效果主要由其对前列腺癌死亡率的影响驱动。

下图给出了 DES 和安慰剂下前列腺癌死亡累积发生率的标准化总差异以及可分离的间接差异,与随机化后时间的关系。可分离的间接效应随着时间的推移而增加,但在整个随访过程中仍然很低。在诊断后的 36 个月 (3 年),前列腺癌死亡的标准化累积发病率的总差异为 7.2%,而间接效应的估计值为 1.1%。这相当于与安慰剂相比, DES 下前列腺癌死亡率的降低是由于 DES 对其他原因死亡率的影响。因此,治疗对前列腺癌死亡率的总影响并不受对其他原因死亡影响的高度影响。

4. 参考资料

  • Syriopoulou E, Mozumder S I, Rutherford M J, et al. Estimating causal effects in the presence of competing events using regression standardisation with the Stata command standsurv[J]. arXiv preprint arXiv:2109.03628, 2021. -PDF-
  • Kaplan E L, Meier P. Nonparametric estimation from incomplete observations[J]. Journal of the American statistical association, 1958, 53(282): 457-481. -PDF-
  • Stensrud M J, Young J G, Didelez V, et al. Separable effects for causal inference in the presence of competing events[J]. Journal of the American Statistical Association, 2020: 1-9. -PDF-
  • Mozumder S I, Rutherford M J, Lambert P C. Estimating restricted mean survival time and expected life-years lost in the presence of competing risks within flexible parametric survival models[J]. BMC Medical Research Methodology, 2021, 21(1): 1-20. -PDF-
  • Young J G, Stensrud M J, Tchetgen Tchetgen E J, et al. A causal framework for classical statistical estimands in failure‐time settings with competing events[J]. Statistics in medicine, 2020, 39(8): 1199-1236. -PDF-

5. 相关推文

Note:产生如下推文列表的 Stata 命令为:
lianxh 生存分析, m
安装最新版 lianxh 命令:
ssc install lianxh, replace

  • 专题:计量专题
    • Stata: 生存分析一文读懂
  • 专题:生存分析
    • mrsprep:边际相对生存分析的数据准备
    • stata生存分析:stset命令详解

连享会 · 文本分析 | 爬虫 | 机器学习

New! Stata 搜索神器:lianxhsongbl  GIF 动图介绍
搜: 推文、数据分享、期刊论文、重现代码 ……
👉 安装:
. ssc install lianxh
. ssc install songbl
👉  使用:
. lianxh DID 倍分法
. songbl all

🍏 关于我们

  • 连享会 ( www.lianxh.cn,推文列表) 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 直通车: 👉【**百度一下:**连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。


您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存