## ----setup, include = FALSE--------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----eval=FALSE--------------------------------------------------------------- # # From CRAN (once published) # install.packages("mhtopt") # # # Development version from GitHub # # remotes::install_github("dviviano/mhtopt", subdir = "r/mhtopt") ## ----------------------------------------------------------------------------- library(mhtopt) ## ----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) ## ----------------------------------------------------------------------------- # Default: cf_share = 0.46, J_bar = 3 (from Sertkaya et al. 2016) mht_critical(J = 5, alpha_bar = 0.05) ## ----------------------------------------------------------------------------- # beta = 0.13, iota = 0.075 (from J-PAL data, Appendix A) mht_critical(J = 5, alpha_bar = 0.05, model = "cobbdouglas") ## ----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) ## ----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")]) ## ----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) ## ----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 ## ----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) ## ----------------------------------------------------------------------------- # 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) ## ----------------------------------------------------------------------------- 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) ## ----------------------------------------------------------------------------- zstats <- c(3.2, 2.1, 1.8, 1.3, 0.8, -0.2) mht_test(z = zstats, alpha_bar = 0.05) ## ----------------------------------------------------------------------------- # 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)) ## ----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) ## ----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) ## ----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))]) ## ----eval=FALSE--------------------------------------------------------------- # mht_est(fit, # vars = c("treat1", "treat2", "treat3"), # alpha_bar = 0.05, # onesided = FALSE) ## ----------------------------------------------------------------------------- # 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)