# Weibull survival model with random intercept

We can define a multilevel survival model as follows,

$$h(t_{ij}) = h_{0}(t_{ij}) \exp(\mathbf{X}_{ij}\mathbf{\beta} + u_{i})$$

where

$$h_{0}(t_ij) = \lambda \gamma t_{ij} ^ {\gamma-1}$$

and

$$u_{i} \sim N(0,\sigma^{2})$$

We have our vector of baseline covariates $\mathbf{X}_{ij}$ and associated conditional log hazard ratios $\mathbf{\beta}$ (conditional on the frailty, $u_{i}$). 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 $X_{1}$ (coded 0/1), and a continuous covariate, $X_{2}$, 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 x1 x2 M1[centre]@1, family(weibull, failure(event)))

Fitting fixed effects model:

Fitting full model:

Iteration 0:   log likelihood = -2093.3598
Iteration 1:   log likelihood = -1967.5735
Iteration 2:   log likelihood = -1963.2057
Iteration 3:   log likelihood = -1963.2003
Iteration 4:   log likelihood = -1963.2003

Mixed effects regression model                  Number of obs     =      6,000
Log likelihood = -1963.2003
------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
stime:       |
x1 |   1.008781   .0363988    27.71   0.000     .9374409    1.080122
x2 |  -.9623624   .0603866   -15.94   0.000    -1.080718   -.8440067
M1[centre] |          1          .        .       .            .           .
_cons |   .1779824   .0976107     1.82   0.068     -.013331    .3692959
log(gamma) |   .6896406    .013228    52.13   0.000     .6637142     .715567
-------------+----------------------------------------------------------------
centre:      |
sd(M1) |   .8972619   .0672296                      .7747128    1.039196
------------------------------------------------------------------------------



which gives us an estimate of $\sigma =$ 0.897 (95% CI: 0.775, 1.039), showing substantial heterogeneity between centres.