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('<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. 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

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

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

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. Vansteelandt 2009 (Stat Sci 24:471) gives 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 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 <- 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

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

Reconciliation: Observational vs MR Mediation

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

Quantitative Thresholds

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

For survival outcomes wrap the LHS as Surv(time, status); HIMA auto-routes to Cox. The old hima_classic() (Zhang 2016 original) is still exported but screens by beta only and misses mediators with strong alpha + weak beta -- prefer the wrapper hima() unless reproducing a 2016-2021 paper.

For highly-correlated mediators (CpG-island clusters, gene-module co-expression): HIMA uses joint significance with BH-FDR on max(p_alpha, p_beta) and handles correlation only weakly. Within hima(), set penalty='MCP' for stronger correlation handling; alternatively pre-reduce the mediator panel by principal components or by clustering correlated mediators and screening the cluster centroid (VanderWeele & Vansteelandt 2014 Epidemiol Methods 2:95).

Time-Varying Mediation Methods

When the mediator is measured at multiple timepoints (or exposure varies over time), natural-effects estimands are not identified; switch to one of:

• g-formula (parametric or Monte Carlo): gfoRmula::gformula_continuous_eof(); CMAverse::cmest(estimation='gformula')
• g-estimation of a structural nested mean model: gesttools::gestSingle() / gestMultiple()
• Marginal structural model with stabilized IPTW: ipw::ipwtm() followed by glm(..., weights=sw)
• Sequential mediation for K timepoints: VanderWeele & Tchetgen Tchetgen 2017 JRSSB 79:917

Choose based on the experimental structure:

• = 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:

library(TwoSampleMR)
exp_E <- extract_instruments('ieu-a-2', clump=TRUE)
m_E <- extract_outcome_data(exp_E$SNP, 'ieu-b-30')
dat_EM <- harmonise_data(exp_E, m_E)
mr_EM <- mr(dat_EM)

exp_M <- extract_instruments('ieu-b-30', clump=TRUE)
exp_M_indep <- exp_M[!exp_M$SNP %in% exp_E$SNP, ]
exp_M_indep <- steiger_filtering(exp_M_indep)
out_M <- extract_outcome_data(exp_M_indep$SNP, 'ieu-a-7')
dat_MY <- harmonise_data(exp_M_indep, out_M)
mr_MY <- mr(dat_MY)

Indirect effect = beta_EM * beta_MY (product of coefficients). CI via delta method or parametric bootstrap of the joint (beta_EM, beta_MY) distribution. Steiger filter on M-instruments is mandatory to ensure the M -> 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.

library(TwoSampleMR); library(MVMR)

total <- mr_ivw(harmonised_data_E_Y)
mvmr_dat <- format_mvmr(BXGs=cbind(beta_E, beta_M),
                        BYG=beta_Y, seBXGs=cbind(se_E, se_M), seBYG=se_Y, RSID=snps)
fstat <- strength_mvmr(mvmr_dat, gencov=0)
mvmr_fit <- ivw_mvmr(mvmr_dat)
direct <- mvmr_fit$coef[1, 'Estimate']
direct_se <- mvmr_fit$coef[1, 'Std. Error']

indirect <- total$b - direct
indirect_se <- sqrt(total$se^2 + direct_se^2)
indirect_ci <- indirect + c(-1.96, 1.96) * indirect_se

Require fstat conditional F > 10 for both E and M independently. If fstat < 10, use Q-statistic-adjusted IVW (qhet_mvmr) or report the result as weak-instrument-limited.

Double-ML Doubly-Robust Mediation

Goal: Avoid model misspecification of both mediator and outcome models via cross-fitted ML nuisance estimators.

Approach: causalweight::medDML uses random forests (or other learners) with sample splitting to estimate nuisance parameters; final estimator is doubly robust.

library(causalweight)

result_dml <- medDML(
  y=dat$outcome, d=dat$treatment, m=dat$mediator,
  x=as.matrix(dat[, covariates]),
  trim=0.05, order=1
)

Reports direct, indirect (via mediator), and total effects with influence-function-based standard errors. Robust to non-linearity and interactions; assumes sequential ignorability still.

Mediational E-Value for Sensitivity

Goal: Report the minimum strength of an unmeasured M-Y confounder required to nullify the observed indirect effect.

Approach: Convert ACME and its CI to a risk-ratio scale, then apply VanderWeele E-value formula.

library(EValue)

acme_rr <- exp(med_result$d0)
acme_lower_rr <- exp(med_result$d0.ci[1])
evalues.RR(acme_rr, lo=acme_lower_rr, hi=NULL)

For binary outcomes, convert ACME on probability scale to RR; for continuous, use evalues.OLS() with the standardized indirect effect. E-value > 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

Common Errors

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 inference; foundation for MR-mediation
• causal-genomics/pleiotropy-detection - MR-mediation instrument validity hinges on pleiotropy diagnostics (MR-Egger intercept, MR-PRESSO)
• causal-genomics/colocalization-analysis - Confirm shared causal variant before causal mediation
• causal-genomics/fine-mapping - Identify the causal variant driving the exposure
• methylation-analysis/differential-cpg-testing - Per-CpG inputs for HIMA EWAS mediation
• differential-expression/deseq2-basics - Expression inputs for eQTL mediation
• multi-omics-integration/mofa-integration - Multi-layer mediator construction
• population-genetics/association-testing - GWAS summary statistics for MR-mediation
• clinical-biostatistics/effect-measures - Risk-ratio / odds-ratio scales for binary outcomes
• machine-learning/model-validation - Cross-fitting and sample splitting for medDML