Non-Standard Families of Tests

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
#> iterations, none of your conclusions will change.  For more details about this
#> change, please read the notes in
#> [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)
#> [1] 0.5454 0.4334 0.4334 0.0001

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