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 (GitHub YinanZheng/HIMA; archived from CRAN 2026-07), 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('<pkg>') 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. On HIMA 2.2.x the classic engine is the top-level hima(X, Y, M, COV.XM=, COV.MY=, Y.family=, penalty=) (positional order X, Y, M) with no formula interface; HIMA 2.2.x is NOT API-compatible with the examples here.
• In 2.3+ the classic engine is hima_classic(X, M, Y, COV.XM=, COV.MY=, Y.type=) (positional order X, M, Y; outcome type via Y.type, not Y.family); use it only to reproduce the original 2016-2021 SIS+penalty pipelines.

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) + de-biased Lasso (DBlasso)	No (linear by default)	Yes	~150	Outcome family limited (gaussian/binomial); HIMA-Cox for survival; HIMA-Pois for count
HILAMA (Wang et al 2025)	High-D mediation with latent confounders	No	Yes (>= 100k)	~500	Newer; benchmarks evolving; requires latent-factor specification
BAMA (Song 2020 Biometrics)	Bayesian high-D continuous shrinkage	No	Yes (~5k)	~200	Slow MCMC; prior sensitivity for very weak mediators
Two-step MR / network MR (Burgess 2015 IJE 44:484)	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 < 10 for either exposure; correlated instruments
medDML (Farbmacher 2022 Econometrics J 25:277)	Double-debiased ML, doubly-robust	Limited (depends on learner)	Moderate (sparsity-friendly)	~500	Severe overlap violations; cross-fitting variance with small n

Methodology evolves; verify against the current CMAverse vignette and the Steen / Vansteelandt natural-effects-model literature before locking analytic choices. Difference-in-coefficients and product-of-coefficients give identical estimates in fully linear-Gaussian models but DIVERGE for any non-linear outcome model (logistic, Cox, Poisson); the counterfactual ACME from mediation::mediate() is the correct quantity for non-linear outcomes.

4-Way Decomposition Framework

VanderWeele 4-way explicitly separates effects from exposure-mediator interaction:

Total Effect = CDE + INTref + INTmed + PIE

Component	Meaning	Active when
CDE	Controlled direct effect (with mediator fixed at reference level)	E directly affects Y
INTref	Interaction-reference -- needs interaction AND exposure	E*M interaction with mediator at reference
INTmed	Mediated interaction -- needs interaction AND exposure AND mediation	E shifts M which then interacts with E
PIE	Pure indirect effect (older "mediation" quantity)	E shifts M which shifts Y additively

Without an exposure-mediator interaction term, INTref = INTmed = 0 and the decomposition collapses to CDE + PIE (= ADE + ACME). With interaction present, traditional ACME mixes PIE and INTmed; the 4-way separation is the only framework that disentangles them. Most epidemiology applications include the interaction term and report all four components (Valeri & VanderWeele 2013 Psychol Methods).

Decision Tree by Scenario

Scenario	Recommended pipeline
Observational, single measured mediator, no plausible E-M interaction, continuous outcome	mediation::mediate() with boot=TRUE, sims=5000; always run medsens()
Observational, single mediator, suspected E-M interaction, any outcome family	CMAverse::cmest(..., EMint=TRUE) -> read CDE, PIE, INTref, INTmed
Observational, BINARY outcome, rare disease (< 10%)	cmest(yreg='logistic', EMint=TRUE, casecontrol=FALSE) -- OR-based 4-way decomposition is valid under rare-disease
Observational, survival outcome	cmest(yreg='coxph') OR HIMA::hima_cox for high-D; report HRs
High-D mediators (EWAS, transcriptome-wide), continuous outcome	HIMA::hima(formula, data.pheno, data.M, mediator.type='gaussian', penalty='DBlasso'); report sigcut (FDR threshold, default 0.05)
High-D mediators with latent confounding (very-high-D EWAS)	HILAMA (2025)
High-D mediators with Bayesian shrinkage (small n, ~5k features)	bama::bama()
Strong genetic IVs for exposure available, single mediator with own IVs	Two-step MR with independent instruments + Steiger filter on mediator
Both E and M have IVs but instruments are weak / correlated	MVMR-mediation with conditional F > 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; |rho_crit| > 0.3 is "reasonably robust" (Imai 2010), |rho_crit| < 0.1 is highly sensitive.
• Mediational E-value (Smith & VanderWeele 2019 Epidemiology 30:835): minimum risk-ratio strength of an unmeasured confounder needed to nullify the observed indirect effect; computed via EValue::evalues.OLS() for linear outcomes or by-hand from ACME risk ratio bounds.
• Reporting BOTH rho-based and E-value sensitivity is standard for high-stakes claims.

Methods-Section Defense of Sequential Ignorability

Template sentence for the methods write-up: "We assumed sequential ignorability conditional on {age, sex, ancestry PCs, cell composition, batch, smoking}. Robustness was assessed via Imai rho_crit at the ACME contrast (medsens, sims = 1000) and the mediational E-value on the risk-ratio scale (Smith & VanderWeele 2019 Epidemiology 30:835)."

Quantitative interpretation thresholds:

• rho_crit: |rho_crit| > 0.3 robust; 0.1-0.3 moderately sensitive; < 0.1 highly sensitive (Imai 2010).
• E-value: E > 2 robust to plausible biological confounding; 1.5-2 moderate; < 1.5 fragile (Smith & VanderWeele 2019).

For high-stakes claims (clinical, drug-target, regulatory submissions) report BOTH rho_crit and the mediational E-value; for exploratory work either alone suffices.

Exposure-induced M-Y confounder

Trigger: A covariate L sits between E and Y, AND is affected by E, AND confounds M-Y.

Mechanism: Standard regression-based mediation cannot adjust for L without blocking part of the indirect effect (collider stratification bias). Adjusting biases CDE; not adjusting biases ACME.

Symptom: Sensitivity to confounder set; ACME flips sign when L is added vs removed.

Operational identification: From the DAG, L is a covariate of M and Y that is also affected by E. VanderWeele TJ, Vansteelandt S & Robins JM 2014 (Epidemiology 25:300) give the criterion: if L is adjusted, part of the indirect E -> 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/interventional indirect effects (Lin SH & VanderWeele TJ 2017 J Causal Inference 5:20150027). 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 <- na.omit(dat[, c('outcome', 'exposure', 'age', 'sex', 'batch', 'pc1', 'pc2')])
dat$batch <- factor(dat$batch)
result <- hima(outcome ~ exposure + age + sex + batch + pc1 + pc2,
               data.pheno=dat, data.M=M_matrix, mediator.type='gaussian')

HIMA mediator-type vs outcome-type mismatch

Trigger: Survival outcome (Surv() on LHS) with mediator.type='compositional' chosen against text mediator panel; or count mediators passed as 'gaussian'.

Mechanism: HIMA v2.3+ auto-detects outcome family from the LHS of formula (continuous, binary, survival, count); mediator.type is set for the mediator data only ('gaussian', 'negbin', 'compositional'). Mismatching mediator-type to the actual mediator distribution biases the screening step.

Symptom: Hazard / rate ratios for indirect effects look implausible; many "significant" mediators fail replication.

Fix: Set mediator.type='gaussian' for continuous (e.g., methylation beta, log-CPM expression), 'negbin' for raw count (RNA-seq), 'compositional' for relative-abundance microbiome. Verify by ?hima in the installed version since the catalogue of mediator types has expanded across releases.

Bootstrap iterations too low

Trigger: sims=100 or sims=500 in early exploration left in for the final report.

Mechanism: ACME CIs from bootstrap have Monte-Carlo error that scales as 1/sqrt(sims); at sims=500 the 95% CI bounds have ~5% MC noise, enough to flip the conclusion at the boundary.

Symptom: Re-running mediate() with a different set.seed() gives substantially different CI bounds.

Fix: sims=1000 minimum for any reported result; sims=5000 for publication; sims=10000 if proximity to zero matters. BCa CIs (boot.ci.type='bca' in mediate()) are slightly more accurate than percentile CIs near zero but require more sims for stability.

Two-step MR instrument independence

Trigger: Same set of SNPs used as instruments for E in step 1 and for M in step 2.

Mechanism: If a SNP affects both E and M, the M-instrument violates exclusion restriction (the SNP-Y association is not exclusively through M). Estimates are biased toward the direct effect.

Symptom: Two-step MR shows large indirect effect; replacing M-instruments with non-overlapping SNPs makes it vanish.

Fix: Apply Steiger filter on the mediator: keep only SNPs where the SNP-M F-statistic exceeds SNP-E F-statistic (or where SNP explains more variance in M than E). For MVMR-mediation, require conditional F > 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 (<= 10%), Valeri & VanderWeele 2013 formulas apply; for common outcomes use risk-ratio scale or marginal effects rather than ORs.

Required Reporting for Publication

Component	Required
ACME estimate + 95% CI (BCa preferred)	Yes
ADE + 95% CI	Yes
Total effect	Yes
Proportion mediated	Yes when total > 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 guideline (Lee H et al 2021 JAMA 326:1045) 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.1, trust MR null
Observational ACME null; MR-mediation significant	Measurement error in M attenuated observational estimate	Trust MR (regression-dilution-free) IF instruments pass Steiger and pleiotropy tests (MR-Egger intercept, MR-PRESSO)
Both significant with same sign	Convergent evidence	High-confidence mediation; report effect size from the more-conservative estimate
Both significant with opposite signs	At least one is biased; revisit confounder structure and IV assumptions	Do not pool; investigate via cross-method sensitivity

Operational rule for high-stakes claims (clinical / drug-target mediation): Require (1) significant observational ACME, (2) Imai 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

Quantitative Thresholds

Quantity	Threshold	Source / Rationale
sims (bootstrap iterations)	>= 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	< 0.1	Same
Mediational E-value -- robust	> 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 < 0.05	Default; report q-values not raw p
MVMR conditional F per exposure	> 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 <= 10%	Valeri & VanderWeele 2013

Working Code Patterns

Single-Mediator Observational with Sensitivity

Goal: Decompose a genotype-disease effect via measured molecular mediator with explicit sensitivity to unmeasured M-Y confounding.

Approach: Fit mediator and outcome models, bootstrap ACME/ADE, run medsens() for Imai rho-based sensitivity.

library(mediation)

med_model <- lm(expression ~ genotype + age + sex + pc1 + pc2 + pc3, data=dat)
out_model <- glm(disease ~ genotype + expression + age + sex + pc1 + pc2 + pc3,
                 data=dat, family=binomial)

med_result <- mediate(med_model, out_model,
                      treat='genotype', mediator='expression',
                      boot=TRUE, sims=5000, boot.ci.type='bca')
summary(med_result)

sens <- medsens(med_result, rho.by=0.05, effect.type='indirect', sims=1000)
summary(sens)

d0, z0, n0, tau.coef slots return ACME, ADE, proportion mediated, total effect.

4-Way Decomposition with Exposure-Mediator Interaction

Goal: Separate CDE, PIE, INTref, INTmed when exposure-mediator interaction is biologically plausible (e.g., gene-environment interaction modifying mediator effect).

Approach: Use CMAverse regression-based estimator with EMint=TRUE; bootstrap CIs.

library(CMAverse)

result_4way <- cmest(
  data=dat, model='rb',
  outcome='disease', exposure='genotype', mediator='expression',
  basec=c('age','sex','pc1','pc2'),
  EMint=TRUE,
  mreg=list('linear'), yreg='logistic',
  astar=0, a=1, mval=list(0),
  estimation='paramfunc', inference='bootstrap', nboot=1000
)
summary(result_4way)

CMAverse reports the 4-way decomposition (Vanderweele 2014): for continuous outcomes the components are cde, intref, intmed, pnie (or pie), te, pm; for non-continuous outcomes (logistic / Cox / Poisson) the ratio versions Rcde, Rpnde, Rtnde, Rpnie, Rtnie are reported. When EMint=TRUE, additional proportion-attributable-to-interaction terms (int, pe) are included. Verify column names with summary(result)$results in the installed CMAverse version, since naming has evolved.

High-Dimensional EWAS Mediation (HIMA2)

Goal: Among thousands of candidate CpG mediators, identify those mediating an exposure-outcome effect with FDR control.

Approach: HIMA v2.3+ uses a formula interface and auto-detects outcome family (Gaussian / binomial / Cox / Poisson). Screening + MCP/DBlasso penalisation + joint significance with BH; the sigcut argument controls the FDR threshold (default 0.05).

library(HIMA)

dat <- na.omit(dat[, c('outcome', 'exposure', 'age', 'sex', 'cell_pc1', 'cell_pc2')])
M_matrix <- as.matrix(beta_values)

result <- hima(
  outcome ~ exposure + age + sex + cell_pc1 + cell_pc2,
  data.pheno=dat,
  data.M=M_matrix,
  mediator.type='gaussian',          # 'negbin' for count, 'compositional' for microbiome
  penalty='DBlasso',                  # default; alternatives 'MCP', 'SCAD', 'lasso'
  scale=TRUE,
  sigcut=0.05,
  parallel=TRUE, ncore=8, verbose=TRUE
)
# result is a data.frame of significant mediators below sigcut