---
title: "Reproducing the Article's Analyses With POEMED"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Reproducing the Article's Analyses With POEMED}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7,
                      fig.height = 4.2)
```

This 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:

```{r load}
library(POEMED)
```

# 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.

```{r who-data}
data(WHO_health_mediation)
dim(WHO_health_mediation)
table(WHO_health_mediation$region)
```

## Section 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.

```{r who-imr}
imr <- WHO_mediation_analysis("imr", groupings = c("global", "region", "income"))
imr[, c("group", "n_countries", "pval_hdmm", "pval_pe", "active_mediators")]
```

```{r who-imr-values, include = FALSE}
p_hdmm <- function(g) formatC(signif(imr$pval_hdmm[imr$group == g], 2),
                              digits = 2, format = "fg", flag = "#")
p4 <- function(g) sprintf("%.4f", imr$pval_hdmm[imr$group == g])
```

Read 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` = `r p_hdmm("ALL")`) while PE is essentially zero.
* **SEAR (South-East Asia):** HDMM sees nothing (`pval_hdmm` =
  `r p_hdmm("SEAR")`) while PE detects external health expenditure per
  capita in constant 2021 US dollars (`ext_usd2021_pc`).
* **WPR (Western Pacific):** HDMM sees nothing (`pval_hdmm` =
  `r p_hdmm("WPR")`) while PE detects compulsory health insurance as a
  share of current health expenditure (`chi_che`).
* **Low income:** HDMM sees nothing (`pval_hdmm` = `r p_hdmm("Low")`)
  while PE detects domestic private health expenditure (`pvtd_usd2021`).
* **EMR (Eastern Mediterranean), EUR (Europe), and high income:** the
  benchmark rejects (`pval_hdmm` = `r p_hdmm("EMR")`, `r p_hdmm("EUR")`,
  and `r p_hdmm("High")`; 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
`r p4("SEAR")` 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):

```{r who-summary}
summary(imr)
```

The 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):

```{r who-full}
WHO_mediation_analysis("imr", groupings = "global", full_table = TRUE)
```

## The 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:

```{r who-others, eval = FALSE}
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:

```{r who-region-income, eval = FALSE}
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:

```{r who-cell}
cell <- WHO_mediation_design("imr", region = "AFR", income = "Upper-middle")
c(n_countries = cell$n_countries, mediators = ncol(cell$M))

keep <- apply(cell$M, 2L, sd) > 0          # FALSE for a constant indicator
cell$indicators[!keep]
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"), ]
```

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,
`r sprintf("%.4f", fit$value[fit$term == "pval_hdmm"])`, 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:

```{r who-rank-deficient, warning = FALSE}
cell <- WHO_mediation_design("imr", region = "EMR", income = "Lower-middle")
range(cell$M[, "dom_che"] + cell$M[, "ext_che"])   # 100 in every row

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
```

The four benchmark p-values,
`r paste(sprintf("%.4f", rd$pval_hdmm), collapse = ", ")`, 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.

```{r who-constant-y}
sapply(c("EUR", "WPR"), function(r)
  unique(WHO_mediation_design("pou", region = r, income = "High")$Y))
```

## Why 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.

```{r who-screen, warning = FALSE}
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)
```

# Reproducibility Notes

* The global, regional, and income analyses above are deterministic.
  Called with no explicit grid arguments, `WHO_mediation_analysis()`
  searches `seq(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

```{r session-info}
sessionInfo()
```