Optimal Multiple Hypothesis Testing

Davide Viviano, Kaspar Wuthrich, Paul Niehaus

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

# From CRAN (once published)
install.packages("mhtopt")

# Development version from GitHub
# remotes::install_github("dviviano/mhtopt", subdir = "r/mhtopt")
library(mhtopt)

Quick Start

# 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

Default calibrations

Linear model (FDA calibration):

# Default: cf_share = 0.46, J_bar = 3 (from Sertkaya et al. 2016)
mht_critical(J = 5, alpha_bar = 0.05)
#> 
#> ------------------------------------------------------------ 
#>   Optimal MHT Critical Values
#>   Viviano, Wuthrich, and Niehaus (2026)
#> ------------------------------------------------------------ 
#> 
#>   Cost model:          Linear (Eq. 27)
#>   Number of hypotheses: 5
#>   Benchmark alpha:      0.0500
#>   Sample size ratio:    1.00
#> 
#> ------------------------------------------------------------ 
#>   Optimal test size:    0.021250
#>   Optimal z-threshold:  2.0286
#> 
#>   Bonferroni size:      0.010000
#>   Bonferroni z-thresh:  2.3263
#>   Sidak size:           0.010206
#>   Sidak z-threshold:    2.3187
#>   Unadjusted size:      0.050000
#> ------------------------------------------------------------

Cobb-Douglas (J-PAL calibration):

# beta = 0.13, iota = 0.075 (from J-PAL data, Appendix A)
mht_critical(J = 5, alpha_bar = 0.05, model = "cobbdouglas")
#> 
#> ------------------------------------------------------------ 
#>   Optimal MHT Critical Values
#>   Viviano, Wuthrich, and Niehaus (2026)
#> ------------------------------------------------------------ 
#> 
#>   Cost model:          Cobb-Douglas (Appendix A)
#>   Number of hypotheses: 5
#>   Benchmark alpha:      0.0500
#>   Sample size ratio:    1.00
#> 
#> ------------------------------------------------------------ 
#>   Optimal test size:    0.012327
#>   Optimal z-threshold:  2.2468
#> 
#>   Bonferroni size:      0.010000
#>   Bonferroni z-thresh:  2.3263
#>   Sidak size:           0.010206
#>   Sidak z-threshold:    2.3187
#>   Unadjusted size:      0.050000
#> ------------------------------------------------------------

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

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

Return value

A data.frame with columns:

Attributes: alpha_opt, alpha_bonf, alpha_bar, J, model, onesided, nm_ratio

Examples

Multiple treatment arms:

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:

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

# 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:

# 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.

# Linear model (default)
mht_critical(J = 5, alpha_bar = 0.05)
#> 
#> ------------------------------------------------------------ 
#>   Optimal MHT Critical Values
#>   Viviano, Wuthrich, and Niehaus (2026)
#> ------------------------------------------------------------ 
#> 
#>   Cost model:          Linear (Eq. 27)
#>   Number of hypotheses: 5
#>   Benchmark alpha:      0.0500
#>   Sample size ratio:    1.00
#> 
#> ------------------------------------------------------------ 
#>   Optimal test size:    0.021250
#>   Optimal z-threshold:  2.0286
#> 
#>   Bonferroni size:      0.010000
#>   Bonferroni z-thresh:  2.3263
#>   Sidak size:           0.010206
#>   Sidak z-threshold:    2.3187
#>   Unadjusted size:      0.050000
#> ------------------------------------------------------------

# Cobb-Douglas
mht_critical(J = 5, alpha_bar = 0.05, model = "cobbdouglas")
#> 
#> ------------------------------------------------------------ 
#>   Optimal MHT Critical Values
#>   Viviano, Wuthrich, and Niehaus (2026)
#> ------------------------------------------------------------ 
#> 
#>   Cost model:          Cobb-Douglas (Appendix A)
#>   Number of hypotheses: 5
#>   Benchmark alpha:      0.0500
#>   Sample size ratio:    1.00
#> 
#> ------------------------------------------------------------ 
#>   Optimal test size:    0.012327
#>   Optimal z-threshold:  2.2468
#> 
#>   Bonferroni size:      0.010000
#>   Bonferroni z-thresh:  2.3263
#>   Sidak size:           0.010206
#>   Sidak z-threshold:    2.3187
#>   Unadjusted size:      0.050000
#> ------------------------------------------------------------

# Larger sample size => less conservative
mht_critical(J = 5, alpha_bar = 0.05, nm_ratio = 2.0)
#> 
#> ------------------------------------------------------------ 
#>   Optimal MHT Critical Values
#>   Viviano, Wuthrich, and Niehaus (2026)
#> ------------------------------------------------------------ 
#> 
#>   Cost model:          Linear (Eq. 27)
#>   Number of hypotheses: 5
#>   Benchmark alpha:      0.0500
#>   Sample size ratio:    2.00
#> 
#> ------------------------------------------------------------ 
#>   Optimal test size:    0.035312
#>   Optimal z-threshold:  1.8079
#> 
#>   Bonferroni size:      0.010000
#>   Bonferroni z-thresh:  2.3263
#>   Sidak size:           0.010206
#>   Sidak z-threshold:    2.3187
#>   Unadjusted size:      0.050000
#> ------------------------------------------------------------

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

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)
#> 
#> ----------------------------------------------------------------- 
#>   Multiple Hypothesis Testing Results
#>   Viviano, Wuthrich, and Niehaus (2026)
#> ----------------------------------------------------------------- 
#> 
#>   Hypotheses tested:    6
#>   Benchmark alpha:      0.0500
#>   Cost model:           Linear
#> 
#> ------------------------------------------------------- 
#>   Procedure                 Test size   Rejections
#> ------------------------------------------------------- 
#>   Optimal (model-based)      0.020052            2
#>   Bonferroni                 0.008333            1
#>   Holm (step-down)           0.0500 *            1
#>   BH (FDR control)           0.0500 *            2
#>   Unadjusted                 0.050000            4
#> ------------------------------------------------------- 
#>   * Step-wise; effective threshold varies by rank
#> 
#>   Detailed results:
#>  p_value reject_optimal reject_bonferroni reject_holm reject_bh
#>    0.003           TRUE              TRUE        TRUE      TRUE
#>    0.015           TRUE             FALSE       FALSE      TRUE
#>    0.030          FALSE             FALSE       FALSE     FALSE
#>    0.048          FALSE             FALSE       FALSE     FALSE
#>    0.120          FALSE             FALSE       FALSE     FALSE
#>    0.500          FALSE             FALSE       FALSE     FALSE
#>  reject_unadjusted
#>               TRUE
#>               TRUE
#>               TRUE
#>               TRUE
#>              FALSE
#>              FALSE

Alternatively, pass z-statistics instead of p-values:

zstats <- c(3.2, 2.1, 1.8, 1.3, 0.8, -0.2)
mht_test(z = zstats, alpha_bar = 0.05)
#> 
#> ----------------------------------------------------------------- 
#>   Multiple Hypothesis Testing Results
#>   Viviano, Wuthrich, and Niehaus (2026)
#> ----------------------------------------------------------------- 
#> 
#>   Hypotheses tested:    6
#>   Benchmark alpha:      0.0500
#>   Cost model:           Linear
#> 
#> ------------------------------------------------------- 
#>   Procedure                 Test size   Rejections
#> ------------------------------------------------------- 
#>   Optimal (model-based)      0.020052            2
#>   Bonferroni                 0.008333            1
#>   Holm (step-down)           0.0500 *            1
#>   BH (FDR control)           0.0500 *            1
#>   Unadjusted                 0.050000            3
#> ------------------------------------------------------- 
#>   * Step-wise; effective threshold varies by rank
#> 
#>   Detailed results:
#>       p_value reject_optimal reject_bonferroni reject_holm reject_bh
#>  0.0006871379           TRUE              TRUE        TRUE      TRUE
#>  0.0178644206           TRUE             FALSE       FALSE     FALSE
#>  0.0359303191          FALSE             FALSE       FALSE     FALSE
#>  0.0968004846          FALSE             FALSE       FALSE     FALSE
#>  0.2118553986          FALSE             FALSE       FALSE     FALSE
#>  0.5792597094          FALSE             FALSE       FALSE     FALSE
#>  reject_unadjusted
#>               TRUE
#>               TRUE
#>               TRUE
#>              FALSE
#>              FALSE
#>              FALSE

mht_table() — Generate Reference Tables

Reproduces Table 1 from the paper (or custom tables).

# Default: reproduces paper Table 1 exactly
# (J = 1..9 plus infinity, multiple alpha_bar and nm_ratio values)
mht_table()
#> 
#> ===================================================================================================================================================== 
#>   Optimal Critical Values (Linear (Eq. 27) model)
#>   Viviano, Wuthrich, and Niehaus (2026)
#> ===================================================================================================================================================== 
#> 
#>       n/m=50%                         n/m=100%                        n/m=150%                        n/m=200%                        Sidak           
#> ----------------------------------------------------------------------------------------------------------------------------------------------------- 
#>   |J| a=0.025 a=0.050 a=0.100 a=0.150 a=0.025 a=0.050 a=0.100 a=0.150 a=0.025 a=0.050 a=0.100 a=0.150 a=0.025 a=0.050 a=0.100 a=0.150 a=0.025 a=0.050 
#> ----------------------------------------------------------------------------------------------------------------------------------------------------- 
#>     1  0.021  0.043  0.086  0.129  0.025  0.050  0.100  0.150  0.029  0.057  0.114  0.171  0.032  0.064  0.128  0.192  0.025  0.050 
#>     2  0.012  0.025  0.050  0.075  0.016  0.032  0.064  0.096  0.020  0.039  0.078  0.117  0.023  0.046  0.092  0.138  0.013  0.025 
#>     3  0.010  0.019  0.038  0.057  0.013  0.026  0.052  0.078  0.017  0.033  0.066  0.099  0.020  0.040  0.080  0.120  0.008  0.017 
#>     4  0.008  0.016  0.032  0.048  0.012  0.023  0.046  0.069  0.015  0.030  0.060  0.090  0.019  0.037  0.074  0.111  0.006  0.013 
#>     5  0.007  0.014  0.028  0.043  0.011  0.021  0.042  0.064  0.014  0.028  0.057  0.085  0.018  0.035  0.071  0.106  0.005  0.010 
#>     6  0.007  0.013  0.026  0.039  0.010  0.020  0.040  0.060  0.014  0.027  0.054  0.081  0.017  0.034  0.068  0.102  0.004  0.009 
#>     7  0.006  0.012  0.024  0.036  0.010  0.019  0.038  0.058  0.013  0.026  0.052  0.079  0.017  0.033  0.067  0.100  0.004  0.007 
#>     8  0.006  0.012  0.023  0.035  0.009  0.019  0.037  0.056  0.013  0.026  0.051  0.077  0.016  0.033  0.065  0.098  0.003  0.006 
#>     9  0.006  0.011  0.022  0.033  0.009  0.018  0.036  0.054  0.013  0.025  0.050  0.075  0.016  0.032  0.064  0.096  0.003  0.006 
#>   Inf  0.004  0.007  0.014  0.021  0.007  0.014  0.028  0.042  0.011  0.021  0.042  0.063  0.014  0.028  0.056  0.084  0.000  0.000 
#> -----------------------------------------------------------------------------------------------------------------------------------------------------

# Custom table
mht_table(alpha_bar = 0.05,
          J_range = 1:10,
          nm_ratios = c(0.5, 1.0, 1.5, 2.0))
#> 
#> ===================================================== 
#>   Optimal Critical Values (Linear (Eq. 27) model)
#>   Viviano, Wuthrich, and Niehaus (2026)
#> ===================================================== 
#> 
#>       n/m=50% n/m=100% n/m=150% n/m=200% Sidak           
#> ----------------------------------------------------- 
#>   |J| a=0.050 a=0.050 a=0.050 a=0.050 a=0.025 a=0.050 
#> ----------------------------------------------------- 
#>     1  0.043  0.050  0.057  0.064  0.025  0.050 
#>     2  0.025  0.032  0.039  0.046  0.013  0.025 
#>     3  0.019  0.026  0.033  0.040  0.008  0.017 
#>     4  0.016  0.023  0.030  0.037  0.006  0.013 
#>     5  0.014  0.021  0.028  0.035  0.005  0.010 
#>     6  0.013  0.020  0.027  0.034  0.004  0.009 
#>     7  0.012  0.019  0.026  0.033  0.004  0.007 
#>     8  0.012  0.019  0.026  0.033  0.003  0.006 
#>     9  0.011  0.018  0.025  0.032  0.003  0.006 
#>    10  0.011  0.018  0.025  0.032  0.003  0.005 
#> -----------------------------------------------------

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

# 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:

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

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

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

# Small study (half the benchmark size)
mht_critical(J = 3, alpha_bar = 0.05, nm_ratio = 0.5)
#> 
#> ------------------------------------------------------------ 
#>   Optimal MHT Critical Values
#>   Viviano, Wuthrich, and Niehaus (2026)
#> ------------------------------------------------------------ 
#> 
#>   Cost model:          Linear (Eq. 27)
#>   Number of hypotheses: 3
#>   Benchmark alpha:      0.0500
#>   Sample size ratio:    0.50
#> 
#> ------------------------------------------------------------ 
#>   Optimal test size:    0.019010
#>   Optimal z-threshold:  2.0746
#> 
#>   Bonferroni size:      0.016667
#>   Bonferroni z-thresh:  2.1280
#>   Sidak size:           0.016952
#>   Sidak z-threshold:    2.1212
#>   Unadjusted size:      0.050000
#> ------------------------------------------------------------

# Benchmark size
mht_critical(J = 3, alpha_bar = 0.05, nm_ratio = 1.0)
#> 
#> ------------------------------------------------------------ 
#>   Optimal MHT Critical Values
#>   Viviano, Wuthrich, and Niehaus (2026)
#> ------------------------------------------------------------ 
#> 
#>   Cost model:          Linear (Eq. 27)
#>   Number of hypotheses: 3
#>   Benchmark alpha:      0.0500
#>   Sample size ratio:    1.00
#> 
#> ------------------------------------------------------------ 
#>   Optimal test size:    0.026042
#>   Optimal z-threshold:  1.9424
#> 
#>   Bonferroni size:      0.016667
#>   Bonferroni z-thresh:  2.1280
#>   Sidak size:           0.016952
#>   Sidak z-threshold:    2.1212
#>   Unadjusted size:      0.050000
#> ------------------------------------------------------------

# Large study (double the benchmark)
mht_critical(J = 3, alpha_bar = 0.05, nm_ratio = 2.0)
#> 
#> ------------------------------------------------------------ 
#>   Optimal MHT Critical Values
#>   Viviano, Wuthrich, and Niehaus (2026)
#> ------------------------------------------------------------ 
#> 
#>   Cost model:          Linear (Eq. 27)
#>   Number of hypotheses: 3
#>   Benchmark alpha:      0.0500
#>   Sample size ratio:    2.00
#> 
#> ------------------------------------------------------------ 
#>   Optimal test size:    0.040104
#>   Optimal z-threshold:  1.7495
#> 
#>   Bonferroni size:      0.016667
#>   Bonferroni z-thresh:  2.1280
#>   Sidak size:           0.016952
#>   Sidak z-threshold:    2.1212
#>   Unadjusted size:      0.050000
#> ------------------------------------------------------------

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.