# 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


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

megenreg (stime x1 x2 M1[centre], family(weibull, failure(event)))

Getting starting values:

Fitting full model:

Iteration 0:   log likelihood = -2090.3163
Iteration 1:   log likelihood = -1967.0645
Iteration 2:   log likelihood = -1963.2032
Iteration 3:   log likelihood = -1963.2003
Iteration 4:   log likelihood = -1963.2003

Log likelihood = -1963.2003                     Number of obs     =      6,000

------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
ln_lambda1   |
x1 |   1.008781   .0363988    27.71   0.000     .9374409    1.080122
x2 |  -.9623624   .0603866   -15.94   0.000    -1.080718   -.8440067
_cons |   .1779824   .0976107     1.82   0.068     -.013331    .3692958
-------------+----------------------------------------------------------------
ln_gamma1    |
_cons |   .6896406    .013228    52.13   0.000     .6637142     .715567
-------------+----------------------------------------------------------------
lns1_1       |
_cons |  -.1084075   .0749276    -1.45   0.148    -.2552628    .0384478
------------------------------------------------------------------------------


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