You can fit parametric survival models in Stata using streg. You can fit multilevel parametric survival models using mestreg. You can now fit Bayesian parametric survival models by simply typing bayes: in front of streg and mestreg!
Let's look at several examples.
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 | |
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. |
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 | |
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 | |
Learn more about the general features of the bayes prefix.
Learn more about Stata's Bayesian analysis and survival-time features.
Read more about the bayes prefix and Bayesian analysis in the Stata Bayesian Analysis Reference Manual.