//==================//
// Exercise 244
// Created May 2019
//==================//

use melanoma if stage==1, clear
drop if year8594==1
drop if age < 15

// merge in the external rates
stset surv_mm, failure(status=1,2) scale(12) id(id) exit(time 120.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,  keep(match master)

stset surv_mm, fail(status==1,2) scale(12)

// New age groups according to the International Cancer Survival Standard (ICSS)
label drop agegrp
drop agegrp
egen agegrp=cut(age), at(0 15 45 55 65 75 200) icodes
label variable agegrp "Age group"
label define agegrp 1 "15-44" 2 "45-54" 3 "55-64" 4 "65-74" 5 "75+" 
label values agegrp agegrp

// Generate dummy variables for agegrp. -stpm2- supports 
// factor variables (i.) for main effects but not time-varying effects
// Compare age distribution to the standard population
tab agegrp, gen(agegrp)

// Create weights to use for age-standardisation
// We use ICSS 2 weights for melanoma 
// 15-44 years: 28%, 
// 45-54 years: 17%, 
// 55-64 years: 21%, 
// 65-74 years: 20%,
// 75+   years: 14%
recode agegrp (1=0.28) (2=0.17) (3=0.21) (4=0.20) (5=0.14), gen(ICSSwt)
local total= _N
bysort agegrp: generate a_age = _N/`total'
generate w = ICSSwt/a_age

// look at the weights
tabstat w, statistics( mean min max) by(agegrp)

stpm2 agegrp2 agegrp3 agegrp4 agegrp5, scale(h) df(4) eform ///
  tvc(agegrp2 agegrp3 agegrp4 agegrp5) dftvc(1) bhazard(rate)

// Create values of follow-up time for which to predict survival
// See: https://pclambert.net/software/stpm2/stpm2_timevar/
range temptime 0 10 101

// marginal (population-averaged) unstandardised survival
predict s_unweighted, meansurv timevar(temptime)

// marginal (population-averaged) survival standardised to ICSS
predict s_weighted, meansurv timevar(temptime) meansurvwt(w) 

list temptime s_unweighted s_weighted if inlist(temptime, 1, 5, 10)


// Can also use standsurv
stpm2_standsurv , timevar(temptime) ///
		at1(.) at2(.,atindweights(w)) ///
		atvar(s_unweighted2 s_weighted2)


twoway 	(line s_weighted temptime , sort lpattern(shortdash) lwidth(medthick) lcolor(black)) ///
		(line s_unweighted temptime , sort lpattern(solid) lwidth(medthick) lcolor(black)) ///
		, legend(label(1 "Marginal net survival (ICSS standard)") label(2 "Marginal net survival (internal standard)") ring(0) pos(7) col(1)) ///
		scheme(sj) name(surv1, replace) ysize(8) xsize(11) ///
		subtitle("`text`i''", size(*1.0)) ytitle("Marginal net survival", size(*1.0)) xtitle("Years since diagnosis", size(*1.0)) ///
		ylabel(0.7 0.8 0.9 1.0, labsize(*1.0) angle(0)) yscale(range(0.7 1)) xlabel(, labsize(*1.0))