A user-defined hazard model: fractional polynomials

This tutorial will illustrate some of the more advanced capabilities of merlin when modelling survival data, but hopefully with a rather simple example. In some previous work, Paul Lambert and I developed stgenreg for modelling survival data with a general user-specified hazard function, which is then estimated using numerical integration. Before continuing with this example, please look through the accompanying material provided with stgenreg; the methods paper can be found here and a software tutorial can be found here.

Consider our standard proportional hazards model,

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

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) = \gamma_{0} + \gamma_{1} t + \gamma_{2} \log(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. I’ll use the catheter dataset to fit some models.

. . webuse catheter, clear
(Kidney data, McGilchrist and Aisbett, Biometrics, 1991)

Our dataset consists of the following,

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

  +----------------------------------------+
  | patient   time   infect   age   female |
  |----------------------------------------|
  |       1     16        1    28        0 |
  |       1      8        1    28        0 |
  |       2     13        0    48        1 |
  |       2     23        1    48        1 |
  |       3     22        1    32        0 |
  |----------------------------------------|
  |       3     28        1    32        0 |
  +----------------------------------------+

with patient our individual patient identifier, time is our time of infection at the catheter insertion point, infect is our event indicator with an event being an infection, age is 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, I 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. This is really easy to do:

. 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)
>         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 within

mata:

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 smoe intermediate vectors/matrices that we’ll need. Explicit declaration of each object’s type is good programming practice.

	real matrix linpred
	real colvector gammas

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

	linpred = merlin_util_xzb(gml,t)

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 as follows,

	gammas = merlin_util_ap(gml,1)\merlin_util_ap(gml,2)

Finally we need to return our hazard function, which is done very simply,

	return(exp(linpred :+ merlin_fp(t,(1,0)) * gammas))

This makes use of the internal merlin_fp() function, which returns fractional polynomials, in this case an FP2 function with powers 1 and 0. 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,

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

Fitting full model:

Iteration 0:   log likelihood =      -7424  
Iteration 1:   log likelihood = -338.98462  
Iteration 2:   log likelihood = -335.05093  
Iteration 3:   log likelihood = -334.39268  
Iteration 4:   log likelihood = -334.39242  
Iteration 5:   log likelihood = -334.39242  

Mixed effects regression model                  Number of obs     =         76
Log likelihood = -334.39242
------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
time:        |            
         age |   .0035436   .0092129     0.38   0.701    -.0145134    .0216005
      female |  -.8639103   .2920078    -2.96   0.003    -1.436235   -.2915854
       _cons |  -4.081803   .6379747    -6.40   0.000     -5.33221   -2.831396
        ap:1 |  -.0310417   .1456653    -0.21   0.831    -.3165404     .254457
        ap:2 |  -.0010321   .0017629    -0.59   0.558    -.0044872    .0024231
------------------------------------------------------------------------------

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 hazfunction(). 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.

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,

. 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 = -333.70064  
Iteration 1:   log likelihood = -330.15029  
Iteration 2:   log likelihood = -329.84764  
Iteration 3:   log likelihood = -329.81924  
Iteration 4:   log likelihood = -329.81923  

Mixed effects regression model                  Number of obs     =         76
Log likelihood = -329.81923
------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
time:        |            
         age |   .0077962   .0139771     0.56   0.577    -.0195984    .0351909
      female |  -1.654983   .5303959    -3.12   0.002    -2.694539   -.6154256
 M1[patient] |          1          .        .       .            .           .
       _cons |  -4.649158   .9124842    -5.10   0.000    -6.437594   -2.860722
        ap:1 |    .214335   .2014649     1.06   0.287    -.1805289    .6091989
        ap:2 |    .000749     .00213     0.35   0.725    -.0034258    .0049237
-------------+----------------------------------------------------------------
patient:     |            
      sd(M1) |   .9374593   .2864192                      .5150949    1.706151
------------------------------------------------------------------------------

Which gives us a standard deviation for the random intercept of $\sigma = $0.94, indicating substantial heterogeneity between patients.

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

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

Fitting full model:

Iteration 0:   log likelihood = -616.27647  (not concave)
Iteration 1:   log likelihood = -393.69721  (not concave)
Iteration 2:   log likelihood = -338.73509  (not concave)
Iteration 3:   log likelihood = -330.69894  (not concave)
Iteration 4:   log likelihood = -326.60214  
Iteration 5:   log likelihood = -311.11739  
Iteration 6:   log likelihood = -288.33066  (not concave)
Iteration 7:   log likelihood = -288.09968  
Iteration 8:   log likelihood = -280.39076  
Iteration 9:   log likelihood = -273.69532  (not concave)
Iteration 10:  log likelihood = -272.56442  
Iteration 11:  log likelihood =  -269.6367  
Iteration 12:  log likelihood = -269.10948  
Iteration 13:  log likelihood = -269.07633  
Iteration 14:  log likelihood = -269.07691  
Iteration 15:  log likelihood = -269.07687  
Iteration 16:  log likelihood = -269.07687  

Mixed effects regression model                  Number of obs     =         76
Log likelihood = -269.07687
------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
time:        |            
         age |   .3845645   .0681016     5.65   0.000     .2510878    .5180412
    age#fp() |  -.0850253   .0148018    -5.74   0.000    -.1140363   -.0560144
      female |  -1.805518   .6942311    -2.60   0.009    -3.166186   -.4448501
 M1[patient] |          1          .        .       .            .           .
       _cons |  -19.06475   3.330444    -5.72   0.000     -25.5923    -12.5372
        ap:1 |   3.816655   .8008698     4.77   0.000      2.24698    5.386331
        ap:2 |  -.0036826   .0028632    -1.29   0.198    -.0092944    .0019292
-------------+----------------------------------------------------------------
patient:     |            
      sd(M1) |   1.289988   .4325322                      .6686187    2.488818
------------------------------------------------------------------------------

I’ve formed an interaction between age and $\log(t)$ by using the # notation, using the fp() element. In this case we find evidence of a time-dependent effect of age. Note I also used the zeros option - if you try it without it it fails to converge, as it didn’t like my starting values. Instead it worked well with just the zero vector.