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.

This example used `merlin`

version 1.15.1.