//==================//
// EXERCISE 140
// REVISED June 2019
//==================//

use colon, clear
drop if stage ==0
gen female = sex==2

/* Cancer mortality */
stset surv_mm, failure(status==1) scale(12) exit(time 120.5)
tab status _d
sts generate surv1_sex = s, by(female) 

/* Other cause mortality */
stset surv_mm, failure(status==2) scale(12) exit(time 120.5)
tab status _d
sts generate surv2_sex = s, by(female) 

/* Plot the complement of the Kaplan-Meier survival estimate for each cause */
gen prob1_sex = 1-surv1_sex
gen prob2_sex = 1-surv2_sex

twoway (line prob1_sex _t if female == 0, sort lcolor(black) lpattern(solid)) ///
	   (line prob2_sex _t if female == 0, sort lcolor(black) lpattern(dash)), ///
	   ytitle("Probability of Death") ///
	   xtitle("Time Since Diagnosis (Years)") ///
	   legend(order(1 "Cancer" 2 "Other Causes")) ///
	   ylabel(0(0.1)0.5, angle(0) format(%3.1f)) name(KM, replace) 

/* Estimate the cumulative incidence function (CIF) */
stset surv_mm, failure(status==1) scale(12) exit(time 120.5)
stcompet CIF_sex=ci, compet1(2) by(sex)
gen CIF_sex_cancer=CIF_sex if status==1
gen CIF_sex_other=CIF_sex if status==2

/* Plot the cumulative incidence function along with the complement of the Kaplan-Meier estimates */
twoway (line prob1_sex _t if female==0, sort lcolor(black) lpattern(solid)) ///
	   (line prob2_sex _t if female==0, sort lcolor(black) lpattern(dash)) ///
	   (line CIF_sex_cancer _t if female==0, sort lcolor(gs10) lpattern(solid)) ///
	   (line CIF_sex_other _t if female==0, sort lcolor(gs10) lpattern(dash)), ///
	   ytitle("Probability of Death") ///
	   xtitle("Time Since Diagnosis (Years)") ///
	   legend(order(1 "Cancer K-M" 2 "Other K-M" 3 "Cancer CIF" 4 "Other CIF")) ///
	   ylabel(0(0.1)0.5, angle(0) format(%3.1f)) name(KMandCIF, replace)
	

/* Estimate CIF by age group */
stset surv_mm, failure(status==1) scale(12) exit(time 120.5)
stcompet CIF_age=ci, compet1(2) by(agegrp)
	 
twoway	(line CIF_age _t if agegrp == 0 & status == 1, sort connect(stepstair)) ///	 
		(line CIF_age _t if agegrp == 1 & status == 1, sort connect(stepstair)) ///	 
		(line CIF_age _t if agegrp == 2 & status == 1, sort connect(stepstair)) ///	 
		(line CIF_age _t if agegrp == 3 & status == 1, sort connect(stepstair)) ///	 
		, legend(order(1 "<45" 2 "45-59" 3 "60-74" 4 "75+") ring(0) pos(5) cols(1)) ///
		xtitle("Years since diagnosis") ///
		ytitle("CIF") ///
		title("Cancer") ///
		name(CIF_age1,replace) 

twoway	(line CIF_age _t if agegrp == 0 & status == 2, sort connect(stepstair)) ///	 
		(line CIF_age _t if agegrp == 1 & status == 2, sort connect(stepstair)) ///	 
		(line CIF_age _t if agegrp == 2 & status == 2, sort connect(stepstair)) ///	 
		(line CIF_age _t if agegrp == 3 & status == 2, sort connect(stepstair)) ///	 
		, legend(order(1 "<45" 2 "45-59" 3 "60-74" 4 "75+") ring(0) pos(11) cols(1)) ///
		xtitle("Years since diagnosis") ///
		ytitle("CIF") ///
		title("Other causes") ///
		name(CIF_age2,replace)
		
graph combine CIF_age1 CIF_age2, nocopies ycommon		
		
/* Estimate CIF by stage */
stcompet CIF_stage=ci, compet1(2) by(stage)
	 
twoway	(line CIF_stage _t if stage == 1 & status == 1, sort connect(stepstair)) ///	 
		(line CIF_stage _t if stage == 2 & status == 1, sort connect(stepstair)) ///	 
		(line CIF_stage _t if stage == 3 & status == 1, sort connect(stepstair)) ///	 
		, legend(order(1 "local" 2 "regional" 3 "distant") ring(0) pos(5) cols(1)) ///
		xtitle("Years since diagnosis") ///
		ytitle("CIF") ///
		title("Cancer") ///
		name(CIF_stage1,replace) 

twoway	(line CIF_stage _t if stage == 1 & status == 2, sort connect(stepstair)) ///	 
		(line CIF_stage _t if stage == 2 & status == 2, sort connect(stepstair)) ///	 
		(line CIF_stage _t if stage == 3 & status == 2, sort connect(stepstair)) ///	 
		, legend(order(1 "local" 2 "regional" 3 "distant") ring(0) pos(1) cols(1)) ///
		xtitle("Years since diagnosis") ///
		ytitle("CIF") ///
		title("Other causes") ///
		name(CIF_stage2,replace) 	

graph combine CIF_stage1 CIF_stage2, nocopies ycommon		

/* Now, estimate CIFs from a Cox model instead */
// Stset the data, specify that death from cancer is the outcome of interest and fit the cox model 
stset surv_mm, failure(status==1) scale(12) exit(time 120.5)
stcox female
// Obtain the cause-specific hazard functions from Cox model for cancer
* Predict baseline hazard
predict h0_cancer, basehc
* Sort time within descending order of death indicator variable 
**	and only keep the baseline hazard for one row in _t.
gsort _t -_d
by _t: replace h0_cancer = . if _n > 1
* Baseline CSH rate for males (female=0).
gen h_cancer_male = h0_cancer
* Baseline CSH rate for females (female=1).
gen h_cancer_female = h0_cancer*exp(_b[female])

// Now do the same with other causes as the outcome of interest
stset surv_mm, failure(status==2) scale(12) exit(time 120.5)
stcox female
* Predict baseline hazard
predict h0_other, basehc
* Sort time within descending order of death indicator variable 
**	and only keep the baseline hazard for one row in _t.
gsort _t -_d
by _t: replace h0_other = . if _n > 1
* Baseline CSH rate for males (female=0).
gen h_other_male = h0_other
* Baseline CSH rate for females (female=1).
gen h_other_female = h0_other*exp(_b[female])

* We want to keep the CSH rates at each time-point and cause
drop if missing(h0_cancer) & missing(h0_other)
foreach i in cancer other {
	replace h0_`i' = 0 if missing(h0_`i')
	replace h_`i'_male = 0 if missing(h_`i'_male)
	replace h_`i'_female = 0 if missing(h_`i'_female)
}
sort _t

// Calculate all-cause survival functions.
* Calculate all-cause survival function for those not on treatment
gen S_male = exp(sum(log(1- h_cancer_male - h_other_male)))
* Calculate all-cause survival function for those on treatment
gen S_female = exp(sum(log(1- h_cancer_female - h_other_female)))
// Calculate cause-specific CIFs
foreach i in cancer other {
	gen cif_male_`i' = sum(S_male[_n-1]*h_`i'_male)
	gen cif_female_`i' = sum(S_female[_n-1]*h_`i'_female)
}

/* Plot stacked CIF's for cancer and other causes for males and females */
gen male_tot1_cox = cif_male_cancer
gen male_tot2_cox = cif_male_cancer + cif_male_other

gen female_tot1_cox = cif_female_cancer
gen female_tot2_cox = cif_female_cancer + cif_female_other


twoway (area male_tot2_cox _t, sort fintensity(100)) ///
	   (area male_tot1_cox _t, sort fintensity(100)), ///
	    ylabel(0(0.2)1, angle(0) format(%3.1f)) ///
		ytitle("Probability of Death") xtitle("Time Since Diagnosis (Years)") ///
		legend(order(2 "Cancer" 1 "Other") size(small)) ///
        title("Males") plotregion(margin(zero)) name(malestack_cox, replace)
		
twoway (area female_tot2_cox _t, sort fintensity(100)) ///
	   (area female_tot1_cox _t, sort fintensity(100)), ///
	    ylabel(0(0.2)1, angle(0) format(%3.1f)) ///
		ytitle("Probability of Death") xtitle("Time Since Diagnosis (Years)") ///
		legend(order(2 "Cancer" 1 "Other") size(small)) ///
        title("Females") plotregion(margin(zero)) name(femalestack_cox, replace)
		
graph combine malestack_cox femalestack_cox, ycommon xcommon nocopies


/* Fit a competing risks model using the flexible parametric approach */

* Re-load data
use colon, clear
drop if stage ==0
gen female = sex==2

* Similar to how we did with the Cox model, we fit separate cause-specific FPMs
* To do so, we stset the data with the outcome of interest each time and this time, we will store the estimates.

// Create dummy variables for age group.
tab agegrp, gen(agegrp)

* Fit cause-specific FPM for cancer (k=1)
stset surv_mm, failure(status==1) scale(12) exit(time 120.5)
stpm2 female agegrp2 agegrp3 agegrp4, scale(hazard) df(4) eform
estimates store cancer

* Fit cause-specific FPM for other causes (k=2)
stset surv_mm, failure(status==2) scale(12) exit(time 120.5)
stpm2 female agegrp2 agegrp3 agegrp4, scale(hazard) df(4) eform
estimates store other

* Time to make predictions at.
range tempt 0 10 101

* Predict conditional estimates using standsurv.
* Don't forget if _n == 1 to predict for one individual and specify covariate pattern.

standsurv if _n==1, timevar(tempt) ci cif ///
	crmodels(cancer other) verbose ///
	atvar(CIF_male CIF_female) ///
	at1(female 0 agegrp2 0 agegrp3 0 agegrp4 1) at2(female 1 agegrp2 0 agegrp3 0 agegrp4 1)

/* Estimate the empirical cumulative incidence function (CIF) */
stset surv_mm, failure(status==1) scale(12) exit(time 120.5)
stcompet CIF_sex=ci if agegrp4 == 1, compet1(2) by(sex)
gen CIF_sex_cancer=CIF_sex if status==1
gen CIF_sex_other=CIF_sex if status==2
	
* Compare with empirical estimates
tw (line CIF_female_cancer CIF_female_other tempt, sort lpattern(solid dash) color(black black)) ///
	(line CIF_sex_cancer _t if female==1 & agegrp4==1, sort lcolor(gs10) lpattern(solid)) ///
	(line CIF_sex_other _t if female==1 & agegrp4==1, sort lcolor(gs10) lpattern(dash)), ///
		name(fpmcsh, replace) ytitle("Probability of Death") ///
	   xtitle("Time Since Diagnosis (Years)") ///
	   legend(order(1 "Cancer FPM" 2 "Other FPM" 3 "Cancer AJ" 4 "Other AJ")) ///
	   ylabel(0(0.1)0.5, angle(0) format(%3.1f)) 
	  
/* Repeat with time-dependent effects for age group and sex */
* Fit cause-specific FPM for cancer (k=1)
stset surv_mm, failure(status==1) scale(12) exit(time 120.5)
stpm2 female agegrp2 agegrp3 agegrp4, scale(hazard) df(4) tvc(sex agegrp2 agegrp3 agegrp4) dftvc(3) eform
estimates store cancer_tde
* Fit cause-specific FPM for other causes (k=2)
stset surv_mm, failure(status==2) scale(12) exit(time 120.5)
stpm2 female agegrp2 agegrp3 agegrp4, scale(hazard) df(4) tvc(sex agegrp2 agegrp3 agegrp4) dftvc(3) eform
estimates store other_tde
* Predict CIFs
standsurv if _n==1, timevar(tempt) ci cif ///
	crmodels(cancer_tde other_tde) verbose ///
	atvar(CIF_male CIF_female) ///
	at1(female 0 agegrp2 0 agegrp3 0 agegrp4 1) at2(female 1 agegrp2 0 agegrp3 0 agegrp4 1)

* Compare with empirical estimates
tw (line CIF_female_cancer_tde CIF_female_other_tde tempt, sort lpattern(solid dash) color(black black)) ///
	(line CIF_sex_cancer _t if female==1 & agegrp4==1, sort lcolor(gs10) lpattern(solid)) ///
	(line CIF_sex_other _t if female==1 & agegrp4==1, sort lcolor(gs10) lpattern(dash)), ///
	name(fpmcsh_tde, replace) ytitle("Probability of Death") ///
	   xtitle("Time Since Diagnosis (Years)") ///
	  legend(order(1 "Cancer FPM (tde)" 2 "Other FPM (tde)" 3 "Cancer AJ" 4 "Other AJ")) ///
	   ylabel(0(0.1)0.5, angle(0) format(%3.1f)) 
	
	
/* Fit a competing risks model using Fine and Gray's method */
stset surv_mm, failure(status==1) scale(12) exit(time 120.5)
stcrreg i.sex, compete(status == 2)	
stpepemori sex, compet(2)

/* Plot the associated cancer-specific cumulative incidence functions */
predict cif_males, basecif
gen cif_females = 1 - (1-cif_males)^exp(_b[2.sex])

graph twoway line cif_males cif_females _t, ///
	sort connect(step step) yscale(range(0 0.6)) ylabel(0(0.2)0.6) ///
	legend(order(1 "Males" 2 "Females")) ytitle(Cause-specific cumulative incidence) xtitle(Time since diagnosis (years)) ///
	title(First principles) name(cif_sex, replace) 

stcurve, cif at1(sex = 1) at2(sex=2) title(Stcurve) name(cif_sex2, replace)

graph combine cif_sex cif_sex2, ycommon 

/* Fit a competing risks model and plot the CIFs for deaths due to other causes than cancer */
stset surv_mm, failure(status==2) scale(12) exit(time 120.5)
stcrreg i.sex, compete(status == 1)
stcurve, cif at1(sex = 1) at2(sex=2) title(Other causes) legend(order(1 "Males" 2 "Females"))