Flexible parametric survival analysis with frailty

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. We’ll use merlin to estimate our model. 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.

The model

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

h_{i j} (t) = h_{0} (t) \exp (X_{i j} β + b_{i})

where

b_{i} ~ N (0, σ^{2})

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

H_{i j} (t) = H_{0} (t) \exp (X_{i j} β + b_{i})

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 $X_{i j}$ and associated conditional log hazard ratios $β$ (conditional on the frailty, $b_{i}$ ). We’ll simulate our own data to try it out on.

The data

The data represent a multi-centre trial scenario, with 100 centres and each centre recruiting 60 patients, resulting in 6000 observations. Each centre gets its own random intercept, $b_{i}$ , drawn from a normal distribution with $σ = 1$ . Each patient has two covariates, a binary covariate $X_{1}$ (coded 0/1), and a continuous covariate, $X_{2}$ , within the range [0,1].

Stata

. clear

. set seed 2468

. set obs 100                             // centres
Number of observations (_N) was 0, now 100.

. gen centre = _n

. gen b = rnormal(0, 1)                   // centre random intercept, sd 1

. expand 60                               // 60 patients in each centre
(5,900 observations created)

. sort centre

. gen x1 = runiform() < 0.5               // binary covariate

. gen x2 = runiform()                     // continuous covariate, on [0,1]

We then simulate the survival times with avalon standard. The baseline is a two-component Weibull mixture, whose hazard falls steeply at first and then rises again, a shape no Weibull can take but a spline can follow: distribution(weibull) mixture asks for it, lambda() and gamma() give each component’s scale and shape, and pmix() the weight on the first. covariates() takes each variable with its log hazard ratio: 1 for $X_{1}$ and −1 for $X_{2}$ , and the centre’s random intercept enters with a coefficient of 1. Follow-up stops at two years, and patients still event-free then are censored there.

Stata

. avalon standard stime event,                    ///
>         distribution(weibull) mixture           /// baseline hazard:
>         lambda(0.2 0.8) gamma(1.5 0.8)          ///   a two-component
>         pmix(0.5)                               ///   Weibull mixture
>         covariates(x1 1 x2 -1 b 1)              /// log hazard ratios
>         maxtime(2)                              //  two years' follow-up

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

Stata

. list centre stime event x1 x2 in 1/5, noobs

  +--------------------------------------------+
  | centre       stime   event   x1         x2 |
  |--------------------------------------------|
  |      1   1.2599356       1    1    .496487 |
  |      1   .44048204       1    0   .5419227 |
  |      1    1.896829       1    0   .4834752 |
  |      1   .16377865       1    1   .0108418 |
  |      1   .29317754       1    1    .188637 |
  +--------------------------------------------+

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.

Fitting it with merlin

Stata

. 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 = -4406.7135
Iteration 1:  Log likelihood = -4289.8175
Iteration 2:  Log likelihood = -4288.6445
Iteration 3:  Log likelihood = -4288.6438
Iteration 4:  Log likelihood = -4288.6438

Mixed effects regression model                           Number of obs = 6,000
Log likelihood = -4288.6438
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
stime:       |
          x1 |   .9805431   .0363849    26.95   0.000       .90923    1.051856
          x2 |  -1.003614   .0608453   -16.49   0.000    -1.122869   -.8843597
  M1[centre] |          1          .        .       .            .           .
       _cons |  -1.046078   .1105226    -9.46   0.000    -1.262698   -.8294578
-------------+----------------------------------------------------------------
centre:      |
      sd(M1) |   1.029415   .0765263                      .8898409    1.190881
------------------------------------------------------------------------------
    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 (its name must be M followed by a positive integer), 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 $σ =$ 1.029 (95% CI: 0.890, 1.191), showing substantial heterogeneity between centres.

Because we simulated the data, we know what the model should find. We set $σ = 1$ . The log hazard ratios for $X_{1}$ and $X_{2}$ were 1 and −1, and are estimated at 0.981 and −1.004. All three are within one standard error of the truth. The baseline was a Weibull mixture, and three spline degrees of freedom were enough to recover the covariate effects.

Working from stset data

If your data are already stset, merlin can use the variables that stset creates. After stset stime, failure(event), the same model is merlin (_t x1 x2 M1[centre]@1, family(rp, df(3) failure(_d))), which gives exactly the estimates above, and with delayed entry you would add ltruncated(_t0) to the family() options.

There are lots of predictions available post-estimation, including marginal ones, which integrate over the frailty distribution; just take a look at help merlin postestimation for more details.

MethodFlexible parametric survival models

Flexible parametric survival analysis with frailty

The model

The data

Fitting it with merlin

Working from stset data

Need this applied to your own data?