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.