Power-Enhanced Test for High-Dimensional Mediation

Tests whether any mediator is active among a large, possibly intercorrelated set of candidate mediators, for a continuous, binary, or count outcome, with the power-enhanced (PE) test of Yu and Kelley (in press). This is the main entry point of POEMED: name the outcome type and it dispatches to the appropriate model, returning one tidy table that reports the benchmark Wald test on the total indirect effect side by side with the PE test, plus the individual mediators the screen identified as active.

Usage

pe_mediation(
  X,
  Y,
  M,
  Z = NULL,
  outcome = c("continuous", "binary", "count"),
  method = c("Bonferroni", "BH", "BY"),
  scale = TRUE,
  error_level = 0.05,
  report_all_methods = FALSE,
  drop_constant = FALSE,
  lambda_grid = seq(0.05, 10, length.out = 100)
)

Arguments

X: Numeric matrix of exposures with $n$ rows and $q$ columns ($q$ is usually 1). A numeric vector is treated as a single exposure.
Y: Numeric outcome vector of length $n$. Continuous for outcome = "continuous", coded 0/1 for "binary", and nonnegative integer counts for "count".
M: Numeric matrix of candidate mediators with $n$ rows and $p$ columns; $p$ may exceed $n$. Column names, when present, label the mediators in the print footer and the mediator table; without them mediators are reported by column position.
Z: Optional numeric matrix of confounders with $n$ rows. Default NULL (no confounders).
outcome: Outcome type: "continuous" (linear model), "binary" (logistic model), or "count" (Poisson log-link model).
method: Multiplicity adjustment used to screen individual mediators in the PE component: "Bonferroni" (familywise error rate control, the default), "BH", or "BY" (false discovery rate control by the procedures of Benjamini and Hochberg, 1995, and of Benjamini and Yekutieli, 2001; "BH" assumes independence or positive dependence among mediators, "BY" is valid under arbitrary dependence). See Details for the screening threshold each one uses.
scale: Logical; if TRUE (default) the exposures, mediators, and confounders are each standardized to mean 0 and standard deviation 1 before fitting, as in the article, and a continuous outcome is centered but not rescaled. The total indirect effect is then the change in the outcome per standard deviation of each exposure: in the units of Y for a continuous outcome, in log odds for a binary outcome, and in log mean count for a count outcome. With FALSE the inputs are used exactly as supplied, and the total indirect effect is per unit of each exposure. The continuous-outcome fit has no intercept, and the default tuning grid was chosen for standardized columns, so FALSE is meant for inputs that are already standardized (and, for a continuous outcome, a Y that is already centered); on raw inputs it fits a different model and can even change the sign of the estimate. Because a continuous Y is not rescaled, the tuning grid is not expressed in the outcome's units, so the active mediators identified, both p-values, and the verdict depend on the units of Y: the same data with Y multiplied by 10 can identify a different set of mediators.
error_level: Target error rate for the mediator-screening step, in $(0, 1)$. Interpreted as the familywise error rate when method = "Bonferroni" and as the false discovery rate when method = "BH" or "BY". Default 0.05. This is distinct from the significance level used to test the global null (which is the user's choice when reading pval_pe). The PE test is asymptotically valid for any fixed value in $(0, 1)$, but in finite samples a larger value lets the screen pass more spurious mediators under the null and raises the PE test's size further above the benchmark's (see Details).
report_all_methods: Logical; if TRUE, the identified active set is computed for all three multiplicity methods (Bonferroni, BH, BY) off the same fit and recorded for comparison, retrievable with pe_selection(). The method argument still drives the primary result. Default FALSE.
drop_constant: Logical; how to handle constant (zero-variance) mediator columns, which cannot be standardized or carry signal. If FALSE (the default) they are an error; if TRUE they are dropped with a warning and the remaining mediators are renumbered back to their original column positions in the identified active set. Useful for small subgroups where some mediators happen to be constant.
lambda_grid: Numeric vector of candidate tuning parameters for the SCAD penalty (Fan and Li, 2001) in the penalized mediator fit, which for a binary outcome is adaptively rescaled (see pe_mediation_logistic()). The fit is repeated at each value and the one minimizing the high-dimensional BIC (HBIC) of Wang, Kim, and Li (2013) is kept: whatever the grid, the value selected is the best of the values offered, since HBIC is minimized over the grid and not over every positive lambda. The default, seq(0.05, 10, length.out = 100) for every outcome family, is written out in the signature so that it can be edited in place. Every value is fit (the search is not adaptive), so the cost of a fit grows with the length of the grid. The grid can be tuned to the situation. Fewer values, or a range narrowed to where earlier fits put the HBIC minimum, make each fit faster, which matters when many fits are run (a simulation study, for example); a wider or finer grid is worth trying when the selected value is the smallest or largest value of the grid, in which case a non-empty fit warns (class poemed_grid_boundary) and the print footer, which reports the full model's choice, says so. A continuous outcome also fits a reduced model, the outcome on the mediators and confounders without the exposure, tuned on the same grid by default (see lambda_grid_reduced in pe_mediation_linear()); its choice is lambda_selected_reduced in the "tuning" attribute, and the warning covers it as well, naming whether the full model, the reduced model, or both sat at the end. For a continuous outcome the grid is fixed while Y is centered but not rescaled, so the value HBIC chooses, and with it the active mediators identified, depends on the units of Y (see scale). The grid searched and the value chosen are reported in the "tuning" attribute and the print footer.

Value

A tidy data.frame of class poemed_tbl. See pe_mediation_linear() for the row schema and the attributes. The active mediators the screen identified (see Details for how they relate to the truly active mediators and to the penalized-selected candidates) are in the "active_mediators" attribute and the print footer; the tuning parameter HBIC chose is in the "tuning" attribute and the footer. The total indirect effect is on the scale described under scale.

Details

Standard high-dimensional mediation tests target 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 benchmark reported here is such a test, the Wald test $S_n$ on $\beta$. For a continuous outcome it is the test of Guo et al. (2022), which extends the partially penalized Wald test of Guo et al. (2023) to observed confounders. For binary and count outcomes it is the test of Guo et al. (2024).

The PE test of Yu and Kelley (in press) follows the power enhancement principle of Fan, Liao, and Yao (2015). It adds to $S_n$ the component (equation 2.9 of the article) $$J_m = \sqrt{p} \sum_{i=1}^{q} \sum_{j \in \hat{S}} \left| \frac{\hat\alpha_{m,j}}{\hat\sigma_{m,j}} \right| \left| \frac{\hat\Gamma_{x,i,j}}{\hat\sigma_{\Gamma,i,j}} \right| \mathbf{1}\!\left\{ \max(p_{1,j}, p_{2,i,j}) < \frac{\alpha_{\mathrm{lvl}}}{s q \log\log n} \right\},$$ a sum over the $s$ penalized-selected candidates $\hat S$ of the product of the (standardized) mediator-on-outcome and exposure-on-mediator statistics, kept only for pairs where both paths are individually significant; $\alpha_{\mathrm{lvl}}$ is error_level. That indicator is the Bonferroni screen used by method = "Bonferroni". Under method = "BH" or "BY" the indicator is instead $\mathbf{1}\{\tilde p_{i,j} < \alpha_{\mathrm{lvl}} / \log\log n\}$, where $\tilde p_{i,j}$ is the Benjamini and Hochberg (1995) or Benjamini and Yekutieli (2001) adjusted value of $\max(p_{1,j}, p_{2,i,j})$ over the $s q$ pairs (equations S.22 to S.25 of the article's supplement). Because $J_m$ accumulates magnitudes, opposite-signed indirect effects reinforce rather than cancel, so the test $M_{PE} = S_n + J_m$ stays powerful under heterogeneous and contrasting mediation. The indicator also identifies individual active mediators, with asymptotic familywise error rate control under method = "Bonferroni" or false discovery rate control under "BH" / "BY". In this documentation an active mediator is one whose effect is nonzero (the article's usage, as in the global null of no active mediator). The screen identifies active mediators, the candidates it finds statistically significant under the multiplicity method, and the fit reports them as its active mediators (the "active_mediators" attribute): the identified set is the method's estimate of the truly active set. Selected refers to the first stage, that is, to the candidates the penalized fit keeps (the penalized-selected candidates $\hat S$) and to the tuning parameter HBIC selects.

Both statistics are referred to a chi-square distribution with $q$ degrees of freedom. Under the global null of no active mediator, $J_m$ is zero with probability tending to one, so $M_{PE}$ has the same chi-square limit as $S_n$ and the PE test is asymptotically valid for any fixed error_level, by Theorem 1 of Yu and Kelley (in press). That guarantee is asymptotic. In finite samples, under a null in which some mediators affect the outcome but the exposure affects none of them, the screen passes one of those mediators with probability of about error_level / log(log(n)), which is 0.030 at $n = 200$ with error_level = 0.05, and $J_m$ is then large enough to reject. The PE test's size can therefore exceed the benchmark's by up to about that amount. For example, pe_power_curve(200, 100, pattern = "homogeneous", c1_grid = 0, n_rep = 2000, seed = 113) gave rejection rates at the 0.05 level of 0.0475 for the benchmark and 0.0620 for the PE test over 2,000 null samples, and 0.0510 for the PE test with error_level = 0.01 (Monte Carlo standard errors about 0.005). A smaller error_level narrows the gap, at some cost in identifying active mediators.

When to Use POEMED

POEMED is for the global question, “is there any active mediator among many candidates?”, followed by the identification of individual active mediators. It assumes a linear or generalized linear outcome model, a sparse set of active mediators, standardized inputs, and complete data. Some practical limits, seen on simulated data:

Sample size. The penalized selection needs $n$ large relative to $\log p$ and to the signal. When it is not, as at $n = 100$ with $p = 2000$ under the contrasting pattern, the selection typically keeps one or two mediators, the screen identifies one of them or none, and the test has little left to draw on.
Strongly correlated mediators. Correlation among the mediators weakens both the selection and the screen: the screen can come up empty, and the identified active set varies from sample to sample. An autoregressive correlation of 0.9 among neighboring mediators is enough to see this; 0.5 is not.
Active mediators identified in finite samples. The screen's control of its error rate is asymptotic, so the identified active set can hold a mediator that is not truly active, the more so the stronger the signal. The rate depends on the tuning grid, and a grid that leaves out the smallest values of lambda lowers it at the cost of recall.
Rare binary outcomes. With only a few percent events the penalized logistic fit that HBIC selects is empty, and the function says so.
Overdispersed counts. The count model is Poisson. Under negative binomial counts its size held up, but a share of the fits came back empty. A count fit can also take several times as long as a continuous fit of the same design.
Missing values are not handled; supply complete cases.

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), the DMAR package is the tool; POEMED's contribution is the high-dimensional global test. The other high-dimensional methods the article compares against are HILMA (Zhou, Wang, & Zhao, 2020), GlobalTest (Djordjilovic et al., 2019), HDMT (Dai, Stanford, & LeBlanc, 2022), and DACT (Liu et al., 2022). POEMED does not depend on them. HDMT is on CRAN, and GlobalTest is the Bioconductor package globaltest. The CRAN package freebird, which implemented HILMA, was archived in July 2026, and the CRAN package named DACT is an unrelated clinical-trials package, not the method of Liu et al. The package implementing the DACT method of Liu et al. is also named DACT, but it is not on CRAN and is instead available through their GitHub repository. The vignette poemed-vs-competitors shows how to call HILMA and GlobalTest beside a POEMED fit.

How to Cite

If you use POEMED in published work, please cite Yu and Kelley (in press), the article that introduces its methods and the package. citation("POEMED") gives the full reference and a BibTeX entry.

References

Yu, X., & Kelley, K. (in press). Power Enhancement in High-Dimensional Heterogeneous Mediation Analysis. Journal of the American Statistical Association.

Fan, J., Liao, Y., & Yao, J. (2015). Power enhancement in high-dimensional cross-sectional tests. Econometrica, 83(4), 1497–1541. doi:10.3982/ECTA12749

Guo, X., Li, R., Liu, J., & Zeng, M. (2022). High-dimensional mediation analysis for selecting DNA methylation loci mediating childhood trauma and cortisol stress reactivity. Journal of the American Statistical Association, 117(539), 1110–1121. doi:10.1080/01621459.2022.2053136

Guo, X., Li, R., Liu, J., & Zeng, M. (2023). Statistical inference for linear mediation models with high-dimensional mediators and application to studying stock reaction to COVID-19 pandemic. Journal of Econometrics, 235(1), 166–179. doi:10.1016/j.jeconom.2022.03.001

Guo, X., Li, R., Liu, J., & Zeng, M. (2024). Estimations and tests for generalized mediation models with high-dimensional potential mediators. Journal of Business & Economic Statistics, 42(1), 243–256. doi:10.1080/07350015.2023.2174548

Fan, J., & Li, R. (2001). Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the American Statistical Association, 96(456), 1348–1360. doi:10.1198/016214501753382273

Wang, L., Kim, Y., & Li, R. (2013). Calibrating nonconvex penalized regression in ultra-high dimension. The Annals of Statistics, 41(5), 2505–2536. doi:10.1214/13-AOS1159

Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate: A practical and powerful approach to multiple testing. Journal of the Royal Statistical Society, Series B, 57(1), 289–300. doi:10.1111/j.2517-6161.1995.tb02031.x

Benjamini, Y., & Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. The Annals of Statistics, 29(4), 1165–1188. doi:10.1214/aos/1013699998

Zhou, R. R., Wang, L., & Zhao, S. D. (2020). Estimation and inference for the indirect effect in high-dimensional linear mediation models. Biometrika, 107(3), 573–589. doi:10.1093/biomet/asaa016

Djordjilovic, V., Page, C. M., Gran, J. M., Nost, T. H., Sandanger, T. M., Veierod, M. B., & Thoresen, M. (2019). Global test for high-dimensional mediation: Testing groups of potential mediators. Statistics in Medicine, 38(18), 3346–3360. doi:10.1002/sim.8199

Dai, J. Y., Stanford, J. L., & LeBlanc, M. (2022). A multiple-testing procedure for high-dimensional mediation hypotheses. Journal of the American Statistical Association, 117(537), 198–213. doi:10.1080/01621459.2020.1765785

Liu, Z., Shen, J., Barfield, R., Schwartz, J., Baccarelli, A. A., & Lin, X. (2022). Large-scale hypothesis testing for causal mediation effects with applications in genome-wide epigenetic studies. Journal of the American Statistical Association, 117(537), 67–81. doi:10.1080/01621459.2021.1914634

Author

Xiufan Yu and Ken Kelley

Examples

set.seed(113)
# Contrasting mediation: active mediators whose effects cancel, so the
# total indirect effect is zero. Compare the two p-values: the benchmark
# Wald test's (pval_hdmm) is 0.19 while the PE test's (pval_pe) is tiny.
d <- simulate_mediation_data(n = 200, p = 60, pattern = "contrasting",
                             outcome = "continuous", c1 = 1)
pe_mediation(d$X, d$Y, d$M, outcome = "continuous")
#> Warning: HBIC selected lambda = 0.05, the smallest value of `lambda_grid`, for the full model and the reduced model of the benchmark test. The criterion is minimized over the grid alone, so its minimum may lie beyond that end or between the end and its neighbor. Consider extending the grid past it and using a finer partition (more values), and compare the active mediators identified and the p-values across grids.
#>  term                  value   
#>  stat_hdmm             1.691   
#>  pval_hdmm             0.1934  
#>  stat_pe               873.9   
#>  j_pe                  872.2   
#>  pval_pe               < 0.0001
#>  total_indirect_effect -0.08409
#>  n_active_mediators    3       
#>  df                    1       
#>  n_candidate_mediators 60      
#>  n_observations        200     
#> 
#> Outcome model: continuous (linear)
#> Active mediators identified (3): 2, 3, 4
#> Tuning parameter (HBIC): lambda = 0.05 from 100 values in [0.05, 10] (the grid's lower end)
# HBIC chose 0.05, the smallest value of the default grid, so the fit
# warned (class poemed_grid_boundary). A grid that extends below it and
# is spaced more finely moves the choice inside the grid, and the screen
# identifies the same three mediators.
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
#> [1] 0.06020101
attr(finer, "active_mediators")
#> [1] 2 3 4