program define mrsprep_example
  version 16.0
  syntax [, EGNUMBER(integer 1) LEAVEDATA]

  if "`leavedata'" != "" {
    qui desc
    if `r(N)'!=0 | `r(k)'!=0 {
      di as error `"To leave the example data in memory,"' _newline ///
                   "You need to start with an empty dataset."
      exit
    }
  }

if `egnumber' == 1 {

if "`leavedata'" == "" {
  qui frame
  local currentframe `r(currentframe)'
  preserve
}

display "\\ Load and stset melanoma data"
display `". use "https://pclambert.net/data/melanoma.dta""'
display ". stset exit, origin(dx) failure(status=1,2) id(id) exit(time dx + 10*365.24) scale(365.24)"
use "https://pclambert.net/data/melanoma.dta"
stset exit, origin(dx) failure(status=1,2) id(id) exit(time dx + 10*365.24) scale(365.24)

display "// Use mrsprep to expand data and calculate weights"
display ///
`". mrsprep using "https://pclambert.net/data/popmort.dta", pmother(sex) agediag(age) datediag(dx) ///"' _newline ///
"                                                           verbose breaks(0(0.2)10) " _newline

mrsprep using "https://pclambert.net/data/popmort.dta", pmother(sex) agediag(age) datediag(dx) ///
                           verbose breaks(0(0.2)10) newframe(, replace)
display _newline
display "// incorprate weights when using stset"
display ". stset tstop [iw=wt], enter(tstart) failure(event==1)"
stset tstop [iw=wt], enter(tstart) failure(event==1)

display "// Fit marginal model using the weighted mean hazard" _newline
display ". stpm2, scale(hazard) df (5) bhazard(meanhazard_wt) vce(cluster id)"
stpm2, scale(hazard) df (5) bhazard(meanhazard_wt) vce(cluster id)

display "// predict marginal relative survival"
display ". range tt 0 10 101"
display ". predict s_mrs, surv timevar(tt) ci"
range tt 0 10 101
predict s_mrs, surv timevar(tt) ci

di ///
". twoway (line s_mrs* tt, lcolor(red..) lpattern(solid dash dash)) ///" _newline ///
"         , legend(off)                                             ///" _newline ///
"           ylabel(0.6(0.1)1, format(%3.1f))                        ///" _newline ///
"           ytitle(Marginal relative survival)                      ///" _newline ///
"   	    xtitle(Years from diagnosis)                               " _newline
twoway (line s_mrs* tt, lcolor(red..) lpattern(solid dash dash)) ///
       , legend(off)                                             ///
	 ylabel(0.6(0.1)1, format(%3.1f))                        ///
         ytitle("Marginal relative survival")                    ///
	 xtitle("Years from diagnosis")

if "`leavedata'" == "" {
  frame change `currentframe'
  frame drop mrs_data
  restore
}
}




  else if `egnumber' == 2 {

if "`leavedata'" == "" {
  qui frame
  local currentframe `r(currentframe)'
  preserve
}


qui frame change default
display "\\ Load and stset melanoma data"
display `". use "https://pclambert.net/data/melanoma.dta""'
display ". stset exit, origin(dx) failure(status=1,2) id(id) exit(time dx + 10*365.24) scale(365.24)"
use "https://pclambert.net/data/melanoma.dta"
stset exit, origin(dx) failure(status=1,2) id(id) exit(time dx + 10*365.24) scale(365.24)

di "// change age groups to those defined in ICSS" _newline
di ". drop agegrp" _newline
di ". egen agegrp=cut(age), at(0 45 55 65 75 200) icodes" _newline
di ". replace agegrp = agegrp + 1" _newline
di `". label variable agegrp "Age group""' _newline
di `". label define agegrplab 1 "0-44" 2 "45-54" 3 "55-64" 4 "65-74" 5 "75+", replace"' _newline
di ". label values agegrp agegrplab" _newline

drop agegrp
egen agegrp=cut(age), at(0 45 55 65 75 200) icodes
replace agegrp = agegrp + 1
label variable agegrp "Age group"
label define agegrplab 1 "0-44" 2 "45-54" 3 "55-64" 4 "65-74" 5 "75+", replace
label values agegrp agegrplab

di "// create weights in reference population (ICSS)" _newline
di "recode agegrp (1=0.28) (2=0.17) (3=0.21) (4=0.20) (5=0.14), gen(ICSSwt)" _newline
recode agegrp (1=0.28) (2=0.17) (3=0.21) (4=0.20) (5=0.14), gen(ICSSwt)

di "// calculate relative weights" _newline
di ". local total= _N" _newline
di ". bysort agegrp: gen a_age = _N/`total'" _newline
di ". gen double wt_age = ICSSwt/a_age" _newline
local total= _N
bysort agegrp: gen a_age = _N/`total'
gen double wt_age = ICSSwt/a_age

di  "// Prepare data for marginal model" _newline
di ///
`". mrsprep using "https://pclambert.net/data/popmort.dta", pmother(sex) agediag(age) datediag(dx) ///"' _newline ///
"                           verbose breaks(0(0.2)10)               ///" _newline ///
"                           indweights(wt_age)                     ///" _newline ///
"                           newframe(mrs_stand, replace)" _newline
mrsprep using "https://pclambert.net/data/popmort.dta", pmother(sex) agediag(age) datediag(dx) ///
                                verbose breaks(0(0.2)10)               ///
                           indweights(wt_age)                     ///
                           newframe(mrs_stand, replace)

di  "// incorprate weights when using stset"
di ". stset tstop [iw=wt], enter(tstart) failure(event==1)" _newline
stset tstop [iw=wt], enter(tstart) failure(event==1)

di "// Fit marginal model using the weighted mean hazard" _newline
di ". stpm2, scale(hazard) df (5) bhazard(meanhazard_wt) vce(cluster id)" _newline
stpm2, scale(hazard) df (5) bhazard(meanhazard_wt) vce(cluster id)

di "// predict externally age standardized marginal relative survival" _newline
di ". range tt 0 10 101" _newline
di ". predict s_mrs, surv timevar(tt) ci" _newline
range tt 0 10 101
predict s_mrs, surv timevar(tt) ci

di "// Plot results"
di ///
". twoway (line s_mrs* tt, lcolor(red..) lpattern(solid dash dash)) ///" _newline ///
"       , legend(off)                                             ///" _newline ///
"	  ylabel(0.6(0.1)1, format(%3.1f))                        ///" _newline ///
"         ytitle(Marginal relative survival)                      ///" _newline ///
"	  xtitle(Years from diagnosis)" _newline

if "`leavedata'" == "" {
  frame change `currentframe'
  frame drop mrs_data
  restore
}

  }
  else if `egnumber' == 3 {

if "`leavedata'" == "" {
  qui frame
  local currentframe `r(currentframe)'
  preserve
}

qui frame change default
display "\\ Load and stset melanoma data"
display `". use "https://pclambert.net/data/melanoma.dta""'
display ". stset exit, origin(dx) failure(status=1,2) id(id) exit(time dx + 10*365.24) scale(365.24)"
use "https://pclambert.net/data/melanoma.dta"
stset exit, origin(dx) failure(status=1,2) id(id) exit(time dx + 10*365.24) scale(365.24)

display "// Use mrsprep to expand data and calculate weights"
display ///
`". mrsprep using "https://pclambert.net/data/popmort.dta", pmother(sex) agediag(age) datediag(dx) ///"' _newline ///
"                                                           verbose breaks(0(0.2)10) " _newline

di "// change age groups to those defined in ICSS" _newline
di ". drop agegrp" _newline
di ". egen agegrp=cut(age), at(0 45 55 65 75 200) icodes" _newline
di ". replace agegrp = agegrp + 1" _newline
di `". label variable agegrp "Age group""' _newline
di `". label define agegrplab 1 "0-44" 2 "45-54" 3 "55-64" 4 "65-74" 5 "75+", replace"' _newline
di ". label values agegrp agegrplab" _newline

drop agegrp
egen agegrp=cut(age), at(0 45 55 65 75 200) icodes
replace agegrp = agegrp + 1
label variable agegrp "Age group"
label define agegrplab 1 "0-44" 2 "45-54" 3 "55-64" 4 "65-74" 5 "75+", replace
label values agegrp agegrplab

di "// create weights in reference population (ICSS)" _newline
di "recode agegrp (1=0.28) (2=0.17) (3=0.21) (4=0.20) (5=0.14), gen(ICSSwt)" _newline
recode agegrp (1=0.28) (2=0.17) (3=0.21) (4=0.20) (5=0.14), gen(ICSSwt)


di "// Proportion within each age group by sex" _newline
di ". gen female = sex == 2" _newline
di ". bysort female: egen totalsex = total(sex)" _newline
di ". bysort agegrp female: gen a_age_sex = _N/totalsex" _newline
di ". gen double wt_age_sex = ICSSwt/a_age_sex" _newline
gen female = sex == 2
bysort female: egen totalsex = total(sex)
bysort agegrp female: gen a_age_sex = _N/totalsex
gen double wt_age_sex = ICSSwt/a_age_sex

di "// Prepare data for marginal model" _newline
di `". mrsprep using "https://pclambert.net/data/popmort.dta", pmother(sex) agediag(age) datediag(dx) ///"' ///
"                           pmmaxyear(2000)                        ///" _newline ///
"                           verbose breaks(0(0.2)10)               ///" _newline ///
"                           indweights(wt_age_sex)                 ///" _newline ///
"                           by(female)                             ///" _newline ///
"                           newframe(mrs_stand, replace)  " _newline
mrsprep using "https://pclambert.net/data/popmort.dta", pmother(sex) agediag(age) datediag(dx) ///
                           pmmaxyear(2000)                        ///
                           verbose breaks(0(0.2)10)               ///
                           indweights(wt_age_sex)                 ///
                           by(female)                             ///
                           newframe(mrs_stand, replace)

di "// incorprate weights when using stset" _newline
di ". stset tstop [iw=wt], enter(tstart) failure(event==1)" _newline
stset tstop [iw=wt], enter(tstart) failure(event==1)

di "// Fit proportional hazards marginal model" _newline
di ". stpm2 female, scale(hazard) df (5) bhazard(meanhazard_wt) vce(cluster id)" _newline
stpm2 female, scale(hazard) df (5) bhazard(meanhazard_wt) vce(cluster id)

di "// Relax proportional hazards assumption" _newline
di ". stpm2 female, scale(hazard) df (5) bhazard(meanhazard_wt) vce(cluster id) ///" _newline ///
"              tvc(female) dftvc(3)" _newline

 // predict externally age standardized marginal relative survival by sex
di "range tt 0 10 101" _newline
di "predict s_mrs_male,   surv timevar(tt) ci at(female 0)" _newline
di "predict s_mrs_female, surv timevar(tt) ci at(female 1)" _newline
range tt 0 10 101
predict s_mrs_male,   surv timevar(tt) ci at(female 0)
predict s_mrs_female, surv timevar(tt) ci at(female 1)

di "//Plot results" _newline
di ". twoway (line s_mrs_male* tt, lcolor(red..) lpattern(solid dash dash)) ///" _newline ///
"       (line s_mrs_female* tt, lcolor(blue..) lpattern(solid dash dash))   ///" _newline ///
"       , legend(off)                                                       ///" _newline ///
"	 ylabel(0.6(0.1)1, format(%3.1f))                                   ///" _newline ///
"         ytitle(Marginal relative survival)                                ///" _newline ///
"	 xtitle(Years from diagnosis)" _newline

twoway (line s_mrs_male* tt, lcolor(red..) lpattern(solid dash dash)) ///
       (line s_mrs_female* tt, lcolor(blue..) lpattern(solid dash dash)) ///
       , legend(off)                                             ///
	 ylabel(0.6(0.1)1, format(%3.1f))                        ///
         ytitle(Marginal relative survival)                    ///
	 xtitle(Years from diagnosis)

if "`leavedata'" == "" {
  frame change `currentframe'
  frame drop mrs_data
  restore
}

  }

end