mrsprep:边际相对生存分析的数据准备

发布时间:2021-05-31 阅读 2001

Stata连享会   主页 || 视频 || 推文 || 知乎

温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。

New! lianxh 命令发布了:
随时搜索推文、Stata 资源。安装命令如下:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh

课程详情 https://gitee.com/lianxh/Course

课程主页 https://gitee.com/lianxh/Course

⛳ Stata 系列推文:

PDF下载 - 推文合集

作者:保瑞 (中南财经政法大学)
邮箱br19965956@163.com


目录


1. 背景简介

mrsprep 是为拟合边际相对生存参数模型而进行的数据准备命令。它可将观测样本外推到观察数据之外的时间进行生存预测,并计算时依权重 (time-dependent weights),即预期生存的倒数。在运行 mrsprep 命令后,我们可以得到每个事件时间点的平均预期风险率,并扩展数据,从而在后续模型拟合中将时依权重纳入其中。通过扩展数据和计算加权预期死亡率和时依权重,就可以使用标准相对生存估计命令 (如 stpm2strcs) 来估计边际相对生存率。

在此基础上,还可以纳入个体权重进行外部标准化,从而实现实现在年龄、性别上的可比。与此同时,在使用 mrsprep 命令时,需要使用 stset 命令对数据结构进行设定。在数据扩展时,id() 是必选项。并且,在使用 mrsprep 命令前,每个 id 只对应一行数据,而在使用命令之后,会生成一个新的数据结构。

2. 概念界定

在使用癌症登记数据量化癌症患者的生存率时,通常需要进行相对生存分析,即估计边际相对生存 (Marginal relative survival)。

  • 相对生存分析:是指通过将患癌病人与没有患癌的一般人群进行比较,从而把因癌症引起的死亡和所有其他原因引起的死亡分开,以研究癌症对患者生存情况的影响的一种分析方法;
  • 相对生存率:是指癌症人群观察生存率与一般人群期望生存率的比值。

3. mrsprep 命令

*命令安装
ssc install mrsprep, replace
*命令语法
 mrsprep using filename[if] [in], agediag(varname) 
         datediag(varname) breaks(numlist) [options]

相关选项介绍:

  • agediag(varname):设定诊断年龄变量;
  • datediag(varname):设定诊断日期变量;
  • breaks(numlist):设定计算时依权重的分割点,并且间隔越窄,精确度越高;
  • by(varlist):按照括号里面的分类变量分别计算平均预期风险率;
  • indweights(varname):个体权重;
  • keep(varlist):保留扩展数据中需要的变量;
  • newframe(framename):给创建的新框架名;
  • pmage(varname):设定给一般人群死亡率文件里的年龄变量,默认名称是 _age,这个变量不能存放在患者数据文件里,必须放在一般人群死亡率文件;
  • pmother(varname):设定一般人群死亡率文件中的其他变量命名,通常包括性别或地区,所有列出的变量应同时放在患者数据集和一般人群死亡率文件里;
  • pmrate(varname):设定一般人群死亡率文件中的比率变量,默认为 rate,比率的单位为 %/每人年。若一般人群死亡率文件中仅有一年的生存概率数据,可以通过 rate = -ln(survprob) 计算得到这一比率,其中 survprob 表示生存概率;
  • pmmaxage(#):设定在一般人群死亡率文件中提供的一般人口死亡率的最大年龄。超过最大年龄则令其等于最大年龄,默认最大年龄为 99 岁;
  • pmmaxyear(#):设定在一般人群死亡率文件中提供的一般人口死亡率的最大年份,超过最大年份则令其等于最大年份;
  • pmyear(varname):设定在一般人群死亡率文件中的年份变量。默认为 _year,该变量不能放在患者数据文件里,而应该放在在一般人群死亡率文件中;
  • verbose:显示计算过程的细节信息。

Note:由于相对生存分析是通过将患者与一般人群进行比较,因此需要注意区分一般人群死亡率数据文件和患者数据文件。

4. Stata 实操

需要注意的是,在运行 mrsprep 之前,需要通过 stset 命令声明数据是生存分析数据。其语法如下:

stset timevar, origin(time exp) failure(failvar[==numlist]) 
      id(idvar) enter(time exp) exit(time exp)

其中,

  • timevar:是事件退出日期的变量;
  • origin():是说明一个样本从什么时候开始,事件可能发生,即样本可能面临发生风险的时间;
  • failure():选项用于指定事件变量;
  • id():是每一个样本的标识,在生存分析的数据中,有时候一个样本有几行观测值,这时候就需要 id 来帮助识别样本;
  • enter():主要说明一个样本从什么时候进入研究;
  • exit():主要说明一个样本从什么时候离开研究。

在运行 stset 命令之后,Stata 将自动生成四个新变量,即:

  • _t0:表示起始时间点,默认为 0,但如果有左截断,则 _t0 不为 0;
  • _st:表示数据的状态是否完好可用;
  • _d:表示事件的结果,是否发生;
  • _t:表示每个样本事件持续的时间。

4.1 数据拓展和时依权重计算

. *调用数据
. use "https://pclambert.net/data/melanoma.dta", clear
(Skin melanoma, diagnosed 1975-94, follow-up to 1995)
. *原始数据结构
. dataex id sex age mmdx yydx surv_mm surv_yy dx exit in 1/5

----------------------- copy starting from the next line --------------------
[CODE]
* Example generated by -dataex-. To install: ssc install dataex
clear
input int id byte(sex age mmdx) int yydx float(surv_mm surv_yy) int(dx exit)
1 2 81 2 1981  26.5  2.5 7703  8510
2 2 75 9 1975  55.5  4.5 5742  7432
3 2 78 2 1978 177.5 14.5 6626 12029
4 2 75 8 1975  29.5  2.5 5715  6613
5 2 81 7 1981  57.5  4.5 7860  9611
end
format %d dx
format %d exit
label values sex sex
label def sex 2 "Female", modify
[/CODE]
------------------ copy up to and including the previous line ---------------
Listed 5 out of 7775 observations
. *设置生存分析数据结构
. stset exit, origin(dx) failure(status=1,2) id(id) ///
>       exit(time dx + 10*365.24) scale(365.24)
. *生成5个新的变量_st、_d、_origin、_t、_t0 
. dataex id sex age mmdx yydx surv_mm surv_yy dx exit ///
>        _st _d _origin _t _t0 in 1/5

----------------------- copy starting from the next line ---------
[CODE]
* Example generated by -dataex-. To install: ssc install dataex
clear
input int id byte(sex age mmdx) int yydx float(surv_mm surv_yy) 
int(dx exit) byte(_st _d) int _origin double _t byte _t0
1 2 81 2 1981  26.5  2.5 7703  8510 1 1 7703 2.2095060781951594 0
2 2 75 9 1975  55.5  4.5 5742  7432 1 1 5742  4.627094513196802 0
3 2 78 2 1978 177.5 14.5 6626 12029 1 0 6626  9.999999999999998 0
4 2 75 8 1975  29.5  2.5 5715  6613 1 1 5715  2.458657321213449 0
5 2 81 7 1981  57.5  4.5 7860  9611 1 1 7860   4.79410798379148 0
end
format %d dx
format %d exit
label values sex sex
label def sex 2 "Female", modify
[/CODE]
------------------ copy up to and including the previous line ----
Listed 5 out of 7775 observations
. *扩展数据和计算权重,以0.2为分割点进行数据扩展
. mrsprep using "https://pclambert.net/data/popmort.dta", ///
>         pmother(sex) agediag(age) datediag(dx) verbose  ///
>                                 breaks(0(0.2)10)
. *展示id=1的患者
. dataex id age dx sex meanhazard_wt tstart tstop event t ///
>        wt if id == 1 

----------------------- copy starting from the next line ------------------
[CODE]
* Example generated by -dataex-. To install: ssc install dataex
clear
input double(id age dx sex meanhazard_wt tstart tstop event t wt)
1 81 7703 2         999   0    .2 0 2.2095060781951594 1.0144102352694586
1 81 7703 2         999  .2    .4 0 2.2095060781951594 1.0438566628056225
1 81 7703 2         999  .4    .6 0 2.2095060781951594  1.073635007266397
1 81 7703 2         999  .6    .8 0 2.2095060781951594 1.1037253339775324
1 81 7703 2         999  .8     1 0 2.2095060781951594 1.1346589898978077
1 81 7703 2         999   1     1.2 0 2.2095060781951594 1.1663819557285744
1 81 7703 2         999 1.2     1.4 0 2.2095060781951594 1.1989120162714493
1 81 7703 2         999 1.4     1.6 0 2.2095060781951594 1.2323493309378348
1 81 7703 2         999 1.6     1.8 0 2.2095060781951594 1.2667192027868366
1 81 7703 2         999 1.8     2 0 2.2095060781951594 1.3032625252158494
1 81 7703 2         999   2     2.2 0 2.2095060781951594 1.3414414994396304
1 81 7703 2        .028 2.2   2.209 1 2.2095060781951594 1.3615266152925936
end
[/CODE]
------------------ copy up to and including the previous line -------------
Listed 12 out of 219460 observations

以为 1 号患者为例,相关描述如下:

  • 1 号个体为女性,诊断年龄为 81 岁,对应 12 个不同区间;
  • tstart 表示区间开始位置,tstop 表示区间结束位置;
  • 在最后一个区间,患者在诊断后 2.209 年后死亡 (event = 1);
  • 每个时间区间都 0.2 年;
  • 每个区间对应一个时依权重 (wt),等于预期生存的倒数。因此随着预期生存时间的缩短,时依权重随时间增加;
  • meanhazard_wt 表示子在每个个体时间点上的加权预期生存率,分割区间对应的加权预期生存率设定为 999;
  • 拟合相对设生存模型时,需要用到 stpm2 / strcs 命令。

在扩展数据和计算时依权重后,需要再次使用 stset 命令设置生存分析数据结构。此时,需要设定区间开始位置 tstop 和区间结束位置 tstart,以及权重 wt

. stset tstop [iw=wt], enter(tstart) failure(event==1)

接着,使用灵活参数生存模型 (Flexible parametric survival models, FPSM) 进行估计。该模型由 Royston 和 Parmar 开发,目的是为了更好地平衡 Cox 模型和参数模型的缺陷。它的优势主要体现在:

  • 利用到了基线风险信息,可以更容易地进行生存预测;
  • 相比于 Cox 模型,它可以外推到观察数据之外的时间进行生存预测;
  • 与 Cox 模型相比,该模型得到的是一个相对平滑的关于基线风险函数的变化趋势,不会出现因过拟合而剧烈波动的情况 (闫沛静等,2020);
  • 相比于其他参数模型,它更加灵活。

具体地,我们使用 stpm2 命令拟合边际模型,并使用时依权重聚类稳健标准误 vce(cluster id)

. *命令安装
. ssc install stpm2, replace
. ssc install rcsgen, replace
. stpm2, scale(hazard) df (5) bhazard(meanhazard_wt) vce(cluster id)

note: delayed entry models are being fitted
Iteration 0:   log pseudolikelihood = -11548.878  
Iteration 1:   log pseudolikelihood = -11339.064  
Iteration 2:   log pseudolikelihood = -11338.253  
Iteration 3:   log pseudolikelihood =  -11338.25  
Iteration 4:   log pseudolikelihood =  -11338.25  
Log pseudolikelihood =  -11338.25               Number of obs     =    219,460
                                 (Std. Err. adjusted for 7,775 clusters in id)
------------------------------------------------------------------------------
             |               Robust
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
xb           |
       _rcs1 |      0.981      0.030    33.24   0.000        0.923       1.039
       _rcs2 |      0.177      0.017    10.19   0.000        0.143       0.211
       _rcs3 |      0.072      0.011     6.51   0.000        0.050       0.093
       _rcs4 |      0.004      0.008     0.54   0.586       -0.011       0.019
       _rcs5 |      0.007      0.006     1.32   0.187       -0.004       0.019
       _cons |     -1.957      0.034   -56.74   0.000       -2.024      -1.889
------------------------------------------------------------------------------
. *-预测边际相对生存并画图
. twoway (line s_mrs* tt, lcolor(red..) lpattern(solid dash dash)), ///
>        legend(off) ylabel(0.6(0.1)1, format(%3.1f))  ///
>        ytitle("Marginal relative survival")          ///
>        xtitle("Years from diagnosis")

结果如下图所示,表明患者的边际相对生存率逐年下降。

4.2 加入个体权重进行外部标准化

由于不同组别的患者年龄分布不同,因此需要分别计算不同年龄组的时依权重,从而对年龄进行标准化。


. use "https://pclambert.net/data/melanoma.dta", clear
(Skin melanoma, diagnosed 1975-94, follow-up to 1995)

. *使用exit()选项来限制样本为10年的跟踪观察样本
. stset exit, origin(dx) failure(status=1,2) id(id) ///
>       exit(time dx + 10*365.24) scale(365.24)

. *根据国际癌症存活标准人口对年龄进行标准化,并生成ICSS权重
. drop agegrp
. egen agegrp=cut(age), at(0 45 55 65 75 200) icodes
. replace agegrp = agegrp + 1
. label variable agegrp "Age group"
. label define agegrplab 1 "0-44" 2 "45-54" 3 "55-64" 4 "65-74" 5 "75+", replace
. label values agegrp agegrplab
. recode agegrp (1=0.28) (2=0.17) (3=0.21) (4=0.20) (5=0.14), gen(ICSSwt)

. *用indweights()选项设定个体层面的权重
. local total= _N
. bysort agegrp: gen a_age = _N/`total'
. gen double wt_age = ICSSwt/a_age
. mrsprep using "https://pclambert.net/data/popmort.dta", ///
>         pmother(sex) agediag(age) datediag(dx) verbose  ///
>         breaks(0(0.2)10) indweights(wt_age)             ///
>         newframe(mrs_stand, replace)
. *预测边际相对生存并画图
. stset tstop [iw=wt], enter(tstart) failure(event==1)
. stpm2, scale(hazard) df (5) bhazard(meanhazard_wt) vce(cluster id)
. range tt 0 10 101
. predict s_mrs, surv timevar(tt) ci
. twoway (line s_mrs* tt, lcolor(red..) lpattern(solid dash dash)),   ///
>        legend(off) ylabel(0.6(0.1)1, format(%3.1f)) ytitle(Marginal ///
>        relative survival) xtitle(Years from diagnosis)

结果如下图所示:

4.3 拟合累积风险函数参数模型

在前文基础上,本节主要介绍如何在相对生存分析中纳入协变量。

. use "https://pclambert.net/data/melanoma.dta", clear
. stset exit, origin(dx) failure(status=1,2)  ///
>       id(id) exit(time dx + 10*365.24) scale(365.24)
. drop agegrp
. egen agegrp=cut(age), at(0 45 55 65 75 200) icodes
. replace agegrp = agegrp + 1
. label variable agegrp "Age group"
. label define agegrplab 1 "0-44" 2 "45-54" 3 "55-64" 4 "65-74" 5 "75+", replace
. label values agegrp agegrplab
. recode agegrp (1=0.28) (2=0.17) (3=0.21) (4=0.20) (5=0.14), gen(ICSSwt)

. *分别计算男性和女性的ICSS权重
. gen female = sex == 2
. bysort female: egen totalsex = total(sex)
. bysort agegrp female: gen a_age_sex = _N/totalsex
. gen double wt_age_sex = ICSSwt/a_age_sex

. *使用by(female)选项分别估计男性和女性样本
. mrsprep using "https://pclambert.net/data/popmort.dta",        ///
>         pmother(sex) agediag(age) datediag(dx) pmmaxyear(2000) ///
>         verbose breaks(0(0.2)10) indweights(wt_age_sex) by(female)
. *预测边际相对生存并画图
. stset tstop [iw=wt], enter(tstart) failure(event==1)
. stpm2 female, scale(hazard) df (5) bhazard(meanhazard_wt) vce(cluster id)
. stpm2 female, scale(hazard) df (5) bhazard(meanhazard_wt)  ///
>       vce(cluster id) tvc(female) dftvc(3)
. range tt 0 10 101
. predict s_mrs_male, surv timevar(tt) ci at(female 0)
. predict s_mrs_female, surv timevar(tt) ci at(female 1)
. twoway (line s_mrs_male* tt, lcolor(red..)                    ///
>        lpattern(solid dash dash)) (line s_mrs_female* tt,     ///
>        lcolor(blue..) lpattern(solid dash dash)), legend(off) ///
>        ylabel(0.6(0.1)1, format(%3.1f)) ytitle(Marginal       ///
>        relative survival) xtitle(Years from diagnosis)    

结果如下图所示:

5. 参考文献

  • Lambert P. STPM2: Stata module to estimate flexible parametric survival models[J]. 2020. -Link-
  • Lambert P C, Syriopoulou E, Rutherford M R. Direct modelling of age standardized marginal relative survival through incorporation of time-dependent weights[J]. BMC medical research methodology, 2021, 21(1): 1-11. -Link-
  • Rutherford M J, Dickman P W, Coviello E, et al. Estimation of age-standardized net survival, even when age-specific data are sparse[J]. Cancer epidemiology, 2020, 67: 101745. -PDF-
  • 闫沛静, 李镜, 冉孟冬, 惠旭, 侯利莎, 杨克虎, 朱彩蓉. 灵活参数生存分析模型简介及应用[J]. 中国卫生统计, 2020, 37(05):776-779. -Link-

6. 相关推文

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

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,400+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxh 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh