Cluster Boostrap with Unequally Sized Clusters

Question

I need to perform a bootstrap for variance estimation on a GEE model for clustered data that I am analyzing. I understand that I need to use a clustered bootstrap for this, which is pretty much the same thing as the usual nonparametric bootstrap, only the clusters are sampled rather than the individual observations within a cluster. I've read a few articles on this and the articles always assume clusters of equal size. Can the clustered bootstrap be modified for clusters of unequal size? If so, how is this done? In the usual bootstrap, resampling with replacement is performed until the dataset is the same size of the original dataset. How do we perform a cluster bootstrap with unequal sized clusters such that the $i$th bootstrap sample is the same size as the original dataset? I could imagine a situation where the cluster sizes are different such that it might be impossible to obtain datasets of the same size as the original dataset if we sample clusters, stead of elements.

Lastly is there a procedure in SAS or R that will perform the cluster bootstrap?

Answer 1

This is explained quite nicely in Sherman and leCessie's paper, "A comparison between bootstrap methods and generalized estimating equations for correlated outcomes in generlized linear models." On page 905, they note:

"If as often may be the case, there are blocks of different sizes, then the algorithm can be modified as follows: let $m_i$, denote the number of blocks of size $i$, $i = 1,. . . ,I$. For each $i$ resample $m_i$ times with replacement from the $m_i$ blocks and compute $\hat{\beta}^*$ from the $n = \sum_{i=1}^Iim_i$ resampled observations. This conditioning on block size guarantees a total resarnple size equal to the original sample size, making the bootstrap replicates "comparable". If, however, $I$ is large it may be more attractive to resample $m$ times from the entire set of blocks. We will call this the "All Block" bootstrap. This algorithm gives a random total sample size, $n^*$, say, which makes the replicates less comparable. A reasonable approach to make them more comparable is to let the replicate be $(n^*/n)^{1/2}\hat{\beta}^*$, as suggested by Efron and Tibshirani (1993, p.101) for a whole block bootstrap in the time series setting."

REFERENCES:

Michael Sherman & Saskia le Cessie (1997) A comparison between bootstrap methods and generalized estimating equations for correlated outcomes in generalized linear models, Communications in Statistics - Simulation and Computation, 26:3, 901-925, DOI: 10.1080/03610919708813417

Answer 2

I wrote something in R for my own use, based on the quote from Sherman and Cessie (1997) in StatsStudent's answer.

It implements bootstrap replicates on clustered data with clusters of different sizes.

It makes sure that clusters sampled more than once (due to replacement) are treated as distinct clusters within bootstrap samples (especially important in estimations involving fixed-effects along the cluster dimension).

Finally, note that this approach allows to compute the one-way cluster bootstrap standard error of any statistic computed via a custom, user-provided function. In this respect, it is more flexible than sandwich::vcovBS or sandwich::vcovCL which accept only some fitted model objects as input.

EDIT (improvement)

The Sherman and Cessie (1997) approach that draws clusters from pools of same-size clusters is incorrect I think. More precisely, the same-size-cluster pools are systematically smaller than the full pool of clusters. Therefore, this approach produces standard errors biased towards zero. The most extreme case being when every cluster is of different size (e.g. clustering on year level and no cross-section has the same size). In this case, the same cluster gets sampled for every size, all bootstrap replicates are thus identical and the standard error is exactly zero. I leave my old code implementing this approach at the end of my answer.

The better solution is to draw in the full pool of clusters such that the replicate data set has, on average, the same size as the original data set.

library(boot)
library(sandwich)

data("PetersenCL", package = "sandwich")

# construct an unbalanced panel with different numbers of firms every year. 
datalist <- list()
for(yr in 1:7){
  datalist[[yr]] <- subset(PetersenCL, (firm %in% c(1:(2*yr)) & year == yr))
}
data <- bind_rows(datalist)

cluster_var <- "year"

cluster_names = unique(data[,cluster_var])

# get the different cluster sizeS. This is necessary to cluster bootstrapping with clusters of different sizes. 
sizes <- table(data[,cluster_var])

# don't bother the size of every cluster. Use the average cluster size only. 
avg_cl_size <- mean(sizes)
# sample size of original data
N <- nrow(data)
# Number of draws, for the bootstrap replicate data sets to be as large as original data, on average. 
n_draws <- N/avg_cl_size

ran.gen_cluster <- function(original_data, arg_list){
  
  cl_boot_dat <- NULL
  
  # for later: 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, N/avg_cl_size
  sample_cl <- sample(x = cluster_names, 
                      size = n_draws, 
                      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 name was sampled bc of replacement
    # vector to select obs. that are within the sampled clusters. 
    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 <- rbind(cl_boot_dat, clda) 
  }
  return(cl_boot_dat)
}

#the returned data ARE NOT the same dimension as input data. It only has the same dimension on average over bootstrap replicates. 
test_boot_d <- ran.gen_cluster(original_data = data)
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")])

# custom estimation function
est_fun <- function(est_data){
  
  est <- lm(as.formula("y ~ x"), est_data)
  
  # statistics we want to evaluate the variance of:
  return(est$coefficients)
}

# Run the bootstrap (see ?boot::boot for more details on the arguments)
set.seed(1234)
boot(data = data, 
     statistic = est_fun, 
     ran.gen = ran.gen_cluster,
     mle = list(), # this argument cannot be left empty 
     sim = "parametric",
     parallel = "no",
     R = 400)

Test that it computes the same standard error as sandwich::vcovBS and compare with the asymptotic solution of sandwich::vcovCL

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

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

OLDER, FLAWED APPROACH

library(boot)
library(sandwich)

## Make some necessary objects
# unbalanced panel data
data("PetersenCL", package = "sandwich")
data <- subset(PetersenCL, !(firm %in% c(1:10) & year == 10))

cluster_var <- "firm"

# get the different cluster sizeS. This is necessary to cluster bootstrapping with clusters of different sizes. 
sizes <- table(data[,cluster_var])
u_sizes <- sort(unique(sizes))

# names and numbers of clusters of every sizes
cl_names <- list()
n_clusters <- list()
for(s in u_sizes){
  cl_names[[s]] <- names(sizes[sizes == s])
  n_clusters[[s]] <- length(cl_names[[s]])
}

par_list <- list(unique_sizes = u_sizes,
                 cluster_names = cl_names,
                 number_clusters = n_clusters)

## Design the bootstrap sampling function. 
# It will tell boot::boot how to sample data at each replicate 
    ran.gen_cluster_bysize <- function(original_data, arg_list){
  
  cl_boot_dat <- NULL

  # 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]) 
  
  for(s in arg_list[["unique_sizes"]]){
    # sample, in the vector of names of clusters of size s, as many draws as there are clusters of that size, with replacement
    sample_cl_s <- sample(arg_list[["cluster_names"]][[s]], 
                          arg_list[["number_clusters"]][[s]], 
                          replace = TRUE) 
    
    # because of replacement, some names are sampled more than once
    sample_cl_s_tab <- table(sample_cl_s)
    
    # 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 2015)    
    
    for(n in 1:max(sample_cl_s_tab)){ # from 1 to the max number of times a name was sampled bc of replacement
      # vector to select obs. that are within the sampled clusters. 
      names_n <- names(sample_cl_s_tab[sample_cl_s_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 <- rbind(cl_boot_dat, clda) 
    }
  }  
  return(cl_boot_dat)
}

#test that the returned data are the same dimension as input
test_boot_d <- ran.gen_cluster_bysize(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")])


## Custom ("black-box") estimation function
est_fun <- function(est_data){
  
  est <- lm(as.formula("y ~ x"), est_data)
  
  # statistics we want to evaluate the variance of:
  return(est$coefficients[2])
}

# see ?boot::boot for more details on these arguments
boot(data = data, 
      statistic = est_fun, 
      ran.gen = ran.gen_cluster_bysize,
      mle = par_list,
      sim = "parametric",
      parallel = "no",
      R = 200)

With this clustering structure, the result is close to the asymptotic solution (available for unbalanced data, but only for lm or glm objects):

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

But it does not replicate the result from vcovBS, even for balanced data:

## compare with sandwich::vcovBS, on balanced panel data
set.seed(1234)
custom_bs <- boot(data = PetersenCL, 
                  statistic = est_fun, 
                  ran.gen = ran.gen_cluster_bysize,
                  mle = par_list,
                  sim = "parametric",
                  parallel = "no",
                  R = 200)
sd(custom_bs$t[,2])

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