//==================// // EXERCISE 251 // REVISED JUNE 2019 //==================// /* (a) Load melanoma data and merge in the expected mortality */ clear all use melanoma, clear stset surv_mm, failure(status=1,2) scale(12) id(id) exit(time 60.5) gen _age = min(int(age + _t),99) gen _year = int(yydx + _t) sort _year sex _age merge m:1 _year sex _age using popmort.dta, keep(match master) tab agegrp, gen(agegrp) stpm2 sex agegrp2-agegrp4, scale(hazard) bhazard(rate) df(5) /// tvc(sex agegrp2-agegrp4) dftvc(2) range temptime 0 5 1000 predict nm1, failure at(sex 1) zeros timevar(temptime) predict nm2, failure at(sex 1 agegrp2 1) zeros timevar(temptime) predict nm3, failure at(sex 1 agegrp3 1) zeros timevar(temptime) predict nm4, failure at(sex 1 agegrp4 1) zeros timevar(temptime) twoway (line nm1 nm2 nm3 nm4 temptime) /// , legend(order(1 "0-44" 2 "45-59" 3 "60-74" 4 "75+") ring(0) pos(11) cols(1)) /// ylabel(0(0.1)0.5, angle(h) format(%3.1f)) name("nm,replace") /// xtitle("Time since diagnosis") ytitle("Net probability of death") /* (b) estimate and plot crude mortality */ mata function calc_allcause(at) return(at[1]+at[2]) range temptime2 0 5 101 gen aged = . gen dated = mdy(1,1,1985) in 1 replace aged = 40 in 1 standsurv if _n==1, at1(sex 1 agegrp2 0 agegrp3 0 agegrp4 0) verbose timevar(temptime2) /// atvar(crprob1) crudeprob stub2(cancer other) /// expsurv(using(popmort.dta) /// datediag(dated) /// agediag(aged) /// pmrate(rate) /// pmage(_age) /// pmyear(_year) /// pmother(sex) /// pmmaxyear(1985) /// at1(sex 1)) /// userfunction(calc_allcause) /// userfunctionvar(allcause1) transform(none) replace aged = 55 in 1 standsurv if _n==1, at1(sex 1 agegrp2 1 agegrp3 0 agegrp4 0) verbose timevar(temptime2) /// atvar(crprob2) crudeprob stub2(cancer other) /// expsurv(using(popmort.dta) /// datediag(dated) /// agediag(aged) /// pmrate(rate) /// pmage(_age) /// pmyear(_year) /// pmother(sex) /// pmmaxyear(1985) /// at1(sex 1)) /// userfunction(calc_allcause) /// userfunctionvar(allcause2) transform(none) replace aged = 70 in 1 standsurv if _n==1, at1(sex 1 agegrp2 0 agegrp3 1 agegrp4 0) verbose timevar(temptime2) /// atvar(crprob3) crudeprob stub2(cancer other) /// expsurv(using(popmort.dta) /// datediag(dated) /// agediag(aged) /// pmrate(rate) /// pmage(_age) /// pmyear(_year) /// pmother(sex) /// pmmaxyear(1985) /// at1(sex 1)) /// userfunction(calc_allcause) /// userfunctionvar(allcause3) transform(none) replace aged = 80 in 1 standsurv if _n==1, at1(sex 1 agegrp2 0 agegrp3 0 agegrp4 1) verbose timevar(temptime2) /// atvar(crprob4) crudeprob stub2(cancer other) /// expsurv(using(popmort.dta) /// datediag(dated) /// agediag(aged) /// pmrate(rate) /// pmage(_age) /// pmyear(_year) /// pmother(sex) /// pmmaxyear(1985) /// at1(sex 1)) /// userfunction(calc_allcause) /// userfunctionvar(allcause4) transform(none) twoway (line crprob?_cancer temptime2) /// , legend(order(1 "40" 2 "55" 3 "70" 4 "80") ring(0) pos(11) cols(1)) /// ylabel(0(0.1)0.5, angle(h) format(%3.1f)) name("cm,replace") /// xtitle("Time since diagnosis") ytitle("Probability of death") title("Cancer") /* (c) crude prob of death due to other causes */ twoway (line crprob?_other temptime2) /// , legend(order(1 "40" 2 "55" 3 "70" 4 "80") ring(0) pos(11) cols(1)) /// ylabel(0(0.1)0.5, angle(h) format(%3.1f)) name("cm,replace") /// xtitle("Time since diagnosis") ytitle("Probability of death") title("Other causes") /* (d) stacked graphs */ local title1 "40" local title2 "55" local title3 "70" local title4 "80" forvalues i = 1/4 { twoway (area crprob`i'_cancer temptime2) /// (rarea allcause`i' crprob`i'_cancer temptime2) /// (area allcause`i' temptime2, base(1)) /// , ylabel(0(0.2)1.0, angle(h) format(%3.1f)) /// xtitle("Time since diagnosis") ytitle("crude probability of death") /// legend(order(1 "P(Dead Cancer)" 2 "P(Dead Other Causes)" 3 "P(Alive)") /// cols(3)) plotregion(margin(zero)) title(`title`i'') /// name(cm_stack`i',replace) } grc1leg cm_stack1 cm_stack2 cm_stack3 cm_stack4, nocopies cols(4) name(crpob_agegrp, replace) /* (e) Compare with net prob */ local title1 "40" local title2 "55" local title3 "70" local title4 "80" forvalues i = 1/4 { twoway (area crprob`i'_cancer temptime2) /// (rarea allcause`i' crprob`i'_cancer temptime2) /// (area allcause`i' temptime2, base(1)) /// (line nm`i' temptime, sort lcolor(black)) /// , ylabel(0(0.2)1.0, angle(h) format(%3.1f)) /// xtitle("Time since diagnosis") ytitle("crude probability of death") /// legend(order(1 "Cancer" 2 "Other Causes" 3 "Alive" 4 "Net") /// cols(4)) plotregion(margin(zero)) title(`title`i'') /// name(cm_stack`i',replace) } grc1leg cm_stack1 cm_stack2 cm_stack3 cm_stack4, nocopies cols(4) name(crpob_age, replace) /* (f) Advanced: splines for age */ rcsgen age, gen(rcsage) df(4) orthog global knots `r(knots)' matrix Rage = r(R) stpm2 sex rcsage1-rcsage4, scale(hazard) df(5) bhazard(rate) /// tvc(sex rcsage1-rcsage4) dftvc(2) range temptime2 0 5 101 gen aged = . gen dated = mdy(1,1,1985) in 1 replace aged = 40 in 1 rcsgen , scalar(40) knots($knots) rmatrix(Rage) gen(c) standsurv if _n==1, at1(sex 1 rcsage1 `=c1' rcsage2 `=c2' rcsage3 `=c3' rcsage4 `=c4') verbose timevar(temptime2) /// atvar(crprob_age40) crudeprob stub2(cancer other) /// expsurv(using(popmort.dta) /// datediag(dated) /// agediag(aged) /// pmrate(rate) /// pmage(_age) /// pmyear(_year) /// pmother(sex) /// pmmaxyear(1985) /// at1(sex 1)) /// userfunction(calc_allcause) /// userfunctionvar(allcause_age40) transform(none) replace aged = 55 in 1 rcsgen , scalar(55) knots($knots) rmatrix(Rage) gen(c) standsurv if _n==1, at1(sex 1 rcsage1 `=c1' rcsage2 `=c2' rcsage3 `=c3' rcsage4 `=c4') verbose timevar(temptime2) /// atvar(crprob_age55) crudeprob stub2(cancer other) /// expsurv(using(popmort.dta) /// datediag(dated) /// agediag(aged) /// pmrate(rate) /// pmage(_age) /// pmyear(_year) /// pmother(sex) /// pmmaxyear(1985) /// at1(sex 1)) /// userfunction(calc_allcause) /// userfunctionvar(allcause_age55) transform(none) replace aged = 70 in 1 rcsgen , scalar(70) knots($knots) rmatrix(Rage) gen(c) standsurv if _n==1, at1(sex 1 rcsage1 `=c1' rcsage2 `=c2' rcsage3 `=c3' rcsage4 `=c4') verbose timevar(temptime2) /// atvar(crprob_age70) crudeprob stub2(cancer other) /// expsurv(using(popmort.dta) /// datediag(dated) /// agediag(aged) /// pmrate(rate) /// pmage(_age) /// pmyear(_year) /// pmother(sex) /// pmmaxyear(1985) /// at1(sex 1)) /// userfunction(calc_allcause) /// userfunctionvar(allcause_age70) transform(none) replace aged = 80 in 1 rcsgen , scalar(80) knots($knots) rmatrix(Rage) gen(c) standsurv if _n==1, at1(sex 1 rcsage1 `=c1' rcsage2 `=c2' rcsage3 `=c3' rcsage4 `=c4') verbose timevar(temptime2) /// atvar(crprob_age80) crudeprob stub2(cancer other) /// expsurv(using(popmort.dta) /// datediag(dated) /// agediag(aged) /// pmrate(rate) /// pmage(_age) /// pmyear(_year) /// pmother(sex) /// pmmaxyear(1985) /// at1(sex 1)) /// userfunction(calc_allcause) /// userfunctionvar(allcause_age80) transform(none) twoway (line crprob?_cancer temptime2) (line crprob_age??_cancer temptime2, lpattern(dash..)) /// , legend(order(1 "40" 2 "55" 3 "70" 4 "80") ring(0) pos(11) cols(1)) /// ylabel(0(0.1)0.5, angle(h) format(%3.1f)) name("cm_d,replace") /// xtitle("Time since diagnosis") ytitle("Probability of death") title("Cancer") twoway (line crprob?_other temptime2) (line crprob_age??_other temptime2, lpattern(dash..)) /// , legend(order(1 "40" 2 "55" 3 "70" 4 "80") ring(0) pos(11) cols(1)) /// ylabel(0(0.1)0.5, angle(h) format(%3.1f)) name("cm_o,replace") /// xtitle("Time since diagnosis") ytitle("Probability of death") title("Other causes") // Calculate CIF differences replace aged = 80 in 1 rcsgen , scalar(80) knots($knots) rmatrix(Rage) gen(c) standsurv if _n==1, at1(sex 1 rcsage1 `=c1' rcsage2 `=c2' rcsage3 `=c3' rcsage4 `=c4') /// at2(sex 2 rcsage1 `=c1' rcsage2 `=c2' rcsage3 `=c3' rcsage4 `=c4') verbose timevar(temptime2) /// atvar(crprob_age80_male crprob_age80_female) crudeprob stub2(cancer other) ci /// expsurv(using(popmort.dta) /// datediag(dated) /// agediag(aged) /// pmrate(rate) /// pmage(_age) /// pmyear(_year) /// pmother(sex) /// pmmaxyear(1985) /// at1(sex 1)) /// contrast(difference) /// userfunctionvar(crprob_sexdiff) transform(none) twoway (line _contrast2_cancer_lci _contrast2_cancer_uci temptime2, lpattern(dash..)) (line _contrast2_cancer temptime2, lpattern(solid)) /// , ylabel(, angle(h) format(%3.2f)) name("cm_d,replace") /// xtitle("Time since diagnosis") ytitle("CIF Difference") title("Cancer") /* (g) Advanced: plot over age */ gen tyr = 5 in 1 gen ageplot = . foreach i in cancer other { gen crprob5yr_`i' = . gen crprob5yr_`i'_lci = . gen crprob5yr_`i'_uci = . } gen allcause5yr = . gen allcause5yr_lci = . gen allcause5yr_uci = . local j = 1 forvalues a = 40/90 { replace aged = `a' in 1 replace ageplot = `a' in `j' rcsgen , scalar(`a') knots($knots) rmatrix(Rage) gen(c) standsurv if _n==1, at1(sex 1 rcsage1 `=c1' rcsage2 `=c2' rcsage3 `=c3' rcsage4 `=c4') verbose timevar(tyr) /// atvar(crprob_del) crudeprob stub2(cancer other) ci /// expsurv(using(popmort.dta) /// datediag(dated) /// agediag(aged) /// pmrate(rate) /// pmage(_age) /// pmyear(_year) /// pmother(sex) /// pmmaxyear(1985) /// at1(sex 1)) /// userfunction(calc_allcause) /// userfunctionvar(allcause_del) transform(none) foreach c in cancer other { replace crprob5yr_`c' = crprob_del_`c'[1] in `j' replace crprob5yr_`c'_lci = crprob_del_`c'_lci[1] in `j' replace crprob5yr_`c'_uci = crprob_del_`c'_uci[1] in `j' replace allcause5yr = allcause_del[1] in `j' replace allcause5yr_lci = allcause_del_lci[1] in `j' replace allcause5yr_uci = allcause_del_uci[1] in `j' } capture drop crprob_del* allcause_del* local j = `j' + 1 } tw (line crprob5yr_cancer crprob5yr_cancer_lci crprob5yr_cancer_uci ageplot, lpattern(solid dash dash) lcolor(black..)), /// legend(order(2 "95% CI") ring(0) pos(11) cols(1)) /// ylabel(0(0.2)1, angle(h) format(%3.1f)) name("cm_d,replace") /// xtitle("Age at diagnosis") ytitle("Probability of death") title("Cancer") tw (line crprob5yr_other crprob5yr_other_lci crprob5yr_other_uci ageplot, lpattern(solid dash dash) lcolor(black..)), /// legend(order(2 "95% CI") ring(0) pos(11) cols(1)) /// ylabel(0(0.2)1, angle(h) format(%3.1f)) name("cm_o,replace") /// xtitle("Age at diagnosis") ytitle("Probability of death") title("Other") tw (line allcause5yr allcause5yr_lci allcause5yr_uci ageplot, lpattern(solid dash dash) lcolor(black..)), /// legend(order(2 "95% CI") ring(0) pos(11) cols(1)) /// ylabel(0(0.2)1, angle(h) format(%3.1f)) name("cm_all,replace") /// xtitle("Age at diagnosis") ytitle("Probability of death") title("All-causes") grc1leg cm_d cm_o cm_all, nocopies ycommon cols(3) title("At 5 years since diagnosis")