In s3alfisc/rwolf: Fast Computation of Romano-Wolf Corrected p-Values for Linear Regression Models

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

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)
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)

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

s3alfisc/rwolf documentation built on Jan. 30, 2024, 10:42 a.m.