// Probabilistic sensitivity analysis and survival models
// Red Door Analytics
// https://reddooranalytics.se/resources/probabilistic-sensitivity-analysis-and-survival-models/
//
// 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.

//@ simulate
clear
set seed 686
set obs 686
gen hormon = runiform() < 0.35          // hormonal therapy
gen cens = runiform(1,7)                // follow-up, 1 to 7 years
avalon standard rectime censrec,        /// time to recurrence, event
        distribution(weibull)           /// Weibull baseline
        lambda(0.1) gamma(1.3)          /// scale and shape
        covariates(hormon -0.4)         /// log hazard ratio
        maxtime(cens)                   //  censored at end of follow-up

//@ stset
stset rectime, failure(censrec)

//@ fit
merlin (_t hormon, family(weibull, failure(_d))), nolog

//@ store
estimates store m1

//@ timevar
gen tvar = 5 in 1

//@ predict
predict s5_fit, survival at(hormon 1) timevar(tvar)

//@ list-predict
list tvar s5_fit in 1

//@ truth
// the true S(5) for a treated patient
display exp(-0.1 * 5^1.3 * exp(-0.4))

//@ copy-b
matrix b1 = e(b)

//@ replace-b
matrix b1[1,1] = -1
matrix list b1

//@ from-b
merlin (_t hormon, family(weibull, failure(_d))), from(b1) iter(0)

//@ list-eb
matrix list e(b)

//@ predict-altered
predict s5_alt, survival at(hormon 1) timevar(tvar)
list tvar s5_alt in 1

//@ psa
local Nsim = 200
if _N < `Nsim' set obs `Nsim'   // a row for each draw's prediction
//to store our predictions
gen s5 = .
//loop over Nsim draws
set seed 13490        // <- don't forget!
forvalues i=1/`Nsim' {
        local draw = rnormal(-1,0.1)
        mat b1[1,1] = `draw'
        qui merlin (_t hormon, family(weibull, failure(_d))),   ///
                from(b1) iter(0)
        //double check our hormon coefficient matches the random draw
        assert [_cmp_1_1_1][_cons]==`draw'
        //predict survival at 5 years
        qui predict s5_draw, survival at(hormon 1) timevar(tvar)
        //store my new estimate
        qui replace s5 = s5_draw[1] in `i'
        drop s5_draw
}

//@ psa-summary
summarize s5

//@ final-draw
display "hormon coefficient in memory: " [_cmp_1_1_1][_cons]
estimates restore m1
display "hormon coefficient in memory: " [_cmp_1_1_1][_cons]

//@ merlin-psa
merlin psa s5_psa, survival at(hormon 1) timevar(tvar)          ///
        reps(1000) seed(4217)                                   ///
        genmean(s5_psa_mean) gensd(s5_psa_sd)
list s5_psa s5_psa_lci s5_psa_uci s5_psa_mean s5_psa_sd in 1, abbrev(12)

//@ delta
predict s5_dm, survival at(hormon 1) timevar(tvar) ci
list s5_dm s5_dm_lci s5_dm_uci in 1, abbrev(10)

//@ micro-simulate
keep in 1                       // one treated patient...
replace hormon = 1
expand 100000                   // ...copied 100,000 times
avalon model stime event, model(m1) maxtime(5) seed(3418736)

//@ micro-list
list hormon stime event in 1/5

//@ micro-count
count if event==0
di "S(5) = " r(N)/_N