A user-defined / custom hazard model

This tutorial will illustrate some of the more advanced capabilities of merlin when modelling survival data, but with the aim of using an accessible example. During my PhD, Paul Lambert and I developed stgenreg in Stata for modelling survival data with a general user-specified hazard function, with the generality achieved by using numerical integration to calculate the cumulative hazard function, i.e. the difficult bit in the likelihood (Crowther and Lambert 2013, 2014).

Defining a new model

Consider our standard proportional hazards model,

h (t) = h_{0} (t) \exp (X β)

Within a general hazard model, we can essentially specify any function for our baseline hazard function, subject to the constraint that $h_{0} (t) > 0$ for all $t > 0$ . The easiest way to do this is to model on the log hazard scale. Let’s model our baseline log hazard function with fractional polynomials, such as,

\log h_{0} (t) = γ_{0} + γ_{1} \log (t) + γ_{2} t

This model can be fitted using stgenreg, but with the introduction of merlin, we can do the same as stgenreg, and a whole lot more.

Coding it up

We’ll simulate some data to fit some models to, shaped like Stata’s catheter dataset (McGilchrist and Aisbett, 1991): each patient has a catheter inserted twice, and we record the time to infection at the insertion point each time. We simulate 300 patients, each with an age at baseline between 10 and 69 and an indicator for being female. Infections cluster within patients, so each patient also gets a frailty, a normally distributed random intercept with standard deviation 1, shared by both of their infection times.

Stata

. clear all

. set seed 1991

. set obs 300                             // patients
Number of observations (_N) was 0, now 300.

. 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
(300 observations created)

. sort patient

We simulate the infection times, in days, from the model we’ll build up to, using avalon user, which takes a user-defined log hazard function written in Mata, with {t} standing for time (hence the colon operators). The baseline is $\log h_{0} (t) = - 8 + 0.5 \log (t) - 0.004 t$ , so $γ_{0} = - 8$ , $γ_{1} = 0.5$ and $γ_{2} = - 0.004$ , a hazard that rises to a peak at 125 days and falls after that. The log hazard ratios are 0.02 for each year of age and −1 for being female, and the frailty enters with a coefficient of 1. The effect of age does not change over time. Follow-up stops at 500 days, and catheters still free of infection then are censored there.

Stata

. 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

Our dataset consists of the following,

Stata

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

  +---------------------------------------------+
  | patient        time   infect   age   female |
  |---------------------------------------------|
  |       1   15.448549        1    44        0 |
  |       1   78.984023        1    44        0 |
  |       2   20.808838        1    38        1 |
  |       2   325.77116        1    38        1 |
  |       3   3.8442837        1    21        0 |
  |---------------------------------------------|
  |       3   8.6844023        1    21        0 |
  +---------------------------------------------+

with patient our individual patient identifier, time our time, in days, to infection at the catheter insertion point, infect our event indicator, with an event being an infection, age patient age at baseline and female a binary indicator variable. We immediately see that patients can experience multiple infections, and so we have events nested within patients. For now, we will ignore this clustering.

To fit our model with fractional polynomials for our baseline log hazard function, we need to write a little Mata function which calculates and returns our hazard function. merlin has the capabilities to allow you to define your own likelihood function, or hazard and/or cumulative hazard function, instead of the inbuilt distributions. This is achieved by providing a suite of utility functions. See help merlin user for all the details and documentation. We’ll get straight to our function:

Stata

. mata:
------------------------------------------------- mata (type end to exit) ------
: 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
--------------------------------------------------------------------------------

Let’s go through it line by line. First thing to note is that we declare a chunk of Mata code between mata: and end.

We need to define a function, called whatever we like, in this case I’ll call it userhaz(), which returns a real matrix, and it’s going to have two inputs. The first is a transmorphic object called gml. This is the internal struct which contains all the information needed by merlin in the background. You shouldn’t attempt to alter its contents. The second argument is a real colvector which I’m calling t. This represents the vector of time points that we wish to calculate our hazard function at. Internally, our function will be called by merlin at both our core time variable, and our quadrature points needed to calculate the cumulative hazard, and therefore this second input is needed.

Next we declare some intermediate vectors/matrices that we’ll need, real matrix linpred and real colvector gammas. Explicit declaration of each object’s type is good programming practice and reduces the likelihood of errors.

Now we call our first utility function merlin_util_xzb(), passing it the gml structure and also our time vector, in linpred = merlin_util_xzb(gml,t). This returns our main complex linear predictor; it’s as simple as that.

Passing t is not optional once the model has anything time-dependent in it: the linear predictor has to be evaluated at the same quadrature points as the hazard. Because we pass it, any time-dependent effects that we specify in our linear predictor, or calls to the elements EV[], dEV[] or iEV[] etc. (see the joint model examples with merlin), are automatically taken care of!

We then have two other ancillary parameters to handle, i.e. the coefficients of the fractional polynomial terms, which we extract using merlin_util_ap(gml,i) where i is the ancillary parameter number. In this case we have two extra parameters to estimate, so we build a column vector called gammas with gammas = merlin_util_ap(gml,1) \ merlin_util_ap(gml,2).

Finally we need to return our hazard function, which is done very simply, with return(exp(linpred :+ merlin_fp(t,(0,1)) * gammas)).

This makes use of the internal merlin_fp() function, which returns fractional polynomials, in this case an FP2 function with powers 0 and 1, i.e. $\log (t)$ and $t$ , so the first ancillary parameter is $γ_{1}$ and the second $γ_{2}$ . In just a few lines of code we have defined our model framework, which can now be used with anything specified in the linear predictor when we fit our merlin models. This provides a very powerful modelling framework.

Let’s now fit a model using our userhaz() function. We can call merlin as follows,

Stata

. merlin (time    age female,                     ///
>                 family(user, hfunc(userhaz)     ///
>                 failure(infect) nap(2)))

Fitting full model:

Iteration 0:  Log likelihood = -160708.35
Iteration 1:  Log likelihood = -3131.2067  (not concave)
Iteration 2:  Log likelihood = -2697.1808
Iteration 3:  Log likelihood =   -2637.83
Iteration 4:  Log likelihood = -2632.5529
Iteration 5:  Log likelihood = -2632.5464
Iteration 6:  Log likelihood = -2632.5464

Fixed effects regression model                             Number of obs = 600
Log likelihood = -2632.5464
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
time:        |
         age |   .0186875     .00304     6.15   0.000     .0127293    .0246457
      female |  -.7994153   .1042907    -7.67   0.000    -1.003821   -.5950094
       _cons |  -7.618134   .3983978   -19.12   0.000    -8.398979   -6.837289
        ap:1 |   .5150055   .1039266     4.96   0.000     .3113131    .7186979
        ap:2 |  -.0069097   .0009158    -7.55   0.000    -.0087046   -.0051148
------------------------------------------------------------------------------

I’m telling merlin that I want to fit a model with a user defined family, and in particular I provide the name of the Mata function through hfunction() (abbreviated to hfunc() in the code). The survival time variable and event indicator are declared as normal. I also tell it that there are 2 ancillary parameters to estimate through nap(2). In my linear predictor I’ve adjusted for age and female.

Endless extensions…

Given that we inevitably have correlation between events suffered by the same patient, we can now add in a random intercept at the patient level to account for this, by adding M1[patient]@1,

Stata

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

Fitting fixed effects model:

Fitting full model:

Iteration 0:  Log likelihood =   -2628.88
Iteration 1:  Log likelihood = -2613.8854
Iteration 2:  Log likelihood = -2613.3848
Iteration 3:  Log likelihood = -2613.3808
Iteration 4:  Log likelihood = -2613.3808

Mixed effects regression model                             Number of obs = 600
Log likelihood = -2613.3808
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
time:        |
         age |   .0247622   .0046504     5.32   0.000     .0156476    .0338768
      female |  -1.134893   .1645169    -6.90   0.000     -1.45734   -.8124458
 M1[patient] |          1          .        .       .            .           .
       _cons |  -8.674229   .5033374   -17.23   0.000    -9.660752   -7.687706
        ap:1 |   .6951221   .1181227     5.88   0.000     .4636059    .9266384
        ap:2 |  -.0060641   .0009568    -6.34   0.000    -.0079394   -.0041888
-------------+----------------------------------------------------------------
patient:     |
      sd(M1) |   .9134782   .1166026                      .7112873    1.173144
------------------------------------------------------------------------------

This gives us a standard deviation for the random intercept of $σ = 0.913$ , indicating substantial heterogeneity between patients. We simulated $σ = 1$ .

This is now the model we simulated from, so we can check it against the truth. The log hazard ratios for age and being female are estimated at 0.0248 and −1.135 (we simulated 0.02 and −1), and $γ_{0}$ , $γ_{1}$ and $γ_{2}$ , the _cons, ap:1 and ap:2 rows, at −8.674, 0.695 and −0.006064 (we simulated −8, 0.5 and −0.004). All but $γ_{2}$ are within two standard errors of the truth, and $γ_{2}$ is only a little over two away, which among six estimates is not unusual.

Now look back at the first model, which ignored the clustering. Its log hazard ratios for age and being female, 0.0187 and −0.799, are pulled towards zero, which is what leaving a frailty out of a proportional hazards model does.

We can investigate non-proportional hazards, for example in the effect of age as follows, remembering to add the timevar() option,

Stata

. merlin (time    age age#fp(time, powers(0))     ///
>                 female M1[patient]@1,           ///
>                 family(user, hfunc(userhaz)     ///
>                 failure(infect) nap(2))         ///
>                 timevar(time))
variables created for model 1, component 2: _cmp_1_2_1 to _cmp_1_2_1

Fitting fixed effects model:

Fitting full model:

Iteration 0:  Log likelihood = -2628.9738
Iteration 1:  Log likelihood = -2613.7621
Iteration 2:  Log likelihood = -2613.3033
Iteration 3:  Log likelihood = -2613.2992
Iteration 4:  Log likelihood = -2613.2993

Mixed effects regression model                             Number of obs = 600
Log likelihood = -2613.2993
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
time:        |
         age |   .0182832   .0165632     1.10   0.270      -.01418    .0507464
    age#fp() |    .001423   .0034961     0.41   0.684    -.0054291    .0082752
      female |  -1.137362   .1648652    -6.90   0.000    -1.460492   -.8142327
 M1[patient] |          1          .        .       .            .           .
       _cons |  -8.388509   .8549694    -9.81   0.000    -10.06422     -6.7128
        ap:1 |   .6316927   .1937822     3.26   0.001     .2518865    1.011499
        ap:2 |  -.0060151   .0009631    -6.25   0.000    -.0079027   -.0041276
-------------+----------------------------------------------------------------
patient:     |
      sd(M1) |   .9168317   .1166251                      .7145178     1.17643
------------------------------------------------------------------------------

I’ve formed a multiplicative interaction between age and $\log (t)$ by using the # notation with the fp() element. Here there is no evidence that the effect of age changes over time (the interaction is 0.0014, p = 0.684), which is as it should be: we simulated an effect of age that is constant over time.

There’s a lot more we can do of course – the use of fractional polynomials is merely an example, since within our Mata function, we can utilise anything we like. This example has hopefully set the scene for revealing some of merlin’s more advanced capabilities.

MethodSurvival analysis

A user-defined / custom hazard model

Defining a new model

Coding it up

Endless extensions…

Need this applied to your own data?