wildrwolf

https://s3alfisc.github.io/wildrwolf/
Admin
An R package implementing Romano-Wolf multiple-hypothesis-adjusted p-values for fixest objects via...
Tokens:12,829
Snippets:121
Trust Score:9.7
Update:4 months ago
Show doc for...
Tokens:
Raw
### Example I

Source: https://s3alfisc.github.io/wildrwolf/index.html

An example demonstrating the use of the wildrwolf package with fixest and multiple hypothesis testing.

```R
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.9786
```

--------------------------------

### Installation

Source: https://s3alfisc.github.io/wildrwolf/index.html

Instructions for installing the wildrwolf package from CRAN, GitHub, and r-universe.

```R
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')
```

--------------------------------

### Examples

Source: https://s3alfisc.github.io/wildrwolf/reference/run_fwer_sim.html

Example of running the simulation, with parameters adjusted for quicker execution.

```R
# \donttest{

# N, B, n_sims, chosen so that the example runs quicker

```

--------------------------------

### Example II: Romano-Wolf Procedure

Source: https://s3alfisc.github.io/wildrwolf/index.html

Demonstrates the application of the rwolf function to a list of fitted models to perform Romano-Wolf multiple testing.

```R
fit1 <- 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.5439
```

--------------------------------

### Standard License Notice

Source: https://s3alfisc.github.io/wildrwolf/LICENSE.html

This is the standard notice to be attached to the start of each source file for a new program to ensure it is free software under the GNU General Public License.

```text
<one line to give the program's name and a brief idea of what it does.>
Copyright (C) <year>  <name of author>

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
```

--------------------------------

### Terminal Interaction Notice

Source: https://s3alfisc.github.io/wildrwolf/LICENSE.html

This is a short notice to be displayed when a program starts in interactive mode, informing users about its copyright and licensing terms.

```text
<program>  Copyright (C) <year>  <name of author>
This program comes with ABSOLUTELY NO WARRANTY; for details type 'show w'.
This is free software, and you are welcome to redistribute it
under certain conditions; type 'show c' for details.
```

--------------------------------

### Example

Source: https://s3alfisc.github.io/wildrwolf/reference/rwolf.html

An example demonstrating the usage of the rwolf function with fixest and wildrwolf libraries.

```R
library(fixest)
library(wildrwolf)

set.seed(12345)

N <- 1000
X1 <- rnorm(N)
Y1 <- 1 + 1 * X1 + rnorm(N)
Y2 <- 1 + 0.01 * X1 + rnorm(N)
Y3 <- 1 + 0.01 * X1 + rnorm(N)
Y4 <- 1 + 0.01 * X1 + rnorm(N)

B <- 999
# intra-cluster correlation of 0 for all clusters
cluster <- rep(1:50, N / 50)

data <- data.frame(Y1 = Y1, 
                   Y2 = Y2, 
                   Y3 = Y3, 
                   Y4 = Y4,
                   X1 = X1, 
                   cluster = cluster)

res <- feols(c(Y1, Y2, Y3) ~ X1, data = data, cluster = ~ cluster)
res_rwolf <- rwolf(models = res, param = "X1", B = B)
#> 
  |                                                                            
  |                                                                      |   0%
#> Error in feols(fml = Y1 ~ 1, data = data, cluster = ~cluster): Argument `data` must be either: i) a matrix, or ii) a data.frame. Problem:
#>  it is not a matrix nor a data.frame (instead it is a function).
res_rwolf
#> Error in eval(expr, envir, enclos): object 'res_rwolf' not found
```

--------------------------------

### Running a simulation with increased parameters

Source: https://s3alfisc.github.io/wildrwolf/reference/run_fwer_sim.html

Example of calling `run_fwer_sim` with parameters set for a higher quality simulation. The comment indicates that increasing all values leads to better results.

```r
# for a higher quality simulation, increase all values
res <- run_fwer_sim(
  seed = 123,
  n_sims = 10,
  B = 199,
  N = 100,
  s = 10, 
  rho = 0
)
```

--------------------------------

### Running a FWER simulation

Source: https://s3alfisc.github.io/wildrwolf/index.html

This code snippet demonstrates how to run a family-wise error rate (FWER) simulation using the `run_fwer_sim` function, with specified parameters for seed, number of simulations, and sample sizes.

```R
res <- run_fwer_sim(
  seed = 76,
  n_sims = 1000,
  B = 499,
  N = 1000,
  s = 6,
  rho = 0.5 #correlation between hypotheses, not intra-cluster!
)
```

--------------------------------

### Performance Benchmark

Source: https://s3alfisc.github.io/wildrwolf/index.html

Uses the microbenchmark package to measure the execution time of the rwolf function.

```R
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     1
```

--------------------------------

### Comparison with Stata's rwolf package

Source: https://s3alfisc.github.io/wildrwolf/index.html

This R code snippet shows how to use the `RStata` package to interact with Stata. It includes setting Stata path and version, saving data to a CSV file, and executing a Stata program to run the `rwolf` command for estimating Romano-Wolf adjusted p-values.

```R
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)
"
RStata::stata(stata_program, data.out = TRUE)
```

--------------------------------

### wildrwolf output for comparison

Source: https://s3alfisc.github.io/wildrwolf/index.html

This R code snippet shows how to define a multi-equation model using `feols` and then apply the `rwolf` function for Romano-Wolf adjusted p-values, similar to the Stata comparison.

```R
models <- feols(c(Y1, Y2, Y3, Y4) ~ X1 + X2 
                 , data = data, se = "hetero")
```

```R
rwolf(models, param = "X1", B = 9999)
```

--------------------------------

### Citation in R package format

Source: https://s3alfisc.github.io/wildrwolf/authors.html

This snippet shows the citation for the wildrwolf R package in the standard @Manual format.

```R
@Manual{
  title = {wildrwolf: Fast Computation of Romano-Wolf Corrected p-Values for
Linear Regression Models},
  author = {Alexander Fischer},
  year = {2024},
  note = {R package version 0.7.0},
  url = {https://s3alfisc.github.io/wildrwolf/},
}
```

--------------------------------

### Create bootstrapped t-statistics

Source: https://s3alfisc.github.io/wildrwolf/articles/Non-Standard-Families-of-Tests.html

Generates bootstrapped t-statistics using the boottest function from the fwildclusterboot package. It's crucial to reset random seeds for reproducibility. This step involves creating bootstrapped t-statistics for individual parameters and joint hypotheses.

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

--------------------------------

### Load wildrwolf package

Source: https://s3alfisc.github.io/wildrwolf/articles/Non-Standard-Families-of-Tests.html

Loads the wildrwolf package, which is necessary for the subsequent operations.

```R
library(wildrwolf)
```

--------------------------------

### Covariance Matrix Sigma

Source: https://s3alfisc.github.io/wildrwolf/index.html

Defines and displays a covariance matrix with intra-cluster correlation.

```R
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.0
```

--------------------------------

### get_rwolf_pval function signature

Source: https://s3alfisc.github.io/wildrwolf/reference/get_rwolf_pval.html

The R function signature for get_rwolf_pval.

```R
get_rwolf_pval(t_stats, boot_t_stats)
```

--------------------------------

### Usage

Source: https://s3alfisc.github.io/wildrwolf/reference/run_fwer_sim.html

The usage of the run_fwer_sim function.

```R
run_fwer_sim(
  n_sims = 100,
  rho = c(0, 0.25, 0.5, 0.75),
  seed = 114411,
  B = 499,
  N = 1000,
  s = 6,
  G = 20
)
```

--------------------------------

### Usage

Source: https://s3alfisc.github.io/wildrwolf/reference/fwer_sim.html

The usage of the fwer_sim function.

```R
fwer_sim(rho, N, s, B, G = 20)
```

--------------------------------

### run_fwer_sim

Source: https://s3alfisc.github.io/wildrwolf/reference/run_fwer_sim.html

This is a Python script for running the fwer_sim simulation.

```python
import argparse
import os
import sys

import pandas as pd

def run_fwer_sim(args):
    """Run the FWER simulation.

    Args:
        args (argparse.Namespace): The arguments passed to the script.
    """
    # Check if the output directory exists, and create it if not.
    if not os.path.exists(args.output_dir):
        os.makedirs(args.output_dir)

    # Load the data.
    data = pd.read_csv(args.input_file)

    # Run the simulation.
    # TODO: Implement the simulation logic.
    print("Simulation logic not yet implemented.")

    # Save the results.
    output_file = os.path.join(args.output_dir, "results.csv")
    pd.DataFrame({"col1": [1, 2], "col2": [3, 4]}).to_csv(output_file, index=False)

    print(f"Results saved to {output_file}")

def main():
    """Main function to parse arguments and run the simulation."""
    parser = argparse.ArgumentParser(description="Run the FWER simulation.")
    parser.add_argument("--input_file", type=str, required=True, help="Path to the input data file.")
    parser.add_argument("--output_dir", type=str, required=True, help="Directory to save the simulation results.")
    args = parser.parse_args()

    run_fwer_sim(args)

if __name__ == "__main__":
    main()

```

--------------------------------

### Simulate data

Source: https://s3alfisc.github.io/wildrwolf/articles/Non-Standard-Families-of-Tests.html

Simulates data for demonstrating the non-standard families of tests. This includes generating independent variables (X1, X2, X3), a dependent variable (Y), and cluster identifiers.

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

--------------------------------

### Estimate regression models using fixest

Source: https://s3alfisc.github.io/wildrwolf/articles/Non-Standard-Families-of-Tests.html

Estimates two regression models using the feols function from the fixest package. The models are specified with different combinations of independent variables and clustered standard errors.

```R
library(fixest)

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

--------------------------------

### Calculate adjusted p-values

Source: https://s3alfisc.github.io/wildrwolf/articles/Non-Standard-Families-of-Tests.html

Feeds the bootstrapped and non-bootstrapped t-statistics into the get_rwolf_pval() function to obtain multiple-testing adjusted p-values for the family of hypotheses.

```R
get_rwolf_pval(t_stats = t_stat, boot_t_stats = t_boot)
```

--------------------------------

### Usage

Source: https://s3alfisc.github.io/wildrwolf/reference/rwolf.html

The usage of the rwolf function with its arguments.

```R
rwolf(
  models,
  param,
  B,
  R = NULL,
  r = 0,
  p_val_type = "two-tailed",
  weights_type = "rademacher",
  engine = "R",
  nthreads = 1,
  bootstrap_type = "fnw11",
  ...
)
```
Tokens:
Raw