The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
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 development version from GitHub with:
# install.packages("devtools")
::install_github("s3alfisc/wildrwolf")
devtools
# 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)
<- 1000
N <- rnorm(N)
X1 <- rnorm(N)
X2 <- 0.5
rho <- matrix(rho, 4, 4); diag(sigma) <- 1
sigma <- MASS::mvrnorm(n = N, mu = rep(0, 4), Sigma = sigma)
u <- 1 + 1 * X1 + X2
Y1 <- 1 + 0.01 * X1 + X2
Y2 <- 1 + 0.4 * X1 + X2
Y3 <- 1 + -0.02 * X1 + X2
Y4 for(x in 1:4){
<- paste0("Y", x)
var_char assign(var_char, get(var_char) + u[,x])
}
<- data.frame(Y1 = Y1,
data Y2 = Y2,
Y3 = Y3,
Y4 = Y4,
X1 = X1,
X2 = X2,
#group_id = group_id,
splitvar = sample(1:2, N, TRUE))
<- feols(c(Y1, Y2, Y3, Y4) ~ csw(X1,X2),
fit data = data,
se = "hetero",
ssc = ssc(cluster.adj = TRUE))
# clean workspace except for res & data
rm(list= ls()[!(ls() %in% c('fit','data'))])
<- wildrwolf::rwolf(
res_rwolf1 models = fit,
param = "X1",
B = 9999
)#> | | | 0% | |========= | 12% | |================== | 25% | |========================== | 38% | |=================================== | 50% | |============================================ | 62% | |==================================================== | 75% | |============================================================= | 88% | |======================================================================| 100%
<- lapply(fit, function(x) pvalue(x)["X1"]) |> unlist()
pvals
# 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.9786
<- feols(Y1 ~ X1 , data = data)
fit1 <- feols(Y1 ~ X1 + X2, data = data)
fit2 <- feols(Y2 ~ X1, data = data)
fit3 <- feols(Y2 ~ X1 + X2, data = data)
fit4
<- rwolf(
res_rwolf2 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.5439
The 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 5.041018 5.041018 5.041018 5.041018 5.041018 5.041018 1
We test \(S=6\) hypotheses and generate data as
\[Y_{i,s,g} = \beta_{0} + \beta_{1,s} D_{i} + u_{i,g} + \epsilon_{i,s} \] where \(D_i = 1(U_i > 0.5)\) and \(U_i\) is drawn from a uniform distribution, \(u_{i,g}\) is a cluster level shock with intra-cluster correlation \(0.5\), and the idiosyncratic error term is drawn from a multivariate random normal distribution with mean \(0_S\) and covariance matrix
<- 6
S <- 0.5
rho <- matrix(rho, 6, 6)
Sigma 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.0
with \(\rho \geq 0\). We assume that \(\beta_{1,s}= 0\) for all \(s\).
This experiment imposes a data generating process as in equation (9) in Clarke, Romano and Wolf, with an additional error term \(u_g\) for \(G=20\) clusters and intra-cluster correlation 0.5 and \(N=1000\) observations.
You can run the simulations via the run_fwer_sim()
function attached in the package.
# note that this will take some time
<- run_fwer_sim(
res 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.5
library(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)
"
::stata(stata_program, data.out = TRUE)
RStata
# 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:
<- feols(c(Y1, Y2, Y3, Y4) ~ X1 + X2
models 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
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.