// Joint longitudinal and competing risks models: simulation, estimation and prediction
// Red Door Analytics
// https://reddooranalytics.se/resources/joint-longitudinal-and-competing-risks-models-simulation-estimation-and-prediction/
//
// Every number and figure 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 7254
set obs 500
gen id = _n

//@ random-effects
gen b0 = rnormal(0,1)
gen b1 = rnormal(0,0.2)
gen trt = runiform()>0.5

//@ simulate
//cause one
local maxt 5
local l1 0.1
local g1 1.2
local l2 0.05
local g2 1.5
avalon msm stime state event,                   /// new variables
        hazard1(                                /// cause 1 (cancer)
                user(`l1':*`g1':*{t}:^(`g1'-1)  /// user-defined function
                     :* exp(0.5                 /// biomarker
                     :* (b0 :+ (0.1:+b1):*{t})  ///   trajectory
                     ))                         ///
                covariates(trt -0.5))           /// treatment effect (PH)
        hazard2(                                /// cause 2 (other)
                user(`l2':*`g2':*{t}:^(`g2'-1)  /// user-defined function
                     :* exp( -0.3 :*            /// biomarker
                     (b0 :+ (0.1:+b1):*{t})     ///   trajectory
                     ))                         ///
                covariates(trt -0.2))           /// treatment effect (PH)
        maxtime(`maxt')                         //  admin. censoring

//@ list
list id trt stime* state* event* if _n<=5

//@ indicators
gen byte cancer = state1==2
gen byte other = state1==3
drop stime0 state0 state1 event1

//@ expand
expand 5
bys id : gen time = _n-1
drop if time>stime1

//@ biomarker
gen xb = b0 + 0.1 * time + b1 * time
gen y = rnormal(xb,0.5)                 //measurement error

//@ wide
bys id (time) : replace stime1 = . if _n>1
bys id (time) : replace cancer = . if _n>1
bys id (time) : replace other = . if _n>1

//@ fit
merlin (y                       /// biomarker outcome
          time                  /// fixed linear time
          time#M2[id]@1         /// random slope on linear time
          M1[id]@1,             /// random intercept
          family(gaussian)      /// distribution
          timevar(time))        /// variable representing time
       (stime1                  /// survival time
          trt                   /// baseline treatment
          EV[y],                /// expected value of biomarker
          family(weibull,       /// distribution
               failure(cancer)) /// cause-specific event indicator
          timevar(stime1))      /// timevar
       (stime1                  /// survival time
          trt                   /// baseline treatment
          EV[y],                /// expected value of biomarker
          family(weibull,       /// distribution
               failure(other))  /// cause-specific event indicator
          timevar(stime1))      //  timevar

//@ timevar
range tvar 0 5 500

//@ cif
predict cif1, cif marginal outcome(2) at(trt 0) timevar(tvar)
predict cif2, cif marginal outcome(2) at(trt 1) timevar(tvar)
predict cif3, cif marginal outcome(3) at(trt 0) timevar(tvar)
predict cif4, cif marginal outcome(3) at(trt 1) timevar(tvar)

//@ stack
gen totalcif1 = cif1 + cif3
gen totalcif2 = cif2 + cif4

//@ plot-stacked
preserve
keep tvar cif1-cif4 totalcif1 totalcif2
keep if tvar < .
expand 2
bys tvar : gen byte group = _n - 1           // 0 placebo, 1 treated
gen top    = cond(group, totalcif2, totalcif1)
gen bottom = cond(group, cif4, cif3)
label define group 0 "Placebo group" 1 "Treated group"
label values group group
twoway (area top tvar) (area bottom tvar),                          ///
        by(group, note("") legend(pos(6)))                           ///
        xtitle("Time since entry") ytitle("Cumulative incidence")    ///
        legend(cols(1) order(1 "Prob. of death due to cancer"        ///
                2 "Prob. of death due to other causes"))             ///
        xlabel(0(1)5) ylabel(0(0.1)1, angle(h) format(%2.1f))

//@ export-stacked
graph export "stacked-cif.svg", replace
restore

//@ cifdiff
predict diffcif, cifdifference marginal outcome(2)       ///
                 at1(trt 1) at2(trt 0) timevar(tvar)     ///
                 ci

//@ plot-diff
twoway  (rarea diffcif_lci diffcif_uci tvar)                    ///
        (line diffcif tvar)                                     ///
        , xtitle("Time since entry")                            ///
        ytitle("Difference in cumulative incidence")            ///
        title("CIF({it:t} | treated) - CIF({it:t} | placebo)")  ///
        ylabel(,angle(h) format(%3.2f))                         ///
        legend(pos(6) rows(1) order(2 "Diff. in CIF" 1 "95% confidence interval"))

//@ export-diff
graph export "cif-difference.svg", replace