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.

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.4974    0.4974   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: imr

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 = 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            1

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

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

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:

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_pc

The 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.5

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.

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   FALSE