Three-level models: IPD meta-analysis of recurrent events

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_{i j k} (t) = h_{0} (t) \exp (X_{i j} β + b_{1 i} + b_{2 i j})

where

b_{1 i} ~ N (0, σ_{1}^{2}) and b_{2 i j} ~ N (0, σ_{2}^{2})

By assuming a random intercept at the trial level (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.

Simulating the data

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,

Stata

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

Stata

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

Stata

. gen id = _n

. gen b2 = rnormal(0,1)

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

Stata

. gen trt = runiform()>0.5

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

Stata

. expand 2
(1,500 observations created)

. sort trial id

Now we can simulate our observation-level survival times, by using avalon, the merlin family’s simulator, and its standard subcommand. I’ll assume a Weibull baseline hazard function with $λ = 0.1$ and $γ = 1.2$ , and a log hazard ratio of $β = - 0.5$ for the treatment effect,

Stata

. avalon standard stime died,             /// new variables
>         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. 3 (trial)
>                    b2 1)                //  random int. lev. 2 (patient)

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. For simplicity, each patient’s two gap times are drawn independently given the random effects, each censored at 5 years, so a second gap time exists even when the first was censored; in a real study it would only start after a first event. The clock-forward approach, i.e. one continuous timescale with delayed entry, is sketched in our joint frailty models tutorial.

Fitting the model

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

Stata

. 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.799
Iteration 1:  Log likelihood =   -4179.08
Iteration 2:  Log likelihood = -4175.5672
Iteration 3:  Log likelihood = -4175.5566
Iteration 4:  Log likelihood = -4175.5564

Mixed effects regression model                           Number of obs = 3,000
Log likelihood = -4175.5564
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
stime:       |
         trt |   -.544609   .0759017    -7.18   0.000    -.6933735   -.3958445
   M1[trial] |          1          .        .       .            .           .
M2[trial>id] |          1          .        .       .            .           .
       _cons |  -1.972639    .211429    -9.33   0.000    -2.387032   -1.558246
  log(gamma) |   .1926321   .0262564     7.34   0.000     .1411705    .2440937
-------------+----------------------------------------------------------------
trial:       |
      sd(M1) |   .8249628   .1522666                      .5745459    1.184524
-------------+----------------------------------------------------------------
id:          |
      sd(M2) |   .9743125   .0592236                      .8648842    1.097586
------------------------------------------------------------------------------

Because we simulated the data, we know what the model should find: a log hazard ratio for treatment of −0.5, a log baseline scale of log 0.1 = −2.303 and log shape of log 1.2 = 0.182, and random intercepts with a standard deviation of 1 at both the trial and the patient level. We estimate −0.545, −1.973 and 0.193, and standard deviations of 0.82 for trials and 0.97 for patients, all within 1.6 standard errors of the truth. The trial-level standard deviation and the log baseline scale are the least precise of these (standard errors 0.15 and 0.21), since only 15 trials inform them.

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_{i j k} (t) = h_{0} (t) \exp (b_{11 i} + (β + b_{12 i}) X_{i j} + b_{2 i j})

where

b_{11 i} ~ N (0, σ_{11}^{2}), b_{12 i} ~ N (0, σ_{12}^{2}) and b_{2 i j} ~ N (0, σ_{2}^{2})

Simulating with a random treatment effect

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

Stata

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

Stata

. expand 100
(1,485 observations created)

. sort trial

and then generate our unique patient id variable, and our patient-level random effect,

Stata

. gen id = _n

. gen b2 = rnormal(0,1)

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

Stata

. gen trt = runiform()>0.5

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

Stata

. 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 the trial-specific treatment effect, multiplied by treatment, using

Stata

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

and then simulate our survival times as follows,

Stata

. avalon standard stime died,             /// new variables
>         distribution(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.

Fitting the random-treatment-effect model

We can fit our data-generating model, i.e. the true model, with merlin as follows. With two random effects at the trial level we ask for 14 integration points, twice the default: at the default 7, the log likelihood moves by 0.12 when the points are doubled, a sign that 7 are too few here,

Stata

. 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
>         , intpoints(14)                 //  twice the default points

Fitting fixed effects model:

Fitting full model:

Iteration 0:  Log likelihood = -4116.2869
Iteration 1:  Log likelihood = -4017.6581
Iteration 2:  Log likelihood = -4008.7275
Iteration 3:  Log likelihood = -4008.6974
Iteration 4:  Log likelihood = -4008.6974

Mixed effects regression model                           Number of obs = 3,000
Log likelihood = -4008.6974
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
stime:       |
         trt |  -.5193968   .3029808    -1.71   0.086    -1.113228    .0744347
trt#M1[tri~] |          1          .        .       .            .           .
   M2[trial] |          1          .        .       .            .           .
M3[trial>id] |          1          .        .       .            .           .
       _cons |  -2.030801   .2268589    -8.95   0.000    -2.475436   -1.586165
  log(gamma) |   .2035926   .0259125     7.86   0.000     .1528051    .2543802
-------------+----------------------------------------------------------------
trial:       |
      sd(M1) |   1.127428    .222247                      .7661146    1.659144
      sd(M2) |     .83477   .1652206                      .5663619    1.230381
-------------+----------------------------------------------------------------
id:          |
      sd(M3) |    1.04181   .0575327                      .9349356      1.1609
------------------------------------------------------------------------------

Here the true values are an average log hazard ratio for treatment of −0.5, varying across trials with a standard deviation of 1 (sd(M1)), random intercepts with a standard deviation of 1 at the trial (sd(M2)) and patient (sd(M3)) levels, and the same baseline as before. We estimate −0.519 for the average treatment effect, standard deviations of 1.13, 0.83 and 1.04, and −2.031 and 0.204 for the log scale and log shape (true values −2.303 and 0.182), all within 1.2 standard errors of the truth. Notice how much less precise the average treatment effect is now, with a standard error of 0.30 against 0.08 in the first example: with a treatment effect that varies between trials, it is an average over just 15 of them.

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.

MethodsSimulation studies Recurrent events

Three-level survival models: IPD meta-analysis of recurrent events

Random intercepts

Simulating the data

Fitting the model

Random trial and random treatment effect

Simulating with a random treatment effect

Fitting the random-treatment-effect model

Need this applied to your own data?