Reproducing the Article's Analyses With POEMED
Source:vignettes/reproducing-the-article.Rmd
reproducing-the-article.RmdThis vignette accompanies the article “Power Enhancement in High-Dimensional Heterogeneous Mediation Analysis” (Yu and Kelley, in press at the Journal of the American Statistical Association). It reruns the article’s empirical data analysis (its Section 5) with POEMED and reads the package’s results against the printed tables, cell by cell. The companion vignette Power Enhancement for High-Dimensional Mediation: A Guided Tour introduces the method itself.
Scope of the comparison. The article benchmarks the power-enhanced tests (PE-HDMM for linear outcomes, PE-HDGMM for generalized outcomes) against several competitors: HDMM, HILMA, and GlobalTest. POEMED implements the proposed PE tests and the HDMM benchmark (the total-indirect-effect Wald test that the PE component augments). Every POEMED result therefore reports the central contrast of the article’s figures and tables, PE versus HDMM; HILMA and GlobalTest are separate methods outside this package.
Load the package before running any of the code below:
Empirical Data Analysis (Article Section 5)
The article asks whether health-care expenditure mediates the
relationship between economic growth and population health across the
world’s countries. POEMED ships that exact data set,
WHO_health_mediation: a country-by-year panel (2000–2021)
with
-
Exposure
X: annual growth in GDP per capita (World Bank), -
Mediators
M: 57 health-expenditure indicators (WHO Global Health Expenditure Database), -
Outcomes
Y: five health outcomes, infant mortality (imr), under-five mortality (u5mr), life expectancy (leb), low birthweight (lbw), and undernourishment (pou),
with WHO region and World Bank income group available as the grouping and confounding variables.
data(WHO_health_mediation)
dim(WHO_health_mediation)
#> [1] 2002 68
table(WHO_health_mediation$region)
#>
#> AFR AMR EMR EUR SEAR WPR
#> 704 616 176 198 110 198Section 5, Table 1: Infant Mortality (IMR)
WHO_mediation_analysis() builds a table in the layout of
the article’s per-outcome tables: it fits the global model (region,
income, and year as confounders), then a separate model within each WHO
region (income and year as confounders) and each income group (region
and year as confounders), and reports the HDMM and PE p-values together
with the active mediators the screen identifies.
The function searches one tuning grid,
seq(0.1, 2, length.out = 100), for the global all-country
model and for every subgroup model alike. It is the grid behind the
article’s real-data tables (the article’s text gives no grid; it states
only that the tuning parameter is chosen by the HBIC criterion), and it
is the default of both lambda_grid, for the subgroup
models, and lambda_grid_global, for the global model, so no
explicit grid arguments are needed. The article’s analysis fit each
model with pe_mediation_linear(), passing this grid as
lambda_grid and again as lambda_grid_reduced,
the grid for the reduced model (the outcome on the mediators and
confounders without the exposure). lambda_grid_reduced
defaults to lambda_grid, so the function’s fits, which
leave it at that default, are the same fits.
imr <- WHO_mediation_analysis("imr", groupings = c("global", "region", "income"))
imr[, c("group", "n_countries", "pval_hdmm", "pval_pe", "active_mediators")]
#> group n_countries pval_hdmm pval_pe active_mediators
#> ALL 91 0.0644 < 0.0001 gge_gdp
#> AFR 32 0.4434 0.4434 none
#> AMR 28 0.4972 0.4972 none
#> EMR 8 0.0261 0.0261 none
#> EUR 9 < 0.0001 < 0.0001 none
#> SEAR 5 0.6483 < 0.0001 ext_usd2021_pc
#> WPR 9 0.8681 < 0.0001 chi_che
#> Low 16 0.3091 < 0.0001 pvtd_usd2021
#> Lower-middle 27 0.9543 0.9543 none
#> Upper-middle 26 0.4481 0.4481 none
#> High 22 0.0389 < 0.0001 oops_che, shi_che
#>
#> Outcome: imrRead this against the article’s Table 1. The highlights:
-
Global (ALL): general government expenditure as a
share of GDP (
gge_gdp) is the active mediator the screen identifies; HDMM is borderline (pval_hdmm= 0.064) while PE is essentially zero. -
SEAR (South-East Asia): HDMM sees nothing
(
pval_hdmm= 0.65) while PE detects external health expenditure per capita in constant 2021 US dollars (ext_usd2021_pc). -
WPR (Western Pacific): HDMM sees nothing
(
pval_hdmm= 0.87) while PE detects compulsory health insurance as a share of current health expenditure (chi_che). -
Low income: HDMM sees nothing
(
pval_hdmm= 0.31) while PE detects domestic private health expenditure (pvtd_usd2021). -
EMR (Eastern Mediterranean), EUR (Europe), and high
income: the benchmark rejects (
pval_hdmm= 0.026, 0.000094, and 0.039; Table 1 prints the EUR value as 0.0001). In the high-income group PE identifies out-of-pocket spending and social health insurance as shares of current health expenditure (oops_che,shi_che); in EMR and EUR no single indicator passes the screen, so the PE p-value equals the benchmark’s. - AFR, AMR, lower-middle income, and upper-middle income: neither test detects mediation.
In the global model, SEAR, WPR, and the low-income group, heterogeneous mediation hides the signal from the conventional test and the PE component recovers it. Among the income groups the strongest signals appear in the low- and high-income groups, as Yu and Kelley (in press) report in Table 1.
The table above agrees with the article’s Table 1 row for row. Every
group’s identified active set is the one the table lists, and every
benchmark p-value matches the printed value to four decimals except
SEAR, which is 0.6483 here where the table prints 0.6482. Where the
table prints a PE p-value below 1e-12, the package prints
< 0.0001, its display floor; the values behind those
entries are below 1e-12 as well.
summary() digests the whole table, highlighting the
groups where the power-enhanced test detects mediation the benchmark
misses (the article’s Table 1 sets those four rows in bold):
summary(imr)
#> POEMED comparison: IMR
#> groups: 11 (11 with data), alpha_level = 0.05
#> detected by benchmark (HDMM): 3 by power-enhanced (PE): 7
#> PE detects mediation in 4 groups the benchmark misses:
#> group n_countries pval_hdmm pval_pe active_mediators
#> ALL 91 0.0644 < 0.0001 gge_gdp
#> SEAR 5 0.6483 < 0.0001 ext_usd2021_pc
#> WPR 9 0.8681 < 0.0001 chi_che
#> Low 16 0.3091 < 0.0001 pvtd_usd2021
#> most-flagged mediators:
#> chi_che 1
#> ext_usd2021_pc 1
#> gge_gdp 1
#> oops_che 1
#> pvtd_usd2021 1
#> shi_che 1The supplement also gives an extended version of
each table (its Table S.8 for IMR) that breaks the power-enhanced result
out by multiplicity method. Pass full_table = TRUE to get
that layout, with the PE p-value and the active mediators identified
reported separately under Bonferroni (FWER), BH, and BY (FDR):
WHO_mediation_analysis("imr", groupings = "global", full_table = TRUE)
#> grouping group n_countries pval_hdmm pval_pe_bonferroni active_bonferroni
#> global ALL 91 0.0644 < 0.0001 gge_gdp
#> pval_pe_bh active_bh pval_pe_by active_by
#> < 0.0001 chi_che, gge_gdp < 0.0001 gge_gdp
#>
#> Outcome: imrThe Other Outcomes, and the Region-by-Income Cells
The remaining article outcomes (undernourishment in the article’s Table 2, and life expectancy, under-five mortality, and low birthweight in the supplement’s Tables S.10 to S.12) run identically; only the outcome code changes:
for (y in c("u5mr", "leb", "lbw", "pou"))
print(WHO_mediation_analysis(y, groupings = c("global", "region", "income")))The article also drills into the region-by-income
cells. WHO_mediation_analysis() returns the
region-by-income block of the table in one call, with NA
for the cells that have no countries or whose fit fails:
imr_cells <- WHO_mediation_analysis("imr", groupings = "region_income")To look inside one cell, build it with
WHO_mediation_design(), which assembles the exposure,
outcome, mediator, and confounder blocks for a chosen subset, and fit it
the way WHO_mediation_analysis() does. A small cell can
contain indicators that never vary inside it. Such a column carries no
signal and cannot be standardized, so the code below drops it by hand
before standardizing, as the function does, and the print footer then
names the active mediators identified by their indicator codes. The code
also passes the function’s default grid, so this fit is the
AFR-Upper-middle row of the block above:
cell <- WHO_mediation_design("imr", region = "AFR", income = "Upper-middle")
c(n_countries = cell$n_countries, mediators = ncol(cell$M))
#> n_countries mediators
#> 4 57
keep <- apply(cell$M, 2L, sd) > 0 # FALSE for a constant indicator
cell$indicators[!keep]
#> [1] "chi_pvt_che"
fit <- pe_mediation(scale(cell$X), cell$Y - mean(cell$Y),
scale(cell$M[, keep, drop = FALSE]),
Z = cell$Z, outcome = "continuous", scale = FALSE,
lambda_grid = seq(0.1, 2, length.out = 100))
fit[fit$term %in% c("pval_hdmm", "pval_pe"), ]
#> term value
#> pval_hdmm 0.0148
#> pval_pe < 0.0001
#>
#> Outcome model: continuous (linear)
#> Active mediators identified (1): pvtd_ncu2021_pc
#> Tuning parameter (HBIC): lambda = 0.138 from 100 values in [0.1, 2]The active mediator identified, domestic private health expenditure
per capita in constant 2021 national currency units
(pvtd_ncu2021_pc), is the one the article’s Table 1 lists
for this cell, and the benchmark p-value, 0.0148, is the value the table
prints.
Six Cells Reported as NA
Across the region-by-income blocks of the article’s tables for its
five outcomes (Table 1, Table 2, and Tables S.8 to S.12), six cells that
the article prints with a value come back from
WHO_mediation_analysis() as NA, each with a
warning that names the cell and quotes the reason. They are of two
kinds.
Four rank-deficient cells. In the EMR-Lower-middle
cells of IMR and U5MR and the AFR-Upper-middle cells of LEB and U5MR,
two of the 57 indicators are exactly collinear: domestic and external
health expenditure as shares of current health expenditure
(dom_che, ext_che) sum to 100 in every
country-year of those cells. The AFR-Upper-middle cells add a second
identity, since compulsory and social health insurance as shares of
current health expenditure (chi_che, shi_che)
coincide there (compulsory private insurance, chi_pvt_che,
is zero throughout, which is why it was dropped as constant above). The
benchmark test inverts the covariance matrix of the mediators the
penalized fit selected, and when that selection keeps both members of a
pair the matrix is singular; the fit stops with the error that the
system is computationally singular, and the function reports the cell as
NA. (The same two identities hold in the IMR cell above,
where the selection keeps ext_che and chi_che
without their partners, so the test goes through.) In three of the four
hand fits below HBIC chooses 0.1, the smallest value of the article’s
grid, for the full model, the reduced model, or both; the fits keep that
grid as the article’s analysis did, so the warning a single
pe_mediation() fit gives at an end of its grid is not shown
in this section. The article’s values come back when the collinearity is
removed before the fit, as the article’s analysis did: it dropped
dom_che and vfa_che (near-exact partners of
ext_che and cfa_che) in the EMR-Lower-middle
cells, and shi_che and ext_che in the
AFR-Upper-middle cells:
cell <- WHO_mediation_design("imr", region = "EMR", income = "Lower-middle")
range(cell$M[, "dom_che"] + cell$M[, "ext_che"]) # 100 in every row
#> [1] 100 100
fit_cell <- function(outcome, region, income, drop) {
cell <- WHO_mediation_design(outcome, region = region, income = income)
keep <- apply(cell$M, 2L, sd) > 0 & !(cell$indicators %in% drop)
fit <- pe_mediation(scale(cell$X), cell$Y - mean(cell$Y),
scale(cell$M[, keep, drop = FALSE]),
Z = cell$Z, outcome = "continuous", scale = FALSE,
lambda_grid = seq(0.1, 2, length.out = 100))
active <- cell$indicators[keep][attr(fit, "active_mediators")]
data.frame(outcome = outcome, cell = paste(region, income, sep = "-"),
dropped = paste(drop, collapse = ", "),
pval_hdmm = fit$value[fit$term == "pval_hdmm"],
active_mediators = if (length(active))
paste(active, collapse = ", ") else "none")
}
rd <- rbind(fit_cell("imr", "EMR", "Lower-middle",
drop = c("dom_che", "vfa_che")),
fit_cell("u5mr", "EMR", "Lower-middle",
drop = c("dom_che", "vfa_che")),
fit_cell("leb", "AFR", "Upper-middle",
drop = c("shi_che", "ext_che")),
fit_cell("u5mr", "AFR", "Upper-middle",
drop = c("shi_che", "ext_che")))
rd
#> outcome cell dropped pval_hdmm active_mediators
#> 1 imr EMR-Lower-middle dom_che, vfa_che 0.08372121 none
#> 2 u5mr EMR-Lower-middle dom_che, vfa_che 0.09144904 none
#> 3 leb AFR-Upper-middle shi_che, ext_che 0.13114741 none
#> 4 u5mr AFR-Upper-middle shi_che, ext_che 0.03894879 pvtd_ncu2021_pcThe four benchmark p-values, 0.0837, 0.0914, 0.1311, 0.0389, are the
values the article prints for these cells (Table 1 for IMR, Table S.10
for LEB, and Table S.11 for U5MR), and in the U5MR AFR-Upper-middle cell
the active mediator identified is the one Table S.11 lists,
pvtd_ncu2021_pc.
Two constant-outcome cells. In the EUR-High and
WPR-High cells of the undernourishment analysis (PoU, the article’s
Table 2), the outcome is 2.5 in every country-year of the cell: the
World Bank publishes a prevalence of undernourishment below 2.5 percent
as 2.5, its lower reporting bound, and the high-income countries of
those two regions are all at that bound. There is no variation in the
outcome to mediate, so pe_mediation() refuses the fit with
the message that Y is constant, and the function reports
the two cells as NA. Table 2 prints a p-value of 1 for
both.
sapply(c("EUR", "WPR"), function(r)
unique(WHO_mediation_design("pou", region = r, income = "High")$Y))
#> EUR WPR
#> 2.5 2.5Why a Mediator? Opening Up the Screening Step
For any fit, pe_mediators() shows the evidence the PE
component acted on: the two path statistics, the screening p-value (the
larger of the two path p-values), and the active flag
(TRUE for a mediator the screen identified as active), for
every candidate the penalized fit selected. This is the per-mediator
detail behind the last column of the article’s tables. In this global
model HBIC picks 0.1, the smallest value of the grid, for the full model
and for the reduced model; across the five outcomes it picks 0.1 for the
full model, the reduced model, or both in 67 of the 128 cells that are
fit. The warning a single fit gives at an end of its grid is not shown
here either.
des <- WHO_mediation_design("imr") # global IMR design
gfit <- 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))
pe_mediators(gfit)
#> mediator name t_outcome t_exposure screen_p active
#> 2 che_pc_usd 8.43 -1.686 0.0918 FALSE
#> 4 gghed 3.027 0.6476 0.5172 FALSE
#> 6 ext 1.141 1.534 0.2540 FALSE
#> 9 pvtd_che 1.058 0.537 0.5913 FALSE
#> 10 oops_che 0.9397 0.9471 0.3474 FALSE
#> 14 gghed_gge -2.405 1.231 0.2183 FALSE
#> 16 pvtd_pc_usd 2.309 -0.8475 0.3967 FALSE
#> 17 oop_pc_usd -8.1 -1.315 0.1886 FALSE
#> 19 cfa_che -2.355 0.6753 0.4995 FALSE
#> 21 chi_che -3.603 3.111 0.0019 FALSE
#> 22 shi_che 2.082 3.703 0.0374 FALSE
#> 25 vhi_che -3.755 -2.494 0.0126 FALSE
#> 26 gge_gdp -5.025 -3.732 0.0002 TRUE
#> 29 gghed_usd -4.748 0.4069 0.6841 FALSE
#> 31 ext_usd 2.95 3.451 0.0032 FALSE
#> 35 ext_ncu_pc -3.875 -1.373 0.1697 FALSE
#> 38 pvtd_ppp_pc -0.5555 -0.2293 0.8187 FALSE
#> 39 ext_ppp_pc -1.532 2.224 0.1256 FALSE
#> 40 pvtd_gdp 2.492 -1.064 0.2875 FALSE
#> 41 ext_gdp -4.303 0.325 0.7452 FALSE
#> 45 ext_ncu2021 -3.454 1.731 0.0834 FALSE
#> 52 pvtd_ncu2021_pc -4.381 0.1229 0.9022 FALSE
#> 53 ext_ncu2021_pc 5.762 -1.532 0.1255 FALSE
#> 57 ext_usd2021_pc 7.531 1.361 0.1735 FALSEReproducibility Notes
- The global, regional, and income analyses above are deterministic.
Called with no explicit grid arguments,
WHO_mediation_analysis()searchesseq(0.1, 2, length.out = 100), the grid behind the article’s real-data tables, for the all-country model and for every subgroup model, and returns the p-values and the active mediators identified that are shown here. - Exact real-data p-values depend on the preprocessing. The article’s
text says that covariates are standardized before fitting.
WHO_mediation_analysis()instead applies the preprocessing of the authors’ analysis code behind the article’s tables: it standardizes the exposure and mediators, centers the outcome, and leaves the confounders (the region and income indicators and the year) unscaled. The region and income factors enter through treatment contrasts with alphabetical reference levels. The guided-tour vignette shows how to do this by hand and what standardizing the confounders as well changes.
References
Yu, X., & Kelley, K. (in press). Power Enhancement in High-Dimensional Heterogeneous Mediation Analysis. Journal of the American Statistical Association.
Session Information
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
#> [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
#> [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
#> [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] POEMED_1.0.0
#>
#> loaded via a namespace (and not attached):
#> [1] cli_3.6.6 knitr_1.52 rlang_1.3.0 xfun_0.61
#> [5] otel_0.2.0 generics_0.1.4 textshaping_1.0.5 jsonlite_2.0.0
#> [9] htmltools_0.5.9 ragg_1.5.2 sass_0.4.10 glmnet_5.1
#> [13] rmarkdown_2.32 grid_4.6.1 evaluate_1.0.5 jquerylib_0.1.4
#> [17] fastmap_1.2.0 foreach_1.5.2 yaml_2.3.12 lifecycle_1.0.5
#> [21] compiler_4.6.1 codetools_0.2-20 fs_2.1.0 Rcpp_1.1.2
#> [25] lattice_0.22-9 systemfonts_1.3.2 digest_0.6.39 R6_2.6.1
#> [29] ncvreg_3.16.0 splines_4.6.1 shape_1.4.6.1 bslib_0.12.0
#> [33] Matrix_1.7-5 tools_4.6.1 iterators_1.0.14 survival_3.8-6
#> [37] pkgdown_2.2.1 cachem_1.1.0 desc_1.4.3