The boottest() function
The fwildclusterboot package consists of one key
function, boottest(). It implements the fast wild bootstrap
and works with regression objects of type lm,
felm and fixest from base R and the
lfe and fixest packages.
To start, we create a random data set with two cluster variables
(group_id1 & group_id2), two fixed effects and a set of covariates.
The icc_ arguments control the cluster variable’s
intra-cluster correlation.
# set seed via dqset.seed for boot_algo = "R" & Rademacher, Webb & Normal weights
dqrng::dqset.seed(2352342)
# set 'familiar' seed for all other algorithms and weight types
set.seed(23325)
# load data set voters included in fwildclusterboot
data(voters)
# estimate the regression model via lm
lm_fit <- lm(proposition_vote ~ treatment + ideology1 + log_income + Q1_immigration , data = voters)
# model with interaction
lm_fit_interact <- lm(proposition_vote ~ treatment + ideology1 + log_income:Q1_immigration , data = voters)The boottest() function has 4 required and several
optional arguments. The required objects are
- object: a regression object of type
lm,fixestorfelm - clustid: a character vector that defines the clustering variables
- param: a character vector of length one - the model parameter to be tested
- B: the number of bootstrap iterations
# boottest on an object of type lm
boot_lm <- boottest(lm_fit, clustid = "group_id1", param = "treatment", B = 9999)To tests for an interaction, it is important to use the coefficient names that are internally created by the modeling function.
#names(coef(lm_fit_interact))
boot_lm_interact <- boottest(lm_fit_interact, clustid = "group_id1", param = "log_income:Q1_immigration1", B = 9999)boottest() further allows for multivariable tests.
Suppose we’re interested in testing the null hypothesis \(0.6*treatment + 0.2*ideology1 = 0.02\). To
test such a hypothesis, one would have to specify the hypothesis via the
param, R and r arguments:
boot_multi <- boottest(lm_fit, clustid = "group_id1", param = c("treatment", "ideology1"), R = c(0.6, 0.2), r = 0.02, B = 9999)To access the estimation results, boottest() comes with
summary(), tidy() and glance()
methods. The tidy() method returns the estimation results
in a data.frame. summary() returns additional information
on top of the test statistics reported by tidy().
Theglance() method enables the use of output formatting
tools from the modelsummary package.
# fwildclusterboot's internal summary() method
summary(boot_lm)
#> boottest.lm(object = lm_fit, param = "treatment", B = 9999, clustid = "group_id1")
#>
#> Hypothesis: 1*treatment = 0
#> Observations: 300
#> Bootstr. Type: rademacher
#> Clustering: 1-way
#> Confidence Sets: 95%
#> Number of Clusters: 40
#>
#> term estimate statistic p.value conf.low conf.high
#> 1 1*treatment = 0 0.073 3.709 0.001 0.033 0.113
summary(boot_multi)
#> boottest.lm(object = lm_fit, param = c("treatment", "ideology1"),
#> B = 9999, clustid = "group_id1", R = c(0.6, 0.2), r = 0.02)
#>
#> Hypothesis: 0.6*treatment+0.2*ideology1 = 0.02
#> Observations: 300
#> Bootstr. Type: rademacher
#> Clustering: 1-way
#> Confidence Sets: 95%
#> Number of Clusters: 40
#>
#> term estimate statistic p.value conf.low
#> 1 0.6*treatment+0.2*ideology1 = 0.02 0.048 2.346 0.025 0.024
#> conf.high
#> 1 0.072
if(requireNamespace("modelsummary")){
# summary via the modelsummary package
library(modelsummary)
msummary(list(boot_lm, boot_lm_interact),
estimate = "{estimate} ({p.value})",
statistic = "[{conf.low}, {conf.high}]")
}
#> Loading required namespace: modelsummary| Model 1 | Model 2 | |
|---|---|---|
| 1*treatment = 0 | 0.073 (0.001) | |
| [0.033, 0.113] | ||
| 1*log_income × Q1_immigration1 = 0 | -0.038 (0.000) | |
| [-0.057, -0.019] | ||
| Num.Obs. | 300 | 300 |
| R2 | 0.316 | 0.339 |
| R2 Adj. | 0.288 | 0.311 |
| AIC | -82.1 | -92.2 |
| BIC | -30.2 | -40.4 |
| Log.Lik. | 55.025 | 60.102 |
A plot() method allows the user to inspect the bootstrap
t-statistics:
plot(boot_lm)
Multiway Clustering
The boottest() function supports clustering of any
dimension. E.g. for two-way clustering, one simply needs to specify the
names of the cluster variables in a character vector.
boot_lm <- boottest(lm_fit, clustid = c("group_id1", "group_id2"), param = "treatment", B = 9999)
summary(boot_lm)
#> boottest.lm(object = lm_fit, param = "treatment", B = 9999, clustid = c("group_id1",
#> "group_id2"))
#>
#> Hypothesis: 1*treatment = 0
#> Observations: 300
#> Bootstr. Type: rademacher
#> Clustering: 2-way
#> Confidence Sets: 95%
#> Number of Clusters: 40 20 253
#>
#> term estimate statistic p.value conf.low conf.high
#> 1 1*treatment = 0 0.073 3.842 0.004 0.031 0.115The Heteroskedastic Bootstrap
If you drop the clustid argument,
boottest() will run a heteroskedasticity robust wild
bootstrap via the ‘R-lean’ algorithm. At the moment, the null hypothesis
is always imposed, only Rademacher and Webb weights are supported, and
no confidence intervals are computed. Further, no regression weights are
supported. As all algorithms in fwildclusterboot, p-values
are calculated based on pivotal t-statistics.
boot_lm <- boottest(lm_fit, param = "treatment", B = 9999)
summary(boot_lm)
#> boottest.lm(object = lm_fit, param = "treatment", B = 9999)
#>
#> Hypothesis: 1*treatment = 0
#> Observations: 300
#> Bootstr. Type: rademacher
#> Clustering: 0-way
#> Confidence Sets: 95%
#>
#> term estimate statistic p.value conf.low conf.high
#> 1 1*treatment = 0 0.073 3.096 0.002 NA NA
boot_lm$boot_algo
#> [1] "R-lean"Choice of Bootstrap Weights
Furthermore, the user can choose among four different weighting
distribution via the type argument: Rademacher, Mammen,
Normal and Webb. By default, boottest() uses the Rademacher
distribution.
boot_lm_rade <- boottest(lm_fit,
clustid = c("group_id1", "group_id2"),
param = "treatment",
B = 999,
type = "rademacher")
boot_lm_webb <- boottest(lm_fit,
clustid = c("group_id1", "group_id2"),
param = "treatment",
B = 999,
type = "webb")
if(requireNamespace("modelsummary")){
library(modelsummary)
msummary(list(boot_lm_rade, boot_lm_webb),
estimate = "{estimate} ({p.value})",
statistic = "[{conf.low}, {conf.high}]")
}| Model 1 | Model 2 | |
|---|---|---|
| 1*treatment = 0 | 0.073 (0.006) | 0.073 (0.006) |
| [0.031, 0.115] | [0.029, 0.117] | |
| Num.Obs. | 300 | 300 |
| R2 | 0.316 | 0.316 |
| R2 Adj. | 0.288 | 0.288 |
| AIC | -82.1 | -82.1 |
| BIC | -30.2 | -30.2 |
| Log.Lik. | 55.025 | 55.025 |
Other function arguments
Via the function argument sign_level, the user can
control the significance level of the test. The default value is
sign_level = 0.05, which corresponds to a 95% confindence interval.
boot_lm_5 <- boottest(lm_fit,
clustid = c("group_id1"),
param = "treatment", B = 9999,
sign_level = 0.05)
boot_lm_10 <- boottest(lm_fit,
clustid = c("group_id1"),
param = "treatment", B = 9999,
sign_level = 0.10)
if(requireNamespace("modelsummary")){
library(modelsummary)
msummary(list(boot_lm_5, boot_lm_10),
estimate = "{estimate} ({p.value})",
statistic = "[{conf.low}, {conf.high}]")
}| Model 1 | Model 2 | |
|---|---|---|
| 1*treatment = 0 | 0.073 (0.000) | 0.073 (0.001) |
| [0.033, 0.113] | [0.040, 0.106] | |
| Num.Obs. | 300 | 300 |
| R2 | 0.316 | 0.316 |
| R2 Adj. | 0.288 | 0.288 |
| AIC | -82.1 | -82.1 |
| BIC | -30.2 | -30.2 |
| Log.Lik. | 55.025 | 55.025 |
In the case of multiway clustering, the user might want to specify
the bootstrap clustering level. By default, boottest chooses the
clustering level with the highest number of clusters as
bootcluster = "max". Other choices are the minimum cluster,
or independent clustering variables.
boot_lm1 <- boottest(lm_fit,
clustid = c("group_id1", "group_id2"),
param = "treatment",
B = 9999,
bootcluster = "min")
boot_lm2 <- boottest(lm_fit,
clustid = c("group_id1", "group_id2"),
param = "treatment",
B = 9999,
bootcluster = "group_id1")
if(requireNamespace("modelsummary")){
library(modelsummary)
msummary(list(boot_lm1, boot_lm2),
estimate = "{estimate} ({p.value})",
statistic = "[{conf.low}, {conf.high}]")
}| Model 1 | Model 2 | |
|---|---|---|
| 1*treatment = 0 | 0.073 (0.005) | 0.073 (0.008) |
| [0.031, 0.111] | [0.027, 0.115] | |
| Num.Obs. | 300 | 300 |
| R2 | 0.316 | 0.316 |
| R2 Adj. | 0.288 | 0.288 |
| AIC | -82.1 | -82.1 |
| BIC | -30.2 | -30.2 |
| Log.Lik. | 55.025 | 55.025 |
Fixed Effects
Last, boottest() supports out-projection of fixed
effects in the estimation stage via lfe::felm() and
fixest::feols(). Within the bootstrap, the user can choose
to project out only one fixed effect, which can be set via the
fe function argument. All other fixed effects specified in
either felm() or feols() are treated as sets
of binary regressors.
if(requireNamespace("fixest")){
# estimate the regression model via feols
library(fixest)
feols_fit <- feols(proposition_vote ~ treatment + ideology1 + log_income | Q1_immigration , data = voters)
boot_feols <- boottest(feols_fit,
clustid = "group_id1",
param = "treatment",
B = 9999,
fe = "Q1_immigration")
}
#> Loading required namespace: fixestThe Subcluster Bootstrap
In the case of few treated clusters, MacKinnon
and Webb (2018) suggest to use subclusters to form the bootstrap
distribution. boottest() allows the user to specify
subclusters via the bootcluster argument.
boot_min <- boottest(lm_fit,
clustid = c("group_id1", "group_id2"),
param = "treatment",
B = 9999,
bootcluster = "min")
boot_var <- boottest(lm_fit,
clustid = c("group_id1", "group_id2"),
param = "treatment",
B = 9999,
bootcluster = "group_id1")
boot_2var <- boottest(lm_fit,
clustid = c("group_id1", "group_id2"),
param = "treatment",
B = 9999,
bootcluster = c("group_id1", "Q1_immigration"))
if(requireNamespace("modelsummary")){
library(modelsummary)
msummary(model = list(boot_min, boot_2var),
estimate = "{estimate} ({p.value})",
statistic = "[{conf.low}, {conf.high}]")
}| Model 1 | Model 2 | |
|---|---|---|
| 1*treatment = 0 | 0.073 (0.003) | 0.073 (0.010) |
| [0.032, 0.110] | [0.027, 0.119] | |
| Num.Obs. | 300 | 300 |
| R2 | 0.316 | 0.316 |
| R2 Adj. | 0.288 | 0.288 |
| AIC | -82.1 | -82.1 |
| BIC | -30.2 | -30.2 |
| Log.Lik. | 55.025 | 55.025 |
Regression Weights / Weighted Least Squares (WLS)
If regression weights are specified in the estimation stage via
lm(), feols() or felm(),
boottest() incorporates the weights into the bootstrap
inference:
# regression with weights / WLS
lm_w_fit <- lm(proposition_vote ~ treatment + ideology1 + log_income, weights = voters$weights, data = voters)
boot_w_lm <- boottest(lm_w_fit,
clustid = "group_id1",
param = "treatment",
B = 9999)Parallel Execution
A major bottleneck for the performance of boottest() is
a large matrix multiplication, which includes the bootstrap weights
matrix on the right. In order to speed up the computation, this
multiplication calls the c++ Eigen library, which allows for
parallelization of dense matrix products. By default,
boottest() uses one thread. Note that there is a cost of
parallelization due to communication overhead. As a rule of thumb, if
boottest() takes more than 10 seconds per execution, using
a second thread might speed up the bootstrap.
The number of threads can be specified via the nthreads
argument of boottest():
boot_lm <- boottest(lm_fit,
clustid = "group_id1",
param = "treatment",
B = 9999,
nthreads = 2)