Should I bootstrap at the cluster level or the individual level?

Question

I have a survival model with patients nested in hospitals that includes a random-effect for the hospitals. The random effect is gamma-distributed, and I am trying to report the 'relevance' of this term on a scale that is easily understood.

I have found the following references which use the Median Hazard Ratio (a bit like the Median Odds Ratio), and calculated this.

Bengtsson T, Dribe M: Historical Methods 43:15, 2010

However, now I wish to report the uncertainty associated with this estimate using the bootstrap. The data is survival data, and hence there are multiple observations per patient, and multiple patients per hospital. It seems obvious that I need to cluster the patient observations when re-sampling. But I don't know if I should cluster the hospitals too (i.e. resample hospitals, rather than patients?

I am wondering if the answer depends on the parameter of interest, and so would be different if the target was something that was relevant at the patient level rather than the hospital level?

I have listed the stata code below in case that helps.

cap program drop est_mhr
program define est_mhr, rclass
stcox patient_var1 patient_var2 ///
    , shared(hospital) ///
    noshow
local twoinvtheta2 = 2 / (e(theta)^2)
local mhr = exp(sqrt(2*e(theta))*invF(`twoinvtheta2',`twoinvtheta2',0.75))
return scalar mhr = `mhr'
end

bootstrap r(mhr), reps(50) cluster(hospital): est_mhr

Answer 1

Imagine that you conducted a study about children educational achievements. You took a random sample of schools from some area and from each school one class was included in the study. You conducted analysis and now want to use bootstrap to obtain confidence intervals for your estimates. How to do it?

First, notice that your data is hierarchical, it has several levels: schools, classes within schools, and students within classes. Since there is only one class per school, so the second level is nonexistent in your data. We can assume that there are some similarities within schools and differences between schools. If there are similarities within schools then if you sampled pupils at random, not taking into consideration their school membership you could possibly destroy the hierarchical structure of your data.

In general, there are several options:

1. sample students with replacement,
2. sample whole schools with replacement,
3. first sample schools with replacement and then sample students (a) with replacement, or (b) without replacement.

It appears that the first approach is the worst one. Recall that bootstrap sampling should somehow imitate the sampling process in your study and you were sampling schools rather than individual students. Choosing between (2) and (3) is more complicated, but hopefully you can find research papers considering this topic (e.g. Rena et al. 2010, Field and Welsh, 2007). Generally options (2) or (3b) are preferable as it seems that including too much levels of sampling with replacement leads to biased results. You can find more information about this topic also in books by Efron and Tibshirani (1994) and Davison and Hinkley (1997). Notice that we have similar problem with bootstrapping time-series data and in this case we also rather sample whole blocks of series (e.g. whole season if we assume seasonality) rather than individual observations because otherwise the time structure would get destroyed. In practice there is no one-size-fits-all solution but with complicated data structures you should choose such bootstrap sampling scheme that best fits your data and your problem and if possible use a simulation study to compare different solutions.

---

Davison, A.C. and Hinkley, D.V. (1997). Bootstrap Methods and Their Application. Cambridge.

Efron, B. and Tibshirani, R.J. (1994). An Introduction to the Bootstrap. CRC Press.

Ren, S., Lai, H., Tong, W., Aminzadeh, M., Hou, X., & Lai, S. (2010). Nonparametric bootstrapping for hierarchical data. Journal of Applied Statistics, 37(9), 1487-1498.

Field, C. A., & Welsh, A. H. (2007). Bootstrapping clustered data. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(3), 369-390.

Answer 2

The answer seems to be that the resampling process needs to take account of the structure of the data. There is a nice explanation here (along with some R code to implement this).

http://biostat.mc.vanderbilt.edu/wiki/Main/HowToBootstrapCorrelatedData

Thanks to the pointer from the UCLA Statistical Consulting Group.

I have written a speedier (but less flexible) version of the code snippet linked to above - check here for updates and details.

rsample2 <- function(data=tdt, id.unit=id.u, id.cluster=id.c) {
require(data.table)

setkeyv(tdt,id.cluster)
# Generate within cluster ID (needed for the sample command)
tdt[, "id.within" := .SD[,.I], by=id.cluster, with=FALSE]

# Random sample of sites
bdt <- data.table(sample(unique(tdt[[id.cluster]]), replace=TRUE))
setnames(bdt,"V1",id.cluster)
setkeyv(bdt,id.cluster)

# Use random sample of sites to select from original data
# then
# within each site sample with replacement using the within site ID
bdt <- tdt[bdt, .SD[sample(.SD$id.within, replace=TRUE)],by=.EACHI]

# return data sampled with replacement respecting clusters
bdt[, id.within := NULL] # drop id.within
return(bdt)
}