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,

Within a general hazard model, we can essentially specify any function for our baseline hazard function, subject to the constraint that for all . 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,

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 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, 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:

```
. 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 some intermediate vectors/matrices that we’ll need. Explicit declaration of each object’s type is good programming practice and reduces the likelihood of errors.

```
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)`

Because we pass our time vector t, 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`

), which may be time dependent, 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`

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.04913
Iteration 3: log likelihood = -334.39269
Iteration 4: log likelihood = -334.39242
Iteration 5: log likelihood = -334.39242
Fixed effects regression model Number of obs = 76
Log likelihood = -334.39242
------------------------------------------------------------------------------
| Coefficient 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 .6379746 -6.40 0.000 -5.33221 -2.831396
ap:1 | -.0310417 .1456652 -0.21 0.831 -.3165403 .2544568
ap:2 | -.0010321 .0017629 -0.59 0.558 -.0044872 .002423
------------------------------------------------------------------------------
```

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()`

. 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`

```
. 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 = -336.03458
Iteration 1: log likelihood = -330.8822
Iteration 2: log likelihood = -330.0283
Iteration 3: log likelihood = -329.82039
Iteration 4: log likelihood = -329.81922
Iteration 5: log likelihood = -329.81923
Mixed effects regression model Number of obs = 76
Log likelihood = -329.81923
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
time: |
age | .0077964 .0139772 0.56 0.577 -.0195984 .0351912
female | -1.655001 .5303999 -3.12 0.002 -2.694565 -.6154359
M1[patient] | 1 . . . . .
_cons | -4.649168 .9124864 -5.10 0.000 -6.437609 -2.860728
ap:1 | .2143389 .2014649 1.06 0.287 -.1805251 .6092029
ap:2 | .000749 .00213 0.35 0.725 -.0034258 .0049238
-------------+----------------------------------------------------------------
patient: |
sd(M1) | .9374698 .28642 .5151034 1.706162
------------------------------------------------------------------------------
```

This gives us a standard deviation for the random intercept of , 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))
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 = -274.64296
Iteration 1: log likelihood = -270.67862 (not concave)
Iteration 2: log likelihood = -270.42118
Iteration 3: log likelihood = -269.14092
Iteration 4: log likelihood = -269.09187
Iteration 5: log likelihood = -269.07631
Iteration 6: log likelihood = -269.07693
Iteration 7: log likelihood = -269.07687
Iteration 8: log likelihood = -269.07687
Mixed effects regression model Number of obs = 76
Log likelihood = -269.07687
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
time: |
age | .3845644 .0681264 5.64 0.000 .2510392 .5180896
age#fp() | -.0850253 .0148071 -5.74 0.000 -.1140466 -.0560039
female | -1.805518 .6942804 -2.60 0.009 -3.166283 -.4447536
M1[patient] | 1 . . . . .
_cons | -19.06474 3.331674 -5.72 0.000 -25.5947 -12.53478
ap:1 | 3.816654 .8011637 4.76 0.000 2.246402 5.386906
ap:2 | -.0036826 .0028632 -1.29 0.198 -.0092944 .0019292
-------------+----------------------------------------------------------------
patient: |
sd(M1) | 1.289988 .4326155 .6685336 2.489132
------------------------------------------------------------------------------
```

I’ve formed an multiplicative interaction between `age`

and $\log(t)$ by using the `#`

notation, using the `fp()`

element.

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.

Subscribe to get new post alerts directly in your inbox: