This example takes a look at incorporating a frailty, or random intercept, into a flexible parametric survival model, and how to fit them in Stata. First we’ll use `merlin`

to estimate our model, and then the more user-friendly wrapper function `stmixed`

. More details on these models can be found in the following papers:

- Crowther MJ, Look MP, Riley RD. Multilevel mixed effects parametric survival models using adaptive Gauss-Hermite quadrature with application to recurrent events and individual participant data meta-analysis.
*Statistics in Medicine*2014;33(22):3844-3858. - Crowther MJ. Multilevel mixed-effects parametric survival analysis: Estimation, simulation, and application.
*The Stata Journal*2019;19(4):931-949.

We can define a multilevel proportional hazards survival model as follows,

where

Now the flexible parametric model specifies the linear predictor on the log cumulative hazard scale, so in actual fact we have,

but with no time-dependent effects, log cumulative hazard ratios can be interpreted as log hazard ratios. So, we have our vector of baseline covariates and associated conditional log hazard ratios (conditional on the frailty, ). We can load an example simulated dataset from the `stmixed`

package

`. use http://fmwww.bc.edu/repec/bocode/s/stmixed_example1, clear`

which represents a multi-centre trial scenario, with 100 centres and each centre recuiting 60 patients, resulting in 6000 observations. Two covariates were collected, a binary covariate (coded 0/1), and a continuous covariate, , within the range [0,1].

Let’s inspect the first few rows of our dataset,

```
. list centre stime event x1 x2 in 1/5, noobs
+-------------------------------------------+
| centre stime event x1 x2 |
|-------------------------------------------|
| 1 .2714585 1 1 .4306063 |
| 1 .7353575 1 0 .5285968 |
| 1 .8213144 1 0 .1442798 |
| 1 .758931 1 1 .8612115 |
| 1 .9547321 0 0 .501694 |
+-------------------------------------------+
```

where `centre`

denotes the centre number for each patient, `stime`

contains our survival time, `event`

our indicator variable for those that died (1) or were censored (0), `x1`

our binary covariate, and `x2`

our continuous covariate. We can fit our model as follows

```
. merlin (stime /// survival time
> x1 x2 /// fixed covariates
> M1[centre]@1 /// random intercept
> , ///
> family(rp, df(3) failure(event))) // baseline distribution
variables created: _rcs1_1 to _rcs1_3
Fitting fixed effects model:
Fitting full model:
Iteration 0: log likelihood = -2195.8633
Iteration 1: log likelihood = -1964.8829
Iteration 2: log likelihood = -1962.755
Iteration 3: log likelihood = -1962.7449
Iteration 4: log likelihood = -1962.7449
Mixed effects regression model Number of obs = 6,000
Log likelihood = -1962.7449
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
stime: |
x1 | 1.008152 .0367602 27.43 0.000 .936103 1.0802
x2 | -.96152 .0605554 -15.88 0.000 -1.080206 -.8428335
M1[centre] | 1 . . . . .
_cons | -1.499116 .0994654 -15.07 0.000 -1.694065 -1.304167
-------------+----------------------------------------------------------------
centre: |
sd(M1) | .8965598 .0673174 .7738693 1.038702
------------------------------------------------------------------------------
Warning: Baseline spline coefficients not shown - use ml display
```

The key part of the syntax is specifying a normally distributed random effect, which we’ve called `M1`

(it must start with a capital letter, then a number), with the variable defining the cluster in square brackets `[centre]`

, and finally, we tell `merlin`

not to estimate a coefficient, but rather constrain it to be 1. Our model converges nicely, and we observe an estimate of 0.897 (95% CI: 0.774, 1.039), showing substantial heterogeneity between centres.

`merlin`

has a challenging syntax because it emcompasses a lot of different modelling frameworks. To make it more user-friendly, we’ve developed a few wrapper functions to make our lives easier. One of these is `stmixed`

which is specifically designed to have more consistent syntax for specifying a random effects model in Stata. It essentially follows the same design as `mestreg`

. It also has the added benefit of requiring the data to be `stset`

before use. So we can fit the same model with,

```
. stset stime, failure(event)
Survival-time data settings
Failure event: event!=0 & event<.
Observed time interval: (0, stime]
Exit on or before: failure
--------------------------------------------------------------------------
6,000 total observations
0 exclusions
--------------------------------------------------------------------------
6,000 observations remaining, representing
3,427 failures in single-record/single-failure data
3,550.738 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 1.985816
. stmixed x1 x2 || centre:, distribution(rp) df(3)
Random effect M1: Intercept at level centre
variables created: _rcs1_1 to _rcs1_3
Fitting fixed effects model:
Fitting full model:
Iteration 0: log likelihood = -2195.8633
Iteration 1: log likelihood = -1964.8829
Iteration 2: log likelihood = -1962.755
Iteration 3: log likelihood = -1962.7449
Iteration 4: log likelihood = -1962.7449
Mixed effects survival model Number of obs = 6,000
Log likelihood = -1962.7449
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_t: |
x1 | 1.008152 .0367602 27.43 0.000 .936103 1.0802
x2 | -.96152 .0605554 -15.88 0.000 -1.080206 -.8428335
M1[centre] | 1 . . . . .
_cons | -1.499116 .0994654 -15.07 0.000 -1.694065 -1.304167
-------------+----------------------------------------------------------------
centre: |
sd(M1) | .8965598 .0673174 .7738693 1.038702
------------------------------------------------------------------------------
Warning: Baseline spline coefficients not shown - use ml display
```

obtaining identical estimates, unsurprisingly. `stmixed`

gives us some helpful output telling us what it’s named the random effect, and what level it has been specified.

There’s lots of predictions available post-estimation, just take a look at `help stmixed postestimation`

for more details.

Subscribe to get new post alerts directly in your inbox: