The wildrwolf package implements Romano-Wolf
multiple-hypothesis-adjusted p-values for objects of type fixest and
fixest_multi from the fixest package via a wild (cluster) bootstrap.
Because the bootstrap-resampling is based on the
fwildclusterboot
package, wildrwolf is usually really fast.
The package is complementary to wildwyoung (still work in progress), which implements the multiple hypothesis adjustment method following Westfall and Young (1993).
Adding support for multi-way clustering is work in progress.
You can install the package from CRAN and the development version from GitHub with:
install.packages("wildrwolf")
# install.packages("devtools")
devtools::install_github("s3alfisc/wildrwolf")
# from r-universe (windows & mac, compiled R > 4.0 required)
install.packages('wildrwolf', repos ='https://s3alfisc.r-universe.dev')library(wildrwolf)
library(fixest)
set.seed(1412)
N <- 1000
X1 <- rnorm(N)
X2 <- rnorm(N)
rho <- 0.5
sigma <- matrix(rho, 4, 4); diag(sigma) <- 1
u <- MASS::mvrnorm(n = N, mu = rep(0, 4), Sigma = sigma)
Y1 <- 1 + 1 * X1 + X2
Y2 <- 1 + 0.01 * X1 + X2
Y3 <- 1 + 0.4 * X1 + X2
Y4 <- 1 + -0.02 * X1 + X2
for(x in 1:4){
var_char <- paste0("Y", x)
assign(var_char, get(var_char) + u[,x])
}
data <- data.frame(Y1 = Y1,
Y2 = Y2,
Y3 = Y3,
Y4 = Y4,
X1 = X1,
X2 = X2,
#group_id = group_id,
splitvar = sample(1:2, N, TRUE))
fit <- feols(c(Y1, Y2, Y3, Y4) ~ csw(X1,X2),
data = data,
se = "hetero",
ssc = ssc(cluster.adj = TRUE))
# clean workspace except for res & data
rm(list= ls()[!(ls() %in% c('fit','data'))])
res_rwolf1 <- wildrwolf::rwolf(
models = fit,
param = "X1",
B = 9999
)
#> | | | 0% | |========= | 12% | |================== | 25% | |========================== | 38% | |=================================== | 50% | |============================================ | 62% | |==================================================== | 75% | |============================================================= | 88% | |======================================================================| 100%
pvals <- lapply(fit, function(x) pvalue(x)["X1"]) |> unlist()
# Romano-Wolf Corrected P-values
res_rwolf1
#> model Estimate Std. Error t value Pr(>|t|) RW Pr(>|t|)
#> 1 1 0.9896609 0.04204902 23.53588 8.811393e-98 0.0001
#> 2 2 0.9713667 0.03201663 30.33945 9.318861e-144 0.0001
#> 3 3 -0.007682607 0.04222391 -0.1819492 0.8556595 0.9786
#> 4 4 -0.02689601 0.03050616 -0.8816584 0.3781741 0.7402
#> 5 5 0.411529 0.04299497 9.571561 7.9842e-21 0.0001
#> 6 6 0.3925661 0.03096423 12.67805 2.946569e-34 0.0001
#> 7 7 0.0206361 0.04405654 0.4684003 0.6396006 0.9112
#> 8 8 0.001657765 0.03337464 0.04967138 0.9603942 0.9786fit1 <- feols(Y1 ~ X1 , data = data)
fit2 <- feols(Y1 ~ X1 + X2, data = data)
fit3 <- feols(Y2 ~ X1, data = data)
fit4 <- feols(Y2 ~ X1 + X2, data = data)
res_rwolf2 <- rwolf(
models = list(fit1, fit2, fit3, fit4),
param = "X1",
B = 9999
)
#> | | | 0% | |================== | 25% | |=================================== | 50% | |==================================================== | 75% | |======================================================================| 100%
res_rwolf2
#> model Estimate Std. Error t value Pr(>|t|) RW Pr(>|t|)
#> 1 1 0.9896609 0.04341633 22.79467 6.356963e-93 0.0001
#> 2 2 0.9713667 0.03186495 30.48386 9.523796e-145 0.0001
#> 3 3 -0.007682607 0.04403736 -0.1744566 0.861542 0.8568
#> 4 4 -0.02689601 0.03130345 -0.8592027 0.3904352 0.5439The above procedure with S=8 hypotheses, N=1000 observations and
k %in% (1,2) parameters finishes in around 5 seconds.
if(requireNamespace("microbenchmark")){
microbenchmark::microbenchmark(
"Romano-Wolf" = wildrwolf::rwolf(
models = fit,
param = "X1",
B = 9999
),
times = 1
)
}
#> | | | 0% | |========= | 12% | |================== | 25% | |========================== | 38% | |=================================== | 50% | |============================================ | 62% | |==================================================== | 75% | |============================================================= | 88% | |======================================================================| 100%
#> Unit: seconds
#> expr min lq mean median uq max neval
#> Romano-Wolf 3.604916 3.604916 3.604916 3.604916 3.604916 3.604916 1We test
S <- 6
rho <- 0.5
Sigma <- matrix(rho, 6, 6)
diag(Sigma) <- 1
Sigma
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 1.0 0.5 0.5 0.5 0.5 0.5
#> [2,] 0.5 1.0 0.5 0.5 0.5 0.5
#> [3,] 0.5 0.5 1.0 0.5 0.5 0.5
#> [4,] 0.5 0.5 0.5 1.0 0.5 0.5
#> [5,] 0.5 0.5 0.5 0.5 1.0 0.5
#> [6,] 0.5 0.5 0.5 0.5 0.5 1.0with
This experiment imposes a data generating process as in equation (9) in
Clarke, Romano and Wolf, with an
additional error term
You can run the simulations via the run_fwer_sim() function attached
in the package.
# note that this will take some time
res <- run_fwer_sim(
seed = 76,
n_sims = 1000,
B = 499,
N = 1000,
s = 6,
rho = 0.5 #correlation between hypotheses, not intra-cluster!
)Both Holm’s method and wildrwolf control the family wise error rates,
at both the 5 and 10% significance level.
res
#> reject_5 reject_10 rho
#> fit_pvalue 0.999 0.999 0.5
#> fit_pvalue_holm 0.000 0.000 0.5
#> fit_padjust_rw 0.000 0.000 0.5library(RStata)
# initiate RStata
options("RStata.StataPath" = "\"C:\\Program Files\\Stata17\\StataBE-64\"")
options("RStata.StataVersion" = 17)
# save the data set so it can be loaded into STATA
write.csv(data, "c:/Users/alexa/Dropbox/rwolf/inst/extdata/readme.csv")
# estimate with stata via Rstata
stata_program <- "
clear
set more off
import delimited c:/Users/alexa/Dropbox/rwolf/inst/data/readme.csv
set seed 1
rwolf y1 y2 y3 y4, indepvar(x1) controls(x2) reps(9999)
"
RStata::stata(stata_program, data.out = TRUE)
# Romano-Wolf step-down adjusted p-values
#
#
# Independent variable: x1
# Outcome variables: y1 y2 y3 y4
# Number of resamples: 9999
#
#
# ------------------------------------------------------------------------------
# Outcome Variable | Model p-value Resample p-value Romano-Wolf p-value
# --------------------+---------------------------------------------------------
# y1 | 0.0000 0.0001 0.0001
# y2 | 0.3904 0.3755 0.6070
# y3 | 0.0000 0.0001 0.0001
# y4 | 0.9586 0.9596 0.9596
# ------------------------------------------------------------------------------For comparison, wildrwolf produces the following output:
models <- feols(c(Y1, Y2, Y3, Y4) ~ X1 + X2
, data = data, se = "hetero")rwolf(models, param = "X1", B = 9999)
#> | | | 0% | |================== | 25% | |=================================== | 50% | |==================================================== | 75% | |======================================================================| 100%
#> model Estimate Std. Error t value Pr(>|t|) RW Pr(>|t|)
#> 1 1 0.9713667 0.03201663 30.33945 9.318861e-144 0.0001
#> 2 2 -0.02689601 0.03050616 -0.8816584 0.3781741 0.5922
#> 3 3 0.3925661 0.03096423 12.67805 2.946569e-34 0.0001
#> 4 4 0.001657765 0.03337464 0.04967138 0.9603942 0.9618