温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装命令如下:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
⛳ Stata 系列推文:
作者:保瑞 (中南财经政法大学)
邮箱:br19965956@163.com
目录
mrsprep
是为拟合边际相对生存参数模型而进行的数据准备命令。它可将观测样本外推到观察数据之外的时间进行生存预测,并计算时依权重 (time-dependent weights),即预期生存的倒数。在运行 mrsprep
命令后,我们可以得到每个事件时间点的平均预期风险率,并扩展数据,从而在后续模型拟合中将时依权重纳入其中。通过扩展数据和计算加权预期死亡率和时依权重,就可以使用标准相对生存估计命令 (如 stpm2
和 strcs
) 来估计边际相对生存率。
在此基础上,还可以纳入个体权重进行外部标准化,从而实现实现在年龄、性别上的可比。与此同时,在使用 mrsprep
命令时,需要使用 stset
命令对数据结构进行设定。在数据扩展时,id()
是必选项。并且,在使用 mrsprep
命令前,每个 id 只对应一行数据,而在使用命令之后,会生成一个新的数据结构。
在使用癌症登记数据量化癌症患者的生存率时,通常需要进行相对生存分析,即估计边际相对生存 (Marginal relative survival)。
*命令安装
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:由于相对生存分析是通过将患者与一般人群进行比较,因此需要注意区分一般人群死亡率数据文件和患者数据文件。
需要注意的是,在运行 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 将自动生成四个新变量,即:
. *调用数据
. 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 号患者为例,相关描述如下:
stpm2
/ strcs
命令。在扩展数据和计算时依权重后,需要再次使用 stset
命令设置生存分析数据结构。此时,需要设定区间开始位置 tstop 和区间结束位置 tstart,以及权重 wt。
. stset tstop [iw=wt], enter(tstart) failure(event==1)
接着,使用灵活参数生存模型 (Flexible parametric survival models, FPSM) 进行估计。该模型由 Royston 和 Parmar 开发,目的是为了更好地平衡 Cox 模型和参数模型的缺陷。它的优势主要体现在:
具体地,我们使用 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")
结果如下图所示,表明患者的边际相对生存率逐年下降。
由于不同组别的患者年龄分布不同,因此需要分别计算不同年龄组的时依权重,从而对年龄进行标准化。
. 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)
结果如下图所示:
在前文基础上,本节主要介绍如何在相对生存分析中纳入协变量。
. 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)
结果如下图所示:
Note:产生如下推文列表的 Stata 命令为:
lianxh 生存分析 dataex, m
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会-常见问题解答:
✨ https://gitee.com/lianxh/Course/wikis
New!
lianxh
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh