# Simulation and estimation of three-level survival models: IPD meta-analysis of recurrent event data

In this example I’ll look at the analysis of clustered survival data with three levels. This kind of data arises in the meta-analysis of recurrent event times, where we have observations (events or censored), k (level 1), nested within patients, j (level 2), nested within trials, i (level 3).

### Random intercepts

The first example will consider a model with a random intercept at level 2 and a random intercept at level 3,

$$h_{ijk}(t) = h_{0}(t)\exp (X_{ij}\beta + b_{1i} + b_{2ij})$$

where

$$b_{1i} \sim N(0,\sigma_{1}^{2}) \quad\text{and}\quad b_{2ij} \sim N(0,\sigma_{2}^{2})$$

By assuming a random intercept at the trial level 3, we are assuming that each trial comes from a distribution of trials, which arguably is a strong and perhaps implausible assumption. We may prefer to stratify by trial, but for the sake of illustration, here I’ll proceed with a random trial effect. A second important note is that our covariate effects, , are interpreted as conditional on the random effects.

I’ll start by showing how to simulate data from such a model. Let’s assume we have data from 15 trials, so we generate a `trial`

variable which will indicate which trial the observations come from, and also our trial level random effect

```
. clear
. set seed 498093
. set obs 15
Number of observations (_N) was 0, now 15.
. gen trial = _n
. gen b1 = rnormal(0,1)
```

Now let’s assume we recruit 100 patients within each trial. We can use `expand`

for this,

```
. expand 100
(1,485 observations created)
. sort trial
```

where `expand`

will replace each existing observation with 100 copies of itself. We then generate a unique patient `id`

variable, and our patient level random effect

```
. gen id = _n
. gen b2 = rnormal(0,1)
```

We also generate a binary treatment variable, at the patient level,

`. gen trt = runiform()>0.5`

Finally, we allow each patient to experience up to two events, and here we can use `expand`

again,

```
. expand 2
(1,500 observations created)
. sort trial id
```

Now we can simulate our observation level survival times, by using `survsim`

. I’ll assume a Weibull baseline hazard function with and , and a log hazard ratio of for the treatment effect,

```
. survsim stime died, distribution(weibull) /// baseline distribution
> lambda(0.1) gamma(1.2) /// baseline parameters
> maxtime(5) /// admin. censoring
> covariates(trt -0.5 /// trt effect
> b1 1 /// random int. lev. 1
> b2 1) // random int. lev. 2
```

By incorporating the `b1`

and `b2`

variables into the linear predictor with a coefficient of 1, we have a very convenient way of including random intercepts into our data-generating process.

Note that I’m simulating the recurrent events under a clock-reset approach, i.e. after each event the timescale is reset to 0. In a further example I will look at the clock-forward approach, i.e. one continuous timescale.

We can fit our data-generating model, i.e. the true model, with `merlin`

as follows,

```
. merlin (stime /// recurrent event times
> trt /// treatment (fixed effect)
> M1[trial]@1 /// random int. at trial level
> M2[trial>id]@1 /// random int. at id level
> , ///
> family(weibull, failure(died))) // distribution & event indicator
Fitting fixed effects model:
Fitting full model:
Iteration 0: log likelihood = -4250.898 (not concave)
Iteration 1: log likelihood = -4183.9242
Iteration 2: log likelihood = -4176.4582
Iteration 3: log likelihood = -4175.4024
Iteration 4: log likelihood = -4175.3996
Iteration 5: log likelihood = -4175.3996
Mixed effects regression model Number of obs = 3,000
Log likelihood = -4175.3996
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
stime: |
trt | -.5449902 .0757995 -7.19 0.000 -.6935545 -.396426
M1[trial] | 1 . . . . .
M2[trial>id] | 1 . . . . .
_cons | -1.99857 .2224338 -8.99 0.000 -2.434532 -1.562608
log(gamma) | .192227 .0259834 7.40 0.000 .1413005 .2431534
-------------+----------------------------------------------------------------
trial: |
sd(M1) | .8217196 .156671 .5654985 1.194031
-------------+----------------------------------------------------------------
id: |
sd(M2) | .971307 .0566649 .8663601 1.088967
------------------------------------------------------------------------------
```

alternatively, we can make our lives easier by using `stmixed`

```
. stset stime, failure(died)
Survival-time data settings
Failure event: died!=0 & died<.
Observed time interval: (0, stime]
Exit on or before: failure
--------------------------------------------------------------------------
3,000 total observations
0 exclusions
--------------------------------------------------------------------------
3,000 observations remaining, representing
1,598 failures in single-record/single-failure data
10,166.366 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 5
. stmixed trt || trial: || id: , distribution(weibull)
Random effect M1: Intercept at level trial
Random effect M2: Intercept at level id
Fitting fixed effects model:
Fitting full model:
Iteration 0: log likelihood = -4250.898 (not concave)
Iteration 1: log likelihood = -4183.9242
Iteration 2: log likelihood = -4176.4582
Iteration 3: log likelihood = -4175.4024
Iteration 4: log likelihood = -4175.3996
Iteration 5: log likelihood = -4175.3996
Mixed effects survival model Number of obs = 3,000
Log likelihood = -4175.3996
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_t: |
trt | -.5449902 .0757995 -7.19 0.000 -.6935545 -.396426
M1[trial] | 1 . . . . .
M2[trial>id] | 1 . . . . .
_cons | -1.99857 .2224338 -8.99 0.000 -2.434532 -1.562608
log(gamma) | .192227 .0259834 7.40 0.000 .1413005 .2431534
-------------+----------------------------------------------------------------
trial: |
sd(M1) | .8217196 .156671 .5654985 1.194031
-------------+----------------------------------------------------------------
id: |
sd(M2) | .971307 .0566649 .8663601 1.088967
------------------------------------------------------------------------------
```

where we reassuringly get the same results. Read more about `stmixed`

in this Stata Journal paper.

### Random trial and random treatment effect

Now we consider a second example, but still sticking within a meta-analysis context. In many situations, we would investigate the presence of heterogeneity, i.e. unobserved variability, in any treatment effects. We can do this by modelling a random treatment effect, with the following model.

$$h_{ijk}(t) = h_{0}(t)\exp (b_{11i} + (\beta + b_{12i})X_{ij} + b_{2ij})$$

where

$$b_{11i} \sim N(0,\sigma_{11}^{2}), \quad b_{12i} \sim N(0,\sigma_{12}^{2}) \quad \text{and}\quad b_{2ij} \sim N(0,\sigma_{2}^{2})$$

Once more let’s assume we have data from 15 trials, so we generate a `trial`

variable which will indicate which trial the observations come from, but this time we have a random intercept and also a random treatment effect, which we assume are independent (we could use `drawnorm`

if we wanted to impose correlated random effects),

```
. clear
. set seed 498093
. set obs 15
Number of observations (_N) was 0, now 15.
. gen trial = _n
. gen b11 = rnormal(0,1)
. gen b12 = rnormal(0,1)
```

Again let’s assume we recruit 100 patients within each trial,

```
. expand 100
(1,485 observations created)
. sort trial
```

and then generate our unique patient `id`

variable, and our patient level random effect

```
. gen id = _n
. gen b2 = rnormal(0,1)
```

We also generate a binary treatment variable, at the patient level,

`. gen trt = runiform()>0.5`

Finally, we allow each patient to experience up to two events,

```
. expand 2
(1,500 observations created)
. sort trial id
```

Now to impose our random treatment effect (varying at the trial level), we can first

generate its overall effect using

`. gen trtb12 = (-0.5 + b12) * trt`

and then simulate our survival times as follows,

```
. survsim stime died, dist(weibull) /// distribution
> lambda(0.1) /// scale
> gamma(1.2) /// shape
> maxtime(5) /// admin. censoring
> covariates( /// linear predictor
> trtb12 1 /// fixed + random treatment effect
> b11 1 /// random intercept
> b2 1) // random intercept
```

where once again we include coefficients of 1 on our random effects. We can fit our data-generating model, i.e. the true model, with `merlin`

as follows,

```
. merlin (stime /// outcome
> trt /// fixed treatment effect
> trt#M1[trial]@1 /// random treatment effect
> M2[trial]@1 /// random int. at trial level
> M3[trial>id]@1 /// random int. at id level
> , ///
> family(weibull, failure(died))) // dist. & event indicator
Fitting fixed effects model:
Fitting full model:
Iteration 0: log likelihood = -4116.4692 (not concave)
Iteration 1: log likelihood = -4016.0111
Iteration 2: log likelihood = -4009.4162
Iteration 3: log likelihood = -4009.1759
Iteration 4: log likelihood = -4009.1752
Iteration 5: log likelihood = -4009.1752
Mixed effects regression model Number of obs = 3,000
Log likelihood = -4009.1752
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
stime: |
trt | -.5156967 .3010874 -1.71 0.087 -1.105817 .0744238
trt#M1[tri~] | 1 . . . . .
M2[trial] | 1 . . . . .
M3[trial>id] | 1 . . . . .
_cons | -2.027292 .2255663 -8.99 0.000 -2.469394 -1.585191
log(gamma) | .2010559 .0258223 7.79 0.000 .1504451 .2516668
-------------+----------------------------------------------------------------
trial: |
sd(M1) | 1.120972 .2208733 .7618621 1.649351
sd(M2) | .8307332 .1644114 .5636366 1.224402
-------------+----------------------------------------------------------------
id: |
sd(M3) | 1.03136 .0559169 .9273872 1.14699
------------------------------------------------------------------------------
```

And of course, we can replicate this with `stmixed`

using,

```
. stset stime, failure(died)
Survival-time data settings
Failure event: died!=0 & died<.
Observed time interval: (0, stime]
Exit on or before: failure
--------------------------------------------------------------------------
3,000 total observations
0 exclusions
--------------------------------------------------------------------------
3,000 observations remaining, representing
1,608 failures in single-record/single-failure data
9,961.735 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 5
. stmixed trt || trial: trt || id: , distribution(weibull)
Random effect M1: Intercept at level trial
Random effect M2: trt at level trial
Random effect M3: Intercept at level id
Fitting fixed effects model:
Fitting full model:
Iteration 0: log likelihood = -4116.4692 (not concave)
Iteration 1: log likelihood = -4016.6173
Iteration 2: log likelihood = -4009.4527
Iteration 3: log likelihood = -4009.1758
Iteration 4: log likelihood = -4009.1752
Iteration 5: log likelihood = -4009.1752
Mixed effects survival model Number of obs = 3,000
Log likelihood = -4009.1752
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
_t: |
trt | -.5156688 .3010636 -1.71 0.087 -1.105743 .0744049
M1[trial] | 1 . . . . .
trt#M2[tri~] | 1 . . . . .
M3[trial>id] | 1 . . . . .
_cons | -2.027279 .2255593 -8.99 0.000 -2.469368 -1.585191
log(gamma) | .2010557 .0258222 7.79 0.000 .150445 .2516663
-------------+----------------------------------------------------------------
trial: |
sd(M1) | .8307552 .1643872 .5636896 1.224352
sd(M2) | 1.121014 .2210486 .7616681 1.649895
-------------+----------------------------------------------------------------
id: |
sd(M3) | 1.031359 .0559165 .9273866 1.146988
------------------------------------------------------------------------------
```

This post gave a little taster of simulating and estimating multilevel survival models – go crazy changing the baseline model, more levels, time-dependent effects…the list goes on.

## Latest Resources

Videos

### State-of-the-art statistical models for modern HTA

Videos