library(wildrwolf)

At the moment, the rwolf() function does not support many different types of families of tests. E.g. if you had estimated the following two regression models

$Y = \beta_0 + \beta_1 X_1 + \beta_2 X2 + u$ and

$Y = \beta_0 + \beta_3 X_3 + u$ and wanted to correct the following family of hypotheses for multiple testing:

$H_{0,A}: \beta_1 = 0 \text{ vs } H_{1,A}: \beta_1 \neq 0$ $H_{0,B}: \beta_2 = 0 \text{ vs } H_{1,B}: \beta_2 \neq 0$

$H_{0,C}: \beta_1 + \beta_2 = 0 \text{ vs } H_{1,C}: \beta_1 + \beta_2\neq 0$ $H_{0,D}: \beta_3 = 0 \text{ vs } H_{1,D}: \beta_3 \neq 0.$

Unfortunately, the current API of rwolf() does not yet support such a family of tests. To apply the Romano Wolf corrections for this family, you would have to follow three steps outlined below. But first, let’s simulate some data.

N <- 1000
X1 <- rnorm(N)
X2 <- rnorm(N)
X3 <- rnorm(N)
Y <- 1 + 0.001 * X1 + 0.001 * X2 + 0.5 *X3 + rnorm(N)
cluster <- sample(1:20, N, TRUE)
df <- data.frame(Y = Y, X1 = X1, X2 = X2, X3 = X3, cluster)

### Step 1: Estimate two regression models via fixest

library(fixest)

fit1 <- feols(Y ~ X1 + X2, data = df, cluster = ~cluster)
fit2 <- feols(Y ~ X3, data = df, cluster = ~cluster)

### Step 2: Create bootstrapped t-statistics via fwildclusterboot::boottest(). Make sure to reset the random seeds, so that all calls to boottest() use the same bootstrap weights

library(fwildclusterboot)

set.seed(123)
boot1 <- boottest(fit1, param = ~X1, B = 9999, clustid = ~cluster)
#> Warning: Please note that the seeding behavior for random number generation for
#> boottest() has changed with fwildclusterboot version 0.13.
#>
#> It will no longer be possible to exactly reproduce results produced by versions
#> lower than 0.13.
#>
#> If your prior results were produced under sufficiently many bootstrap
#> [news.md](https://cran.r-project.org/web/packages/fwildclusterboot/news/news.html).
#> This warning is displayed once per session.
#> Too guarantee reproducibility, don't forget to set a global random seed
#> **both** via set.seed() and dqrng::dqset.seed().
#> This message is displayed once every 8 hours.
set.seed(123)
boot2 <- boottest(fit1, param = ~X2, B = 9999, clustid = ~cluster)
set.seed(123)
boot3 <- boottest(fit1, param = ~X1 + X2, B = 9999, clustid = ~cluster)
set.seed(123)
boot4 <- boottest(fit2, param = ~ X3, B = 9999, clustid = ~cluster)

# get the bootstrapped t-stats from boottest
t_boot <- lapply(list(boot1, boot2, boot3, boot4),
function(x) x[["t_boot"]])
t_boot <- Reduce("cbind",t_boot)

# get the non-bootstrap t-stats from boottest
t_stat <-  lapply(list(boot1, boot2, boot3, boot4),
function(x) teststat(x))
t_stat <- Reduce("cbind",t_stat)

### Step 3: feed the bootstrapped and non-bootstrapped t-statistics into the get_rwolf_pval() function

get_rwolf_pval(t_stats = t_stat, boot_t_stats = t_boot)
#>  0.9957 0.9957 0.9957 0.0001

This returns a vector of multiple-testing adjusted pvalues for all hypotheses.