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

## ----load---------------------------------------------------------------------
library(POEMED)

## ----who-data-----------------------------------------------------------------
data(WHO_health_mediation)
dim(WHO_health_mediation)
table(WHO_health_mediation$region)

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

## ----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])

## ----who-summary--------------------------------------------------------------
summary(imr)

## ----who-full-----------------------------------------------------------------
WHO_mediation_analysis("imr", groupings = "global", full_table = TRUE)

## ----who-others, eval = FALSE-------------------------------------------------
# for (y in c("u5mr", "leb", "lbw", "pou"))
#   print(WHO_mediation_analysis(y, groupings = c("global", "region", "income")))

## ----who-region-income, eval = FALSE------------------------------------------
# imr_cells <- WHO_mediation_analysis("imr", groupings = "region_income")

## ----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"), ]

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

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

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

## ----session-info-------------------------------------------------------------
sessionInfo()