--- title: "Optimal Multiple Hypothesis Testing" author: "Davide Viviano, Kaspar Wuthrich, Paul Niehaus" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Optimal Multiple Hypothesis Testing} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} bibliography: references.bib --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ## Overview The `mhtopt` package implements the economically optimal multiple hypothesis testing (MHT) framework from Viviano, Wuthrich, and Niehaus (2026). Unlike ad hoc procedures (Bonferroni, Holm, Benjamini-Hochberg), this package derives the *correct* adjustment from the structure of research costs. **Key insight:** The optimal per-test significance level is (Proposition 4.1): $$ \alpha^*(J) = \frac{C(J, n)}{b \cdot \bar{\omega}(J)} $$ where $C(J, n)$ is the research cost function, $b$ is the researcher's benefit from a rejection, and $\bar{\omega}(J)$ is the sum of hypothesis weights. This simplifies to closed-form expressions depending on the cost structure. ### Package functions | Function | Purpose | |----------|---------| | `mht_est()` | Main function: test J treatment coefficients from a regression | | `mht_critical()` | Compute optimal critical values for given J | | `mht_test()` | Apply adjustment to a vector of p-values | | `mht_cost_estimate()` | Estimate cost function parameters from data | | `mht_table()` | Generate reference tables of critical values | ## Installation ```{r eval=FALSE} # From CRAN (once published) install.packages("mhtopt") # Development version from GitHub # remotes::install_github("dviviano/mhtopt", subdir = "r/mhtopt") ``` ```{r} library(mhtopt) ``` ## Quick Start ```{r eval=FALSE} # After running a regression with multiple treatments fit <- lm(outcome ~ treat1 + treat2 + treat3 + controls, data = mydata) # Test all treatment coefficients with optimal MHT adjustment result <- mht_est(fit, vars = c("treat1", "treat2", "treat3"), alpha_bar = 0.05) # View results: shows rejections under 5 procedures print(result) ``` The function returns rejection decisions under five procedures: - **Optimal** — derived from cost structure (default: FDA calibration) - **Bonferroni** — traditional FWER control - **Holm** — step-down Bonferroni - **Benjamini-Hochberg** — FDR control - **Unadjusted** — no MHT correction ## When to Adjust? The paper's central insight: the need for MHT adjustment depends on the **research cost function** — specifically, whether costs exhibit returns to scale with respect to the number of hypotheses. ### Cost structure determines the optimal correction - **Fixed costs dominate** (`cf_share` near 1, or `beta` near 0): Adding more hypotheses barely increases costs → researchers have incentive to "test many and get lucky" → **optimal correction approaches Bonferroni** - **Variable costs dominate** (`cf_share` near 0, or `beta` near 1): Each hypothesis costs proportionally more → cost itself deters excessive testing → **little/no adjustment needed** - **Intermediate** (empirically relevant): Optimal threshold lies between Bonferroni and unadjusted, determined by cost parameters ### Default calibrations **Linear model (FDA calibration):** ```{r} # Default: cf_share = 0.46, J_bar = 3 (from Sertkaya et al. 2016) mht_critical(J = 5, alpha_bar = 0.05) ``` **Cobb-Douglas (J-PAL calibration):** ```{r} # beta = 0.13, iota = 0.075 (from J-PAL data, Appendix A) mht_critical(J = 5, alpha_bar = 0.05, model = "cobbdouglas") ``` Notice: Cobb-Douglas is *more conservative* because beta ≈ 0 implies near-constant marginal cost (fixed-cost-like). ## Main Function: `mht_est()` Tests J coefficients from a fitted regression model. ### Basic usage ```{r eval=FALSE} mht_est(fit, vars = NULL, # coefficient names (default: all non-intercept) alpha_bar, # benchmark significance level onesided = TRUE, # one-sided (positive) vs two-sided model = "linear", # "linear" or "cobbdouglas" cf_share = 0.46, # fixed cost share (linear model) J_bar = 3, # avg arms per trial (linear model) nm_ratio = 1.0, # sample size ratio n/m mbar = NULL, # benchmark arm size (overrides nm_ratio) beta = 0.13, # arms elasticity (Cobb-Douglas) iota = 0.075) # sample elasticity (Cobb-Douglas) ``` ### Supported model classes - `lm` — base R linear models - `glm` — generalized linear models - `lm_robust`, `iv_robust` — from `estimatr` (cluster/HC standard errors) - `fixest` — from `fixest` (`feols`, `fepois`, etc.) - Any class with a `coef(summary(fit))` method ### Return value A `data.frame` with columns: - `term` — coefficient name - `estimate`, `std_error`, `t_stat`, `p_value` — from the model - `reject_optimal`, `reject_bonferroni`, `reject_holm`, `reject_bh`, `reject_unadjusted` — rejection indicators Attributes: `alpha_opt`, `alpha_bonf`, `alpha_bar`, `J`, `model`, `onesided`, `nm_ratio` ### Examples **Multiple treatment arms:** ```{r eval=FALSE} library(estimatr) fit <- lm_robust(outcome ~ treat_A + treat_B + treat_C + controls, fixed_effects = ~village, clusters = ~hh_id, data = df) result <- mht_est(fit, vars = c("treat_A", "treat_B", "treat_C"), alpha_bar = 0.05) # Extract optimal alpha used attr(result, "alpha_opt") # How many hypotheses rejected under each procedure? colSums(result[, c("reject_optimal", "reject_bonferroni", "reject_holm", "reject_bh", "reject_unadjusted")]) ``` **Two-sided test:** ```{r eval=FALSE} fit_logit <- glm(uptake ~ treat1 + treat2 + controls, family = binomial, data = df) mht_est(fit_logit, vars = c("treat1", "treat2"), alpha_bar = 0.05, onesided = FALSE) ``` **Using benchmark sample size (`mbar`):** ```{r eval=FALSE} # If you know typical arm size is 500 observations # and your study has N = 1500 with J = 3 arms # then nm_ratio = (1500/3) / 500 = 1.0 fit <- lm(outcome ~ treat1 + treat2 + treat3, data = df) mht_est(fit, vars = c("treat1", "treat2", "treat3"), alpha_bar = 0.05, mbar = 500) # automatically computes nm_ratio ``` **Custom cost parameters:** ```{r eval=FALSE} # Suppose you estimated cost parameters from your data est <- mht_cost_estimate(cost, arms, sample_size, alpha_bar = 0.05) mht_est(fit, vars = c("treat1", "treat2", "treat3"), alpha_bar = 0.05, model = "cobbdouglas", beta = est$beta, iota = est$iota) ``` ## Supporting Functions ### `mht_critical()` — Compute Optimal Critical Values The core formula. Computes $\alpha^*$ for given J and cost model. ```{r} # Linear model (default) mht_critical(J = 5, alpha_bar = 0.05) # Cobb-Douglas mht_critical(J = 5, alpha_bar = 0.05, model = "cobbdouglas") # Larger sample size => less conservative mht_critical(J = 5, alpha_bar = 0.05, nm_ratio = 2.0) ``` **Cost model formulas:** Linear (Eq. 26 of the paper): $$ \alpha^* = \alpha_{\text{bar}} \left[ \frac{1 + r/J}{1 + r} + \frac{n/m - 1}{1 + r} \right] $$ where $r = \frac{c_f \cdot \bar{J}}{c_v \cdot m} = \frac{\text{cf\_share} \cdot J_{\text{bar}}}{1 - \text{cf\_share}}$ Cobb-Douglas: $$ \alpha^* = \alpha_{\text{bar}} \cdot J^{\beta - 1} \cdot (n/m)^\iota $$ ### `mht_test()` — Test a Vector of P-Values Use when you have pre-computed p-values (not from a regression). ```{r} pvals <- c(0.003, 0.015, 0.030, 0.048, 0.120, 0.500) result <- mht_test(p = pvals, alpha_bar = 0.05) print(result) ``` Alternatively, pass z-statistics instead of p-values: ```{r} zstats <- c(3.2, 2.1, 1.8, 1.3, 0.8, -0.2) mht_test(z = zstats, alpha_bar = 0.05) ``` ### `mht_table()` — Generate Reference Tables Reproduces Table 1 from the paper (or custom tables). ```{r} # Default: reproduces paper Table 1 exactly # (J = 1..9 plus infinity, multiple alpha_bar and nm_ratio values) mht_table() # Custom table mht_table(alpha_bar = 0.05, J_range = 1:10, nm_ratios = c(0.5, 1.0, 1.5, 2.0)) ``` Returns a formatted table and stores all values in `attr(result, "values")`. ## Estimating Cost Functions ### `mht_cost_estimate()` — Estimate Parameters from Data If you have data on research project costs, you can estimate cost function parameters. **Cobb-Douglas model (default):** Estimates: $\log(C) = \text{const} + \beta \log(J) + \iota \log(n) + \text{controls}$ Returns $\beta$ and $\iota$ with hypothesis tests: - $H_0: \beta = 0$ ⇒ Bonferroni is optimal - $H_0: \beta = 1$ ⇒ no adjustment needed ```{r eval=FALSE} # Example data: cost, number of arms, sample size projects <- data.frame( cost = c(500000, 750000, 1200000, 600000, 900000), arms = c(2, 3, 5, 2, 4), n = c(1000, 1500, 2000, 800, 1200) ) est <- mht_cost_estimate( cost = projects$cost, arms = projects$arms, sample_size = projects$n, alpha_bar = 0.05, model = "cobbdouglas" ) # View estimates print(est) # Use in mht_est() fit <- lm(outcome ~ treat1 + treat2 + treat3, data = mydata) mht_est(fit, vars = c("treat1", "treat2", "treat3"), alpha_bar = 0.05, model = "cobbdouglas", beta = est$beta, iota = est$iota) ``` **Linear model:** ```{r eval=FALSE} est_linear <- mht_cost_estimate( cost = projects$cost, arms = projects$arms, sample_size = projects$n, alpha_bar = 0.05, model = "linear_share" ) # Returns cf_share (fixed cost share) mht_est(fit, vars = c("treat1", "treat2", "treat3"), alpha_bar = 0.05, cf_share = est_linear$cf_share) ``` ## Worked Example ```{r eval=FALSE} # Simulate a simple experiment set.seed(123) n <- 300 df <- data.frame( treat1 = rbinom(n, 1, 0.33), treat2 = rbinom(n, 1, 0.33), treat3 = rbinom(n, 1, 0.33) ) df$outcome <- 0.3 * df$treat1 + 0.2 * df$treat2 + 0.05 * df$treat3 + rnorm(n) # Estimate effects fit <- lm(outcome ~ treat1 + treat2 + treat3, data = df) summary(fit) # Apply optimal MHT adjustment result <- mht_est(fit, vars = c("treat1", "treat2", "treat3"), alpha_bar = 0.05) print(result) # Compare: how many rejections under each procedure? colSums(result[, grep("^reject_", names(result))]) ``` ## One-Sided vs Two-Sided Tests Default behavior (`onesided = TRUE`): tests $H_0: \beta_j \leq 0$ vs $H_1: \beta_j > 0$. **Justification** (Remark 4.1, paper v10): When the status quo is "no intervention," a negative effect means the status quo holds (no social value lost by failing to reject). One-sided tests are appropriate. For two-sided tests (Appendix B.1): ```{r eval=FALSE} mht_est(fit, vars = c("treat1", "treat2", "treat3"), alpha_bar = 0.05, onesided = FALSE) ``` P-values are multiplied by 2 and compared against the same $\alpha^*$. ## Sample Size Adjustment The `nm_ratio` parameter captures how your study's sample size compares to a benchmark: $$ \text{nm\_ratio} = \frac{n/J}{m} $$ where: - $n$ = total sample size in your study - $J$ = number of hypotheses/arms - $m$ = benchmark per-arm sample size **Larger studies are more expensive** → optimal $\alpha^*$ is *less* conservative (encourages research): ```{r} # Small study (half the benchmark size) mht_critical(J = 3, alpha_bar = 0.05, nm_ratio = 0.5) # Benchmark size mht_critical(J = 3, alpha_bar = 0.05, nm_ratio = 1.0) # Large study (double the benchmark) mht_critical(J = 3, alpha_bar = 0.05, nm_ratio = 2.0) ``` The `mbar` argument automates this: if you know the typical per-arm sample size in your field, pass it via `mbar`, and the function computes `nm_ratio` automatically from the model's sample size. ## References Viviano, D., Wuthrich, K., & Niehaus, P. (2026). A Model of Multiple Hypothesis Testing. *arXiv:2104.13367v10*. Sertkaya, A., Wong, H. H., Jessup, A., & Beleche, T. (2016). Key cost drivers of pharmaceutical clinical trials in the United States. *Clinical Trials*, 13(2), 117-126.