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

发布时间：2021-05-31 阅读 2211
分享到微信 分享到微博 收藏文章

 

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 命令后，我们可以得到每个事件时间点的平均预期风险率，并扩展数据，从而在后续模型拟合中将时依权重纳入其中。通过扩展数据和计算加权预期死亡率和时依权重，就可以使用标准相对生存估计命令 (如 stpm2 和 strcs) 来估计边际相对生存率。

在此基础上，还可以纳入个体权重进行外部标准化，从而实现实现在年龄、性别上的可比。与此同时，在使用 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