Bayesian parametric survival models

Let's see it work

Let's look at several examples.

• Parametric survival models

• Multilevel survival models

Parametric survival models

Consider a dataset in which we model the time until hip fracture as a function of age and whether the patient wears a hip-protective device (variable protect). Let's fit a Bayesian Weibull model to these data and compare the results with the classical analysis.

First, we declare our survival data.

. stset time1, id(id) failure(fracture) time0(time0)

Then, we fit a Weibull survival model using streg.

. streg protect age, distribution(weibull)

         failure _d:  fracture
   analysis time _t:  time1
                 id:  id

Weibull PH regression

No. of subjects =          148                  Number of obs    =         206
No. of failures =           37
Time at risk    =         1703
                                                LR chi2(2)       =       49.97
Log likelihood  =   -77.446477                  Prob > chi2      =      0.0000



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

     protect 
 
   .0922046   .0321722    -6.83   0.000     .0465318    .1827072
         age 
 
   1.101041    .038173     2.78   0.005     1.028709    1.178459
       _cons 
 
    .000024   .0000624    -4.09   0.000     1.48e-07    .0039042

       /ln_p 
 
   .4513032   .1265975     3.56   0.000     .2031767    .6994297

           p 
 
   1.570357   .1988033                      1.225289    2.012605
         1/p 
 
   .6367977    .080617                      .4968686     .816134
Note: Estimates are transformed only in the first equation.
Note: _cons estimates baseline hazard.

Finally, to fit a Bayesian survival model, we simply prefix the above streg command with bayes:.

. bayes: streg protect age, distribution(weibull)

Model summary

 
Likelihood:
  _t ~ streg_weibull(xb__t,{ln_p})

Priors:
  {_t:protect age _cons} ~ normal(0,10000)                                 (1)
                  {ln_p} ~ normal(0,10000)                                    
 
(1) Parameters are elements of the linear form xb__t.

Bayesian Weibull PH regression                   MCMC iterations  =     12,500
Random-walk Metropolis-Hastings sampling         Burn-in          =      2,500
                                                 MCMC sample size =     10,000
No. of subjects =        148                     Number of obs    =        206
No. of failures =         37
No. at risk     =       1703
                                                 Acceptance rate  =       .368
                                                 Efficiency:  min =     .05571
                                                              avg =     .09994
Log marginal-likelihood = -107.88854                          max =      .1767



             
 
                                                Equal-tailed    
             
 
Haz. ratio   Std. dev.     MCSE     Median  [95% cred. interval]

_t           
 
 
     protect 
 
  .0956023   .0338626   .001435   .0899154   .0463754   .1787249
         age 
 
  1.103866   .0379671   .001313   1.102685   1.033111   1.180283
       _cons 
 
  .0075815   .0411427   .000979    .000567   4.02e-06   .0560771

        ln_p 
 
  .4473869   .1285796   .004443   .4493192   .1866153   .6912467
Note: Estimates are transformed only in the first equation.
Note: _cons estimates baseline hazard.
Note: Default priors are used for model parameters.

Because the default priors used are noninformative for these data, the above results are similar to those obtained from streg. Instead of the default priors, you can specify your own; see Custom priors.

The hazard ratios are reported by default, but you can use the nohr option with bayes, during estimation or on replay, to report coefficients. Alternatively, you can specify this option with streg during estimation.

. bayes, nohr



             
 
                                                Equal-tailed    
             
 
      Mean   Std. dev.     MCSE     Median  [95% cred. interval]

_t           
 
 
    protect 
 
 -2.407909   .3482806   .015077  -2.408886  -3.070986  -1.721908
        age 
 
  .0982285   .0343418   .001189   .0977484   .0325748    .165754
      _cons 
 
 -7.561389   2.474563   .084712  -7.475201  -12.42343  -2.881028

       ln_p 
 
  .4473869   .1285796   .004443   .4493192   .1866153   .6912467

Unlike streg, bayes: streg reports only the log of the shape parameter. We can use the bayesstats summary command ([BAYES] bayesstats summary) to obtain the estimates of the shape parameter and its reciprocal.

. bayesstats summary (p: exp({ln_p})) (sigma: 1/exp({ln_p}))

Posterior summary statistics                      MCMC sample size =    10,000

           p : exp({ln_p})
       sigma : 1/exp({ln_p})



             
 
                                                Equal-tailed    
             
 
      Mean   Std. dev.     MCSE     Median  [95% cred. interval]

          p 
 
  1.577122    .201685   .006993   1.567245   1.205164   1.996203
      sigma 
 
  .6446338   .0839366   .002879   .6380624   .5009511   .8297629