Multilevel survival models

Let's see it work

Suppose we are interested in modeling the effects of laparoscopic surgery vs. invasive surgery (lap_surg) and patient's age (age) on length of hospital stay for adult patients with appendicitis. We believe that doctors affect the length of a patient's stay, so we include a random effect for doctor (doctor).

Before we fit the model, we use stset to specify that we have survival-time data. Here, los represents the patient's length of stay in hours. We observe the length of stay for all patients; there is no censoring. Therefore, we simply identify our survival-time variable by typing

. stset los

We fit a Weibull model for length of stay:

. mestreg lap_surg age || doctor:, distribution(weibull)

        Failure _d: 1 (meaning all fail)
  Analysis time _t: los

(output omitted)

Mixed-effects Weibull PH regression             Number of obs     =        457
Group variable: doctor                          Number of groups  =        185

                                                Obs per group:
                                                              min =          1
                                                              avg =        2.5
                                                              max =          8

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(2)      =      81.66
Log likelihood = -1712.8779                     Prob > chi2       =     0.0000



          _t 
 
 Haz. ratio   Std. err.      z    P>|z|     [95% conf. interval]

    lap_surg 
 
   7.606505    1.81996     8.48   0.000     4.759086    12.15757
         age 
 
   .9565034   .0149671    -2.84   0.004     .9276137    .9862928
       _cons 
 
   2.48e-21   5.79e-21   -20.32   0.000     2.55e-23    2.41e-19

       /ln_p 
 
    2.37308   .0490428                      2.276958    2.469202

doctor       
 

  var(_cons) 
 
   1.479602   .3157434                      .9738693    2.247963

Note: Estimates are transformed only in the first equation to hazard ratios.
Note: _cons estimates baseline hazard (conditional on zero random effects).
LR test vs. Weibull model: chibar2(01) = 80.96        Prob >= chibar2 = 0.0000

Remember that "failure" is the happy event of departing the hospital. With a hazard ratio greater than 1, laparoscopic surgery is associated with shorter hospital stays. Increased age is associated with longer hospital stays. We also see a large variation across doctors that might have contaminated our results had we not taken them into account.

We can plot the marginal survivor function for surgery methods:

After about two days (48 hours), the marginal probability of remaining in the hospital begins to decrease quickly for laparoscopic surgery patients (lap_surg=0). Those having more invasive surgery (lap_surg=1) can expect to remain in the hospital longer.

We could fit a similar model using streg with shared frailties, but streg assumes the frailties follow a gamma distribution. mestreg makes the often more plausible assumption that random effects are normally distributed, meaning frailties are lognormal. mestreg also lets us extend the types of models that we can fit beyond two-level models with random intercepts.

Suppose we also believe that unobserved hospital characteristics such as discharge procedures may also impact the length of hospital stay. We will assume that doctors practice in only one hospital, so patients are nested within doctors and doctors nested within hospitals. We fit a three-level model with hospital-level (hospital) and doctor-level (doctor) random effects.

.  mestreg lap_surg age || hospital: || doctor:, distribution(weibull)

        Failure _d: 1 (meaning all fail)
  Analysis time _t: los

(output omitted)

Mixed-effects Weibull PH regression               Number of obs     =        457 

        Grouping information



                
 
     No. of       Observations per group
 Group variable 
 
     groups    Minimum    Average    Maximum
 
 
 
       hospital 
 
         12          1       38.1         73
         doctor 
 
        185          1        2.5          8



Integration method: mvaghermite                    Integration pts.  =          7

                                                   Wald chi2(2)      =      99.67
Log likelihood = -1676.8122                        Prob > chi2       =     0.0000


             _t 
 
 Haz. ratio   Std. err.      z    P>|z|     [95% conf. interval]

       lap_surg 
 
   7.025538   1.410615     9.71   0.000      4.73993    10.41327
            age 
 
   .9716887   .0127353    -2.19   0.028     .9470459    .9969728
          _cons 
 
   3.34e-21   7.60e-21   -20.70   0.000     3.84e-23    2.89e-19

          /ln_p 
 
   2.369251   .0477096                      2.275742     2.46276

hospital        
 
                                                                
 
     var(_cons) 
 
   .8649538   .4363381                      .3218047    2.324842

hospital>doctor 
 
                                                                
     var(_cons) 
 
   .7086615   .1753561                      .4363265    1.150975

Note: Estimates are transformed only in the first equation to hazard ratios.
Note: _cons estimates baseline hazard (conditional on zero random effects).
LR test vs. Weibull model: chi2(2) = 153.09               Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.

We draw the same conclusions with regard to age and having laparoscopic surgery as we drew with the previous model. We now find that there is large variation across both doctors and hospitals.