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.

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)
}