// Graph settings grstyle init grstyle set plain, horizontal grid grstyle set legend 4, nobox grstyle set color 538 grstyle set horizontal grstyle set margin "0mm 0mm 0mm 0mm": twoway // Format data use https://www.pclambert.net/data/colonsim_stage, clear // Format datediag to display as a date format datediag %td // Restrict analysis to patients aged 18-99 at diagnosis keep if agediag>=18 & agediag<=99 // Recode the sex variable to match the popmort file replace sex = sex+1 label define label_sex 1 "Male" 2 "Female" label values sex label_sex // Prepare data for model fitting stset t, failure(dead=1) id(id) exit(time 5) // Attained age gen age = min(floor(agediag + _t),100) // Attained calendar year gen year = min(floor(yeardiag + _t),2016) // Merge in life tables merge m:1 age year dep sex using https://www.pclambert.net/data/popmort_uk_2017, /// keep(match master) keepusing(rate) // Non-linear function for age rcsgen agediag, gen(agercs) df(3) orthog global ageknots `r(knots)' matrix Rage =r(R) // Dummy variables tab stage, gen(stage) gen female = sex == 2 gen dep5 = dep == 5 // Interaction terms between stage and deprivation forvalues i = 2/4 { gen stage`i'dep5 = stage`i'*dep5 } // Fit the model stpm2 agercs* female dep5 stage2 stage3 stage4 stage?dep5, scale(hazard) df(5) /// tvc(agercs* dep5 stage2 stage3 stage4) dftvc(3) bhazard(rate) // Marginal range t5 0 5 51 standsurv, verbose /// Display output at1(.) /// Use observed covariate values crudeprob /// Crude probabilities of death timevar(t5) /// Time points used for predictions atvar(marg) /// New variable containing the predictions expsurv(using(https://www.pclambert.net/data/popmort_uk_2017) /// Popmort file agediag(agediag) /// Age at diagnosis in the dataset datediag(datediag) /// Date of diagnosis in the dataset pmage(age) /// Age variable in the popmort file pmyear(year) /// Year variable in the popmort file pmother(dep sex) /// Other variables included in the popmort file pmrate(rate) /// Rate variable in the popmort file pmmaxage(100) /// Maximum age in the popmort file pmmaxyear(2016) /// Maximum year in the popmort file ) // 2 separate graphs line marg_disease t5, sort col("3 144 214") xtitle("Years since Diagnosis") /// ytitle("Probability of Death") ylabel(,format(%3.1f)) title("Cancer") name(marg_disease,replace) line marg_other t5, sort col("254 48 11") xtitle("Years since Diagnosis") /// ytitle("Probability of Death") ylabel(,format(%3.1f)) title("Other Causes") name(marg_other,replace) ylabel(0(0.1)0.4) graph combine marg_disease marg_other // Stacked graph gen alive = 1 gen marg_all = marg_disease + marg_other twoway(area alive t5, sort col("120 172 68") fintensity(30)) /// (area marg_all t5, sort fintensity(30) col("3 144 214")) /// (area marg_other t5, sort fintensity(30) col("254 48 11") /// legend(order(1 "Alive" 2 "Death from Cancer" 3 "Death from Other Causes") /// rows(1) ring(1) pos(6)) xtitle("Years since Diagnosis") /// ytitle("Probability") ylabel(,format(%3.1f))) // Deprivation standsurv, verbose /// Display output at1(dep5 0 stage2dep5 0 stage3dep5 0 stage4dep5 0 ) /// Least deprived at2(dep5 1 stage2dep5=stage2 stage3dep5=stage3 stage4dep5=stage4) /// Most deprived atvar(dep_1 dep_5) /// New variables containing the predictions contrast(difference) /// Calculate the difference between groups contrastvar(diff_dep) /// New variables containing the difference ci /// Calculate confidence intervals crudeprob /// Crude probabilities of death timevar(t5) /// Time points used for predictions expsurv(using(https://www.pclambert.net/data/popmort_uk_2017) /// Popmort file agediag(agediag) /// Age at diagnosis in the dataset datediag(datediag) /// Date of diagnosis in the dataset pmage(age) /// Age variable in the popmort file pmyear(year) /// Year variable in the popmort file pmother(dep sex) /// Other variables included in the popmort file pmrate(rate) /// Rate variable in the popmort file pmmaxage(100) /// Maximum age in the popmort file pmmaxyear(2016) /// Maximum year in the popmort file at1(dep 1) /// Use expected rates for least deprived at2(dep 5) /// Use expected rates for most deprived ) // Crude probabilities for each deprivation group line dep_1_disease t5, sort col("3 144 214") || /// line dep_5_disease t5, sort col("3 144 214") lpattern(-) || /// line dep_1_other t5, sort col("254 48 11") || /// line dep_5_other t5, sort col("254 48 11") lpattern(-) xtitle("Years since Diagnosis") /// ytitle("Probability of Death") ylabel(,format(%3.1f)) /// legend(order(1 "Cancer: Least Deprived" 2 "Cancer: Most Deprived" /// 3 "Other Causes: Least Deprived" 4 "Other Causes: Most Deprived") ring(0) pos(11)) /// ysc(r(0 0.5)) ylabel(0(0.1)0.5) // Graph of the differences twoway(rarea diff_dep_disease_lci diff_dep_disease_uci t5, sort col("3 144 214") fintensity(30)) /// (line diff_dep_disease t5, sort col("3 144 214") legend(off) xtitle("Years since Diagnosis") /// ytitle("Difference in Probability of Death") ylabel(0(0.01)0.05) ysc(r(-0.005 0.05)) /// title("Cancer") ylabel(,format(%3.2f)) name(diff_cancer,replace)) twoway(rarea diff_dep_other_lci diff_dep_other_uci t5, sort col("254 48 11") fintensity(30)) /// (line diff_dep_other t5, sort col("254 48 11") legend(off) xtitle("Years since Diagnosis") /// ytitle("Difference in Probability of Death") ylabel(0(0.01)0.05) ysc(r(-0.005 0.05)) /// title("Other Causes") ylabel(,format(%3.2f)) name(diff_other,replace)) graph combine diff_cancer diff_other // Patient 2510 standsurv if id==2510, /// Only include Patient 2510 verbose /// Display output at1(.) /// Use observed covariate values crudeprob /// Crude probabilities of death timevar(t5) /// Time points used for predictions atvar(id2510) /// New variable containing the predictions expsurv(using(https://www.pclambert.net/data/popmort_uk_2017) /// Popmort file agediag(agediag) /// Age at diagnosis in the dataset datediag(datediag) /// Date of diagnosis in the dataset pmage(age) /// Age variable in the popmort file pmyear(year) /// Year variable in the popmort file pmother(dep sex) /// Other variables included in the popmort file pmrate(rate) /// Rate variable in the popmort file pmmaxage(100) /// Maximum age in the popmort file pmmaxyear(2016) /// Maximum year in the popmort file ) gen all_id2510 = id2510_disease + id2510_other twoway(area alive t5, sort col("120 172 68") fintensity(30)) /// (area all_id2510 t5, sort fintensity(30) col("3 144 214")) /// (area id2510_other t5, sort fintensity(30) col("254 48 11") /// legend(order(1 "Alive" 2 "Death from Cancer" 3 "Death from Other Causes") /// rows(1) ring(1) pos(6)) xtitle("Years since Diagnosis") ytitle("Probability") /// title("Patient 2510") subtitle("Male, Stage 2, Deprivation Group 1, Diagnosed Aged 85 in 2011") /// ylabel(,format(%3.1f))) // Effect of age gen temp_agediag = . gen temp_datediag = mdy(1,1,2011) foreach age in 60 70 80 90 { replace temp_agediag = `age' rcsgen, scalar(`age') knots($ageknots) rmatrix(Rage) gen(c) standsurv if _n==1, verbose /// Display output at1(agercs1 `=c1' agercs2 `=c2' agercs3 `=c3' /// Specify covariate pattern female 0 dep5 0 stage2 1 stage3 0 stage4 0 /// stage2dep5 0 stage3dep5 0 stage4dep5 0) /// crudeprob /// Crude probabilities of death timevar(t5) /// Time points used for predictions atvar(age`age') /// New variable containing the predictions expsurv(using(https://www.pclambert.net/data/popmort_uk_2017) /// Popmort file agediag(temp_agediag) /// Temporary age at diagnosis variable datediag(temp_datediag) /// Temporary date of diagnosis variable pmage(age) /// Age variable in the popmort file pmyear(year) /// Year variable in the popmort file pmother(dep sex) /// Other variables included in the popmort file pmrate(rate) /// Rate variable in the popmort file pmmaxage(100) /// Maximum age in the popmort file pmmaxyear(2016) /// Maximum year in the popmort file at1(dep 1 sex 1) /// Use expected rates for least deprived and male ) // All-cause probability of death gen all_age`age' = age`age'_disease + age`age'_other twoway(area alive t5, sort col("120 172 68") fintensity(30)) /// (area all_age`age' t5, sort fintensity(30) col("3 144 214")) /// (area age`age'_other t5, sort fintensity(30) col("254 48 11") /// legend(order(1 "Alive" 2 "Death from Cancer" 3 "Death from Other Causes") /// rows(1) ring(1) pos(6))xtitle("Years since Diagnosis") ytitle("Probability") /// title("Age `age'") ylabel(,format(%3.1f)) name(age`age',replace)) } grc1leg age60 age70 age80 age90, title("Male, Stage 2, Deprivation Group 1, Diagnosed in 2011") // Effect of stage replace temp_agediag = 75 replace temp_datediag = mdy(1,1,2011) rcsgen, scalar(75) knots($ageknots) rmatrix(Rage) gen(c) standsurv if _n==1, verbose /// at1(agercs1 `=c1' agercs2 `=c2' agercs3 `=c3' female 1 dep5 1 /// Stage 1 stage2 0 stage3 0 stage4 0 stage2dep5 0 stage3dep5 0 stage4dep5 0) /// at2(agercs1 `=c1' agercs2 `=c2' agercs3 `=c3' female 1 dep5 1 /// Stage 2 stage2 1 stage3 0 stage4 0 stage2dep5 1 stage3dep5 0 stage4dep5 0) /// at3(agercs1 `=c1' agercs2 `=c2' agercs3 `=c3' female 1 dep5 1 /// Stage 3 stage2 0 stage3 1 stage4 0 stage2dep5 0 stage3dep5 1 stage4dep5 0) /// at4(agercs1 `=c1' agercs2 `=c2' agercs3 `=c3' female 1 dep5 1 /// Stage 4 stage2 0 stage3 0 stage4 1 stage2dep5 0 stage3dep5 0 stage4dep5 1) /// crudeprob /// Crude probabilities of death timevar(t5) /// Time points used for predictions atvar(stage_1 stage_2 stage_3 stage_4) /// New variables containing the predictions expsurv(using(https://www.pclambert.net/data/popmort_uk_2017) /// Popmort file agediag(temp_agediag) /// Temporary age at diagnosis variable datediag(temp_datediag) /// Temporary date of diagnosis variable pmage(age) /// Age variable in the popmort file pmyear(year) /// Year variable in the popmort file pmother(dep sex) /// Other variables included in the popmort file pmrate(rate) /// Rate variable in the popmort file pmmaxage(100) /// Maximum age in the popmort file pmmaxyear(2016) /// Maximum year in the popmort file at1(dep 5 sex 2) /// Use expected rates for most deprived and female at2(dep 5 sex 2) /// Use expected rates for most deprived and female at3(dep 5 sex 2) /// Use expected rates for most deprived and female at4(dep 5 sex 2) /// Use expected rates for most deprived and female ) forvalues stage = 1/4 { gen all_stage_`stage' = stage_`stage'_disease + stage_`stage'_other twoway(area alive t5, sort col("120 172 68") fintensity(30)) /// (area all_stage_`stage' t5, sort fintensity(30) col("3 144 214")) /// (area stage_`stage'_other t5, sort fintensity(30) col("254 48 11") /// legend(order(1 "Alive" 2 "Death from Cancer" 3 "Death from Other Causes") /// rows(1) ring(1) pos(6)) xtitle("Years since Diagnosis") ytitle("Probability") /// title("Stage `stage'") ylabel(,format(%3.1f)) name(stage`stage',replace)) } grc1leg stage1 stage2 stage3 stage4, title("Female, Age 75, Deprivation Group 5, Diagnosed in 2011")