import os
os.chdir('/Applications/stata/utilities/')
from pystata import config
config.init('mp')
___ ____ ____ ____ ____ ®
/__ / ____/ / ____/ 17.0
___/ / /___/ / /___/ MP—Parallel Edition
Statistics and Data Science Copyright 1985-2021 StataCorp LLC
StataCorp
4905 Lakeway Drive
College Station, Texas 77845 USA
800-STATA-PC https://www.stata.com
979-696-4600 stata@stata.com
Stata license: Single-user 4-core perpetual
Serial number: 501706307451
Licensed to: MALF
IBS.GRANADA
Notes:
1. Unicode is supported; see help unicode_advice.
2. More than 2 billion observations are allowed; see help obs_advice.
3. Maximum number of variables is set to 5,000; see help set_maxvar.
%%stata
cd /Users/MALF/Dropbox/CANSURV-slides-course/Granada/
use melanoma
stset surv_mm, failure(status=1 2) scale(12) exit(time 120.5) id(id)
gen _age = min(int(age + _t),99)
gen _year = int(yydx + _t)
sort _year sex _age
merge m:1 _year sex _age using popmort, keep(match master)
tab agegrp, gen(agegrp)
gen female = sex == 2
/* Fit initial model */
stpm2, df(3) scale(hazard) bhazard(rate)
predict h1, hazard per (1000) ci
predict s1, survival ci
. cd /Users/MALF/Dropbox/CANSURV-slides-course/Granada/
/Users/MALF/Dropbox/CANSURV-slides-course/Granada
. use melanoma
(Skin melanoma, diagnosed 1975-94, follow-up to 1995)
. stset surv_mm, failure(status=1 2) scale(12) exit(time 120.5) id(id)
Survival-time data settings
ID variable: id
Failure event: status==1 2
Observed time interval: (surv_mm[_n-1], surv_mm]
Exit on or before: time 120.5
Time for analysis: time/12
--------------------------------------------------------------------------
7,775 total observations
0 exclusions
--------------------------------------------------------------------------
7,775 observations remaining, representing
7,775 subjects
2,777 failures in single-failure-per-subject data
43,384.625 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 10.04167
. gen _age = min(int(age + _t),99)
. gen _year = int(yydx + _t)
. sort _year sex _age
. merge m:1 _year sex _age using popmort, keep(match master)
Result Number of obs
-----------------------------------------
Not matched 0
Matched 7,775 (_merge==3)
-----------------------------------------
. tab agegrp, gen(agegrp)
Age in 4 |
categories | Freq. Percent Cum.
------------+-----------------------------------
0-44 | 2,046 26.32 26.32
45-59 | 2,238 28.78 55.10
60-74 | 2,280 29.32 84.42
75+ | 1,211 15.58 100.00
------------+-----------------------------------
Total | 7,775 100.00
. gen female = sex == 2
.
. /* Fit initial model */
.
. stpm2, df(3) scale(hazard) bhazard(rate)
Iteration 0: log likelihood = -8731.8053
Iteration 1: log likelihood = -8590.4493
Iteration 2: log likelihood = -8590.0259
Iteration 3: log likelihood = -8590.0249
Iteration 4: log likelihood = -8590.0249
Log likelihood = -8590.0249 Number of obs = 7,775
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
xb |
_rcs1 | .8252308 .0249859 33.03 0.000 .7762595 .8742022
_rcs2 | .2110309 .0235833 8.95 0.000 .1648085 .2572534
_rcs3 | .0631928 .0109672 5.76 0.000 .0416974 .0846881
_cons | -1.813097 .0314253 -57.70 0.000 -1.87469 -1.751505
------------------------------------------------------------------------------
. predict h1, hazard per (1000) ci
. predict s1, survival ci
.
%%stata
twoway (rarea h1_lci h1_uci _t, sort pstyle(ci)) ///
(line h1 _t, sort), name(h1,replace) legend(off) ///
xtitle("time since diagnosis") ///
ytitle("excess mortality rate (per 1000 py's)")
twoway (rarea s1_lci s1_uci _t, sort pstyle(ci)) ///
(line s1 _t, sort), name(s1,replace) legend(off) ///
xtitle("time since diagnosis") ///
ytitle("relative survival")
.
. twoway (rarea h1_lci h1_uci _t, sort pstyle(ci)) ///
> (line h1 _t, sort), name(h1,replace) legend(off) ///
> xtitle("time since diagnosis") ///
> ytitle("excess mortality rate (per 1000 py's)")
.
. twoway (rarea s1_lci s1_uci _t, sort pstyle(ci)) ///
> (line s1 _t, sort), name(s1,replace) legend(off) ///
> xtitle("time since diagnosis") ///
> ytitle("relative survival")
.