## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7,
                      fig.height = 4.2)
library(POEMED)
# Rejection rates quoted in the prose, as percentages. With 40 replications
# a rate such as 0.225 falls on a half percent (22.5), so the half is kept
# instead of being rounded to a whole number in one direction or the other.
pct <- function(rate) format(100 * rate, digits = 3)

## ----sim----------------------------------------------------------------------
set.seed(113)
d <- simulate_mediation_data(n = 200, p = 60, outcome = "continuous",
                             pattern = "contrasting", c1 = 1)

## ----sim-pieces---------------------------------------------------------------
dim(d$X)   # exposure: 200 observations, 1 column
length(d$Y)   # outcome: 200 values
dim(d$M)   # mediators: 200 x 60
d$active_mediators   # which mediators are truly active (the ground truth)
d$beta               # the total indirect effect: zero, up to rounding error

## ----first-fit----------------------------------------------------------------
fit <- pe_mediation(d$X, d$Y, d$M, outcome = "continuous")
fit

## ----finer-grid---------------------------------------------------------------
finer <- pe_mediation(d$X, d$Y, d$M, outcome = "continuous",
                      lambda_grid = seq(0.01, 10, length.out = 200))
attr(finer, "tuning")$lambda_selected
attr(finer, "active_mediators")

## ----read---------------------------------------------------------------------
fit$value[fit$term == "pval_pe"]   # the exact p-value
attr(fit, "active_mediators")      # indices of the active mediators identified

## ----methods, warning = FALSE-------------------------------------------------
des_idx <- function(f) attr(f, "active_mediators")
des_idx(pe_mediation(d$X, d$Y, d$M, outcome = "continuous", method = "Bonferroni"))
des_idx(pe_mediation(d$X, d$Y, d$M, outcome = "continuous", method = "BH"))

## ----error-level, warning = FALSE---------------------------------------------
strict <- pe_mediation(d$X, d$Y, d$M, outcome = "continuous",
                       error_level = 0.01)
strict[strict$term %in% c("pval_hdmm", "j_pe", "pval_pe"), ]
des_idx(pe_mediation(d$X, d$Y, d$M, outcome = "continuous", method = "BH",
                     error_level = 0.01))

## ----glm, warning = FALSE-----------------------------------------------------
set.seed(113)
db <- simulate_mediation_data(n = 250, p = 50, outcome = "binary",
                              pattern = "contrasting", c1 = 1, c2 = 1)
pe_mediation(db$X, db$Y, db$M, outcome = "binary")

## ----power, fig.alt = "Rejection-rate curves for the benchmark and power-enhanced tests"----
set.seed(113)
pc <- pe_power_curve(n = 150, p = 50, outcome = "continuous",
                     pattern = "contrasting", c1_grid = c(0, 0.5, 1),
                     n_rep = 40, lambda_grid = seq(0.05, 2, length.out = 20))
pc

plot(pc, main = "Contrasting Mediation: PE vs. Benchmark")

## ----who-data-----------------------------------------------------------------
data(WHO_health_mediation)
dim(WHO_health_mediation)
WHO_health_mediation[1:3, c("country", "region", "income", "year",
                            "gdp_growth", "imr", "gge_gdp")]

## ----codebook-----------------------------------------------------------------
head(WHO_indicator_codebook, 4)

## ----who-design, warning = FALSE----------------------------------------------
des <- WHO_mediation_design("imr")          # infant mortality rate
c(observations = des$n, mediators = ncol(des$M), confounders = ncol(des$Z))

fit <- pe_mediation(scale(des$X), des$Y - mean(des$Y), scale(des$M),
                    Z = des$Z, outcome = "continuous", scale = FALSE,
                    lambda_grid = seq(0.1, 2, length.out = 100))
identified <- des$indicators[attr(fit, "active_mediators")]
WHO_indicator_codebook$description[
  match(identified, WHO_indicator_codebook$indicator)]

## ----who-analysis-------------------------------------------------------------
res <- WHO_mediation_analysis("imr", groupings = c("global", "region"))
res[, c("group", "n_countries", "pval_hdmm", "pval_pe", "active_mediators")]

## ----who-summary--------------------------------------------------------------
summary(res)

## ----who-raw, warning = FALSE-------------------------------------------------
raw <- pe_mediation(des$X, des$Y, des$M, Z = des$Z, outcome = "continuous",
                    lambda_grid = seq(0.1, 2, length.out = 100))
c(confounders_unscaled = fit$value[fit$term == "pval_hdmm"],
  all_standardized = raw$value[raw$term == "pval_hdmm"])
des$indicators[attr(raw, "active_mediators")]

## ----formula, warning = FALSE-------------------------------------------------
set.seed(113)
d <- simulate_mediation_data(n = 200, p = 60, outcome = "continuous",
                             pattern = "contrasting", c1 = 1)
df <- data.frame(y = d$Y, x = d$X[, 1], d$M)
meds <- grep("^X", names(df), value = TRUE)   # the mediator columns
fit <- pe_mediate(y ~ x, data = df, mediators = meds, outcome = "continuous")

## ----mediators----------------------------------------------------------------
pe_mediators(fit)

## ----plot-mediators, fig.alt = "Per-mediator screening evidence for the candidate mediators the penalized fit kept"----
plot(fit)

## ----broom--------------------------------------------------------------------
glance(fit)

## ----selection, warning = FALSE-----------------------------------------------
fit_all <- pe_mediation(d$X, d$Y, d$M, outcome = "continuous",
                        report_all_methods = TRUE)
pe_selection(fit_all)

## ----ssplan, eval = FALSE-----------------------------------------------------
# plan <- ss_power_pe_mediation(outcome = "continuous", pattern = "contrasting",
#                               p = 50, n_grid = c(100, 150, 200, 250),
#                               target_power = 0.8, n_rep = 100, seed = 113)
# plan
# attr(plan, "recommended_n")

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