Bio Causal Genomics Mediation Analysis

Version Compatibility

Reference examples tested with: R 4.3+, mediation 4.5.0+, CMAverse 0.1.0+ (GitHub BS1125/CMAverse), HIMA >= 2.3.0 (CRAN), bama 1.3+, causalweight 1.0.5+ (medDML), MVMR 0.4+, TwoSampleMR 0.6+, EValue 4.1+, gesttools 1.3+, ipw 1.0.11+.

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('') then ?function_name to verify parameters
• HIMA must be pinned at >= 2.3.0 for the code patterns below; the formula interface hima(formula, data.pheno, data.M, mediator.type, penalty, ...) was introduced in 2.3.0. Users on HIMA 2.2.x must instead call hima_classic() with positional arguments (X=, Y=, M=, COV.XM=, Y.family=) and lose the formula interface; HIMA 2.2.x is NOT API-compatible with the examples here.
• hima_classic() is still available in 2.3+ for the original Zhang 2016 SIS+MCP pipeline; use it only to reproduce 2016-2021 papers.

If code throws an error, introspect the installed package (?hima, args(cmest)) and adapt the example to match the actual API rather than retrying.

Mediation Analysis

"Does expression of GENE_X mediate the SNP-to-disease effect?" -> Decompose the total effect of a treatment (genotype, exposure) on an outcome into direct and indirect paths through one or more mediators, with explicit handling of exposure-mediator interaction, sensitivity to unmeasured confounding, and high-dimensional mediator screening.

• R (single-mediator, observational, sequential-ignorability assumed): mediation::mediate(med_model, out_model, treat='X', mediator='M', boot=TRUE, sims=5000)
• R (4-way decomposition with exposure-mediator interaction): CMAverse::cmest(...EMint=TRUE, estimation='paramfunc', inference='bootstrap', nboot=1000)
• R (high-dimensional / EWAS mediators): HIMA::hima(Y ~ X + covariates, data.pheno, data.M, mediator.type='gaussian', penalty='DBlasso') (modern v2.3+ formula interface)
• R (MR-based mediation): two-step TwoSampleMR with independent instruments OR MVMR::ivw_mvmr for joint direct effect
• R (doubly-robust double-ML): causalweight::medDML(y, d, m, x)

Sequential ignorability (no unmeasured confounder of treatment-mediator, mediator-outcome, treatment-outcome) is the single load-bearing assumption of observational mediation and is fundamentally untestable. Every report should include a sensitivity result (Imai's rho via medsens() or a mediational E-value).

Algorithmic Taxonomy

| Method | Framework | Handles E-M interaction | High-D mediators | Min n | Fails when | |--------|-----------|--------------------------|------------------|-------|------------| | Baron-Kenny (1986) | Additive regression-based product/difference | No | No | ~100 | Any non-linearity, interaction, or binary outcome; deprecated for causal inference | | Imai mediation R (Imai 2010 Psychol Methods) | Counterfactual ACME/ADE with bootstrap | Yes (via interaction term in outcome model) | No | ~200 | Sequential ignorability violated; exposure-induced M-Y confounder; rare binary outcome with logistic outcome model | | VanderWeele 4-way (VanderWeele 2014 Epidemiology 25:749; 2015 OUP) | CDE + PIE + INTref + INTmed decomposition | Native | No | ~300 | Without interaction term reduces to standard mediation; binary outcome needs rare-disease assumption | | CMAverse (Shi 2021 Epidemiology 32:e20) | 6 estimators: regression (rb), weighting (wb), IORW (iorw), natural effect models (ne), MSM (msm), g-formula (gformula) | Yes | No (single M, or M-vector) | ~300 | Estimator-specific; wb fails with rare exposure; msm needs censoring weights for survival | | HIMA1 / hima_classic (Zhang 2016 Bioinformatics 32:3150) | SIS screen by beta (M->Y) + MCP penalty | No | Yes (up to ~10k) | ~150 + p>>n | Misses mediators with strong alpha and weak beta; screening-step false-negatives | | HIMA2 / hima (Perera 2022 BMC Bioinformatics 23:296) | SIS screen by alpha*beta (indirect effect) + MCP | No (linear by default) | Yes | ~150 | Outcome family limited (gaussian/binomial); HIMA-Cox for survival; HIMA-Pois for count | | HILAMA (Zhang 2025) | High-D mediation with latent confounders | No | Yes (>= 100k) | ~500 | Newer; benchmarks evolving; requires latent-factor specification | | BAMA (Song 2020 Biostatistics) | Bayesian high-D continuous shrinkage | No | Yes (~5k) | ~200 | Slow MCMC; prior sensitivity for very weak mediators | | Two-step MR / network MR (Burgess 2017 Eur J Epidemiol 32:377) | IV-based at each step with INDEPENDENT instruments | Implicit (no interaction modeling) | One mediator at a time | Large summary-stat samples | Same SNP used for E and M (violates exclusion); horizontal pleiotropy; Steiger reversal of M->E direction | | MVMR-mediation (Carter & Sanderson 2021 Eur J Epidemiol 36:465-478) | Total minus direct via MVMR | Implicit | Single mediator | Large GWAS samples for both E and M | Conditional F read CDE, PIE, INTref, INTmed | | Observational, BINARY outcome, rare disease ( 10 each | | Observational with rich confounder set, want doubly-robust estimate | causalweight::medDML (double-debiased ML) | | Longitudinal with time-varying confounding | g-formula via CMAverse::cmest(estimation='gformula') OR gfoRmula package | | Exposure-induced confounder of M-Y exists | Interventional indirect effects (Vansteelandt & Daniel 2017); CMAverse::cmest(estimation='msm') |

Sequential Ignorability and Why It Always Needs Sensitivity

Observational mediation requires three no-unmeasured-confounding assumptions. The third (M-Y unmeasured confounder, after conditioning on E) is the most common violator in genomic mediation because biological confounders (cell composition, batch effects, technical mediators) frequently affect both M and Y.

Sequential ignorability untestable

Trigger: Always, by design.

Mechanism: No statistical test can detect an unmeasured confounder of M-Y. Bootstrap CIs assume the assumption holds; they do NOT propagate uncertainty about it.

Symptom: Significant ACME with no sensitivity reported -> reviewer rejects.

Fix: Report at least one of:

• Imai's rho sensitivity: medsens(med_result, rho.by=0.05, sims=1000); the critical rho where ACME crosses 0; |rhocrit| > 0.3 is "reasonably robust" (Imai 2010), |rhocrit| 0.3 robust; 0.1-0.3 moderately sensitive; 2 robust to plausible biological confounding; 1.5-2 moderate; L -> M -> Y pathway is blocked; if L is not adjusted, L confounds the M-Y leg. Both are wrong under natural-effects; the natural indirect effect is simply not identified.

Fix: Switch to interventional indirect effects (Vansteelandt & Daniel 2017 Epidemiology 28:258), NOT natural indirect effects. Use CMAverse::cmest(estimation='msm') with stabilized inverse-probability weights (yields the randomized-interventional analogue), gfoRmula (parametric g-formula), or randomized indirect effects (Lin 2017 Biometrics 73:1109). The interventional indirect is identified under weaker assumptions than the natural indirect.

HIMA covariate or data.pheno error

Trigger: data.pheno contains factor columns with NA, or formula references columns missing from data.pheno.

Mechanism: HIMA v2.3+ uses a formula interface and constructs the design matrix internally from data.pheno; missing values or unparseable formulas surface as cryptic glmnet errors.

Symptom: Pipeline fails inside hima() with a non-obvious storage.mode or model.matrix error.

Fix: Pre-clean data.pheno (drop NA rows for the variables in the formula; convert factors with factor(); ensure all RHS variables exist as columns). Example:

dat  10 for each exposure independently (Sanderson 2019 IJE 48:713).

### Difference vs product of coefficients diverge for non-linear outcomes

**Trigger:** Binary or survival outcome modeled with logistic / Cox.

**Mechanism:** Difference = total - direct; product = alpha * beta. Equivalent under linear-Gaussian; diverge under any link function. Counterfactual ACME from `mediation::mediate()` is the correct quantity; hand-computed product-of-coefficients on logistic output is biased except under rare-disease.

**Fix:** Report only counterfactual ACME (Imai or CMAverse). For OR-based 4-way decomposition on rare outcomes ( effect-size threshold |
| Bootstrap method + sims | percentile / BCa; min 1000, recommend 5000 |
| Sequential ignorability sensitivity | rho_crit (medsens) OR mediational E-value |
| Exposure-mediator interaction test | Coefficient + p; 4-way decomposition if significant |
| Confounder set justification | DAG description |
| Mediator measurement reliability | Cite |
| Sample size + missing-data handling | Yes |
| Mediator / exposure scale | Standardized? log? raw? |

Reference: AGReMA-Mediation guideline (Lee 2021 BMJ 372:n122) and MacKinnon 2008 Introduction to Statistical Mediation Analysis.

## Reconciliation: Observational vs MR Mediation

| Pattern | Likely cause | Action |
|---------|--------------|--------|
| Observational ACME significant; MR-mediation null | Unmeasured M-Y confounding inflated observational estimate; OR weak IVs in MR | Re-run observational with `medsens()`; if rho_crit  0.2 OR mediational E-value > 1.5, (3) directionally consistent MR-mediation result OR documented absence of valid instruments. Single-method mediation claims should be reported as exploratory.

## Anticipated Reviewer Pushback

| Pushback | Standard response |
|----------|-------------------|
| "Sequential ignorability?" | Imai rho_crit reported via `medsens`; mediational E-value reported on the risk-ratio scale |
| "Exposure-induced confounder of M-Y?" | DAG drawn; if L present, switch to `CMAverse::cmest(estimation='msm')` for interventional indirect effect (Vansteelandt & Daniel 2017) |
| "Why this bootstrap method?" | BCa with sims=5000 for publication; percentile fallback when BCa fails to converge (acceleration estimate unstable at boundary) |
| "Why was MR-mediation not done?" | If valid IVs for E and M exist: two-step MR or MVMR-mediation done (see code below); if not, documented absence of trans-instruments |
| "Mediator measured with error?" | Regression calibration (Carroll 2006 Measurement Error in Nonlinear Models) OR sensitivity analysis assuming reliability r = 0.7 (Valeri & VanderWeele 2014) |
| "Why HIMA2 not BAMA?" | HIMA2 = frequentist + FDR control + faster; BAMA = Bayesian when prior information is available; sample-size justification given against simulation rule-of-thumb |
| "Proportion mediated unstable?" | When |total| = 1000 exploratory, >= 5000 publication | Imai 2010; MC error scales 1/sqrt(sims) |
| Proportion mediated -- meaningful | > 0.2 | Convention; weak guideline only -- effect size in absolute terms matters more (MacKinnon 2008) |
| Proportion mediated -- "most of the effect" | > 0.5-0.8 | Convention |
| Imai rho_crit -- robust | > 0.3 | Imai 2010 Psychol Methods 15:309 |
| Imai rho_crit -- sensitive |  2.0 (working convention; the original Smith & VanderWeele 2019 E-value framework does not prescribe a specific cutoff -- magnitude is context-dependent) | Smith & VanderWeele 2019 Epidemiology 30:835 |
| HIMA FDR cutoff | BH FDR  10 each | Sanderson 2019 IJE 48:713 |
| Two-step MR -- F for both stages | > 10 each | Burgess weak-instrument convention |
| Sample size -- single-mediator (Imai) | >= 200 for stable bootstrap | Simulation rule-of-thumb |
| Sample size -- HIMA EWAS | >= 150 with p_mediators up to ~10k | Zhang 2016 simulations |
| Rare-outcome cutoff for OR-based 4-way | outcome prevalence = 3 timepoints required for g-methods to identify time-varying indirect effects
- Longitudinal mediator measurement at EACH timepoint is required (not just baseline)
- MSM is preferred when treatment is binary and time-varying; g-formula when continuous
- Sequential mediation when the causal ordering of multiple mediators is known and stable across time

### MR-Mediation: Two-Step vs MVMR-Mediation

Decision tree:
- Independent instrument sets available for E and M -> two-step MR (Burgess 2017 Eur J Epidemiol 32:377)
- E and M share instruments (common in cis-eQTL / cis-pQTL mediator cases) -> MVMR-mediation (Carter & Sanderson 2021 Eur J Epidemiol 36:465)
- Both feasible -> report both (triangulation)

Two-step code sketch:

```r
library(TwoSampleMR)
exp_E  Y direction (not Y -> M).

### MR-Mediation: Total Minus Direct via MVMR

**Goal:** Estimate the proportion of a genetic-instrument-identified causal effect that flows through a mediator, using independent IVs for E and (E + M).

**Approach:** Univariable MR for total E->Y; MVMR for direct E->Y conditional on M; indirect = total - direct via delta-method CI.

```r
library(TwoSampleMR); library(MVMR)

total  10 for both E and M independently. If `fstat  2 indicates a confounder would need >2-fold associations with both M and Y to nullify the indirect effect (Smith & VanderWeele 2019).

## Tool Install Notes

| Package | Source | Notes |
|---------|--------|-------|
| mediation | CRAN | `install.packages('mediation')`; actively maintained (Imai group) |
| CMAverse | GitHub | `remotes::install_github('BS1125/CMAverse')`; NOT on CRAN; 6 estimators in one interface |
| HIMA | CRAN | `install.packages('HIMA')`; v2.x renamed `hima()` to HIMA2 -- verify with `?hima` |
| bama | CRAN | `install.packages('bama')`; Bayesian; slow MCMC |
| causalweight | CRAN | `install.packages('causalweight')`; medDML for double-ML mediation |
| EValue | CRAN | `install.packages('EValue')`; for mediational E-values |
| TwoSampleMR | r-universe | See causal-genomics/mendelian-randomization for setup |
| MVMR | r-universe | `remotes::install_github('WSpiller/MVMR')`; for MVMR-mediation |
| gfoRmula | CRAN | For longitudinal / time-varying confounders |

## Common Errors

| Error / symptom | Cause | Solution |
|-----------------|-------|----------|
| `Error in storage.mode(x)  proportion ill-defined | Report ACME and TE separately; pm is unstable when |TE| is small |
| Different ACME between `mediation` and CMAverse `rb` | Default `astar/a` levels differ; binary mediator handled differently | Set `astar=0, a=1` explicitly; for binary mediator pass `mval=list(0)` |
| HIMA returns zero significant mediators | Screening too aggressive; or no true mediators | Try `topN=2*sqrt(n)` instead of default; verify with permutation null |
| Two-step MR shows indirect > total | Steiger reversal: M actually causes E; or pleiotropic SNPs | Run MR-Steiger filter; use MR-PRESSO for pleiotropy |
| `medDML` trim removes most data | Severe positivity violation -- few units with overlapping treatment/mediator distributions | Tighten covariate set; check propensity score distributions |

## References

- Baron RM, Kenny DA 1986 J Pers Soc Psychol 51:1173 (original product-of-coefficients)
- Imai K, Keele L, Tingley D 2010 Psychol Methods 15:309 (counterfactual mediation, sequential ignorability)
- VanderWeele TJ 2014 Epidemiology 25:749 (4-way decomposition)
- VanderWeele TJ 2015 Explanation in Causal Inference (OUP) -- canonical textbook
- Valeri L, VanderWeele TJ 2013 Psychol Methods 18:137 (binary outcomes; rare-disease 4-way)
- Vansteelandt S, Daniel RM 2017 Epidemiology 28:258 (interventional / randomized indirect effects)
- Shi B et al 2021 Epidemiology 32:e20 (CMAverse package; 6 estimators)
- Zhang H et al 2016 Bioinformatics 32:3150 (HIMA original)
- Perera C et al 2022 BMC Bioinformatics 23:296 (HIMA2 alpha-beta screening)
- Song Y et al 2020 Biostatistics 21:806 (BAMA Bayesian high-D mediation)
- Burgess S et al 2017 Eur J Epidemiol 32:377 (two-step / network MR)
- Sanderson E et al 2019 IJE 48:713 (MVMR conditional F-statistic)
- Sanderson E et al 2021 IJE 50:1651 (MVMR-mediation tutorial)
- Farbmacher H et al 2022 Econometrics J 25:277 (medDML / double-ML mediation)
- Smith LH, VanderWeele TJ 2019 Epidemiology 30:835 (mediational E-value)

## Related Skills

- causal-genomics/mendelian-randomization - IV-based causal infe

…

## Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

- **Author:** [GPTomics](https://github.com/GPTomics)
- **Source:** [GPTomics/bioSkills](https://github.com/GPTomics/bioSkills)
- **License:** MIT

Install and usage instructions live in the source repository linked above.