Bootstrapping hierarchical/multilevel data (resampling clusters)

Question

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

Following the Davidson and Hinkley textbook [1] 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?

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

Answer 1

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.

Answer 2

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

Answer 3

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")
labs = unique(dat$lab)
b.sample <- replicate(B,
  dat[
    .(sample(labs, replace = TRUE)),
    glm(y ~ x)$coefficients
  ]
)

Answer 4

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.

Answer 5

See there for an answer that you may find useful: Cluster Boostrap with Unequally Sized Clusters It basically proposes code to make bootstrap inference on a custom statistic (not only from a fitted model), with unbalanced clusters, and possibly fixed effects in the same dimension as the clusters.

If your clusters are of equal size (balanced), the code below will be faster:

library(boot)
library(sandwich)

## Make some necessary objects
# unbalanced panel data
data("PetersenCL", package = "sandwich")
data <- PetersenCL

# list of parameters related to the dataset, and the clustering variable
cluster_var <- "firm"

# names and numbers of clusters of size s
par_list <- list(cluster_variable = cluster_var, 
                 cluster_names = unique(data[,cluster_var]),
                 number_clusters = length(unique(data[,cluster_var])))

ran.gen_cluster_blc <- function(original_data, arg_list){
  # to store 
  cl_boot_dat <- list()
  
  # such that we don't have to call it from the list every time
  cluster_var <- arg_list[["cluster_variable"]]
  
  # non-unique names of clusters (repeated when there is more than one obs. in a cluster) 
  nu_cl_names <- as.character(original_data[,cluster_var]) 
  
  # sample, in the vector of names of clusters as many draws as there are clusters, with replacement
  sample_cl <- sample(arg_list[["cluster_names"]], 
                      arg_list[["number_clusters"]], 
                      replace = TRUE) 
  
  # because of replacement, some names are sampled more than once
  # we need to give them a new cluster identifier, otherwise a cluster sampled more than once 
  # will be "incorrectly treated as one large cluster rather than two distinct clusters" (by the fixed effects) (Cameron and Miller, 2015)    
  sample_cl_tab <- table(sample_cl)
  
  for(n in 1:max(sample_cl_tab)){ # from 1 to the max number of times a cluster was sampled bc of replacement
    # vector to select obs. that are within the clusters sampled n times. 
    # seems slightly faster to construct the names_n object beforehand 
    names_n <- names(sample_cl_tab[sample_cl_tab == n])
    sel <- nu_cl_names %in% names_n

    # select data accordingly to the cluster sampling (duplicating n times observations from clusters sampled n times)
    clda <- original_data[sel,][rep(seq_len(sum(sel)), n), ]
    
    #identify row names without periods, and add ".0" 
    row.names(clda)[grep("\\.", row.names(clda), invert = TRUE)] <- paste0(grep("\\.", row.names(clda), invert = TRUE, value = TRUE),".0")
    
    # add the suffix due to the repetition after the existing cluster identifier. 
    clda[,cluster_var] <- paste0(clda[,cluster_var], sub(".*\\.","_",row.names(clda)))
    
    # stack the bootstrap samples iteratively 
    cl_boot_dat[[n]] <- clda
  }
  return(bind_rows(cl_boot_dat))
}

#test that the returned data are the same dimension as input
test_boot_d <- ran.gen_cluster_blc(original_data = data,
                                   arg_list = par_list)

dim(test_boot_d)
dim(data)
# test new clusters are not duplicated (correct if anyDuplicated returns 0)
base::anyDuplicated(test_boot_d[,c("firm","year")])


# test that it computes the same standard error as sandwich::vcovBS, for statistic being a regression coefficient 
est_fun <- function(est_data){
  
  est <- lm(as.formula("y ~ x"), est_data)
  
  # statistics we want to evaluate the variance of:
  return(est$coefficients)
}
# see ?boot::boot for more details on these arguments
set.seed(1234)
boot(data = data, 
     statistic = est_fun, 
     ran.gen = ran.gen_cluster_blc,
     mle = par_list,
     sim = "parametric",
     parallel = "no",
     R = 400)

set.seed(1234)
sdw_bs <- vcovBS(lm(as.formula("y ~ x"), PetersenCL), cluster = ~firm, R=400)#
sqrt(sdw_bs["x","x"])

sdw_cl <- vcovCL(lm(as.formula("y ~ x"), data), cluster = ~firm)
sqrt(sdw_cl["x","x"])

Answer 6

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[1]
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[1])^2 + (attributes(variance$sujet.L)$stddev[2])^2
    var.sujet <- (attributes(variance$sujet.L)$stddev[1])^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(data$juge.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)