This is the first release of POEMED (Power Enhancement for High-Dimensional Mediation Analysis), the companion package to Yu and Kelley (in press), “Power Enhancement in High-Dimensional Heterogeneous Mediation Analysis,” in the Journal of the American Statistical Association. The name comes from POwer-Enhanced MEDiation and is pronounced “POE-med”.
POEMED asks whether any of a large set of candidate mediators carries the effect of an exposure to an outcome and, when the answer is yes, which ones do. The benchmark Wald test of Guo et al. (2022, 2023, 2024) tests the total indirect effect, the sum of the individual indirect effects, so it loses its power when active mediators act in opposite directions and their effects cancel. The power-enhanced test of Yu and Kelley (in press), which follows the power enhancement principle of Fan, Liao, and Yao (2015), adds a component that accumulates the magnitudes of the marginal signals of the mediators that pass a significance screen, so it detects active mediators whether their effects are homogeneous, heterogeneous, or exactly contrasting. The package works for continuous, binary, and count outcomes.
pe_mediation() is the front end. Give it the exposure,
the outcome, the candidate mediators, and any confounders, name the
outcome type ("continuous", "binary", or
"count"), and it reports the benchmark Wald test and the
power-enhanced test side by side, the estimated total indirect effect,
and the active mediators the screen identified. One or several exposures
may be analyzed.pe_mediation_linear(),
pe_mediation_logistic(), and
pe_mediation_poisson() are the outcome-specific workers
that pe_mediation() dispatches to, and they can be called
directly.pe_mediate() is a formula and data-frame front end for
those who prefer naming columns to building matrices. It analyzes the
complete cases and reports how many rows it dropped and which columns
held the missing values; pe_mediation() and its workers
take complete data.?pe_mediation_logistic).?pe_mediation
states the method’s assumptions, describes its practical limits as seen
on simulated data (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), and points to other tools for
a single mediator or a few. The same page states that the power-enhanced
test controls its size asymptotically and that, in finite samples, it
can reject somewhat more often than the benchmark, by up to about
error_level / log(log(n)).example_continuous, example_binary, and
example_count are small data sets simulated with
simulate_mediation_data(), one for each outcome type, each
with three confounders, stored with the truth that generated it (the
truly active mediators and every coefficient, the confounders’
included), and ready for a first fit (see
?example_data).error_level, with
familywise error rate control (method = "Bonferroni", the
default) or false discovery rate control (method = "BH" for
Benjamini-Hochberg or method = "BY" for
Benjamini-Yekutieli). Like the size of the test, this control holds
asymptotically."active_mediators" attribute, counts them in the row
n_active_mediators (which glance() returns as
a column), and lists them in the print footer, as “Active mediators
identified (k): …” or “Active mediators identified: none”.pe_mediators() returns the detail behind a fit for each
mediator the penalized fit selected: its two path statistics (mediator
to outcome and exposure to mediator), its screening p-value, and, in the
logical column active, whether the screen identified it as
active.pe_selection() reports the power-enhanced statistic,
its enhancement component J_m, the global p-value, and the
identified active set (the columns n_active and
active_mediators) under the method the fit used, or under
all three multiplicity methods side by side when the fit was made with
report_all_methods = TRUE.seq(0.05, 10, length.out = 100) for every outcome family,
is the default value of lambda_grid, and any other grid can
be passed in its place. HBIC keeps the best value the grid offers, so
the grid can be suited to the task: shortened to save time when many
fits are run, extended or refined when the chosen value sits at an end
of the grid.pe_mediation_linear() takes the grid for that model as
lambda_grid_reduced, which defaults to
lambda_grid and is validated in the same way. The argument
belongs to that worker alone: pe_mediation() does not
accept it, so call pe_mediation_linear() directly to set
it."tuning" attribute, a
list with lambda_grid (the grid searched),
lambda_selected (the value HBIC selected), and
at_lower_end and at_upper_end (whether that
value is the smallest or the largest of the grid); a continuous fit also
records lambda_grid_reduced and
lambda_selected_reduced for the reduced model. The printed
fit reports the HBIC-selected value and the grid it came from.poemed_grid_boundary). For a continuous
outcome the warning says which model’s choice sat there, the full
model’s, the reduced model’s, or both; the reduced model’s choice is
tested against the ends of its own grid, and the warning names
lambda_grid_reduced when that grid is not the full model’s.
The warning suggests extending the grid past that end, using a finer
partition, and comparing the active mediators identified and the
p-values across grids; when only the reduced model’s choice sits at an
end of a separate lambda_grid_reduced, it says in place of
the comparison that the choice is recorded in the "tuning"
attribute and changes no reported statistic.poemed_empty_fit), says so when it prints, and is reported
as a non-rejection (both tests 0, with p-values of 1), which is expected
when no mediator is active but is not evidence that none is.simulate_mediation_data() generates data under the
homogeneous and contrasting mediation patterns of the article for all
three outcome types, and mediation_ar1_cov() builds the
autoregressive mediator covariance the simulations use.tau is 0.1, 0.2, 0.3, 0.4,
and 0.5 on the first five mediators and noise on the rest, drawn with
standard deviation 0.1 for a continuous outcome and 0.01 for a binary or
count outcome. With confounders (d > 0), their
coefficients on the outcome, alpha_z, and on the mediators,
Gamma_z, are drawn from U(-0.5, 0.5) unless supplied;
Gamma_z is drawn as d values recycled across
the p mediator columns, so each confounder loads equally on
every mediator. For a count outcome, when the largest value of the
log-mean exceeds 5, the whole log-mean vector is mapped affinely onto
the interval from -5 to 5 before the counts are drawn.alpha_z and Gamma_z among
them (NULL when d = 0), the total indirect
effect, and the truly active mediators
(active_mediators).simulate_guo_mediation(), with the
guo_calibration constants, generates the article’s
real-data-motivated heterogeneous setting: a continuous outcome with
eleven active mediators of mixed sign, calibrated to the DNA-methylation
study of Guo et al. (2022).pe_power_curve() runs the Monte Carlo size and power
study across a grid of signal strengths, and
pe_simulation_study() runs it across mediation
patterns.pe_identification_study() estimates how well the screen
recovers the truly active mediators: the familywise error rate, false
discovery rate, precision, and recall of the identified active set under
each multiplicity method.ss_power_pe_mediation() plans the sample size that
reaches a target power by simulation.cores
on Unix-like systems such as macOS and Linux (serially on Windows),
report the number of empty fits at each design point, and take a
seed that reproduces a run exactly at the same value of
cores while leaving the caller’s random-number state as it
was.WHO_health_mediation holds the article’s empirical
data, a panel of 91 World Health Organization members from 2000 to 2021
with economic growth as the exposure, 57 health-expenditure indicators
as candidate mediators, and five health outcomes;
WHO_indicator_codebook defines the indicators.WHO_mediation_design() assembles the exposure, outcome,
mediator, and confounder matrices for any outcome and any region or
income group.WHO_mediation_analysis() runs the article’s subgroup
analyses (global, by region, by income group, and by region and income
group) with the preprocessing and the tuning grid,
seq(0.1, 2, length.out = 100), behind the article’s Tables
1 and 2 and supplement Tables S.8 to S.12. It reproduces those tables
apart from a few cells that ?WHO_mediation_analysis lists
and explains. Six cells that the tables fill come back as
NA: the reproduction vignette shows how to fit by hand the
four in which two indicators are exactly collinear, and why the two in
which the outcome is constant cannot be fit. Each row gives the number
of active mediators identified and their names (the columns
n_active and active_mediators); with
full_table = TRUE it gives the power-enhanced p-value and
the active mediators identified under each multiplicity method (the
columns active_bonferroni, active_bh, and
active_by).pe_mediation(), its workers, or pe_mediate()
is a long table with a term column and a numeric
value column; the tables from pe_mediators(),
pe_selection(), the Monte Carlo and planning functions, and
WHO_mediation_analysis() are wide, with one row per
mediator, multiplicity method, design point, sample size, or group.
Values are stored at full precision and rounded only for display: whole
numbers print without decimals, p-values print to four decimals with
< 0.0001 below that, and other values print to four
significant figures (options(poemed.digits = 6) shows
more).M, the
printed fit and the mediator tables name the mediators, and with several
named exposures the per-exposure rows carry the exposures’ names.plot() draws the power curves of a
pe_power_curve() result, the benchmark and power-enhanced
rejection rates against signal strength, and, for a fit, the screening
evidence for each mediator the penalized fit selected, with the ones
identified as active highlighted (the legend reads “active” and
“screened out”). summary() of a grouped table such as the
output of WHO_mediation_analysis() counts the groups in
which each test detects mediation, lists those in which the
power-enhanced test detects mediation that the benchmark misses, and
names the mediators identified as active most often. The broom verbs
tidy() and glance() are re-exported so they
work after library(POEMED) alone.vignette("POEMED"); a head-to-head comparison of the
power-enhanced test with the benchmark Wald test across mediation
patterns, outcome types, and the WHO data, closing with where to find
HILMA, GlobalTest, HDMT, and DACT,
vignette("poemed-vs-competitors"); and a guide to rerunning
the article’s empirical analysis,
vignette("reproducing-the-article").
browseVignettes("POEMED") lists them.remotes::install_github("yelleKneK/POEMED", build_vignettes = TRUE);
without them, vignette() finds nothing to open. All three
can also be read on the package website, https://yelleknek.github.io/POEMED/, where the guided
tour is the “Get started” page and the other two are under
Articles.seed argument reproduces its results and leaves the
caller’s random-number state as it was.