# Bootstrapping hierarchical/multilevel data (resampling clusters)

I am producing a script for creating bootstrap samples from the cats dataset (from the -MASS- package).

Following the Davidson and Hinkley textbook  I ran a simple linear regression and adopted a fundamental non-parametric procedure for bootstrapping from iid observations, namely pairs resampling.

The original sample is in the form:

Bwt   Hwt

2.0   7.0
2.1   7.2

...

1.9    6.8


Through an univariate linear model we want to explain cats hearth weight through their brain weight.

The code is:

library(MASS)
library(boot)

##################
#   CATS MODEL   #
##################

cats.lm <- glm(Hwt ~ Bwt, data=cats)
cats.diag <- glm.diag.plots(cats.lm, ret=T)

#######################
#   CASE resampling   #
#######################

cats.fit <- function(data) coef(glm(data$Hwt ~ data$Bwt))
statistic.coef <- function(data, i) cats.fit(data[i,])

bootl <- boot(data=cats, statistic=statistic.coef, R=999)


Suppose now that there exists a clustering variable cluster = 1, 2,..., 24 (for instance, each cat belongs to a given litter). For simplicity, suppose that data are balanced: we have 6 observations for each cluster. Hence, each of the 24 litters is made up of 6 cats (i.e. n_cluster = 6 and n = 144).

It is possible to create a fake cluster variable through:

q <- rep(1:24, times=6)
cluster <- sample(q)
c.data <- cbind(cats, cluster)


I have two related questions:

How to simulate samples in accordance with the (clustered) dataset strucure? That is, how to resample at the cluster level? I would like to sample the clusters with replacement and to set the observations within each selected cluster as in the original dataset (i.e. sampling with replacenment the clusters and without replacement the observations within each cluster).

This is the strategy proposed by Davidson (p. 100). Suppose we draw B = 100 samples. Each of them should be composed by 24 possibly recurrent clusters (e.g. cluster = 3, 3, 1, 4, 12, 11, 12, 5, 6, 8, 17, 19, 10, 9, 7, 7, 16, 18, 24, 23, 11, 15, 20, 1), and each cluster should contain the same 6 observations of the original dataset. How to do that in R? (either with or without the -boot- package.) Do you have alternative suggestions for proceeding?

The second question concerns the initial regression model. Suppose I adopt a fixed-effects model, with cluster-level intercepts. Does it change the resampling procedure adopted?

 Davidson, A. C., Hinkley, D. V. (1997). Bootstrap methods and their applications. Cambridge University press.

Resampling the whole clusters has been known in survey statistics for as long as any resampling methods have been used there at all (which is, since mid 1960s), so it is a well established method. See my collection of links at http://www.citeulike.org/user/ctacmo/tag/survey_resampling. Whether boot can do this or not, I don't know; I use survey package when I need to work with survey bootstraps, although the last time I checked, it did not have all the functionality I needed (like some small sample corrections, as far as I can recall).

I don't think applying a particular model such as fixed effects changes things much, but, IMO, the residual bootstrap makes a lot of strong assumptions (the residuals are i.i.d., the model is correctly specified). Every one of them is easily broken, and the cluster structure surely breaks the i.i.d. assumption.

There's been some econometrics literature on wild cluster bootstrap. They pretended they worked in vacuum without all those fifty years of survey statistics research into the topic, so I am not sure as to what to make of it.

• My main doubt about making fixed-effects at cluster level is that in some simulated samples it can happen that we have not selected some of the initial clusters, so that the related cluster-specific intercepts cannot be identified. If you have a look at the code I posted, it should not be a problem from a "mechanical" point of view (at each iteration we can fit a different FE model with just the sampled clusters' intercepts). I was wondering whether there is a "statistical" issue in all of this – Stefano Lombardi Jan 3 '13 at 2:06

I tried to solve the problem myself, and I produced the following code.

Although it works, it could probably be improved in terms of speed. Also, if possible I would have preferred to find a way for using the -boot- package, as it allows to automatically compute a number of bootstrapped confidence intervals through boot.ci...

For simplicity, the starting dataset consists in 18 cats (the "lower-level" observations) nested in 6 laboratories (the clustering variable). The dataset is balanced (n_cluster = 3 for each cluster). We have one regressor, x, for explaining y.

The fake dataset and the matrix where to store results are:

  # fake sample
dat <- expand.grid(cat=factor(1:3), lab=factor(1:6))
dat <- cbind(dat, x=runif(18), y=runif(18, 2, 5))

# empty matrix for storing coefficients estimates and standard errors of x
B <- 50 # number of bootstrap samples
b.sample <- matrix(nrow=B, ncol=3, dimnames=list(c(), c("sim", "b_x", "se_x")))
b.sample[,1] <- rep(1:B)


At each of the B iterations, the following loop samples 6 clusters with replacement, each composed by 3 cats sampled without replacement (i.e. the clusters' internal composition is maintained unaltered). The estimates of the regressor coefficient and of its standard error are stored in the previously created matrix:

  ####################################
#   loop through "b.sample" rows   #
####################################

for (i in seq(1:B)) {

###   sampling with replacement from the clustering variable

# sampling with replacement from "cluster"
cls <- sample(unique(dat$lab), replace=TRUE) cls.col <- data.frame(lab=cls) # reconstructing the overall simulated sample cls.resample <- merge(cls.col, dat, by="lab") ### fitting linear model to simulated data # model fit mod.fit <- function(data) glm(data$y ~ data$x) # estimated coefficients and standard errors b_x <- summary(mod.fit(data=cls.resample))$coefficients[2,1]
se_x <- summary(mod.fit(data=cls.resample))$coefficients[2,2] b.sample[i,2] <- b_x b.sample[i,3] <- se_x }  Hope this helps, Lando • using a for loop must be dominated by using replicate; as a bonus it automatically returns the b.sample array for you. Also, with all the merging here, you're almost certainly better off using data.table and resampling by key. I may contribute an answer when I get to a computer... Question: why are you keeping track of the coefficients' standard errors? – MichaelChirico Sep 27 '15 at 12:09 • Thanks @MichaelChirico, I agree. If I remember well I was saving standard errors for plotting confidence intervals later. – Stefano Lombardi Sep 27 '15 at 12:17 • shouldn't confidence intervals just be the quantiles of the distribution of bootstrap coefficients? i.e. for a 95% confidence interval, quantile(b.sample[,2], c(.025, .975)) – MichaelChirico Sep 27 '15 at 13:13 I had to do this recently and used dplyr. The solution is not as elegant as with data.table, but: library(dplyr) replicate(B, { cluster_sample <- data.frame(cluster = sample(dat$cluster, replace = TRUE))
dat_sample <- dat %>% inner_join(cluster_sample, by = 'cluster')
coef(lm(y ~ x, data = dat_sample))
})


The inner_join repeats every row having a given value x of cluster by the number of times that x appears in cluster_sample.

Here is a much simpler (and almost undoubtedly faster) way to do the bootstrapping using data.table (on @lando.carlissian 's data):

library(data.table)
setDT(dat, key = "lab")
b.sample <-
replicate(B, dat[.(sample(unique(lab), replace = T)),
glm(y ~ x)\$coefficients])


Hi a very simple solution based on split and lapply, no need specific package except "boot", exemple with an estimation of ICC based on nagakawa procedure :

# FIRST FUNCTION : "parameter assesment"
nagakawa <- function(dataICC){
#dataICC <- dbICC
modele <- lmer(indicateur.L ~ 1 + (1|sujet.L) + (1|injection.L) + experience.L, data = dataICC)
variance <- get_variance(modele)
var.fixed <- variance$$var.fixed var.random <- variance$$var.random
var.sujet <- variance$$var.intercept var.resid <- variance$$var.residual
icc.juge1 <- var.random / (var.random + var.fixed + var.resid)

modele <- lmer(indicateur.L ~ 1 + (1 + injection.L|sujet.L) + experience.L, data = dataICC)
variance <- VarCorr(modele)
var.fixed <- get_variance_fixed(modele)
var.random <- (attributes(variance$$sujet.L)$$stddev)^2 + (attributes(variance$$sujet.L)$$stddev)^2
var.sujet <- (attributes(variance$$sujet.L)$$stddev)^2
var.resid <- (attributes(variance)$$sc)^2 icc.juge2 <- var.random / (var.random + var.fixed + var.resid) return(c(as.numeric(icc.juge1),as.numeric(icc.juge2))) }  #SECOND FONCTION : bootstrap function, split on the hirarchical level as you want  nagakawa.boot <- function(data,x){ list.ICC <- split(x = data, f = paste(datajuge.L,data$$injection.L,sep = "_"))
list.BOOT <- lapply(X = list.ICC, FUN = function(y){
y[x,]
})
db.BOOT <- do.call(what = "rbind", args = list.BOOT)
nagakawa(dataICC = db.BOOT)
}


# THIRD : bootstrap execution

ICC.BOOT <- boot(data = dbICC, statistic = nagakawa.boot, R = 1000)