/*************************************************************************** This code available at: http://pauldickman.com/talk/proportional-hazards-cbb-march2022/melanoma.do Code used in the seminar: http://pauldickman.com/talk/proportional-hazards-cbb-march2022/ Cause-specific survival among patients diagnosed with melanoma - Fit main effects model to study the effect of sex (adjusted for age, period, stage) - Assess proportional hazards assumption (test and Schoenfeld residual plots) - fit model where PH assumption is relaxed for sex - use stpm2 to estimate time-varying effect of sex and plot it Paul Dickman, March 2022 ***************************************************************************/ // Read data for patients diagnosed with melanoma // exclude unknown stage (stage==0) use https://pauldickman.com/data/melanoma.dta if stage>0, clear // rename to a more descriptive name rename year8594 period // create dummy variable for sex generate male=(sex==1) label define male 0 "Female" 1 "Male" label values male male // stset for cause-specific survival (status==1 is death due to melanoma) stset surv_mm, fail(status==1) scale(12) exit(time 120) // estimate effect of sex stcox i.male i.period i.agegrp i.stage, noshow nolog estimates store ph // test proportional hazards assumption // this is a test of nonzero slope in the linear regression of the // scaled Schoenfeld residuals on time estat phtest, detail // plot of scaled Schoenfeld residuals on time for sex estat phtest, plot(1.male) name(male_ph, replace) scheme(plotplain) graph export "melanoma_male_ph_plot.pdf", replace // plot of scaled Schoenfeld residuals on time for distant (v localised) estat phtest, plot(3.stage) name(stage_ph, replace) scheme(plotplain) graph export "melanoma_stage_ph_plot.pdf", replace // relax the proportional hazards assumptiom for sex // assume log(HR) varies as a linear function of time stcox i.male i.period i.agegrp i.stage, tvc(i.male) texp(_t) noshow nolog estimates store nonph // compare the parameter estimates from PH nand non-PH models estimates table ph nonph, equations(1) eform b(%9.6f) modelwidth(12) // Test PH assumption using likelihood ratio test lrtest nonph ph // estimate the time-varying hazard ratio using flexible parametric models stpm2 male i.period i.agegrp i.stage, scale(h) df(5) eform tvc(male) dftvc(3) estimates store stpm2_nonph // compare the parareter estimates for Cox and stpm2 (they will be very similar) estimates table ph nonph stpm2_nonph, equations(1) eform b(%9.6f) modelwidth(12) // predict the time-varying HR for 51 values of time between 0 and 10 range temptime 0 10 51 predict hr, hrnumerator(male 1) ci timevar(temptime) // plot the time-varying HR twoway (rarea hr_lci hr_uci temptime, color(red%25)) /// (line hr temptime, sort lcolor(red)) /// , legend(off) ysize(8) xsize(11) scheme(plotplain) yscale(log) /// ytitle("Hazard ratio (male/female)") name("hrtvc", replace) /// xtitle("Years since diagnosis") graph export "melanoma_male_hr_tvc.pdf", replace // empirical hazards as a function of follow-up time sts graph, hazard by(sex) kernel(epan2) name(h, replace) scheme(plotplain) // Wald test for time-varying effect of sex testparm _rcs_male*