// Competing risks in merlin
// Red Door Analytics
// https://reddooranalytics.se/resources/competing-risks-in-merlin/
//
// Every number on that page comes from running this file, top to bottom.
// It simulates its own data, so it needs nothing but Stata and two merlin-family packages:
//     merlin   https://reddooranalytics.se/software/merlin/
//     avalon   https://reddooranalytics.se/software/avalon/
// The page says which versions produced it, and whether they are released yet.
// Lines beginning //@ mark the sections the page shows.

//@ setup
clear
set seed 1906
set obs 400
gen tumsize = runiform(1, 7)            // tumour diameter, cm
gen byte hypoxic = runiform() < 0.4     // hypoxic tumour
gen byte chemo = runiform() < 0.5       // chemoradiotherapy

//@ simulate
avalon msm dftime state event,                           /// new variables
        hazard1(dist(weibull) lambda(0.02) gamma(1.2)    /// cause 1: pelvic
            covariates(tumsize 0.3 hypoxic 0.5 chemo -0.5)) ///
        hazard2(dist(weibull) lambda(0.01) gamma(1.5)    /// cause 2: distant
            covariates(tumsize 0.1 hypoxic 0.8 chemo -0.2)) ///
        maxtime(5)                                       //  5-year follow-up

//@ indicators
rename dftime1 dftime
gen byte failtype = state1 - 1
label define failtype 0 "none" 1 "pelvic" 2 "distant"
label values failtype failtype
gen byte pelvic  = failtype==1
gen byte distant = failtype==2

//@ tabulate
tab failtype

//@ fit-pelvic
merlin (dftime tumsize hypoxic chemo, family(weibull, failure(pelvic)))
scalar ll_pelvic = e(ll)
matrix b_pelvic = e(b)

//@ fit-distant
merlin (dftime tumsize hypoxic chemo, family(weibull, failure(distant)))
scalar ll_distant = e(ll)
matrix b_distant = e(b)

//@ fit-joint
merlin (dftime tumsize hypoxic chemo, family(weibull, failure(pelvic)))  ///
       (dftime tumsize hypoxic chemo, family(weibull, failure(distant)))
scalar ll_joint = e(ll)
matrix b_joint = e(b)

//@ compare-ll
display "pelvic alone   = " %11.6f ll_pelvic
display "distant alone  = " %11.6f ll_distant
display "their sum      = " %11.6f ll_pelvic + ll_distant
display "both together  = " %11.6f ll_joint

//@ compare-b
matrix b_separate = b_pelvic, b_distant
mata: d = st_matrix("b_joint") :- st_matrix("b_separate")
mata: d'
mata: max(abs(d))