## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7,
                      fig.height = 4.2)
library(POEMED)
set.seed(113)
# Every study below passes this tuning grid: a simulation study fits
# hundreds of models, and 20 values are several times faster per fit than
# the package default of 100 values from 0.05 to 10 (see ?pe_mediation).
grid20 <- seq(0.05, 2, length.out = 20)

## ----two-patterns-------------------------------------------------------------
lead <- pe_simulation_study(n = 200, p = 60, outcome = "continuous",
                            c1_grid = 1, n_rep = 20, lambda_grid = grid20)
lead

## ----size---------------------------------------------------------------------
size <- pe_power_curve(n = 150, p = 50, outcome = "continuous",
                       pattern = "contrasting", c1_grid = 0, n_rep = 40,
                       lambda_grid = grid20)
size[, c("rejection_hdmm", "rejection_pe")]

## ----draw-helper--------------------------------------------------------------
draw <- function(pc, main) {
  plot(pc$c1, pc$rejection_pe, type = "b", pch = 19, ylim = c(0, 1),
       xlab = expression(c[1]), ylab = "rejection rate", main = main)
  lines(pc$c1, pc$rejection_hdmm, type = "b", pch = 1, lty = 2)
  abline(h = 0.05, col = "grey60", lty = 3)
  legend("topleft", c("PE", "HDMM"), pch = c(19, 1), lty = c(1, 2), bty = "n")
}
grid <- c(0, 0.5, 1)

## ----dominate, fig.height = 4, fig.alt = "Two panels of rejection-rate curves against c1, for fully canceling and partially canceling heterogeneous mediation. PE is at or above HDMM throughout and far above it at c1 = 1 in both panels."----
p <- 50
# Mixed-sign mediators that only partially cancel: a small but nonzero total
# indirect effect, where the benchmark is weak and the enhancement is not.
alpha_mix <- c(1, -0.9, 0.8, -0.7, 0.6, -0.5, rep(0, p - 6))
tau_mix   <- c(rep(0.3, 6), rep(0, p - 6))

pc_cancel <- pe_power_curve(n = 150, p = p, outcome = "continuous",
                            pattern = "contrasting", c1_grid = grid, n_rep = 20,
                            lambda_grid = grid20)
pc_mixed  <- pe_power_curve(n = 150, p = p, outcome = "continuous",
                            c1_grid = grid, n_rep = 20, lambda_grid = grid20,
                            outcome_args = list(alpha_m = alpha_mix, tau = tau_mix))

op <- par(mfrow = c(1, 2))
draw(pc_cancel, "fully canceling")
draw(pc_mixed, "partially canceling")
par(op)

## ----homo-curve, fig.height = 4, fig.width = 4.2, fig.alt = "Rejection-rate curves under homogeneous mediation. Both curves rise with c1, and PE is at or above HDMM throughout, with the widest gap at the middle value, c1 = 0.5."----
pc_homo <- pe_power_curve(n = 150, p = p, outcome = "continuous",
                          pattern = "homogeneous", c1_grid = grid, n_rep = 20,
                          lambda_grid = grid20)
draw(pc_homo, "homogeneous")

## ----cross-outcome------------------------------------------------------------
at_c1 <- function(outcome, c2) {
  pc <- pe_power_curve(n = 180, p = 50, outcome = outcome,
                       pattern = "contrasting", c1_grid = 1, c2 = c2,
                       n_rep = 20, lambda_grid = grid20)
  c(HDMM = pc$rejection_hdmm, PE = pc$rejection_pe)
}
round(rbind(continuous = at_c1("continuous", 0.5),
            binary     = at_c1("binary", 1),
            count      = at_c1("count", 0.4)), 3)

## ----identification-----------------------------------------------------------
id <- pe_identification_study(n = 150, p = 60, outcome = "continuous",
                              pattern = "contrasting", c1_grid = c(0, 1),
                              n_rep = 20, methods = "Bonferroni",
                              lambda_grid = grid20)
id

## ----identification-long, eval = FALSE----------------------------------------
# # The longer runs quoted above (a seeded run with cores above 1 repeats
# # exactly at the same number of cores).
# long <- lapply(list(grid20, seq(0.2, 0.39, length.out = 20)), function(g)
#   pe_identification_study(n = 150, p = 60, outcome = "continuous",
#                           pattern = "contrasting", c1_grid = c(0, 1),
#                           n_rep = 400, methods = "Bonferroni",
#                           lambda_grid = g, cores = 4, seed = 113))

## ----who----------------------------------------------------------------------
imr <- WHO_mediation_analysis("imr", groupings = c("global", "region", "income"))
imr_summary <- summary(imr)
imr_summary

## ----external, eval = FALSE---------------------------------------------------
# d <- simulate_mediation_data(n = 300, p = 500, outcome = "continuous",
#                              pattern = "contrasting", c1 = 0.5)
# pe <- pe_mediation(d$X, d$Y, d$M, outcome = "continuous")
# 
# # HILMA (package 'freebird'): a debiased-lasso total-mediation-effect test.
# hilma_p <- freebird::hilma(d$Y, d$M, d$X)$pvalue_beta_hat
# 
# # GlobalTest (Bioconductor package 'globaltest'): a score test of the
# # mediator block against the outcome.
# gt_p <- globaltest::p.value(globaltest::gt(d$Y, d$M))
# 
# data.frame(method = c("HDMM", "PE-HDMM", "HILMA", "GlobalTest"),
#            pval = c(pe$value[pe$term == "pval_hdmm"],
#                     pe$value[pe$term == "pval_pe"], hilma_p, gt_p))

## ----session-info-------------------------------------------------------------
sessionInfo()