I have a data set (patients (15000) nested in hospitals (50)), with a number of covariates. I built a proportional hazards model in Stata, and entered site as a frailty term. I would like to report how survival varies between hospitals. I can extract the frailties and report these, but this has no straightforward, face-value interpretation.
I have experimented with the following:
- Extract the (estimated) baseline hazard, and multiply this for each site and plot the different curves. There is no way to quantify the certainty of the frailty with this, and it seems to use a lot of ink to describe the variation in one parameter.
- Report a median (or other percentile) survival for each site, and do so ideally using something like a funnel plot. I have done this by extracting the (estimated) baseline hazard after centering the covariates, and then multiplying this by the random effect, and then finding the greatest reported survival time below my chosen centile. For now, I have plotted this against the sample size (for that hospital). There are clearly problems with bounding on the y-axis, and I do not know how I might construct control limits.
- This reference recommends Silcocks P. Hazard ratio funnel plots for survival comparisons. J Epidemiol Community Health. 2009;63:856-861., and uses a fixed-effects model, eventually reporting relative, centred hazard ratios.
Does anyone have any recommendations or comments?
NB: I am using Stata if this is relevant.