
POEMED implements powerful global tests for high-dimensional mediation analysis, for continuous, binary, and count outcomes. It answers the question “is there any active mediator among a large set of candidate mediators?” and, when the answer is yes, reports which individual mediators its screen identifies as active. The name comes from POwer-Enhanced MEDiation and is pronounced “POE-med”.
The conventional high-dimensional mediation test targets the
total indirect effect
beta = Gamma_x' alpha_m, the sum of every mediator’s
individual indirect effect. When some indirect effects are positive and
others negative they can cancel, making beta = 0 even
though active mediators exist, and a test built on beta is
then powerless. The power-enhanced (PE) test adds a component
J_m that accumulates the marginal signal from each
individual mediator. Because J_m sums magnitudes,
opposite-signed effects reinforce rather than cancel, so the test stays
powerful under heterogeneous and contrasting mediation. Under the global
null the enhancement is built to vanish: as the sample size and the
number of candidate mediators grow, the test’s Type I error rate
approaches the nominal level, as the benchmark’s does, and in finite
samples it rejects slightly more often than the benchmark.
POEMED answers the global question, “is there any active mediator
among many candidates?”, and then names the mediators it identifies
as active. Reach for it when the candidate mediators number in the
dozens or more (p may exceed n), the outcome
is continuous, binary, or a count, and the individual mediation effects
may differ in sign. For a single mediator, or a few mediators fit as one
structural model (indirect effects with confidence intervals, likelihood
ratio tests of arbitrary indirect effects, moderated mediation), use DMAR instead. Among
high-dimensional global tests, POEMED implements the power-enhanced
tests of Yu and Kelley (in press) and the HDMM benchmark they extend;
HILMA, GlobalTest, HDMT, and DACT are separate methods in their own
packages. The vignette poemed-vs-competitors shows how to
run HILMA and GlobalTest beside a POEMED fit and where to find HDMT and
DACT. The guided-tour vignette and ?pe_mediation state the
method’s practical limits (a sample size small relative to
log(p) and to the signal, strongly correlated mediators,
the screen’s asymptotic error rate control in finite samples, rare
binary outcomes, overdispersed counts, and missing values).
Once POEMED is on CRAN:
install.packages("POEMED")The development version, from GitHub:
# install.packages("remotes")
remotes::install_github("yelleKneK/POEMED", build_vignettes = TRUE)POEMED imports glmnet and ncvreg;
knitr and rmarkdown build the vignettes. After
installing, browseVignettes("POEMED") opens the guides.
library(POEMED)
# A contrasting setting: active mediators whose indirect effects cancel,
# so the total indirect effect is zero.
set.seed(113)
d <- simulate_mediation_data(n = 200, p = 60, outcome = "continuous",
pattern = "contrasting", c1 = 1)
# One call reports the benchmark Wald test and the power-enhanced test
# side by side, and names the active mediators the screen identified.
pe_mediation(d$X, d$Y, d$M, outcome = "continuous")The benchmark test (pval_hdmm) does not reject the
(zero) total indirect effect; the power-enhanced test
(pval_pe) rejects decisively, and the screen identifies
three of the four truly active mediators as active (the footer reads
“Active mediators identified (3): 2, 3, 4”). The footer also reports the
tuning parameter the high-dimensional BIC chose and the grid it was
chosen from, seq(0.05, 10, length.out = 100) by default for
every outcome type; the grid can be shortened to save time or extended
when the value chosen sits at an end of it, which the fit warns about
(see ?pe_mediation). Run as printed, this fit chooses 0.05,
the smallest value, and warns; a grid extended below it and spaced more
finely (seq(0.01, 10, length.out = 200)) chooses an
interior value and identifies the same three mediators as active.
All data needed to rerun the article’s analyses ship with the
package; nothing is downloaded at run time. The vignette
reproducing-the-article reruns the
empirical analysis cell by cell; the table below maps each article
result to the function that runs its analysis.
| Article result | Function |
|---|---|
| Section 4.1, Figure 3: linear-model size and power | pe_power_curve(outcome = "continuous", pattern = "homogeneous" / "contrasting") |
| Section 4.2, Figures 4 and 5: logistic and Poisson size and power | pe_power_curve(outcome = "binary", c2 = 1) and
pe_power_curve(outcome = "count", c2 = 0.4) |
| Supplement Section S.4.1, Tables S.1 to S.3: FWER/FDR mediator-identification study | pe_identification_study() |
| Supplement Section S.4.4: real-data-motivated heterogeneous setting | simulate_guo_mediation() with the shipped
guo_calibration constants |
| Section 5, Tables 1 (IMR) and 2 (PoU); supplement Tables S.10 (LEB), S.11 (U5MR), and S.12 (LBW) | WHO_mediation_analysis("imr"), "pou",
"leb", "u5mr", "lbw" |
| Supplement Section S.5, Tables S.8 to S.12: PE, PE-BH, and PE-BY columns | WHO_mediation_analysis(..., full_table = TRUE) |
| Where PE detects mediation the benchmark misses | summary() of a WHO_mediation_analysis()
result |
The Monte Carlo studies use the article’s n = 300,
p = 500, and 1,000 replications (the supplement’s
real-data-motivated setting its case-study n = 85 and
p = 1008), and each panel takes hours; set
cores > 1 to parallelize (near-linear). The real-data
analyses are deterministic: under its default tuning grid,
WHO_mediation_analysis() reproduces the active mediators
identified in Tables 1, 2, and S.8 to S.12 (every identified set but one
BH set) and their benchmark p-values to the printed four decimals, apart
from a few cells that differ in the fourth decimal or below 0.0001 and
six cells it does not fit: four in which two indicators are exactly
collinear and two in which the outcome is constant. The reproduction
vignette names those six cells, shows how to fit the four collinear ones
by hand, and shows why the two constant-outcome cells cannot be fit. The
package’s test suite pins every global, region, and income row of Tables
1 and 2 (tests/testthat/test-benchmark.R) and checks every
populated cell of Tables 1, 2, and S.8 to S.12 against the printed
values (tests/testthat/test-who_article_tables.R).
# The global, region, and income models of the article's infant-mortality
# analysis (Table 1); add "region_income" to groupings for the
# region-by-income cells.
WHO_mediation_analysis("imr", groupings = c("global", "region", "income"))browseVignettes("POEMED")POEMED: a guided tour of the method
and the package.reproducing-the-article: the empirical
analysis rerun cell by cell.poemed-vs-competitors: a head-to-head
against the benchmark test.| Function | Purpose |
|---|---|
pe_mediation() |
Front end; dispatches on outcome type (continuous / binary / count) |
pe_mediation_linear() / _logistic() /
_poisson() |
Outcome-specific workers |
pe_mediate() |
Data-frame / formula front end |
pe_mediators() / pe_selection() |
Per-mediator screen detail; PE statistic and identified active set per multiplicity method |
plot() / tidy() / glance() /
summary() |
Plot, broom verbs, and a cross-group comparison digest |
simulate_mediation_data() |
Generate data under homogeneous / contrasting patterns |
simulate_guo_mediation() /
guo_calibration |
The real-data-motivated heterogeneous setting and its constants |
pe_power_curve() /
pe_simulation_study() |
Monte Carlo size and power study |
pe_identification_study() |
Monte Carlo FWER / FDR / precision / recall of mediator identification |
ss_power_pe_mediation() |
Sample-size planning for a target power |
WHO_mediation_design() /
WHO_mediation_analysis() |
Build / run the benchmark health-expenditure analysis |
Every test reports the benchmark Wald test and the power-enhanced
test side by side, plus the estimated total indirect effect. No
confidence interval is reported: the article reports none, and the Wald
interval built on the penalized estimate under-covers in finite samples.
report_all_methods = TRUE adds the identified active set
and the PE p-value under all three multiplicity methods,
drop_constant = TRUE tolerates zero-variance mediators, and
cores parallelizes the Monte Carlo and benchmark routines.
For a continuous outcome, pe_mediation_linear() also takes
lambda_grid_reduced, the tuning grid for the reduced model
it fits beside the full one (by default the same grid as
lambda_grid).
WHO_health_mediation: the empirical panel analyzed in
the article: 91 World Health Organization member states over 2000–2021,
with GDP per capita growth as the exposure, 57 WHO Global Health
Expenditure Database indicators as candidate mediators, and five World
Bank health outcomes (infant and under-five mortality, life expectancy,
low birthweight, undernourishment).WHO_indicator_codebook: definitions of the 57
indicators.guo_calibration: constants for the real-data-motivated
heterogeneous simulation.example_continuous, example_binary,
example_count: small simulated data sets, one per outcome
type, each with three confounders and stored with the truth that
generated it (the truly active mediators and every coefficient).Yu, X., & Kelley, K. (in press). Power Enhancement in High-Dimensional Heterogeneous Mediation Analysis. Journal of the American Statistical Association.
If you use POEMED in published work, please cite the article. Run
citation("POEMED") for the BibTeX entry.
GPL (>= 3).