//
// Example from https://pclambert.net/software/standsurv/standardized_cif/
//
// Standardized cause-specific cumulative functions

// Load Data
use https://www.pclambert.net/data/mgus2, clear
drop if mspike == .

// How many events
tab event

// PCM model
stset survtime, failure(event=1)
stpm2 age mspike male, scale(h) df(4)
estimates store pcm

// Death Model
stset survtime, failure(event=2)
stpm2 age mspike male, scale(h) df(4) tvc(age) dftvc(2)
estimates store death

// times to predict at
range tt 0 30 31

// marginal CIFs
standsurv, crmodels(pcm death) cif ci timevar(tt) atvar(F) verbose

// plot results
twoway  (rarea F_pcm_lci F_pcm_uci tt, color(red%30)) ///
        (line F_pcm tt, color(red)) ///
        (rarea F_death_lci F_death_uci tt, color(blue%30)) ///
        (line F_death tt, color(blue)) ///
		, legend(order(2 "PCM" 4 "Death") cols(1) ring(0) pos(11)) ///
		ylabel(,angle(h) format(%3.2f)) ///
		xtitle("Time from diagnosis (years)") ytitle("cause-specific CIF") ///
		name(cifs, replace)

// Contrasts using standardization
standsurv, crmodels(pcm death) cif ci timevar(tt) verbose ///
    at1(male 1) at2(male 0) atvar(F_male F_female) ///
    contrast(difference) contrastvar(cif_diff)

// Plot PCM standardized curves
twoway  (rarea F_male_pcm_lci F_male_pcm_uci tt, color(red%30)) ///
        (line F_male_pcm tt, color(red)) ///
        (rarea F_female_pcm_lci F_female_pcm_uci tt, color(blue%30)) ///
        (line F_female_pcm tt, color(blue)) ///
		, legend(order(2 "Males" 4 "Females") cols(1) ring(0) pos(11)) ///
		ylabel(, angle(h) format(%3.2f)) ///
		xtitle("Time from diagnosis (years)") ytitle("cause-specific CIF") ///
		title("PCM") ///
		name(pcm, replace)

// Plot Death standardized curves
twoway  (rarea F_male_death_lci F_male_death_uci tt, color(red%30)) ///
        (line F_male_death tt, color(red)) ///
        (rarea F_female_death_lci F_female_death_uci tt, color(blue%30)) ///
        (line F_female_death tt, color(blue)) ///
		, legend(order(2 "Males" 4 "Females") cols(1) ring(0) pos(11)) ///
		ylabel(, angle(h) format(%3.2f)) ///
		xtitle("Time from diagnosis (years)") ytitle("cause-specific CIF") ///
		title("Death") ///
		name(death, replace)

graph combine pcm death, nocopies ycommon


// Plot Contrasts
twoway  (rarea cif_diff_pcm_lci cif_diff_pcm_uci tt, color(red%30)) ///
        (line cif_diff_pcm tt, color(red)) ///
		, legend(off) ///
		ylabel(, angle(h) format(%3.2f)) ///
		xtitle("Time from diagnosis (years)") ytitle("cause-specific CIF") ///
		title("PCM") ///
		name(pcm_diff, replace)

twoway  (rarea cif_diff_death_lci cif_diff_death_uci tt, color(red%30)) ///
        (line cif_diff_death tt, color(red)) ///
		, legend(off) ///
		ylabel(, angle(h) format(%3.2f)) ///
		xtitle("Time from diagnosis (years)") ytitle("cause-specific CIF") ///
		title("Death") ///
		name(death_diff, replace)

	graph combine pcm_diff death_diff, nocopies ycommon

	// Repeat but using different survival models

	// Weibull models with shape parameter a function of age
	stset survtime, failure(event=1)
	streg age mspike male, dist(weibull) //anc(age)
	estimates store pcm2

	// Accelerated failure time log logistic model
	stset survtime, failure(event=2)
	streg age mspike male, dist(llogistic)
	estimates store death2

	// Call standsurv with new models
	standsurv, crmodels(pcm2 death2) cif ci timevar(tt) verbose ///
		at1(male 1) at2(male 0) atvar(F2_male F2_female) ///
		contrast(difference) contrastvar(cif2_diff)

	// Plot results
	twoway  (rarea F2_male_pcm2_lci F2_male_pcm2_uci tt, color(red%30)) ///
			(line F2_male_pcm2 tt, color(red)) ///
			(rarea F2_female_pcm2_lci F2_female_pcm2_uci tt, color(blue%30)) ///
			(line F2_female_pcm2 tt, color(blue)) ///
			, legend(order(2 "Males" 4 "Females") cols(1) ring(0) pos(11)) ///
			ylabel(0(0.2)0.8, angle(h) format(%3.2f)) ///
			xtitle("Time from diagnosis (years)") ytitle("cause-specific CIF") ///
			title("PCM") ///
			name(pcm2, replace)

	twoway  (rarea F2_male_death2_lci F2_male_death2_uci tt, color(red%30)) ///
			(line F2_male_death2 tt, color(red)) ///
			(rarea F2_female_death2_lci F2_female_death2_uci tt, color(blue%30)) ///
			(line F2_female_death2 tt, color(blue)) ///
			, legend(order(2 "Males" 4 "Females") cols(1) ring(0) pos(11)) ///
			ylabel(0(0.2)0.8, angle(h) format(%3.2f)) ///
			xtitle("Time from diagnosis (years)") ytitle("cause-specific CIF") ///
			title("Death") ///
			name(death2, replace)

	twoway  (rarea cif2_diff_pcm2_lci cif2_diff_pcm2_uci tt, color(red%30)) ///
			(line cif2_diff_pcm2 tt, color(red)) ///
			, legend(off) ///
			ylabel(-0.15(0.05)0.05, angle(h) format(%3.2f)) ///
			xtitle("Time from diagnosis (years)") ytitle("cause-specific CIF") ///
			title("PCM") ///
			name(pcm2_diff, replace)

	twoway  (rarea cif2_diff_death2_lci cif2_diff_death2_uci tt, color(red%30)) ///
			(line cif2_diff_death2 tt, color(red)) ///
			, legend(off) ///
			ylabel(-0.15(0.05)0.05, angle(h) format(%3.2f)) ///
			xtitle("Time from diagnosis (years)") ytitle("cause-specific CIF") ///
			title("Death") ///
			name(death2_diff, replace)

	graph combine pcm2 death2 pcm2_diff death2_diff, nocopies