// A user-defined / custom hazard model
// Red Door Analytics
// https://reddooranalytics.se/resources/a-user-defined-custom-hazard-model/
//
// 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 all
set seed 1991
set obs 300                             // patients
gen patient = _n
gen age = runiformint(10, 69)           // age at baseline
gen female = runiform() < 0.5
gen b = rnormal(0, 1)                   // patient frailty, sd 1
expand 2                                // two catheter insertions each
sort patient

//@ simulate
avalon user time infect,                                ///
        loghazard(-8 :+ 0.5:*log({t}) :- 0.004:*{t})    /// log h0(t)
        covariates(age 0.02 female -1 b 1)              /// log HRs
        maxtime(500)                                    //  days
drop b

//@ list
list patient time infect age female in 1/6, noobs

//@ userhaz
mata:
real matrix userhaz(transmorphic gml, real colvector t)
{
        real matrix linpred
        real colvector gammas

        linpred = merlin_util_xzb(gml,t)
        gammas = merlin_util_ap(gml,1)\merlin_util_ap(gml,2)
        return(exp(linpred :+ merlin_fp(t,(0,1)) * gammas))
}
end

//@ fit-fixed
merlin (time    age female,                     ///
                family(user, hfunc(userhaz)     ///
                failure(infect) nap(2)))

//@ fit-frailty
merlin (time    age female                      ///
                M1[patient]@1,                  ///
                family(user, hfunc(userhaz)     ///
                failure(infect) nap(2)))

//@ fit-tde
merlin (time    age age#fp(time, powers(0))     ///
                female M1[patient]@1,           ///
                family(user, hfunc(userhaz)     ///
                failure(infect) nap(2))         ///
                timevar(time))