(Clustered) Bootstrap Covariance Matrix Estimation
Description
Objectoriented estimation of basic bootstrap covariances, using
simple (clustered) casebased resampling, plus more refined methods
for lm
and glm
models.
Usage
vcovBS(x, ...)
vcovBS(x, cluster = NULL, R = 250, start = FALSE, type = "xy", ...,
fix = FALSE, use = "pairwise.complete.obs", applyfun = NULL, cores = NULL,
center = "mean")
vcovBS(x, cluster = NULL, R = 250, type = "xy", ...,
fix = FALSE, use = "pairwise.complete.obs", applyfun = NULL, cores = NULL,
qrjoint = FALSE, center = "mean")
vcovBS(x, cluster = NULL, R = 250, start = FALSE, type = "xy", ...,
fix = FALSE, use = "pairwise.complete.obs", applyfun = NULL, cores = NULL,
center = "mean")
Arguments
x 
a fitted model object.

cluster 
a variable indicating the clustering of observations,
a list (or data.frame ) thereof, or a formula specifying
which variables from the fitted model should be used (see examples).
By default (cluster = NULL ), either attr(x, "cluster") is used
(if any) or otherwise every observation is assumed to be its own cluster.

R 
integer. Number of bootstrap replications.

start 
logical. Should coef(x) be passed as start
to the update(x, subset = ...) call? In case the model x
is computed by some numeric iteration, this may speed up the bootstrapping.

type 
character (or function). The character string specifies the type of
bootstrap to use: In the default and glm method the three types
"xy" , "fractional" , and "jackknife" are available.
In the lm method there are additionally "residual" , "wild"
(or equivalently: "wildrademacher" or "rademacher" ),
"mammen" (or "wildmammen" ), "norm"
(or "wildnorm" ), "webb" (or "wildwebb" ).
Finally, for the lm method type can be a function(n)
for drawing wild bootstrap factors.

... 
arguments passed to methods. For the default method, this is
passed to update , and for the lm method to lm.fit .

fix 
logical. Should the covariance matrix be fixed to be
positive semidefinite in case it is not?

use 
character. Specification passed to cov for
handling missing coefficients/parameters.

applyfun 
an optional lapply style function with arguments
function(X, FUN, ...) . It is used for refitting the model to the
bootstrap samples. The default is to use the basic lapply
function unless the cores argument is specified (see below).

cores 
numeric. If set to an integer the applyfun is set to
mclapply with the desired number of cores ,
except on Windows where parLapply with
makeCluster(cores) is used.

center 
character. For type = "jackknife" the coefficients from
all jacknife samples (each dropping one observational unit/cluster) can be
centered by their "mean" (default) or by the original fullsample
"estimate" .

qrjoint 
logical. For residualbased and wild boostrap (i.e.,
type != "xy" ), should the bootstrap sample the dependent variable
and then apply the QR decomposition jointly only once? If FALSE ,
the boostrap applies the QR decomposition separately in each iteration
and samples coefficients directly. If the sample size (and the number of
coefficients) is large, then qrjoint = TRUE maybe significantly
faster while requiring much more memory.

Details
Clustered sandwich estimators are used to adjust inference when errors
are correlated within (but not between) clusters. See the documentation for vcovCL
for specifics about covariance clustering. This function allows
for clustering in arbitrarily many cluster dimensions (e.g., firm, time, industry), given all
dimensions have enough clusters (for more details, see Cameron et al. 2011).
Unlike vcovCL
, vcovBS
uses a bootstrap rather than an asymptotic solution.
Basic (clustered) bootstrap covariance matrix estimation is provided by
the default vcovBS
method. It samples clusters (where each observation
is its own cluster by default), i.e., using casebased resampling. For obtaining
a covariance matrix estimate it is assumed that an update
of the model with the resampled subset
can be obtained, the coef
extracted, and finally the covariance computed with cov
.
The update
model is evaluated in the environment(terms(x))
(if available).
To speed up computations two further arguments can be leveraged.

Instead of lapply
a parallelized function such as
parLapply
or mclapply
can be specified to iterate over the bootstrap replications. For the latter,
specifying cores = ...
is a convenience shortcut.

When specifying start = TRUE
, the coef(x)
are passed to
update
as start = coef(x)
. This may not be supported by all
model fitting functions and is hence not turned on by default.
The “xy” or “pairs” bootstrap is consistent for heteroscedasticity and clustered errors,
and converges to the asymptotic solution used in vcovCL
as R
, $n$
, and $g$
become large ($n$
and $g$
are the number of
observations and the number of clusters, respectively; see Efron 1979, or Mammen 1992, for a
discussion of bootstrap asymptotics). For small $g$
– particularly under 30 groups – the
bootstrap will converge to a slightly different value than the asymptotic method, due to
the limited number of distinct bootstrap replications possible (see Webb 2014 for a discussion
of this phenomonon). The bootstrap will not necessarily converge to an asymptotic estimate
that has been corrected for small samples.
The xy approach to bootstrapping is generally only of interest to the
practitioner when the asymptotic solution is unavailable (this can happen when using
estimators that have no estfun
function, for example). The residual bootstrap,
by contrast, is rarely of practical interest, because while it provides consistent
inference for clustered standard errors, it is not robust to heteroscedasticity.
More generally, bootstrapping is useful when the bootstrap makes different assumptions than the asymptotic
estimator, in particular when the number of clusters is small and large $n$
or
$g$
assumptions are unreasonable. Bootstrapping is also often effective for nonlinear models,
particularly in smaller samples, where asymptotic approaches often perform relatively poorly.
See Cameron and Miller (2015) for further discussion of bootstrap techniques in practical applications,
and Zeileis et al. (2020) show simulations comparing vcovBS
to vcovCL
in several
settings.
The jackknife approach is of particular interest in practice because it can be shown to be
exactly equivalent to the HC3 (without cluster adjustment, also known as CV3)
covariance matrix estimator in linear models (see MacKinnon,
Nielsen, Webb 2022). If the number of observations per cluster is large it may become
impossible to compute this estimator via vcovCL
while using the jackknife
approach will still be feasible. In nonlinear models (including nonGaussian GLMs) the
jackknife and the HC3 estimator do not coincide but the jackknife might still be a useful
alternative when the HC3 cannot be computed. A convenience interface vcovJK
is provided whose default method simply calls vcovBS(..., type = "jackknife")
.
The fractionalrandomweight bootstrap (see Xu et al. 2020), first introduced by
Rubin (1981) as Bayesian bootstrap, is an alternative to the xy bootstrap when it is
computationally challenging or even impractical to reestimate the model on subsets, e.g.,
when "successes" in binary responses are rare or when the number of parameters is close
to the sample size. In these situations excluding some observations completely is the
source of the problems, i.e., giving some observations zero weight while others receive
integer weights of one ore more. The fractional bootstrap mitigates this by giving
every observation a positive fractional weight, drawn from a Dirichlet distribution.
These may become close to zero but never exclude an observation completly, thus stabilizing
the computation of the reweighted models.
The glm
method works essentially like the default method but calls
glm.fit
instead of update
.
The lm
method provides additional bootstrapping type
s
and computes the bootstrapped coefficient estimates somewhat more efficiently using
lm.fit
(for casebased resampling) or qr.coef
rather than update
. The default type
is casebased resampling
(type = "xy"
) as in the default method. Alternative type
specifications are:

"residual"
. The residual cluster bootstrap resamples the residuals (as above,
by cluster) which are subsequently added to the fitted values to obtain the bootstrapped
response variable: $y^{*} = \hat{y} + e^{*}$
.
Coefficients can then be estimated using qr.coef()
, reusing the
QR decomposition from the original fit. As Cameron et al. (2008) point out,
the residual cluster bootstrap is not welldefined when the clusters are unbalanced as
residuals from one cluster cannot be easily assigned to another cluster with different size.
Hence a warning is issued in that case.

"wild"
(or equivalently "wildrademacher"
or "rademacher"
).
The wild cluster bootstrap does not actually resample the residuals but instead reforms the
dependent variable by multiplying the residual by a randomly drawn value and adding the
result to the fitted value: $y^{*} = \hat{y} + e \cdot w$
(see Cameron et al. 2008). By default, the factors are drawn from the Rademacher distribution:
function(n) sample(c(1, 1), n, replace = TRUE)
.

"mammen"
(or "wildmammen"
). This draws the wild bootstrap factors as
suggested by Mammen (1993):
sample(c(1, 1) * (sqrt(5) + c(1, 1))/2, n, replace = TRUE, prob = (sqrt(5) + c(1, 1))/(2 * sqrt(5)))
.

"webb"
(or "wildwebb"
). This implements the sixpoint distribution
suggested by Webb (2014), which may improve inference when the number of clusters is small:
sample(c(sqrt((3:1)/2), sqrt((1:3)/2)), n, replace = TRUE)
.

"norm"
(or "wildnorm"
). The standard normal/Gaussian distribution
is used for drawing the wild bootstrap factors: function(n) rnorm(n)
.

Userdefined function. This needs of the form as above, i.e., a function(n)
returning a vector of random wild bootstrap factors of corresponding length.
Value
A matrix containing the covariance matrix estimate.
References
Cameron AC, Gelbach JB, Miller DL (2008).
“BootstrapBased Improvements for Inference with Clustered Errors”,
The Review of Economics and Statistics, 90(3), 414–427.
doi:10.3386/t0344
Cameron AC, Gelbach JB, Miller DL (2011).
“Robust Inference with Multiway Clustering”,
Journal of Business & Economic Statistics, 29(2), 238–249.
doi:10.1198/jbes.2010.07136
Cameron AC, Miller DL (2015).
“A Practitioner's Guide to ClusterRobust Inference”,
Journal of Human Resources, 50(2), 317–372.
doi:10.3368/jhr.50.2.317
Efron B (1979).
“Bootstrap Methods: Another Look at the Jackknife”,
The Annals of Statistics, 7(1), 1–26.
doi:10.1214/aos/1176344552
MacKinnon JG, Nielsen MØ, Webb MD (2022).
“ClusterRobust Inference: A Guide to Empirical Practice”,
Journal of Econometrics, Forthcoming.
doi:10.1016/j.jeconom.2022.04.001
Mammen E (1992).
“When Does Bootstrap Work?: Asymptotic Results and Simulations”,
Lecture Notes in Statistics, 77.
Springer Science & Business Media.
Mammen E (1993).
“Bootstrap and Wild Bootstrap for High Dimensional Linear Models”,
The Annals of Statistics, 21(1), 255–285.
doi:10.1214/aos/1176349025
Rubin DB (1981).
“The Bayesian Bootstrap”,
The Annals of Statistics, 9(1), 130–134.
doi:10.1214/aos/1176345338
Webb MD (2014).
“Reworking Wild Bootstrap Based Inference for Clustered Errors”,
Working Paper 1315, Queen's Economics Department.
https://www.econ.queensu.ca/sites/econ.queensu.ca/files/qed_wp_1315.pdf.
Xu L, Gotwalt C, Hong Y, King CB, Meeker WQ (2020).
“Applications of the FractionalRandomWeight Bootstrap”,
The American Statistician, 74(4), 345–358.
doi:10.1080/00031305.2020.1731599
Zeileis A, Köll S, Graham N (2020).
“Various Versatile Variances: An ObjectOriented Implementation of Clustered Covariances in R.”
Journal of Statistical Software, 95(1), 1–36.
doi:10.18637/jss.v095.i01
See Also
vcovCL
, vcovJK
Examples
data("PetersenCL", package = "sandwich")
m < lm(y ~ x, data = PetersenCL)
suppressWarnings(RNGversion("3.5.0"))
set.seed(1)
cbind(
"classical" = sqrt(diag(vcov(m))),
"HCcluster" = sqrt(diag(vcovCL(m, cluster = ~ firm))),
"BScluster" = sqrt(diag(vcovBS(m, cluster = ~ firm))),
"FWcluster" = sqrt(diag(vcovBS(m, cluster = ~ firm, type = "fractional")))
)
vcovBS(m, cluster = ~ firm + year, type = "wildmammen")
all.equal(
vcovBS(m, cluster = ~ firm, type = "jackknife"),
vcovCL(m, cluster = ~ firm, type = "HC3", cadjust = FALSE),
tolerance = 1e7
)