Decoding Causal Structure: End-to-End Mediation Pathways Inference
Authors: Yulong Li
Xiwei Liu
Feilong Tang
Ming Hu
Jionglong Su
Zongyuan Ge
Imran Razzak
Eran Segal
Venue: NeurIPS 2025 poster ()
Keywords: Mediation Pathways Inference
Primary area: machine_learning_for_sciences

Abstract: Causal mediation analysis is crucial for deconstructing complex mechanisms of action. However, in current mediation analysis, complex structures derived from causal discovery lack direct interpretation of mediation pathways, while traditional mediation analysis and effect estimation are limited by the reliance on pre-specified pathways, leading to a disconnection between structure discovery and causal mechanism understanding. Therefore, a unified framework integrating structure discovery, pathway identification, and effect estimation systematically quantifies mediation pathways under structural uncertainty, enabling automated identification and inference of mediation pathways. To this end, we propose Structure-Informed Guided Mediation Analysis (SIGMA), which guides automated mediation pathway identification through probabilistic causal structure discovery and uncertainty quantification, enabling end-to-end propagation of structural uncertainty from structure learning to effect estimation. Specifically, SIGMA employs differentiable Flow-Structural Equation Models to learn structural posteriors, generating diverse Directed Acyclic Graphs (DAGs) to quantify structural uncertainty. Based on these DAGs, we introduce the Path Stability Score to evaluate the marginal probability of pathways, identifying high-confidence mediation paths. For identified mediation pathways, we integrate Efficient Influence Functions with Bayesian model averaging to fuse within-structure estimation uncertainty and between-structure effect variation, propagating uncertainty to the final effect estimates. In synthetic data experiments, SIGMA achieves state-of-the-art performance in pathway identification accuracy and effect quantification precision under structures uncertainty, concurrent multiple pathways, and nonlinear scenarios. In real-world applications using Human Phenotype Project data, SIGMA identifies mediation effects of sleep quality on cardiovascular health through inflammatory and metabolic pathways, uncovering previously unspecified multiple mediation paths.

Decoding Causal Structure: End-to-End Mediation Pathways Inference
Yulong Li Xiwei Liub*, Feilong TanglAr, Ming Hu
Jionglong Su Zongyuan Ge Imran Razzak
Eran Segal 1.2.t
Mohamed bin Zayed University of Artificial Intelligence Weizmann Institute of Science
Xi' an
University
Monash University
EEuul contribution
iCorresponding authors
Abstract
Causal mediation analysis is crucial for deconstructing complex mechanisms of action. However, in current mediation analysis, complex structures derived from causal discovery lack direct interpretation of mediation pathways, while traditional mediation analysis and effect estimation are limited by the reliance on pre-specified pathways, leading mechanism understanding. Therefore, unified framework integrating structure discovery, pathway identification, and effect estimation systematically quantifies mediation pathways under structural uncertainty, enabling automated identification and inference of mediation pathways. To this end, we propose Structure Informedd Guided Mediation Analysis (SIGMA), which guides automated mediation path- way identification through probabilistic causal structure discovery and uncertainty quantification, enabling end-to-end propagation of structural uncertainty from structure learning to effect estimation. Specifically, SIGMA employs differentiable Flow-Structural Equation Models to learn structural posteriors, generating diverse Directed Acyelic Graphs (DAGss to quantify structural uncertainty. Based on these DAGs, introduce the Path Stability Score to evaluate the marginal probability pathways, identifying high-confidence mediation paths. For identified media- tion pathways, we integrate Efficient Influence Functions with Bayesian model averaging to fuse struitrstructure estimation uncertainty and between-structure effect variation, propagating uncertainty to the final effect estimates. In synthetic data experiments, SIGMA achieves state-of-the-ar performance in pathway identi- fication accuracy and effect quantification precision under structural uncertainty, concurrent multiple pathways, and nonlinear scenarios. In real-world applications using Human Phenotype Project data, SIGMA identifies mediation effects of sleep quality on cardiovascular health through inflammatory and metabolic pathways, uncovering previously unspecified multiple mediation paths.
disconnection between structure discovery and causal
Introduction
Causal Mediation Analysis (CMA) has emerged as critical framework bridging causal inference and machine learning. Recent studies [47 (7] emphasize the evolution of CMA from identifying causal effects to elucidating underlying mechanisms, transition essential for developing machine learning systems capable of intervention reasoning, cross-domain knowledge transfer, and interpretable expla- nations [10] 2 (22. In particular, analyzing causal mediation pathways involving multiple mediators arranged in parallel, sequential, or networked structures [58] provides deeper mechanistic
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

insights in bigh-dimensional contexts commonly encountered in biology, social networks, and behav- ioral sciences [2T 26 40 For instance, elucidating how sleep quality impacts cardiovascular health in Human Phenotype Projects (HPP) 36 requires modeling its interplay with multiple mediators (e.g., inflammation, metabolic) across interrelated pathways 30 61 77. Therefore, quantifying the
effects of mediation pathways under causal structural uncertainty remains a challenge.
Although analyzing complex mediation pathways is both theoretically and practically significant, current methods face several limitations: (i) Structure Discovery (SD): While mainstream causal discovery algorithms [56] 55 72 can automatically learn causal structures, they lack principled mechanisms for quantifying and propagating structural uncertainty, undermining the reliability of downstream analyses. (ii) Pathway Identification (PI): Current structure discovery methods [ 29] fail to explicitly identify specific mediation pathways. Meanwhile, classical mediation analysis approaches 3 27 rely on pre-specified paths, which become infeasible in high-dimensional data due to the exponential growth in the number of candidate paths. (ii)) Effect Estimation (EE): Semi- parametric methods like Efficient Influence Function (EIF) 6 offer theoretical guarantees (e.g.. multiple robustness) [59] 71 41 but typically assume that the mediation structure is known and correctly specified, which creates a contradiction with the structural uncertainty faced in the previous two stages. Consequently, a gap between current methods becomes evident in high-dimennsioal settings like HPP, where discovery yields causal structure without pathway interpretation. while traditional mediation is confined to predefined paths, missing underlying causal mechanisms. In this work, we argue that a unified framework integrating SD, PI, and EE can quantify the mediating pathways under structural uncertainty in real-world settings.
To this end, we propose an end-to- and-to-end Structure-Iforrmed Guided Mediation Analysis (SIGMA) framework, which integrates probabilistic causal structure discovery with uncertainty quantification to guide automatic mediation pathway identification, and propagates structural uncertainty to ro- bust effect estimation, thus enabling integrated inference of causal structures and mediation effects. Specifically, differentiable Flow-Structural Equation Models (44] are leveraged to learn edge param- eters encoding structural uncertainty, from which an ensemble of directed acyclic graphs (DAGs) is subsequently sampled. Building on this DAG ensemble, we introduce the Path Stability Score which performs Bayesian inference by evaluating the marginal probability paths across DAGs to select high- cogh-coofidence mediation pathways, thereby circumventing manual path pre-specification. For each identified mediation pathway, SIGMA employs Efficient Influence Functions estimate pathway-specific mediation effects, using nuisance functions tailored to each DAG structure. Sub- sequently, these pathway-specins estimates aggregated across the DAG ensemble via Bayesian model averaging (BMA) [48), propagating structural uncertainty into the final causal effect estimates. Furthermore, SIGMA integrates variational autoencoder-base imputation module explicitly- model latent dependencies among variables, enabling accurate imputation of missing values. We empirically demonstrate SIGMA's effectiveness on nigh-dimensional. heterogeneous HPP dat:" identifying and quantifying mediation pathways connecting sleep quality to cardiovascular health through inflammatory and metabolic mediators. The contributions are summarized as follows:
We propose an end-to-end framework, SIGMA, unifying probabilistic structure discovery with uncertainty quantification, high confidence path identification, and effect estimation incorporating structural uncertainty, to address mediation analysis under uncertainty structures.
We introduce the Path Stability Score to identify mediation by quantifying their stability across posterior directed acyclic graph samples, thus eliminating manual pre- pre-specification.
We develop a cructure-guided estimation approach, employing Efficient Influence Functions (EIF), combined with Bayesian model averaging, enabling propagation of structural uncertainty and inference by leveraging EIF properties.
Extensive experiments validate SIGMA on real-world HPP data, revealing complex and previously unspecified mediation pathways linking sleep quality to cardiovascular health through inflammatory and metabolic, demonstrating its interpretability in real-world.
2 Dataset
Human Phenotype Project (HPP)'. We conduct our real-world evaluation using the HPP dataset, deeply phenotyped cohort comprising over 6,000 individuals with multi-night home sleep apnea
https::/knooxleeggbbaee. pheno.

---

Sleep Characteristics
Body Mass Index
Restorative sleep ratio
height
Sieep duration
weight
Skeep continulty
Appendicular lean
Sleep efficiency
Visceral fat
Snoring ratio
A/G retio
Snoring ratio
OSA
Cardiovascular System
HRV (during night)
Oxyyen desaturation index
Vescular Health
Oxyyen burcen
Carotid ultresound IMT
Carotid ultresound plaque
Blood pressures
Metabolic Pathways
Pulse rate
GLtMB butyrate ato producers
GLdMB proplonabe producar
Inflemmation
sleop charasteristies
veriables
informmation veriebles
GuSMI vitamin producers
GudMB probiotic farmenters
GutMB-tma producers
GuthiB-ips producers
metabolic variables
cardlovasullrr ables
GLZMB akkormansia
GutMi--uufaa reducers
Exieting pathways
Discovered pathways
GUdMB shannon
Guthia-roteeolyiic bacteria
14 days after visit
days
Measure Collected:
3-nights sleep monitoring
Stool sample collection
Blood tests from Maintenance Organication (HMD)
Biood preszure
sample
Elechrcardiogram (ECG)
Body compositio
Continuous meesures over 14-day post-visit pertod: continuous glucase monitoring Food logging
Carolid
Arthropometrics
(a) Illustration of five body system-level categories from the HPP used in this work. (b) Timeline of data acquisition in the HPP cohort. (c) Conceptual schematic of the SIGMA framework, showing its capacity to identify both known and previously unrecognized mediation pathways under structural uncertainty. (d) Distributions and dependencies among variables in a representative synthetic dataset (MidDim-S), illustrating the indirect effects of mediators M₁, M2, M3 on outcome
Figure 1:
Y under different treatment A and the direct effect of A on Y. Refer Appendix E.3 for details.
testing (HSAT) data across 16,000+ nights. Beyond sleep, HPP includes multimodal data covering 16 other body systems that represent a wide array of physiological functions and environmental exposures 34 This coordinated collection process, as shown in Figure provides temporally aligned, longitudinal measurements across multiple physiological systems. The integration of sleep monitoring, clinical biomarkers, imaging, and omics data within #6-month window requires considerable logistical effort and sustained participant engagement. These characteristics make the HPP dataset particularly valuable testbed for evaluating SIGMA's ability recover complex mediation pathways under structural uncertainty. In total, we utilize data from 11 body systems, covering over 500 raw clinical and physiological attributes. To support tractable and clinically meaningful mediation analysis, we aggregated these into 62 related variables across systems. Each feature was designed to reflect coherent physiological process suitable for pathway-level causal interpretation. Among these, five representative categories analyzed in our work: Sleep Characteristics, Body Mass Index (BMI), Cardiovascular System, Inflammation and Metabolic Pathways, are illustrated in Figure[ Detailed data descriptions can be found in the Appendix D As depicted in Figure Te]. SIGMA framework is capable of validating previously reported mediation pathways [9] while also uncovering novel, clinically meaningful causal mediation mechanisms spanning different systems. We provide results to demonstrate this in Section 5.3 where we quantify these pathways and assess their natural effect.
features by grouping
Synthetic Data. To evaluate SIGMA's ability to identify mediation pathways under controlled yet realistic conditions, we generate synthetic datasets that mirror key statistical and structural charac- teristics of real-world data (Figure d)). Specifically, we simulate missing data patterns with both missing completely at random (MCAR) [50] and missing at random (MAR) mechanisms to reflect non-random measurement gaps common in clinical studies, and introduce heterogeneous variable types (continuous, binary, categorical) and nonlinear relationships to mimic biological complexity. Each synthetic dataset is constructed from randomly sampled DAG using an Erdös-Rénn (ER) (0) model, followed by structural equation modeling (SEM) [46] with configurable nonlinearity and noise. We explicitly embed mediation structures parallel, chain, or bybrid-bbtwer a designated treatment and outcome, with tunable effeett strengths to ensure identifiable direct and indirect effects.
3

---

Functional mechanisms include both linear and nonlinear mappings (e.g., polynomial, sinusoidal), and exogenous noise is drawn from mixed distributions induce non- non-Gaussianity. Natural Direct/Indirect are estimated via Monte Carlo simulations to facilitate quantitative benchmarking. This synthetic setup enables systematic evaluation of SIGMA's robustness across varied structural regimes, pathway configurations, and data quality levels. See Appendix [E]] for detailed generation procedures and E.2 for embedded mediation pathway structures.
Effects (NDE/NIE) 45 49
3 Causal Mediation: Definitions, Notation, and Identification
We consider the problem of causal mediation analysis in observational data. Let A denote the treatment variable, the outcome variable, M = (M....., } a set of k potential mediators, and C' set of covariates. The variables may be continuous, discrete, mixed-type. Let V {A, Y U M UC be the full set of observed variables and assume the underlying causal structure is DAG G (V, E), where Pag(V] denotes the parent set of node V,. We assume that the true graph G satisfies the Causal Markov and Faithfulness assumptions, which enable conditional independencies to be inferred from the graph. Appendix provides detailed discussion of these assumptions, including Causal Sufficiency.
We adopt the Potential Outcomes Framework 50 57 to define causal mediation effects. Let (a) be the potential mediator value under treatment a, and Y(a, m) the potential outcome if A and M
We focus two core mediation effects comparing treatment level to reference level a": (1) The NDE [45] is defined as NDE(a, a²) = EE (a. Al (a")) - Y (a", . ([)]))] quantifying the effect of changing A from a" to on Y while holding M fixed at its natural level under A = a". (2) The NIE [49] is defined as NIE a²) = Y M (a)) - Y (a, M(a))), capturing the effect on Y mediated through M when A change from a" to When the true causal structure is unknown, our goal is to estimate these effects under structural uncertainty. Under certain assumptions [45. the Total Effect. TE(a )=) = (a) - Y(a)) can be decomposed as TE NDE + NIE.
We use Il to denote conditional independence. Appendix A.3 discuss more about following assumptions, in- cluding strategies for practical consideration and the theoretical implications of violation. Throughout, identification denotes recovery of NDE and NIE under ignorability, consistency, and positivity; path identification denotes that a path attains high posterior support as quantified by the PSS. Estimation and plug-in estimators under these assumptions.
Identifying NDE and NIE from observational data requires several key assumptions 45
then proceeds via ceffuuence-fucction--e
Consistency: if A a', then M = (a'. and if A a' and M = m, then Y Y m), for all relevant treatment levels a' (e.g., a' E {aa a" }))) and all m in the support of M.
Ignorability: The following conditional independence relations hold almost surely for all relevant values a. a", m. m': (i) A 11 Y(a", m) 11 M(a) C.
C: (ii) M I1 Y(a, m') | A = a, C; and (iii)
Positivity: For all a and in the support of A and C, the treatment propensity score T(a | C)= = Pr(A = a C = c) (2), for some constants 0 < C1 < e2 < 1. The conditional density function for discrete M = m given A and C = c is bounded between [pI, for some constants 0 P1 P2 < oo, almost surely for all relevant a, c, and m.
These assumptions collectively allow identification of NDE and NIE from observational data, forming the foundation of our mediation analysis under unknown causal graph structure.
4
Methodology
4.1
SIGMA Framework Overview
We propose Structure-Informer Guided Mediation Analysis (SIGMA), a unified framework for causal mediation analysis under structural uncertainty. SIGMA operates in an end-to-end fashion, integrating causal discovery, pathway identification, and effect estimation within probabilistic graphical modeling paradigm. Given observational data over treatment, outcome, mediators, and confounders, SIGMA performs probabilistic structure discovery via a based structural equation model (Flow- SEM), generating posterior ensemble of DAGS that capture causal uncertainty (Section 4.2). then identifies stable mediation pathways using a path stability score (PSS), which quantifies the marginal support for each path across sampled DAGS (Section 4.3) For each retained pathway, SIGMA estimates natural direct and indirect effects using structure-guided nuisance estimation and

---

hybrid influence-funciion/plaee function/plugin strategy, with final estimates aggregated via Bayesian model averaging to propagate structural uncertainty (Section 4.4). This unified design enables automated, uncertainty-aware mediation analysis in complex high-dimensional systems.
4.2 Probabilistic Structure Discovery
To model structural uncertainty in causal mediation, SIGMA employs Flow-SEM to learn distribution over DAGs from observational data. Let V (V...... VeF denote the set of variables, the goal is to infer posterior distribution P(G' D) over DAGs G EV that encode the conditional dependencies among these variables. The validity of this structure discovery step relies several core assumptions, including the Causal Markov, Faithfulness, and Causal Sufficiency assumptions. Refer to Appendix A for a detailed discussion.
4.2.1 Flow-SEM Model Formulation
In Flow-SEM, each variable V; V is modeled as a stochastic function of its parents in candidate DAG G, with the structural equation:
V2 = (Paa) Φ.) + a. (P (V.); U,, U, ~ NV00 1),
(1)
where Pag (V:) denotes the parent set of V, and fi a; are flexible neural network functions parame- terized by Φi. The collection { defines the global functional parameters. This formulation accommodates nonlincarities and heteroskedasticity in the data-generee generating process.
To enable differentiable structure learning, SIGMA introduces continuous weight matrix Wo RPXP, where each entry [Wo]j: reflects the strength or existence
directed edge from V, to V,. zero-diagonal [Wo] i) = O is imposed prevent self-loops. The parent Pag(V, is inferred from the sparsity pattern of We. Under this, the joint distribution P(V) Wo, 9, b) can be derived via a change-of-variables formulation using normalizing flows. Given an observed data sample VO), its log-likelihood computed as:
[[V))/W,93 log P(V Wa T) = Σ Pog (1())) log (Pag (V:1));
where (1) (v.6) (Pac o,)) (Pa 0.) represents the standardized noise for
variable v, This formulation enables likelihood-bart training of both structural and functional components via stochastic optimization.
4.2.2 Model Learning and Optimization
Model parameters (We, 4) are optimized by minimizing regularized loss that balances data fit, structural sparsity, and acyclicity:
L(Wo, :; D)
n log P(V6)/(1 We, +) + Adag (Wo) + a W
(3)
where h ((W)) = tr(exp(W C We)) - serves as continuous relaxation for enforcing acyclic- ity. exp(.) is the matrix exponential and denotes element-wise multiplication: ||W Lisj [olsi encourages sparsity. Further technical details and acyclicity properties are discussed in Appendix A Gradient- optimization is used to train the model. Once trained, We captures both the functional and structural dependencies among variables.
4.2.3 Approximate Structural Posterior and DAG Sampling
To approximate the structural posterior P(G) D), SIGMA interprets each entry in Wo as defining a marginal edge inclusion probability via:
Pij = o Welssd), with pa; = 0,
(4)
where (-) is the sigmoid function. Candidate DAGs {G,) N generated by sampling each edge independently as Bernoulli(p.)). followed by cycle removal to ensure graph validity.
This ensemble of sampled DAGs provides a Monte Carlo approximation of the posterior over causal structures, capturing both legel and pathway-level uncertainty. These samples are used
downstream for pathway selection and mediation effect estimation (see Sections 4.3 and 4.4).
5

---

4.3 Automated Mediation Pathway Identification
Given a posterior ensemble of sampled DAGs {G, = (V, - SIGMA identifies mediation pathways that are robust to structural uncertainty by quantifying the posterior support for each candidate path. This is achieved through the Path Stability Score (PSS), which estimates the marginal probability of a directed path being present in the true causal structure.
Formally, let я = (Vo- -> 1/2 - ... - V&)) denote a directed acyclic path from treatment node Vi = A to outcome node Vi = Y. The PSS of path is defined as:
NaaG
PSS(
Σ 1(TT C.),,
(5)
Nvalid
where I(TT S G.) = 1 if all directed edges in и are present in graph G., and otherwise. The term Nvalid denotes the number of acyclic graphs retained after cycle removal. As shown in Appendix B.1
PSS is an unbiased and consistent estimator of the posterior inclusion probability P( S G D).
To construct the candidate set of mediation paths, SIGMA enumerates all simple directed paths from A to Y in each valid DAG G., subject to a maximum length constraint for computational tractability. These paths are aggregated across DAGs, and their frequencies normalized to obtain the empirical PSS scores. Paths with stability above a user-specified threshold E [0. 1] are retained as high- confidence mediation candidates, forming the set: Ilsblle = {ir PSS(T) 2 7). This selection process effectively propagates structural uncertainty from edge-level variation to the pathway level, allowing downstream effect estimation to focus statistically supported mediation mechanisms. A theoretical justification for using PSS as a Monte Carlo estimator and its asymptotic properties is provided in Appendix B.1 and B.2
4.4 Mediation Effect Estimation
Once a set of stable mediation pathways aabbee is identified from the DAG ensemble, SIGMA proceeds to estimate path-specific NDE and NIE for each retained pathway E Iuabll. Given that each pathway is supported by a subset of DAGs {G, क S G.). effect estimation must account for both within-structure statistical uncertainty and netween-structure causal ambiguity. SIGMA achieves this by combining structure-specit averaging (BMA). Appendix A.33 details the assumptions required for valid estimation, including ignorability, consistency, and positivity conditions.
nuisance estimation with hybrid estimators and Bayesian model
4.4.1 Structure-Guided Nuisance Estimation
For each stable pathway π E Ilstables SIGMA estimates the nuisance functions fi. by leveraging structure-specifie information from each DAG G, in the ensemble. The conditioning sets for all components of f), are defined strictly by the parent sets Pag. (-) within the corresponding DAG. Neural networks are used to fit the conditional relationships implied by these parent sets. To reduce estimation bias and enhance robustness, we use K-fold cross- cross-fitting fittinn to compute out-of-sample predictions across the dataset. This procedure yields all necessary components for downstream effect estimation via both the EIF and plugin estimators, including the required counterfactual outcome predictions.
4.4.2 DAG Effect Estimation
After obtaining nuisance function estimates v) for specific DAG G, (as detailed in Section 4.4. 1). the SIGMA framework calculates path-specific NDE and NIE for a path π. These constitute the
single-DAG effect vector ò, (T) = (m), (T), SIGMA employs a hybrid estimation strategy: for paths of length L(T) (number of nodes) equal (i.e., A + M + Y), an influence function based estimation method is utilized; for paths with L(x)) > 3, a plugin estimator is adopted. All effect quantifications are based on a comparison between a pre-specified exposure level a and a reference level a". For a binary exposure A, these levels are 1 and 0, respectively. For a continuous exposure A, a" is its estimated mean f and a PA k 5 (with k = in this study), where ea and ã. are estimated from the first K-Fold training set for the specific (path, DAG) combination.
Influence function-based Estimation (L(TT) = 3): The point estimate (,(TT) is the sample mean of the corresponding influence function scores vh (Ze; ne. Go), i.c., (,())) 1. E'. (Z n. G.).

---

For binary exposure A, the influence functions VNNDE Z3 17)) and I'NIE (Z: (.) are constructed within the SIGMA framework by combining cross-fitted estimates of several nuisance functions. These include the propensity score P(A, = a a(CA,,), outputs from the mediator model MMIAA.Car which characterizes the conditional behavior of M, given A, and covariates CM,., and the conditional expectation of the outcome E|Y; efficient, bias-corrected estimates for the target counterfactual expectations (c.g.. E[ [a(((',, (a)))) that comprise NDE and NIE. The specific mathematical expressions for these influence functions as utilized within SIGMA are detailed in Appendix A.3.4
a, Md, Cy These IFs are formulated to provide
For continuous exposure A, SIGMA constructs and applies specifie influence function. This func- tion aims to correct plugin estimates of counterfactual expectations (e.g., Yaa- = EE M (aa))) which form the NDE and NIE, by integrating cross-fitted nuisance components using the data-driven intervention levels (a, a") These key nuisance components include: (i) the conditional expectation of A, (C A (ii) estimate of the global homoscedastic residual variance of A, à ACCs.resi (iii)
the conditional expectation of the outcome of the target counterfactual expectation " and their sample mean Fao-. The influence function score Wran. (Zi) used in SIGMA to estimate a counterfactual expectation (such Taa+) ) defined as:
(i())
Ay|A=level,MMa Cy. and (iv) individual plugin predictions
A. -PAA((CA.) (PViA=a.Mas..o - -aaar) (fan* -aa=))
Ve (Z)
(6)
3 AlCaites
The overall influence functions for NDE and NIE derived from linear combinations of these component scores. The theoretical motivation and property analysis of this SIGMA-specific influence function are detailed in Appendix A.3.5 For all influence function-based estimates, the single-DAG variance (T)) is estimated by / ET
x²)2, where vi (T)).
Plugin Estimation (L(m)= > 3): For paths longer than three nodes, SIGMA utilizes a plugin estimator. NDE and NIE are computed by substituting the relevant cross cross-fitted counterfactual outcomes (i.e., indi- vidualized ny (a, M (a))) predictions from Section 4.4.1) into their defining expressions, followed by averaging over the sample. Specifically, NIE, (T) E (iL)) (a, M()) - p (a, M(a")) Its
variance (0,()) is calculated as level counterfactual difference.
( Eh - Act)²) where A denotes the individual-
4.4.3 Aggregation and Uncertainty Propagation
To generate final mediation effect estimates that account for structural uncertainty, SIGMA aggregates single-DAG estimates (T) and their variances (f)( for path æ across the Nputh supporting DAGs. The Bayesian Model Averaged point estimate is θBMA (T) Npan INya 0, (T)) using uniform
weights for each supporting DAG G5- The total BMA variance, ((BMAA(T) (T)), combines the average within-DAG variance (reflecting estimation uncertainty for each structure) and the between-DAG variance (capturing structural uncertainty across different structures) via the law of total variance:
N pash
Nysh
V ((BEA(TT)
1 Npath
((0(T))
1
(Da(T) (BMA(TT)
Npath
Here, (TT( is the Aingle-DAG variance estimate from Section 4.4.2 Confidence intervals
are then constructed using the standard normal approximation: 0BMA (T) + 21-a/2 v (θBMa
This aggregation framework ensures that both estimation uncertainty (within-DAG) and structural uncertainty (between-DAG) are rigorously propagated into the final mediation effect estimates.
5 Experiments
We evaluate SIGMA across two key tasks: (1) causal structure discovery and (2) causal mediation analysis, covering both synthetic datasets with known ground-truth effects and real-world medical cohort from the HPP. For structure discovery, assess the accuracy of learned causal graphs in
7

---

various settings with different graph topologies, nonlinearity levels, and data types. For mediation analysis, we SIGMA's ability to estimate direct and indirect effects of identified mediation pathways. The synthetic data allows controlled henchmarking, while the HPP cohort provides complex, high-dimensional setting with clinically relevant mediation structures. full list of results and details of the experimental setup are in Appendix F.1
Experimental setup Our experimental evaluation comprises two primary tasks: (1) performance assessment of causal structure discovery, and (2) end-to-end causal mediation analysis evaluation. For synthetic data experiments, we construct 10 benchmark datasets that systematically vary feature dimensionality, degree of nonlinearity, and pathway structural complexity assess SIGMA's robust- ness across different complexity regimes. Each dataset contains 6,000 samples with mixed variable types (60% continuous, 30% binary, 10% categorical) and 10% injected missingness to simulate realistic scenarios.
The SIGMA framework is implemented through three-stage pipeline: (1) learning the structural posterior via Flow-SEM and sampling 1,000 DAGs; (2) identifying high-confidence mediation pathways based on Path Stability Score (threshold 7 0.15); (3) estimating natural direct and indirect effects for identified pathways, with estimates aggregated across multiple DAGs via Bayesian model averaging to propagate structural uncertainty. Effect estimation employs a 5-fold cross-fitting strategy, utilizing EIF for paths of length and plug-in estimators for longer paths.
5.1 Causal Discovery Evaluation
Datasets. We generate 10 synthetic datasets with varying characteristics, including different graph sizes from low to high dimensional p E {20, 50, 100, 200), nonlinearity level (Promiin {0, 0.5)), graph structures (parallel e {2, 3, 6} and chain é {1, 2}). and variable types variables (continuous, binary, categorical). Sptesct-specific configurations are provided in Appendix E.3
Figure shows that SIGMA demonstrates consistently strong performance across datasets. On datasets (HighDim-D/P), DECI has a slight advantage in adjacency F1 (0.3814/0. 14/80484669 vs. 0.3741/0.4754), potentially due differences between SIGMA's Flow-SEM continuous optimization framework in capturing weak connections under high-dimensional sparse conditions and DECI's global variational inference approach. In contrast, SIGMA shows comprehensive leadership in orien- tation determination (average improvement of 42. 1%). a capability crucial for accurate identification of mediation pathways. Experimental results confirm that SIGMA under diverse conditions, providing robust foundation for downstream causal effect estimation.
directed structures
O
LowDisp
LowOmID
LowOmn
LonDin
LavD PP
roCCans McDIino
020
(a)
HgDin-
H gDie
(b)
DgDkm
HyDkm
(c)
HgBino
HggnP
PC
Notears
ICALINGAM
Grandag
DECI
SIGMA
Figure 2: Causal discovery on benchmark datasets. We compare SIGMA against PC [33] (lincar)
Notcars (72). ICALiNGAM (55), Grandag (20] and DECI [22] ten synthetic Datasets. (a) shows
adjacency Fl 233. (b) reports orientation Fl (60 and (c) evaluates causal accuracy [5] Refer
Table in Appendix F for detailed results.
5.2 End-fo-end Causal Mediation Analysis
We evaluate SIGMA's effect estimation performance on synthetic data by comparing estimates against ground truth values generated via structural equation models (see Appendix F2) For comparison,

---

use the Lasso, random forest (RF), gradient boosted machine (GBM) and DeepMed [65] Two high- confidence mediation pathways were selected from the synthetic dataset LowDim-L for evaluation. As shown in Table T SIGMA consistently produces low bias and narrow confidence intervals across both natural direct and indirect effects. These results indicate that by sampling multiple DAGs from the structural posterior and applying BMA, SIGMA effectively incorporates both within-graph estimation uncertainty and across-graph structural variation. This leads to robust and stable inference of mediation effects under structural uncertainty.
Table 1: The empirical biases, empirical standard errors (SE) and the lower and upper bounds of the corresponding 95% confidence intervals (CI) of the estimated NDE, NIE and TE across two mediation pathways on synthetic dataset LowDim-L. (Pathway I: V5-V V16; Pathway 2: V12-V10-VV6. Bias is the Monte-Carlo mean deviation from the ground truth.
NDB
NIE
TE
Metbed
DeepMeed
Lasc
GBM
SIGMA
DonpMied
Laso
RF
GBM
SIGMA
DeepMed
Liw
GBM
SIGMA
Emptrial Blas
000667
-0.33 00.3956
-0.2303
0.0713
-0.0090
0.977
01030
0.4347
-6,300
0.0811
0.3416
0.3164
-01300
-0.131
0,0766
Potrway
SE
Upper
D.0077
0:055
000139
03:7552
013345
110833
-00677
0/005
0.2063
0,0113
0,0755
0.880
0.311
00111
..1001
0,0095
0.0152
0.8064
0.187
0.0078
0.006
0.407
010m
0.006
0.066
0.3964
0.0145
0.1873
Imppricall Brac
BO185
-00114
.00653
80158
0
0.397
0.64
0.1522
0.0422
40631
0,158
0.630
6,889
-00513
-00228
Patiway
SE
Uppe
8.000
00877
0.9494
0.0081
0
00061
80871
0.9229
0,0877
03387
0.050
0.000
0.320
0.001
0.448
00128
6.3822
60159
00864
0.1128
0.200
0.003
6.2199
68123
0.1629
0.6995
0.3491
06112
Table 2 reports the performance metrics of SIGMA across different dimensionalities. As feature dimenstonality increases from 500 to 2000, all three core evaluation metrics exhibit declining trend: Adjacency F1 decreases by 22.5% (from 0.603 to 0.467), Orientation F1 decreases by 25.9% (from 0.541 to 0.401). and Path F1 decreases by 26.9% (from 0.587 to 0.429). Despite this performance degradation, SIGMA maintains an edge recovery Fl of 0.467 and pathway identification FI of 0.429 in the 4000-dimensio configuration, with all standard errors controlled within 0.083. These results confirm the scalability of the SIGMA framework for high-dimensional mediation analysis.
Table 2: SIGMA high-dimensional scalability analysis. Setup: three synthetic configurations with p E (000, 1000, 2000]; nonlinearity degree = 0.5: parallel (n 2) and chain (L 2) structures; 6000; 10 independent runs. Metrics are reported as mean (SE).
Dimension (p) Adj Fl (Mean)
Adj Fl (SE)
Orient Fl (Mean)
Orient Fl (SF)
Path Fl (Mean)
Path (SF)
500
0,603
0,044
0.541
0,051
0.587
0.047
1000
0.548
0.058
0.479
0.066
0.521
0.062
2000
0.467
0.074
0.401
0.083
0.429
0.079
5.3 Real Data HPP Validation
Validation the HPP dataset shows that the SIGMA framework successfully identified multiple physiologically meaningful mediation pathways (Figure3), Specifically, the pathway (sleep HRV index + pulse wave velocity - carotid IMT) verified known mechanism whereby sleep quality affects cardiovascular health via vascular function. Previous studies have reported associations between reduced sleep quality and autonomic dysfunction, indicated by decreased heart rate variability (HRV), which is linked to increased arterial stiffness (PWV) and subsequently greater carotid IMT 32 51I. Additionally, SIGMA also identified pathway where sleep snoring severity affects carotid IMT via body fat percentage and blood pressure regulation. This suggest an indirect effect of sleep quality on vascular health mediated through metabolism and blood pressure. Existing clinical studies have linked severe snoring and obstructive sleep apnea (OSA) to obesity via endocrine disturbances 241. Obesity, particularly increased visceral f't, has been linked to elevated blood pressure 773. and sustained high blood pressure is associated with increased carotid IMT [37]. Moreover, as illustrated in Figure the analysis reveals associations among nocturnal hypoxic burden, OSA, body fat composition, carotid IMT, and gut microbiota abundance. This newly identified mediation pathway suggests connections among sleep-disordens cardiovascular health, and gut microbiota ecology. This indicates that abnormal body fat distribution induced by sleep disturbances may further affect cardiovascular health through elevated blood pressure and autonomic dysfunction. Classical medical regression mediation analysis further validated these pathways, with an overall TE error of 0.2734 and consistent NIE and NDE estimation directions. Numerical results and visualization details are provided in the Appendix F.4
breathing, metabolism,
9

---

120
V16
-0.2
0.
Sleep hear rate
41
Body fat composition
182
204
ow
Sleep quality
Carotid IMT
322
V9
85A
OSA
Vascular health
Blood lying
128
10
15A
Blood pressure orthcstatic
Sleep quality
V26
42A
Vascular health
V78
F5A
Sleep hear rate
Carofid IMT
Body fat composition
551
v10
Blood pressure resting
433
Nocturnal hypoxia burden
Carotid IMT
(c)
Blood pressure standing
199
Gut Microbiome ebundence
100
V70
Sleep OSA
Body composition
Figure Correlations of key sleep-derived features with features derived from other body systems. See Table 3 for the corresponding features description of V6-V61. (b) The mediation pathways reproduce a previously reported mechanism; (c) Novely discovered pathways by SIGMA, showing potential links across sleep, metabolic state, vascular function, and gut microbiome composition.
Related Work
CMA has evolved from estimating total treatment effects to decomposing pathways into direct and indirect components [27 45) 49). Semiparametric methods leveraging EIFs have been widely adopted for this purpose, offering multiple robustness and asymptotic efficiency (59 (71). Recent work such as DeepMed (65] advances this line by using deep neural networks estimate high-dimensional nuisance functions without sparsity constraints, achieving semiparametric efficiency under weaker smoothness assumptions. However, these methods generally assume known mediation structure. Parallelly, causal discovery frameworks like NOTEARS 722, GES [56). and DAG-GNN [69] aim to recover causal graphs from observational data but do not target mediation
pathways. DECI (22) proposes variational approach that jointly infers causal structure and structural equations, supporting downstream inference via interventional queries. While powerful, DECI and similar methods [1] 60 focus recovering global structure and interventional distributions rather than quantifying -sweifffcscciicccccc mediation effects. Other works sigh-dimensionll mediation [8] have explored multi-mediator identification but typically rely predefined structures linear assumptions. This disconnect between causal discovery and effect estimation limits their applicability in complex, uncertain systems. Please refer Appendix for detailed literature review.
7 Conclusion, Limitations and Future
The SIGMA framework we proposed provides an end-to-end solution for causal mediation analysis under structural uncertainty by unifying structure discovery, pathway identification, and effect estima- tion. The core innovation of this framework lies in: quantifying causal graph uncertainty through probabilistic structural posterior learning, automatically identifying high-confidence mediation path- ways via the Path Stability Score, and systematically propagating structural uncertainty to final effect estimates through Bayesian model averaging. In validation experiments HPP data, SIGMA confirmed known aleep-cardivvaseula mediation relationships and discovered previously unspecified cross-system mediation pathways, revealing how sleep quality affects cardiovascular health through inflammatory and metabolic pathways, providing new targets for clinical intervention. However, as observational research method, SIGMA remains limited by unmeasured confounders and selection bias. To address these challenges, we plan to conduct clinical intervention experiments, validating the inflammatory and metabolic mediation mechanisms discovered by SIGMA through sleep quality interventions (e.g., CPAP treatment, sleep hygiene improvement), providing stronger causal evidence, and developing sensitivity analysis tools systematically assess the robustness of inference results, thereby facilitating translation from computational discoveries to clinical applications.
10

---

References
[1] Chen Avin, llya Shpitser, and Judea Pearl. Identifiaillity path-specinn effects. In Proceedings of the 19th International Joint Conference Artificial Intelligence (IJCAI), pages 357-363, 2005.
[2] Elias Bareinboim and Judea Pearl. Causal inference and the data-fusion problem. Proceedings of the National Academy Sciences, 113(27):7345-3552 2016.
[3] Reuben M. Baron and David A. Kenny. The modderator-mediatee variable distinction in so- cial psychological research: Conceptual, strategic, and statistical considerations. Journal of Personality and Social Psychology. (1(6):1173-11
1986.
[4] Francesco Beghini, Lauren Mclver, Aitor Blanco- Maguue, Leonard Dubois, Francesco Asnicar, Sagun Maharjan, Ana Mailyan, Paolo Manghi, Matthias Scholz, Andrew Maltez Thomas, et al. Integrating taxonomic, functional, and strain-level profiling of diverse microbial communities with biobakery 3. elife, 10:e65088, 2021.
[5] Ioana Bica, Ahmed Alaa, and Mihaela Van Der Schaar. Time series deconfounder: Estimating treatment effects over time in the presence of hidden confounders. In International conference on machine learning, pages 884 895. PMLR, 2020.
[6] Peter J. Bickel, Chris A. J. Klaassen, Ya acov Ritov, and Jon A. Wellner. Efficient and Adaptive Estimation for Semiparametric Models. Johns Hopkins Series in the Mathematical Sciences. Springer-Verlag, New York, 1993.
[7] Michaël G B Blum, Linda Valeri, Olivier François, Solène Cadiou, Valérie Siroux, Johanna Lepeule, and Rémy Slama. Challenges raised by mediation analysis in high-dimension setting. Environmental Health Perspectives, 128(5),55001, 2020.
[8] Böhnke. Explanation in causal inference: Methods for mediation and interaction. Quarterly Journal of Experimental Psychology, 6946):1243-124, 2016. Epub 2016 Feb 16.
[9] Robert M. Bond, Christopher J. Fariss, Jason J. Jones, Adam D. I. Kramer, Cameron Marlow, Jaime E. Settle, and James H. Fowler. A 61-million personnnnneeonnnnooooo experiment in social influence and political mobilization. Nature, 489(7415):295-299) 2012.
[10] Léon Bottou, Jonas Peters, Joaquin Quiñoneer-CCnnela Denis X. Charles, D. Max Chickering, Elon Portugaly, Dipankar Ray, Patrice Simard, and Ed Snelson. Counterfactual reasoning and learning systems: The example of computational advertising. Journal of Machine Learning Research, 14(1):207-3260, 2013.
[11] Nicolas Bourdillon, Fabien Jeanneret, Masih Nilchian, Manon Albertoni, Peter Ha, and Gré- goire P Millet. Sleep deprivation deteriorates heart rate variability and photoplethysmography Frontiers in Neuroscience, 15:642548, 2021.
[12] Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexandre Drouin. Differentiable causal discovery from interventional data. Advances in Neural Information Processing Systems, 33:21865-22 865-21877 2020.
[13] Lu Cheng, Ruocheng Guo, and Huan Liu. Causal mediation analysis with hidden confounders. In Proceedings of the Fifteenth ACM International Conference on Web Search and Data Mining, pages 113-122. 2022.
[14] Silvia Chiappa. Path-specific counterfactual faimess. In Proceedings of the AAAI conference on artificial intelligence, volume 33, pages 7801-7808, 2019.
[15] Tom Claassen and Tom Heskes. A bayesian approach to constraint-based causal inference. In Proceedings of the 28th Conference on Uncertainty in Artificial Intelligence (UAI), pages 207-216, 2012.
[16] Chris Cundy, Aditya Grover, and Stefano Ermon. Bed nets: Scalable variational approaches for bayesian causal discovery. Advances in Neural Information Processing Systems, 34:7095-7110, 10, 2021.
11

---

[17] Alicia Curth and Mihaela Van der Schaar. On inductive biases for heterogeneous treatment effect estimation. Advances in Neural Information Processing Systems, 15815483-15894, 2021.
[18] Rhian M. Daniel, Bianca L. De Stavola, Simon N. Cousens, and Stijn Vansteelandt. Causal mediation analysis with multiple mediators. Biometricx, 71(1):1-14, 2015.
[19] Alex Dregan, Lauren Rayner, Katrina AS Davis, Ioannis Bakolis, Jorge Arias De La Torre, Jayati Das-Munshi, Stephani L Hatch, Robert Stewart, and Matthew Hotopf. Associations between depression, arterial stiffness, and metabolic syndrome among adults in the uk biobank population study: mediation analysis. JAMA psychiatry, 77(6):598-600 2020.
[20] Erdós and A. Rényi. On random graphs i. Publicationes Muthematicae Debrecen, 6(3-4):290- 297, 1959.
[21] Victoria Garfield. Sleep duration: A review of genome-wide association studies (gwas) in adults from 2007 to 2020. Sleep Medicine Reviews, 56: 101413, 2021.
[22] Tomas Geffner, Javier Antoran, Adam Foster, Wenbo Gong, Chao Ma, Emre Kiciman, Amit Sharma, Angus Lamb, Martin Kukla, Nick Pawlowski, et al. Deep end-to-end causal inference. arXiv preprini urXiv:2202.02195. 2022.
[23] Clark Glymour, Kun Zhang, and Peter Spirtes. Review of causal discovery methods based on graphical models. Frontiers in Genetics, 10:524, 2019.
[24] Igor Alexander Harsch, Peter Christian Konturek, Corinna Koebnick, Peter P Kuehnlein, Fabian S Fuchs. Sayeh Pour Schahin, Gerald H Wiest, Eckhart G Hahn, Tobias Lohmann, and Joachim H Ficker. Leptin and ghrelin levels in patients with obstructive sleep apnoea: effect of cpap treatment. European Respiratory Journal, 22(2):251-257, 2003.
[25] Jason Hartford, Greg Lewis, Kevin Leyton-Brown, and Matt Taddy. Deep iv: A flexible approach for counterfactual prediction. In International Conference on Machine Learning, pages 1414-1423. PMLR, 2017.
[26] Gökhan S. Hotamisligil. Inflammation and metabolic disorders. Nature, 444(7111)8800-86) 2006.
[27] Kosuke Imai, Luke Keele, and Dustin Tingley. A general approach to causal mediation analysis. Psychological Methods, 35(4)330--334, 2010.
[28] Kosuke Imai, Luke Keele, Dustin Tingley, and Teppei Yamamoto. Unpacking the black box: Learning about causal mechanisms from experimental and observational studies. American Political Science Review, 105(4):765-789, 2011.
[29] Kosuke Imai and Teppei Yamamoto. Identification and sensitivity analysis for multiple causal mechanisms: revisiting evidence from framing experiments. Political Analysis, 1112)1141-171, 2013.
[30] Licheng Jiao, Yubin Wang, Xinyu Liu, Ling Li, Feng Liu, Wei Ma, Yao Guo, Peng Chen, Shuai Yang, and Bo Hou. Causal inference meets deep learning: A comprehensive survey, Research, 7(6):0467, 2024.
[31] Fredrik Johansson, Uri Shalit, and David Sontag. Learning representations for counterfactual inference. In International conference on machine learning, pages 3020-3029, PMLR, 2016.
[32] Masahiko Kadoya and Hiroshi Koyama. Sleep, autonomic nervous function and atherosclerosis. International Journal of Molecular Sciences, 20(4):794, 2019.
[33] Markus Kalisch and Peter Bühlman. Estimating high-dimensional directed acyclic graphs with the pc-algorithm. Journal of Machine Learning Research, 8(3), 2007.
[34] Sarah Kohn, Alon Diament, Anastasia Godneva, Raja Dhir, Adina Weinberger, Yotam Reisner, Hagai Rossman, and Eran Segal. Phenome-wide associations of sleep characteristies in the human phenotype project. Nature Medicine, pages 1-12, 2025.
12

---

[35] Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste- Gradient- based neural dag learning, 2020.
[36] Zachary Levine, Iris Kalka, Dmitry Kolobkov, Hagai Rossman, Anastasia Godneva, Smadar Shilo, Ayya Keshet, Daphna Weissglas- Tal Shor, Alon Diament, Yeela Talmor-Barkan, Yaron Aviv, Tom Sharon, Adina Weinberger, and Eran Segal. Genome wide association studies and polygenic risk score phenome-wide association studies across complex phenotypes in the human phenotype project. Med, 5(1):90-101.cc) January 2024.
[37] Sun Min Lim, Hyeon Chang Kim, Hoon Sang Lee, Joo Young Lee, Mina Suh, and Song Vogue Ahn. Association between blood pressure and carotid intima-media thickness. The Journal of Pediatrics, 154(5):667-611 667)))1.116666,,,, 2009.
[38] Phillip Lippe, Taco Cohen, and Efstratios Gavves. Efficient neural causal discovery without acyclicity constraints. arXiv preprint 1XXi::207.11483, 2021.
[39] Christos Louizos, Uri Shalit, Joris M Mooij, David Sontag, Richard Zemel, and Max Welling. Causal effect inference with deep latent-variable models. Advances in neural information processing systems, 30, 2017.
[40] Yunan Luo, Jian Peng, and Jianzhu Ma. When causal inference meets deep learning. Nature Machine Intelligence, 2(8):426-427, 2020.
[41] Marloes H. Maathuis, Markus Kalisch, and Peter B"uhlmann. Estimating high-dimensional intervention effects from observational data. The Annals of Statisties, 37(6A):3133 3164, 2009.
[42] Dominique Makowski, Tam Pham, Zen Lau, Jan C Brammer, François Lespinasse, Hung Pham, Christopher Schölzel, and SH Annabel Chen. Neurokit2: A python toolbox for neuro- physiological signal processing. Behavior research methods, pages 1-8, 2021.
[43] Razieh Nabi and Ilya Shpitser. Fair inference on outcomes. In Proceedings of the AAAI Conference Artificial Intelligence, volume 32, 2018
[44] Achille Nazaret, Justin Hong, Elham Azizi, and David Blei. Stable differentiable causal discovery. arXiv preprint arXiv:2311.10263, 2023.
[45] Judea Pearl. Direct and indirect effects. In Proceedings of the Seventeenth Conference on Uncertainty in Artificial Intelligence, pages 411-420 Morgan Kaufmann, 2001.
[46] Judea Pearl. Causal inference in statistics: An overview. 2009.
[47] Judea Pearl. The causal mediation formula guide to the assessment of pathways and mechanisms. Prevention Science, 13(4):426 436, 2012.
[48] Adrian E. Raftery, David Madigan, and Jennifer A. Hoeting. Bayesian model averaging for linear regression models. Journal of the American Statistical Association, 92(437):17--19,, 1997.
[49] James M. Robins and Sander Greenland. Identifiability and exchangeability for direct and indirect effects. Epidemiology, 3(2):143-155, 1992.
[50] Donald B. Rubin. Estimating causal effects of treatments in randomized and nonrandomized studies. Journal of Educational Psychology, 66(5):688-70) 1974.
[51] Alicia Saz-Lara, Maribel José Francisco López-Gil, and Iván and sleep quality with arterial stiffness: A systematic review and meta-analysis. Sleep Health, 8(6):663-670, 2022.
Linion-LLcas--oooeee Arthur E Mesas, Blanca Notario-Pachest Association between sleep duration
[52] Bernhard Schölkopf, Francesco Locatello, Stefan Bauer, Nan Rosemary Ke, Nikolay Kalchbren- ner, Ambarish Goyal, and Yoshua Bengio. Toward causal representation learning. Proceedings of the IEEE, 10165)6112-634 2021.
13

---

[53] Uri Shalit. Fredrik D Johansson, and David Sontag. Estimating individual treatment effect: generalization bounds and algorithms. In International conference on machine learning, pages 3076-3085. PMLR, 2017.
[54] Claudia Shi, David Blei, and Victor Veitch. Adapting neural networks for the estimation of treatment effects. Advances in neural information processing systems, 32, 2019.
[55] Shohei Shimizu, Patrik O. Hoyer, Aapo arvnnri and Antti Kerminen. A linear non- gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7:2003- 2030, 2006.
[56] Peter Spirtes, Clark N. Glymour, and Richard Scheines. Causation, Prediction, and Search (2nd Ed.). MIT Press, 2000.
[57] Jerzy Splawa-Neyman, Dorota Dabrowska, and T. Speed. On the application of probability theory agricultural experiments. essay on principles. section 9. Statist. Sci., 11 1990.
[58] Daniel L. Taylor, Alisa U. Jackson, Naruemon Narisu, and et al. Integrative analysis of gene expression, dna methylation, physiological traits, and genetic variation in human skeletal muscle. Proceedings of the National Academy of Sciences, 116(22):0088311000)) 10883-10888, 2019.
[59] Eric J. Tchetgen Tchetgen and Ilya Shpitser. Semiparametric theory for causal mediation analysis: Efficiency bounds, multiple robustness, and sensitivity analysis. The Annals of Statistics, 4063):1816-1845, 2012.
[60] Ruibo Tu, Cheng Zhang, Paul Ackermann, Karthika Mohan, Hedvig Kjellström, and Kun Zhang. Causal discovery in the presence of missing data. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 1762-1770. Pmlr, 2019.
[61] Tyler WanderWeele and Stijn Vansteelandt. Mediation analysis with multiple mediators. Epi- demiologic Methods, 2(1):95-115, 2014.
[62] Tyler J. VanderWeele and Peng Ding. Sensitivity analysis in observational research: Introducing the e-value. Annals of Internal Medicine, 167(4):268-277 Aug 2017. Epub 2017 Jul 11.
[63] Victor Veitch, Dhanya Sridhar, and David Blei. Adapting text embeddings for causal inference. In Conference on Uncertainty in Artificial Intelligence, pages 919-928. PMLR, 2020.
[64] Manuel C. Voelkle, Johan H. Oud, Eldad Davidov, and Peter Schmidt. An sem approach to continuous time modeling of panel data: relating authoritarianism and anomia. Psychological Methods, 17(2):176-19) Jun 2012. Erratum in: Psychological Methods. 2012 173:37(3):34.
[65] Siqi Xu, Lin Liu, and Zhonghua Liu. Deepmed: Semiparametric causal mediation analysis with debiased deep learning. Advances in Neural Information Processing Systems, 35:28238-28255 2022.
[66] Ziqi Xu, Debo Cheng, Jiuyong Li, Jixue Liu, Lin Liu, and Ke Wang. Disentangled representation for causal mediation analysis. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 37, pages 10666-10674, 2023.
[67] Liuyi Yao, Sheng Li, Yaliang Li. Mengdi Huai, Jing Gao, and Aidong Zhang. Representation learning for treatment effect estimation from observational data. Advances in neural information processing systems, 31, 2018
[68] Jinsung Yoon, James Jordon, and Mihaela Van Der Schaar. Ganite: Estimation of individualized treatment effects using generative adversarial nets. In International conference on learning representations, 2018
[69] Yue Yu, Jie Chen, Tian Gao, and Mo Yu. Dag-gnn: gnn: Dag structure learning with graph neural networks. In International conference on machine learning, pages 7154-7163. PMLR, 2019.
[70] Weijia Zhang, Lin Liu, and Jiuyong Li. Treatment effect estimation with disentangled latent factors. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 10923-10930, 2021.
14

---

[71] Wenjing Zheng and Mark J. van der Laan. Targeted maximum likelihood estimation of natural direct effects. International Journal of Biostatistics, 8(1):Article 28, 2012.
[72] Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. DAGs with NO TEARS: Continuous optimization for structure learning. In Advances in Neural Information Processing Systems (NeurIPS), volume 31, 2018.
[73] H. Zhou, Li, Li, D. Zheng. J. Yang, and X. Zhuang. Association visceral adiposity index with hypertension (NHANES 2003-2018). Frontiers in Cardiovascular Medicine, 11:1341229, may 2024.
[74] Yukun Zhou, Siegfried K Wagner, Mark A Chia, An Zhao, Moucheng Xu, Robbert Struyven, Daniel C Alexander, Pearse A Keane, al. Automorph: automated retinal vascular morphology quantification via a deep learning pipeline. Translational vision science & technology, 11(7):12- 12, 2022.
[75] Shengyu Zhu, Ignavier Ng, and Zhitang Chen. Causal discovery with reinforcement learning. arXiv preprint 1rXiv:1906.04477, 2019.
15

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
The abstract and introduction of this paper precisely articulate the core contri- butions and intended scope of our proposed SIGMA framework. We clearly define SIGMA as an innovative end-to-end mediation pathway inference method that unifies structure discovery, pathway identification, and effect estimation, with strong emphasis its unique capability to systematically handle structural uncertainty. These claims are fully consistent with the advanced methodology and comprehensive experimental results detailed in the subsequent sections, accurately reflecting the value of our work.
Justification:
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: We discuss the potential limitations of our work in Appendix A provide in-depth exploration of the specific considerations and strategies for addressing core causal assumptions (e.g., Causal Markov, Faithfulness, and Causal Sufficiency) within the SIGMA framework in Appendix A.2 "Robustness Considerations and Sensitivity Analysis SIGMA." " This demonstrates our careful consideration of the applicability boundaries of our method and directions for future improvements.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
This paper provides solid theoretical foundation for all theoretical results. We systematically elaborate on the definitions, notation, and key identification assumptions (including consistency, ignorability, positivity) for causal mediation analysis.
Justification:
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Answer: [Yes]
Justification: We are committed to ensuring high reproducibility of our experimental results. We believe this sufficient and transparent information provides a solid foundation for other researchers reproduce our main experimental results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: We have thoroughly described the generation methods and parame- ters for the synthetic data in Appendix, ensuring the reproducibility of these experi- ments. For the HPP dataset, we have provided link to its official knowledge base (htps://knolledgbbbeeeeeeeee
for reference.
6. Experimental setting/details
Question:
Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
16

---

Justification: We have provided sufficient experimental setting details to ensure readers can understand and evaluate our results.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: We have given due attention the statistical significance of our experimental results in the paper. Specifically, we clearly reports the bias, empirical standard error (SE), root mean squared error (RMSE), and the coverage probability (CP) of their corresponding 95% confidence intervals for the estimated Natural Direct Effect (NDE) and Natural Indirect Effect (NIE). These metrics provide robust statistical evidence for assessing the accuracy and stability of our method's s estimates.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: We discuss in Appendix.
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics http:///wwrrpp.ccccppp/////ceeee
Answer: [Yes]
Justification: We have strictly adhered to the NeurIPS Code Ethics in conducting this research. This study utilizes data from the Human Phenotype Project (HPP), which received formal approval from the Institutional Review Board (IRB) of the Weizmann Institute of Science. Furthermore, all participants provided written informed consent. This ensures that our research aligns with recognized ethical standards in all respects.
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: We provide a preliminary discussion of the potential impacts of our work.As foundational research framework designed uncover causal mechanisms in complex systems, SIGMA holds broad positive application potential fields such as biomedicine, public health, and social sciences, for instance, by aiding in the understanding of disease mechanisms or evaluating the effects of policy interventions. While our primary focus is methodological innovation, we also recognize that any powerful analytical tool must be used responsibly.
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
The SIGMA framework proposed in this paper a causal inference methodol- ogy and, as such, does not constitute a model type with a high risk of misuse (e.g., large-scale language models or image generation models). The HPP dataset we utilized is biomedical research data that has undergone rigorous ethical review, not dataset scraped from the internet or known to have high inherent risks. Therefore, specific safeguards for such models or data are not applicable to this research. account and make best faith effort.
Justification:
12. Licenses for existing assets
17

---

Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We strictly adhere to academic standards by appropriately crediting and citing the creators original owners of all existing assets (including algorithms, datasets, and software tools) used this paper, as detailed in our References section. For instance, for the HPP dataset, we have clearly indicated its origin and cited relevant primary research literature.
13. New assets
Question: Are assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: This paper itself (including the main text and uppendices) serves as com- prehensive and detailed documentation for our innovative SIGMA framework. We have systematically elucidated all components of the SIGMA framework, its solid theoretical underpinnings, a clear algorithmic workflow, and thorough experimental validation within the manuscript. This content provides a complete guide for understanding and applying the SIGMA framework.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: This research did not involve any crowdsourcing experiments. We utilized the existing HPP dataset, which includes human subjects. As detailed in Appendix, the HPP study itself adhered to strict ethical procedures, including written informed consent from all participants. As a secondary analysis of the HPP data, this paper did not directly recruit or interact with subjects, and therefore, providing additional instructions or compensation to participants was not part of our research protocol. Relevant information pertaining to the original HPP study is detailed its respective protocols.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [Yes]
Justification: We place high priority on research ethics. The Human Phenotype Project (HPP) study protocol, from which our data is derived, received formal approval from the Institutional Review Board (IRB) of the Weizmann Institute of Science. Furthermore, all participants in the HPP study provided written informed consent, implying that potential risks of the study were fully disclosed to the subjects, and their consent was obtained. This paper strictly conducts data analysis within this ethical framework.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods this research? Note that if the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, or originality the research, declaration is not required.
Answer: [NA]
Justification: Large Language Models (LLMs) were not utilized as an important, original, or non-standard component in the core methodology of this research. Our research focuses the development and validation of the SIGMA framework, which is based causal
18

---

inference techniques such as structural equation models, path stability scoring, and Bayesian model averaging. LLMs were not employed in the core algorithm design or experimental analysis of this study.
19

---

A Detailed Discussion of Causal Assumptions
A.1 Core Assumptions for Causal Structure Discovery and Identification
The validity of causal structures inferred by Flow-SEM-and by extension, the SIGMA frame- work-reliees on set of core assumptions that connect observed statistical dependencies to underly- ing causal mechanisms. In particular, the following three foundational assumptions enable reliable identification and interpretation of causal graphs from observational data.
A.1.1 The Causal Markov Assumption
Formal Definition. A probability distribution P(V) is said to satisfy the Causal Markov Condi- tion with respect DAG G (V, E) each variable X; E V is conditionally independent of its non-descendants, given its parents in the graph. Formally, for all X, E V:
X NonDese Paa(X.)) (X)().),
(8)
where (X,( denotes the set of parent nodes of X,, and NonDesca(Xi) denotes its non- descendants in G. For DAGs, this local condition is equivalent the global Markov condition defined via d-separation: for disjoint subsets X, Y, Z C V. if Z d-separates X and Y in G, then X TT Y Z in P(V).
Role in SIGMA. The Causal Markov Assumption justifies the structural factorization used by Flow-SEM in SIGMA. Specifically, enables the modeling of the joint distribution P(V We. .) as the product of conditionals:
p II P(V Pa)(V))
P(V Wo. 2)
(9)
which underpins the Structural Equation Model (SEM) formulation:
V = fi + (Pag U/
(10)
where St and a, are flexible functions (e.g., neural networks), and U, NoO 1) represents exogenous noise. This factorization supports the likelihood- training objective used Flow-SEM (see Section 4.2.1). In particular, Equation 3 optimizes the negative log-likelihood while enforcing acyclicity and sparsity constraints. Without the Causal Markov assumption, the inferred structure may not reflect the true causal dependencies, when the model achieves high likelihood under the observed data distribution.
A.1.2 The Faithfulness Assumption
Formal Definition. A distribution P(V) is faithful to DAG G (V,E) if every conditional independence observed in P(V) corresponds to a d-separation G. Formally, for any disjoint subsets X, Y. CV:
X I Z in P(V) + d-separates X from Y in G.
(11)
The forward implication (=) corresponds to the Causal Markov Assumption, while the reverse (+) constitutes the Faithfulness Assumption.
Role in SIGMA. Faithfulness is critical to ensuring that Flow-SEM accurately infers the presence or absence of causal edges. Violations occur when distinct causal pathways counteract each other such that their net effect cancels out. For instance, in competing mediation paths:
BAM, BM, = -BAM, BM,Y,
(12)
which can induce conditional independencies not reflected by any d-separation in the true DAG. While such exact cancellations are measure-7ero events in parameter space, they can still arise in practical scenarios.
20

---

Violations of faithfulness in SIGMA result in biased structure learning. Specifically, the learned posterior P(G, D) may diverge from the true posterior over the causal graph G*:
P(G, D) + P(G") (),
(13)
leading Flow-SEM favor incorrect structural hypotheses. This has downstream effects on mediation pathway identification, particularly through the computation of the Path Stability Score (PSS), defined as:
NOAG Σ I(TT C.).
PSS(
(14)
N
where ı denotes a candidate mediation path, and the indicator function returns 1 all directed edges in π exist in the sampled DAG G,, As a result, true causal pathways may be assigned artificially low PSS values due systematically missing edges, reducing their chances of being selected for downstream effect estimation.
A.1.3 The Causal Sufficiency Assumption
Formal Definition. A variable set Vobs is causally sufficient if for every pair X, Vobss all common causes are also in Vobs. Equivalently, no unmeasured variable U & Vobs is a direct cause of two more variables Vabs.
Role in SIGMA. Causal sufficiency plays a pivotal role across all components of SIGMA. In the structure discovery phase, unmeasured confounders can distort edge estimation, leading to divergence between the learned and true graph posteriors:
P(G D) P(G" [D),
(15)
where G denotes the true underlying causal graph.
In the effect estimation phase (Section 4.4), causal sufficiency affects the correctness of nuisance function estimation. These functions, which are essential for computing NDE and NIE, depend on the appropriate conditioning sets defined by the graph:
i (P)..
(16)
When G, is misspecified due unmeasured confounding, the adjustment sets are invalid, thereby violating key identification assumptions such as:
(a, m), M(a) C,G. and M IL (aa, A= a, C.G.
(17)
SIGMA assumes that the full variable set V = (A,Y)U M U is causally sufficient, or that any substantial unobserved confounding accounted for within C. When this assumption is questionable, violations can be partially mitigated through sensitivity analysis.
A.1.4 SIGMA Framework to Challenges in Core Causal Assumptions
The efficacy inferring causal structures G = (V, E) and quantifying mediational pathway effects from observational data D is intrinsically linked to a set of foundational causal assumptions. The SIGMA framework incorporates specific mechanisms in its design to consider and respond the challenges these assumptions may face in complex real-world systems.
Considerations regarding the Causal Sufficiency Assumption
The Causal Sufficiency assumption, requiring that all common causes U affecting any two variables V., V2 E V within the system are included in the observed set V (i.e., U C V), is an ideal yet often challenging prerequisite for unbiased causal effect estimation. When this condition is potentially violated, implying the existence of unmeasured confounders U & V. SIGMA offers a more robust analytical pathway through the following design considerations:
Structured Integration of High-Dimensiona. Covariates: SIGMA mandates the explicit inclusion
of set of covariates C' within the observed variable set V {A, YJ M U.C. (see Section 3 of the
21

---

main text). The Flow-SEM module (Section 4.2.1 of the main text) incorporates these covariates C into the structural learning process, treating them equivalently to other variables (A, Mr E M,, Y). This allows for the data-driven discovery of relationships between C' and A, M. Y. By employing flexible functions (Paa (V.); Φb) (Equation of the main text), SIGMA can capture complex, non-linear confounding patterns mediated by the observed covariates C, thereby adjusting for such confounding at a structural level, which surpasses traditional methods that often treat C' merely as linear control variables.
Exploration of Structural Posterior Distribution: Rather than relying on a single, potentially misspecified, predefined causal graph, SIGMA leverages Flow-SEM to learn parameterized posterior distribution P, We. 2) over causal structures. From this distribution, an ensemble of DAGs 9= (G. } Nuag sampled (Section 4.2.3 of the main text). This ensemble 9 represents a multiplicity s=1
of causal structure hypotheses compatible with the observed data D and the given variable set V. This strategy acknowledges that, even conditional on V, the true causal structure may be uncertain. Different G2 E G may account for complex dependencies arising from C' in various ways, thus providing a more comprehensive examination of confounding effects mediated by observed variables.
Guided Bayesian Model Averaging: For each candidate pathway т identified via the Path Stability Score (PSS) (Section 4.3 of the main text), SIGMA estimates sathway-speciffe effects (T) under each supporting DAG G2 E Gr C G. The nuisance functions Ats are parameterized based on the parent sets Pac, (-) (including relevant covariates from C') defined by that specific G, of the main text). The final effect estimate (RMA(T) and its variance (OnOA (=))
(Sections 4.4.1 4.4.2 are obtained by BMA over Gx (Section 4.4.3) Equation of the main text). This aggregation ensures that the final inference incorporates not only the statistical uncertainty (()())) within each specific structure but also the structural uncertainty, Varg, (+)), arising from the heterogeneity of causal hypotheses consistent with the observed data V.
Although SIGMA cannot intrinsically detect or fully neutralize biases from truly unmeasured con- founders (i.e., those not included in V), its design aims to maximize the utility of all available information within V (especially covariates C). Through structured modeling of observed covariates, explicit modeling of structural uncertainty, and aggregation over multiple structural hypotheses, SIGMA provides a more systematic and data-driven framework than analyses relying on single, fixed causal model to mitigate biases stemming from confounding by observed variables. For po- tential unmeasured confounding beyond the scope of V, sensitivity analysis, as recommended in Appendix 1.3 of the main paper, remains an important complementary tool for assessing the robustness of conclusions.
Considerations regarding the Causal Faithfulness Assumption
The Causal Faithfulness assumption posits that all conditional independence relations present in the probability distribution P(V) are entailed by d-separations in the true causal graph G", and vice-versa (Sections 3 A.1.2 of the main text). This implies the absence of statistical independencies arising from exact cancellations of path effects (e.g., Exe effect(T)) z O while IIA-.Y * 9). where underlying causal pathways exist. SIGMA's mechanisms offer the following considerations when facing challenges of approximate or exact unfaithfulness:
Continuous Edge Weight Learning and Probabilistic Representation: Flow-SEM (Section 4.2.3
of the main text) learns a continuous adjacency weight matrix We. where entries [Weliil reflect the strength of potential edge V,, + V,, These weights are transformed into edge inclusion probabilities Psi = P(e;; E E)D) via af ([]ss))) (Section 4.2.3 of the main text). In scenarios of approximate unfaithfulness, where the net effect of path is close zero, diminishing the signal for its constituent edges, the corresponding Py: values may be low but not necessarily zero (unless strong sparsity priors are imposed). This allows edges with weak signals a non-zero probability of being included in sampled DAGs, thus mitigating premature exclusion based on hard independence thresholds.
Statistical Filtering via Path Stability Score (PSS): The PSS( (Section 4.3 Equation of the main text) is defined as the frequency of path I's appearance across Nyalid sampled DAGs, Gs, serving as a Monte Carlo estimate, EG-P(G)D p) [L(TT C G)]. of its marginal inclusion probability under the learned posterior D)) If true path's 's signal is attenuated due approximate unfaithfulness affecting one of its constituent edges ejs (resulting in a low p,s), the joint probability of observing ' (and thus its PSS( (T)) will be reduced. By thresholding PSS(TT) with r, SIGMA selects paths that
22

---

exhibit sufficient statistical consistency across the ensemble, effectively focusing on pathways whose structural evidence remains discernible despite data noise and potential complexities in parameter configurations.
Diversity of Structural Ensemble and BMA Weighting: SIGMA generates diverse ensemble of DAGs G. Even if certain parameter configurations lead to the apparent weakness or absence of a specific path in some G2 G due to approximate unfaithfulness, such a path might still manifest in other G, corresponding to different plausible parameterizations, provided some evidential support exists in D. While the BMA procedure in this work employs uniform weights for selected paths, the PSS acts as a preceding filter. Paths consistently undermined by strong approximate unfaithfulness, resulting in PSS values below т, are excluded from the final effect aggregation. This prioritizes mediational mechanisms with more robust and identifiable signals in the observed data.
Consequently, while SIGMA cannot resolve fundamental unidentiffability stemming from perfect parameter cancellations, its probabilistic and tikemble-based framework allows for more nuanced, data-driven assessment of pathways whose signals might be weakened by approximate unfaithfulness. It quantifies the evidential strength for each path via PSS, thereby focusing inferential resources on mediational mechanisms that are statistically stable and discernible under the learned posterior distribution of structures.
Considerations regarding the Causal Markov Condition
The Causal Markov Condition, a cornerstone of graphical causal models, states that any variable V, is conditionally independent of its non-descendants (excluding its parents) given its direct causes Pac(Vi) in DAG G' (Section A.1.3 of the main text). SIGMA's design inherently relies on and leverages this assumption:
The Flow-SEM module (Section 4.2.1 of the main text), central to SIGMA's structure discovery, models inter-variable relationships using Directed Acyclic Graph (DAG). The structural equations V; = f.Paa (Vi); ø,) + a, (Pag Φ;) U, (Equation 1 of the main text) explicitly encode that the distribution of V is determined solely by its parents Pac and an exogenous noise term U..
DAG-based Parametrization:
The lielliiel function log P(V() Wo: + in Flow- SEM (Equation 2 of the main text) is expressed as sum over terms corresponding each variable log o, (Pac (V: o:)) This
Factorization of Likelihood Function:
V; conditioned on its parents (Va i.c., E (oogpu(u)))
factorization is direct mathematical consequence of the Causal Markov Condition, rendering the learning of high-dimensional joint distributions computationally tractable.
Flexible Functional Forms and Structural Exploration: By employing neural networks for fa and oi SIGMA accommodates complex non-linear dependencies. The learning of Wa and subsequent sampling of the DAG ensemble G allow SIGMA explore a rich space of DAG structures that are compatible with the data and inherently satisfy the Markov condition. This combination enables SIGMA to seek the best description of observed data dependencies within the class of Markovian models representable by DAGs.
A.2 Robustness Considerations and Sensitivity Analysis of SIGMA
The theoretical guarantees for SIGMA framework validity rely series of core causal assumptions, including the Causal Markov condition, Faithfulness, and Causal Sufficiency defined in Appen- dices as well as effect identification conditions such Ignorability and Positivity specified in Section 3 of the main text. In practical application contexts, these assumptions may be violated to varying degrees, threatening the validity and reliability of the estimated mediational pathway effects. Therefore, systematic analysis of the robustness and sensitivity of results obtained under the SIGMA framework becomes particularly necessary. Next, we discuss potential methodological frameworks for sensitivity analysis and mathematical prospects for future research.
A.2.1 Sensitivity Analysis under Deviations from Causal Sufficiency
Consider an unobserved potential confounding variable U that simultaneously affects exposure A, the set of mediator variables M, and the outcome variable Y. the presence of such unobserved confounding, we define confounding sensitivity analysis framework modulated by U as follows:
23

---

Assume the confounding variable U is simultaneously associated with exposure and mediators, as well as outcomes and mediators, with their association strengths described by the following two sensitivity parameters I AA Tuy:
sup,, Pr(A a u. C) inf, Pr = a C) TUy = sup, E| m, -
T AU
inf, E E M = m, U -
Here TAU.Tur 2 1, with larger values indicating stronger influence of the unobserved confounding variable. To quantify the robustness of observed effects to the strength of the above unobserved confounding, the E-value [6 can be used, defined as:
value((nma = AC, 2 11:0% 0 E Clr (T)))(T)))
where Au Cu.y (#)) represents the confidence interval for the effect estimate adjusted for U given the sensitivity parameters. A larger E-value indicates that the conclusion is more robust to unobserved confounding; conversely, smaller E-value indicates weaker robustness. Future research could consider further integrating the above sensitivity analysis automatically into the SIGMA process, using simulation studies to quantify its sensitivity performance under various data-generating mechanisms.
A.2.2 Robustness Analysis under Deviations from Causal Faithfulness
Suppose the data-generating mechanism has approximate effect cancellations, where the true effect of some path in the path set A- is ((T) x 0, despite this path existing in the true causal structure (i.c.. approximate unfaithfulness). This scenario leads to the estimated Path Stability Score (PSS, Section 4.3) tending toward lower values. Sensitivity analysis considers path stability threshold T E [0, I]. defining stable path set as:
II, { PSS(#) 2 +}
By studying the impact of changes in T on the final Bayesian Model Averaging (BMA) estimate (FMA( we can examine the sensitivity of (TT7( to different
(0 Ör
π IT,
In applications, if low sensitivity is found (the above derivative is small), indicates that the SIGMA estimation conclusion is relatively robust; conversely, if high sensitivity is found, it may suggest that approximate unfaithfulness significantly affects the analysis, requiring cautious interpretation. Future research could deeply explore the properties of this derivative and its theoretical and empirical performance under different structural assumptions.
A.2.3 Robustness Analysis of Model Specification Errors in Effect Estimation
When estimating path-specific effects δEMA (T). SIGMA involves multiple nuisance functions (e.g.. propensity scores and conditional expectation functions). Setting the bias of these nuisance function estimates as 8, the following form of bias sensitivity analysis can be considered:
Define the sensitivity function of estimator bias nuisance function estimation bias:
30 a8
Lower sensitivity values imply higher robustness of the estimator to model specification errors, and vice versa. The Cross-Fitting strategy adopted in SIGMA provides a certain degree of robustness guarantee, but future research should systematically compare the sensitivity function characteristics of this estimation strategy with other semiparametric doubly robust estimators.

---

A.2.4 Interplay
The three assumptions Causal Markov, Faithfulness, and Causal sufficiency-joiinne define the identification conditions for reliable causal discovery in SIGMA. The Causal Markov condition links the graphical structure to conditional independence relations. The Faithfulness assumption ensures that all observed independencies reflect actual graphical separations rather than coincidental parameter cancellations. Causal Sufficiency guarantees that all relevant common causes are observed, preventing bias from unmeasured confounding. Together, these assumptions enable SIGMA to learn meaningful causal structures from observational data and propagate this structure into consistent mediation effect estimation.
A.3 Identification Assumptions for NDE and NIE
TThe nonparametric identification of Natural Direct Effects (NDE) and Natural Indirect Effects (NIE) requires set of conditional independence assumptions that account for confounding. These effects are defined in the potential outcomes framework (see Section 3 as: NDE(a, a²) B[Y (a, M(a))
Y(a", M(aa( (a")) and NIE(a, a") = (a ((a) Ya, M M(a") These assumptions ensure that observed statistical associations admit valid causal interpretation. In the context of SIGMA, this is
particularly important because the underlying causal graph G is unknown and approximated via a sampled ensemble {G.).
A.3.1 Ignorability Assumptions for NDE and NIE Identification
common set of sufficient conditions for the identification of NDE and NIE, often referred to as "sequential ignorability" or "no unmeasured confounding" assumptions, can be stated conditional on a set pre-treatment covariates C 45
For treatment levels a and a*:
1. Conditional independence of treatment and potential outcomes/mediitooss (Controlling for A-M and A-Y confounding): There is no unmeasured confounding of the effect of A on M, and no unmeasured confounding of the effect of A on Y that is not mediated by M. Formally:
Va' E :: A M(a') C, Ya", Vm M: IL Yaa,, m) C.
(18)
Together, these imply that C' blocks all back-door paths from A to M, and from A to Y that do not go through M. A common shorthand formulation is:
(a m), M()) L LAC.
(19)
Conditional independence of mediator and potential outcomes, given treatment and covari- ates (Controlling for M-Y confounding): This assumption requires that, conditional on treatment A and covariates C, there is no unmeasured confounding between the mediator M and the outcome Y:
Va E A, Vm' E M:: M IL Y m') A = a.
(20)
Under this assumption, among units with the same treatment and covariates, the mediator behaves as if it were randomly assigned with respect to the counterfactual outcome Y(a, m²).
3. Conditional independence of treatment and potential outcomes, given mediator and covari- ates (No treatment-induced M-Y confounding that also depends on A directly): This assumption rules out any variable that is affected by treatment and also confounds the relationship between M and Y. In some formulations, this condition is written as:
A IL (a", m) M m, C.
(21)
However, this expression involves cross-world counterfactuals which are often avoided in practice. Instead, most identification strategies rely on morle-world assumptions derived from the causal graph.
For example, the mediation formula expression:
45
provides means to identify NIE using the following
(a, M(a)))
E A a, M m. = d) P(M = m A -== c) P(C = c).
(22)
25

---

This formula is valid under Assumptions I and 2, as well as the Causal Markov and Faithfulness assumptions for the underlying DAG encoding the relationships among A, Y. Within the SIGMA framework, pathway-specife effect estimation relies on Efficient Influence Functions (EIFs), which are computed relative to each sampled graph Gx. For example, to estimate the effect along a mediation pathway A > M, - Y, SIGMA computes the following nuisance quantities:
Y A, Mj. Pac, (Y) (A, C) P(M, A, Pac, (M,) \ {A), C), P(A Pac,(A), C).
(23)
(24)
These estimates rely on graph-structune adjustment sets and implicitly assume that the correspond- ing back-door paths are blocked. As such, the validity of based identification hinges on the correctness of the local parent sets derived from G.. In cases where no variable is jointly affected by A and also acts as a confounder of the M + Y link, the stronger cross-world independence condition Y(a", m) II M(a) C' can be avoided 49
Role of Graphical Structure G, in Satisfying Ignorability within SIGMA. In SIGMA, the Flow-SEM module generates an ensemble of candidate DAGS {Ga]. For each graph G. = (V, Es), the estimation of nuisance functions i), (see Section 4.4.1 such as the conditional expectation BYY A, M, Zy] and the conditional probability P(M A, 7M), relies on structure- structure-specific conditioning sets. These sets constructed from the local parent sets Pag. (-) defined by each DAG, together with the baseline covariates C. Specifically, we define:
CUU (Y) (A, ZM = Pag, (M) {A).
(25)
This formulation reflects the assumption that causal relationships are encoded in G., and that proper conditioning on parent sets and covariates is sufficient to block non-causal pathways. The theoretical justification rests on Pearl's d-separation criterion and the back-door adjustment principle [46].
Under ignorability Condition 1, which requires IL M(a)) | C and A Yaa", m) | |C, set of variables-denoted C" Cust block all back-door paths from A to M, and from A to Y that are not mediated through M. These paths must begin with an arrow into A, and blocking them ensures that the necessary conditional independencies hold in Gg. For Condition 2, which assumes M IL Y(a, m') A = a, C, the relevant conditioning set must block all back-door paths between M and Y, conditional where conditioning on A = accounts for any paths passing through A as either collider or a non-collider.
A This is achieved by the set U PaG, (Y) \ (A, M) Pag, (M) \ {A),
The central premise is that each sampled DAG G2 is a candidate for the true underlying causal graph G*, and that the observed variable set Vots {A, Y) UM U C satisfies causal sufficiency. When these conditions hold, and if G, accurately reflects the structural relationships encoded in G*, then using the parent sets Pag. (-) in conjunction with C provides a principled approach for selecting conditioning variables based on the d-separation criterion. For instance, if Pag, (M) includes all common causes of M and Y that not affected by A and are not in C, then conditioning on A, C, Pag. (M) \ (A} is sufficient to satisfy the M-Y ignorability condition:
M 1I (a, m/) A = 0, CU PaG, (4).
(26)
However, since each G2 is only a hypothesis, ignorability is not guaranteed unless G, = G* and the covariate set C' is truly sufficient. If either the graph is misspecified or unmeasured confounders are omitted from both C and the parent sets, the resulting conditioning set C'PPa, (-) may be inadequate. This misalignment can lead to biased nuisance function estimates and consequently biased causal effect estimates (T)) While the BMA procedure in SIGMA propagates structural uncertainty by averaging across sampled DAGs, it does not inherently resolve violations of ignorability that are present across all candidate structures. Therefore, rigorous selection of covariates C' remains essential for ensuring the validity of effect estimation under structural uncertainty.
A.3.2 Consistency Assumption
The consistency assumption establishes the fundamental link between potential outcomes (defined under hypothetical interventions) and observed outcomes in the data. It posits that for any individual,
26

---

the observed outcome corresponds precisely to the potential outcome associated with the treatment and mediator values that the individual actually received.
Formally, for all individuals i, the assumption requires:
A, = a M, Mi(a), and (A: = a, -)) Y; = Yi(a,
(27)
which means that if an individual received treatment a, the observed mediator equals the potential mediator under a; and if the individual received treatment a and mediator value m, the observed outcome equals the potential outcome under those values.
These conditions can be restated as follows. For all individuals i, and any treatment level a E A. A; = a, then M, = (a).). This expresses that the observed mediator is what would have been realized under the assigned treatment. Similarly, for all E A and all m E M, if A, M; = m, then Y; = (a m). This implies that the observed outcome reflects the potential outcome under the observed treatment and mediator values. These requirements must hold for all relevant treatment levels (e.g., a, a in the definition of NDE and NIE) and all values of m within the support of M.
and
Consistency is essential for expressing counterfactual estimands such as RI (a, M(a))) in terms of observable quantities. Under this assumption, conditional expectations from the observed distribution P(Y, A, M, C') can be equated to expectations over potential outcomes. For example, the observed expectation E|| A a, M m, C c] corresponds to [E[Y (a, m) A a, M m, C c]. assuming that ignorability holds for the conditioning set.
In the SIGMA framework, consistency enables the interpretation of structure-speciffe nuisance functions as counterfactual components. These include, for each sampled DAG G,:
[ A, M,, Pag. (Y) \ (4. C), P(M PaG, (M) \ C)).C),
(28)
which are estimated from observational data and treated as valid components of ElF-based iden- tification strategies (see Section 4.4.1). This directly affects single-DAG DAG effect estimates (z(
(Section 4 and, through Bayesian model averaging, the final estimates (TM
Violations of consistency occur when the treatment A a or mediator M = m is ambiguously defined in the data. For example, if multiple unobserved variants of A a exist (e.g.. "exercise intervention" encompassing both "light jogging" and "intense weightlifting"), and these variants have distinct effects on downstream variables, then consistency does not hold, even though the data records only a single treatment label. The same logic applies if M is coarsely defined.
SIGMA implicitly assumes that treatment A and mediator M are sufficiently well-defined and measured with adequate granularity such that their causal interpretation is meaningful. Ensuring consistency is not typically addressed at the algorithmic level, but rather requires careful study design or strong domain knowledge to define, operationalize, and measure treatments and mediators in a causally coherent way.
A.3.3 Positivity Assumption
The positivity assumption, sometimes referred to as overlap or the experimental treatment assignment (ETA) assumption, is essential for the nonparametric or semiparametric identification of causal effects, including NDE and NIE. It states that, for any combination of observed pre- trettretment covariates C c, each treatment level A = a and each mediator value . = m must occur with non-zero probability. This ensures adequate empirical support across all necessary strata, which is required to estimate causal contrasts without extrapolation.
Let A denote the set of treatment values, M the support of the mediator(s) M, and C the support of covariates C. The positivity condition for NDE and NIE identification can be stated as:
Ya E A, Vc E C s.t. P(C c) > 0: P(A = C = c) > eA Ya A, Vm E M, Vc C s.t. P(A a, C = c) 0: p(M m A= CM (30)
(29)
where e > and 6M > 0 are small constants ensuring strict positivity. For continuous mediators, the density p(M = m A = a, C = e) must be bounded away from zero over the relevant support of M. In practical applications, especially when using estimators involving inverse probability weighting
27

---

(e.g., ElF-based estimators), stronger bounds such as CΛ S P(A= C= 5 c) CA are often imposed ensure numerical stability.
In SIGMA, positivity required for the estimation of nuisance functions nJ. (Section 4.4.1 such
as the propensity score P(A Pac, (A), C) and the conditional mediator distribution P P(M; A, Pac. (Ayy) \ {A], C"), for each sampled DAG G.. When Equation 29 or Equation 30 violated in practice-i.e. -ic some combinations of A, M, C' occur with near-zero probability the corresponding weights in EIFs or related estimators become excessively large. This can lead to inflated variance or unstable effect estimates (T).
For example, if pMM covariate space, estimating expectations such as E[ A a, M = m.... requires extrapolation into areas with limited or no observed data. This weakens the reliability of the learned nuisance models and introduces instability in the resulting causal effect estimates.
A = a, Pac, (M) \ (4), C c) is close to zero in some region of the
These issues propagate to the final pathway-specific estimates through their effect on the single-DAG estimators (T)), and subsequently affect the BMA-aggregated estimates θEMA (17) Although SIGMA employs flexible nuisance function models (e.g.. neural networks), which may tolerate near-violations better than simpler parametric approaches, such models can still produce biased high-variance estimates when positivity is violated. This underscores the importance of verifying covariate support overlap in practice and ensuring estimation occurs within sufficient empirical support.
Interaction with Other Assumptions. Positivity is intrinsically linked to other core assumptions:
Ignorability: Ignorability requires adjustment for a sufficient set of confounders C' (e.g.. CUPac, (-)). Positivity ensures that such adjustment is empirically feasible by guaranteeing that all levels of A (and M) are present across all strata of C². Without Positivity, even if C" theoretically suffices for Ignorability. the lack of common support can make the adjustment ineffective or impossible (Appendix (A.3.I).
Causal Sufficiency: If Causal Sufficiency is violated, the set C may be inadequate. Even if Positivity holds for the observed C, the inability to adjust for unobserved confounders means Ignorability may still be violated, leading biased estimates irrespective of data overlap on C' (Appendix A.1.3.
Practical Diagnostics for Positivity Violations. While SIGMA does not include automated pro- cedures for testing positivity, several diagnostic practices are recommended to identify potential violations:
Examine Propensity Score Distributions: Plot histograms or density plots of estimated propensity scores P(A a | adjustment set) for each treatment group. Lack of overlap or scores concentrated near 0 or 1 indicates potential Positivity issues.
Assess Covariate Balance: After stratification or weighting by propensity scores (if ap- plicable a diagnostic step), check for balance in covariate distributions across treatment groups using metrics like Standardized Mean Differences (SMDs). Large SMDs (c.g.. > 0.1-0.25) can suggest poor overlap.
Inspect Weight Distributions: For estimators that explicitly inverse probability weights (even internally, like some EIFs), examine the distribution of these weights. A few extremely large weights can dominate the estimate and inflate variance, signaling Positivity problems.
Check for Sparse Strata: For conditional mediator distributions P(M A, adjustment set), particularly with discrete mediators or high-dimensional adjustment sets, tabulate or examine the number of observations in critical (A, M, adjustment set) strata to identify sparse data regions.
Sensitivity Analysis: Evaluate how NDE/NIE estimates change if observations with extreme propensity scores are trimmed or if the analysis is restricted to regions of better covariate overlap (e.g., trimming based on propensity score tails).
Addressing severe positivity violations may require restricting the analysis to subpopulations with sufficient empirical support or applying estimation techniques designed for limited overlap (e.g.. overlap weights or covariate balancing methods). These approaches can help reduce estimation
28

---

variance but may also introduce new assumptions or trade-offs, such as reduced generalizability or increased bias.
A.3.4 Single-DAG Effect Estimation via Influence Functions and Plugin Estimators
This appendix details the specific construction of Efficient Influence Functions (EIFs) for estimating the Natural Direct Effect (NDE) and Natural Indirect Effect (NIE) within the SIGMA framework when the exposure variable A is binary (taking values 0 or 1) and the path length is (i.c.. A - M - Y). These EIFs depend on series of nuisance function estimates obtained through K-Fold cross-fitting (see Section 4.4.1 in the main text for details). We define the following core nuisance functions, where the superscript (s) indicates that these functions are estimated for specific DAG Gg. but will be omitted in the following for brevity:
á (C = P(A = A]Ca): The propensity score for exposure A, i.e., the estimate of the conditional probability of A 1 given covariates C'A (parent nodes of A in G.). For the mediator M (which can be continuous, binary, or categorical):
M Cm) É[M A a, CM]: When M is continuous or binary, the estimate of its conditional expectation given A and covariates Cm (parent nodes of M in G.. excluding A).
far(ma, CM) P(M = A = a, Cm): When M is discrete categorical, the estimate of its conditional probability mass function of taking value given A and CM. For unified notation, we sometimes still use (Mg A a, C'r] to denote the expectation of some function ((M) of M, which can obtained by summing/integrating over fM-
juy (a, m, Cy) E| a, .M m, Cy]: The estimate of the conditional expectation of outcome Y given A = a, M = m, and covariates Cy (parent nodes of Y in G., excluding M))
Furthermore, define the plug-in estimator for counterfactual expectations (obtained via G- computation using the aforementioned nuisance functions and averaged over the sample):
"e'a" = (a")) = - Ei= fay m, dd represents summation over all possible values of the mediator M (discrete M) integration (continuous M), and (fMM(maaa Car,a) represents the estimated conditional distribution based on f C.s,s) or fmmaa""" Ca,,s).
m,
CM.i), where the integral sign
NDE and NIE can be expressed as differences of these counterfactual expectations:
NDE
- (0, M())) Y10 Yoo
NIE = (1, M(1) (1, (()) = 711 7110
Therefore, the influence functions for NDE and NIE are linear combinations of the influence functions of their constituent counterfactual expectation terms. We present the influence function Vi (Z;) for E|Y(a' M(a))) ))] (i.e., "))a"). For observed data Z; = (A;, M,, Y,, CA.,, CM,., Cy,s), the typical form of Wiattn (Zi) is as follows:
T(A, d') SM CM..) (V: - (ay Mj. Cv.a))
P(A, C'(CA,s) FM(Ma
I(A = aa)
(py(d', M,, Cv,,) - [ax m,, Cy.)d (mla", CM,s)
P(A: =a"CAa)
(firld m, Cy.))dF (m))a" (mlla,CC..) -Rara")
(31)
Here, P(A, CAA)). (CA.)) when a' 1, and 1 =CCA)) when a' When M is continuous, far (mla, CMr) is the conditional probability density function; when M is discrete, is the conditional probability mass function. If the mediator model directly estimates the conditional expectation P (a, CM) (e.g.. M is binary or continuous, and some terms in the EIF are simplified for linearity in M), then the parts of the formula involving far are adjusted accordingly. For example, when M is binary (0, 1). fMr (M CM.)) (a Cm.)) M: pus (o, I
29

---

Based on Equation (31. the influence functions for NDE and NIE are respectively:
VNDE Ze; 17s) = 4/330 (2; 1s) 1-o0 (Z;; 17s)
(32)
NIE Ze5 17.) ) y/y (1, 1/y1 (Zg5 (1))
(33)
In its implementation, the SIGMA framework computes the V/NDE( (4) and N/NIE Z4)) scores for each sample according to the above structure, using nuisance function estimates fis obtained via cross-fitting.
A.3.5 Definition and Analysis of the Influence Function for Continuous Exposure in SIGMA
For mediation paths of length ((T) (i.e., A - M + Y) involving a continuous exposure variable A, the SIGMA framework defines and utilizes the influence function described herein to estimate the Natural Direct Effect (NDE) and Natural Indirect Effect (NIE). This section details the mathematical construction of this influence function and the nuisance components it relies upon.
Nuisance Function Components for the AIF
The construction of this influence function predicated on several key nuisance function components, estimated via K-Fold fross-fitting (as detailed in Section 4.4.1) for each (path म, DAG G.) pair. Let Z; (A,, Mj. Y,, CA.s. a Cy.,) denote the data for the i-th observation, where Cx,s represents the values of the parent set (covariates) of variable X under G, for sample i.
Conditional Expectation Model Predictions for Exposure A: The cross-fitted prediction for the conditional expectation of A given its covariates CA,i. obtained from a regression model (a multi-layer perceptron in SIGMA), denoted as:
CA = ETA,/CA.]-
Global Estimate of Conditional Residual Variance of Exposure A: SIGMA assumes homoscedas- ticity for Var(A,(CA.)) positing Var(A. = a Ca.res for all This global residual variance,
denoted 0² AACA is obtained by pooling the cross-fitted residuals e; A; (Ca,s, from all samples and calculating their sample variance (defined as E/= (e, - 2)² or - Ej E)?. corresponding to the specific implementation).
Conditional Expectation Model Predictions for Outcome Y under Intervention on A with Observed M, Cy: For a given intervention level a on A (e.g., AA or μA + kō A as defined in the main text), the cross-fitted prediction of the conditional expectation E|Yi[A, = a, M Cy, is denoted as:
Mas
Individualized Plugin Predictions and Sample Means of Target Counterfactual Expectations: For any counterfactual expectation composing NDE and NIE, such as E[ (a', M (a"))) (denoted ")
The plugin prediction of the counterfactual outcome (a Mi(a")) for sample denoted as "a'a" .i. This prediction is obtained via G-computation, recursively combining the nuisance functions (1)- (3) above and outputs from the mediator model MMIJA.C (i.e., E[M]A a" Carl or P(M m|AA A = a", C A)d the sample mean of these individualized plugin predictions is denoted as "a'a" n E, fa
1
Here, aa", represent specific exposure intervention level combinations required for defining NDE or NIE.
Mathematical Definition of the Influence Function
Within SIGMA, the influence function score Vhea. (Z:) for estimating a target counterfactual ex- pectation (e.g., Naa- ( (a, M(a))) is defined as (corresponding to Equation (6) the main text):
A, (CAs) (.000
Year (Z
(Ki9A - (maa)) (frar,: - faa-). (34)
02 AlCeees
30

---

This influence function comprises two main parts:
textbfFirst Part (Exposure Model Correction Term): A, A.-Aa(Cas) (p A Mass Cy -Paa-) The
term Ar-DA(CA.)) is the score function of the log-conditional likelihood of A|CA with respect 8 Arc
to its mean parameter i (CA,s), under the assumption that A|Ca follows Normal distribution ((AA(CA,)) & AeC), This score function is multiplied by a term related to the outcome model,
(p May.C Pag)), ). which reflects the deviation of the conditional outcome expectation (under intervened A a, observed M, Cy) from the mean target counterfactual expectation.
Second Part (Centered Plugin Term): (fae- 9aa-). This term is the centralization of the individualized plugin estimate You+, , of the target parameter with respect to its sample mean Pea-
The overall influence functions for NDE and NIE are then derived from linear combinations of these influence function scores for their constituent counterfactual expectations (i.e., 710, 7000 M1, where the numerals O and represent the reference level a* and target level a, respectively), according to the definitions of NDE and NIE (NDE 710 700. NIE 711 710). Specifically:
VNDE Z; = 77 Z;; 1.) - 10000 (i,)
(35)
VNNIE Zz5 1).) ) = 4/ (7;; 1i.) 1/-10 (7,; m)s))
(36)
In these combinations, each y (Z₃; fta) term on the right-hand side adopts the structure of Equa- tion (34), utilizing the specific nuisance function components corresponding to its target counterfactual expectation (710. Too, or 71).
Construction Rationale and Discussion
The influence function for continuous exposure defined within the SIGMA framework (Equation 34 is motivated by the goal of providing estimation method that offers first-order bias correction over simple plugin estimators while maintaining computational tractability. Its structure incorporates elements common in influence function theory, such as leveraging score function information from the exposure model and centering plugin estimates.
The definition of this influence function does not explicitly include direct correction terms related to the mediator model P(M A, Cx) beyond its implicit role in the G-computation of fea-,.. In semiparametric theory, fully efficient influence functions for mediation effects often contain more complex terms explicitly dependent on the parameters or score functions of the mediator model(s) to ensure properties like double robustness against misspecification of either the outcome or mediator models. The influence function employed by SIGMA, while incorporating corrections related to the exposure model and centering the outcome predictions, may differ in its statistical properties, particularly double robustness and semiparametric efficiency, compared to such fully specified EIFs. Its performance is thus more reliant on the accurate estimation of all nuisance models, especially the conditional mean of the exposure A (CA) and the conditional expectations related to the outcome j Mar., The adoption of global homoscedastic approximation for the conditional variance of the exposure, 02 Ca.res a choice made for implementational simplicity.
Despite these specific implementational choices, by incorporating corrections related the exposure model, this influence function is expected to offer more robust performance in finite samples compared uncorrected plugin estimators, particularly when the exposure model is well-estimated. A detailed theoretical analysis of this influence function under various data-generating mechanisms and an in-depth study of its finite-sample behavior constitute important directions for future research.
B VAE-based Imputation Module
To address missing data, SG-RMA incorporates a Variational Autoencoder (VAE) specifically adapted for imputing mixed-type tabular data. Formally, consider the raw data matrix X NxP, which may contain missing entries, along with a binary mask matrix M E {0, NxP indicating observed (My = 0) or missing (My 1) values.
Data Preprocessing. Raw features are preprocessed according to their data types: continuous variables undergo standardization, whereas categorical variables are transformed via one-hot encoding,
31

---

expanding the feature dimensionality from P to P'. The resulting processed dataset is denoted as X' E N Missing entries within X' are initially filled with zeros to facilitate VAE training, and the corresponding mask M is concatenated to guide the imputation process explicitly.
VAE Model Specification. The VAE consists of an encoder network (= M) and a decoder network Po (X Given an input vector M;), the encoder produces parameters of a Gaussian posterior distribution in the latent space: mean vector Mız; and variance vector a The latent representation 2; is then sampled via the reparameterization trick & = Mz, + oz @ 6, where ~ <0, 7). The decoder reconstructs data from latent representations, employing distinct output layers for each original feature type: linear layers for continuous features, logistic sigmoid activation for binary features, and softmax logits for categorical features.
The VAE optimization objective maximizes the evidence lower bound (ELBO), defined explicitly for observed entries as:
LvAe(ò, c, 0; m,) (xollt)) [log po (z
(D lIp(=)),
(37)
where oho denotes the observed elements of x/e The reconstruction likelihood log Po 24) aggregates squared errors for continuous features, binary cross- cross-entropy losses for binary features, and categorical cross-entropy losses for categorical features. The prior distribution p(=) standard normal, and /3 is hyperparameter balancing reconstruction quality and latent regularization. Training employs the Adam optimizer, adaptive learning rate schedules, and early stopping criteria based on validation performance.
Imputation Procedure. After model training. missing data imputation proceeds as follows: For each sample i, the trained encoder maps the filled input (X iile M,) to the posterior mean latent representation 1'2.- Passing this latent vector through the decoder yields reconstructed outputs ii-- Reconstruction outputs are then transformed back into the original data space, applying inverse scaling to continuous variables, thresholding via logistic sigmoid for binary variables, and category assignment via argmax operation for categorical variables. Finally, missing entries within the original dataset X are replaced by their respective reconstructions, while observed values remain unchanged. The imputation procedure is summarized in Algorithm
B.1
Path Stability Score: Theoretical Properties
B.1.1
PSS as Posterior Probability Estimator
Section 4.3 of the main text introduces the Path Stability Score (PSS(T)) as a measure to quantify
the evidential support for mediation pathway æ. It is defined PSS(T) N.nd
- ENaad I( C G.).
where {G,) W is an ensemble of Nvalia valid Directed Acyclic Graphs (DAGs) sampled based on the output of the Flow-SEM stage (ectioo|4.2.3// This section provides a rigorous justification for the interpretation of PSS(m (T) Monte Carlo estimate of P( C D), the posterior probability that path π is present a DAG sampled from the (approximate) posterior distribution over graph structures learned by SIGMA given data D.
Let P)G denote this SIGMA-approximatt posterior distribution over the space of DAGs G. derived from the Flow-SEM model and the subsequent DAG sampling mechanism described in Section 4.2.3 The ensemble {G,} Nuaii consists of Nvalid valid DAGS treated as independent and identically distributed (i.i.d.) samples from D).
Consider a specific path ग. We define an indicator random variable (+(G) = I(TT C G), which equals if path is subgraph of G, and otherwise. The posterior probability of path æ under D) is the expectation of this indicator variable:
P(T C GID) EG-
(38)
1. EN'r X,, (G.). is the sample mean of X, (G) based on The Path Stability Score, PSS(T) = NNau these Naald i.i.d. samples.
Theorem (Statistical Properties of PSS). The Path Stability Score PSS( is:
1. An unbiased estimator of P(T C G)D), i.e., E[PSS(m) = P( C G(D).
32

---

Algorithm VAE-based Missing Value Imputation for Mixed- Mixed-Type Data
Require: Trained VAE (encoder lø. decoder po), Dataset X TR with missing values, Feature types ) where Tj E {continuous, binary, categorical). Preprocessing parameters (scaling
factors (j. for continuous features, one-hot encoders for categorical features).
Ensure: Fully imputed dataset Ximputed E TR
1: M +- CreateMissinggMas (a)k
2: X Pe-eproocessFeatuures// (X, (T₃))
3: Xilled + X Fill missing entries in Xiited (identified by M) with 0
4: Ximputed X
5: for
to N in batches do
6:
lo, Encoder ((x (lldss M,)
7:
24 f Us.
8:
X + (ecoders(2:)
9:
Initialize X; R. for j = to do
10:
if T; = continuous then
12:
X, f XLy -0, - 145
13:
else if T binary then
X: -- [sigmoid {1)) > 0.5)
15:
else if T categorical with K2 categories then
Let {X k be the logits for all categories of feature j
17:
X + arg maxke (1.... KY) softmax ( }),
18:
end
19:
if Mij = 1 then
(X (Ximuutedss + Xis
21: end if
22: end for
23: end for
24: if HasRemminnnggMisingg (Ximpuu then
25:
warning "Imputation resulted in remaining NaNs."
26: end if
27:
return Ximputed
2. A consistent estimator of P(T C G)D), ie., PSS(T) 4 P(T C G\D) Nyali cc.
Proof. (1) Unbiasedness: Let {G.) Nuga be the of Nyatia valid DAGs sampled i.idd from P(G
Nusus
[[[SS(=))
X=.),)
N
Nand
((x)((,))
(by linearity of expectation).
N
Since each G, ~ P)G (i.i.d.). E ((x)(G.)) = Bo-P(a))(X+(C)) = P(r S G).. Thus,
E[PSS(=))
P(T P(T G)D).
Naalid
(2) Consistency: The variables (x,),) are i.i.d. Bernoulli random variables with mean p P( S D) and finite variance p(1 - p). By the Weak Law of Large Numbers (and Strong Law for almost sure convergence). their sample mean PSS(TT) converges in probability to as Neali - co.
Therefore, PSS(T) provides statistically sound Monte Carlo approximation the posterior proba- bility of path under the posterior distribution D) approximated by SIGMA. The accuracy of
33

---

this approximation depends on two primary factors: (i) the quality of D) as determined by the Flow-SEM and the DAG sampling mechanism, and (ii) the number of valid DAGs Naaidd
B.1.2 Asymptotic Properties of PSS: Convergence Analysis with Increasing Sample Size
Theorem D establishes that the Path Stability Score, PSS( serves as an unbiased and consistent estimator for p. = P(T C G|D). The consisteney, derived from the Law of Large Numbers, ensures that with sufficiently large ensemble of Nvatia valid DAGS, (SS(=)) converges to Px. This section further elaborates convergence and distributional characteristics, which are fundamental for understanding its reliability and for guiding the choice of Nvalid.
the asymptotic properties of PSS(#) as Nvalid increases, particularly its of
Asymptotic Normality via the Central Limit Theorem (CLT). Since PSS( is a sample mean
of Valid independent and identically distributed (i.i.d)) Bernoulli random variables X,, X,,G.)) = [l( S
G.) (each with mean pa and finite variance Pa (1 Px)), the Central Limit Theorem (CLT) applies.
The CLT states that for sufficiently large Nvalid:
Naaata (PSS(() - p+) 4 (O. PI - Pm)),
(39)
where S denotes convergence in distribution. Equivalently, for large Nvatia, PSS(T)) itself is approxi- mately normally distributed:
(rr pn)
PSS(T) approx N (Pr
Pr
(40)
Nvalid
Variance Analysis and Rate of Convergence.
The exact variance of the PSS(TT) estimator is:
Var(PSS(T)) = P=(1 Po))
(41)
Nealid
This shows that the variance is inversely proportional to Nvalid- The standard error of PSS( is thus SE(PSS(T)) v (1= decreasing at a rate of (////Naaia). This implies that to halve the standard error, Naata must be quadrupled. The variance is maximized when pr = 0.5, indicating that more samples are needed for precise estimation when the true path probability is near 0.5.
Confidence Interval Construction and Practical Considerations for Nvalid- The asymptotic normality (Equation 40 allows for the construction of approximate Wald-type confidence intervals for Pa:
PSS(T) (1S(#)())PSS(T))
PSS(T) 21-0/21
(42)
Nvalid
where 21-0/2 is the (1-a) (1-a/2) 22) quantile of the standard normal distribution, and PSS(TT) is used to estimate Ps the standard error term. This interval provides a range of plausible values for the true posterior probability Px. The reliability of this interval and the PSS estimate itself depends on Nvalid- For probabilities p= very close to 0 or 1, fewer samples may suffice for given precision.
Connection to Bayesian Model Averaging (BMA). The statistical properties of PSS, particularly its interpretation posterior probability estimate, directly support the subsequent Bayesian Model Averaging (BMA) stage of SIGMA (Section 4.4.).. By using a threshold 7 on PSS( (T) to identify stable paths, SIGMA focuses computational resources pathways with substantial evidential support under D). The ensemble ) Nats used to compute PSS is the same basis for BMA, ensuring coherent propagation of structural uncertainty. For stable path the BMA estimate OBMA (T)) and its variance (O (i)) represent the model-averaged effect and its associated uncertainty. This averaging is performed over the plausible graph structures within the ensemble that support path effectively integrating over the structural uncertainty concerning these specific pathways as captured by P))
34

---

B.2 Discussion on path stability score threshold
The Path Stability Score (PSS) is defined the frequency of path in the posterior DAG sample set:
N Σ T( G.), PSS( = N
where и represents the candidate mediational path, G, is the s-th DAG sample drawn from the structural posterior distribution, .N is the total number of DAG samples, and II(.) is the indicator function. Given a threshold E [0, 1] for the path stability score, we define the stable path set as:
Issaa (T) = {TT PSS(T) 2T}:
The selection of threshold 7 directly determines the size and composition of Taaalles thereby further influencing the Bayesian Model Averaging (BMA) estimate of path-specific effects:
(TMA(T)T)= = (T) GRE9.(T)
Σ
(T(TT,
where (x(7) = (G,+ threshold T. and (T) is the mediation effect estimated based on the DAG sample G..
PSS(T) represents the subset DAGs containing path и given
A smaller value of r may lead to the inclusion of more low-frequency paths, increasing the risk of false path discoveries (false positives); a larger value T may exclude paths with real effeets but weaker statistical support (increasing false negatives). Therefore, in practice, the recommended strategy to select a moderate intermediate value (e.g.. t = 0.15), and then make fine adjustments based the specific characteristics of the data and prior knowledge in the relevant domain. We explicitly adopted this strategy in the numerical experiments and practical algorithm implementation in this paper, and observed that this approach can achieve reasonable balance between robustness and sensitivity in path selection.
C Related Work
Causal Structure Learning. Differentiable neural approaches have advanced causal DAG discovery by casting it as continuous optimization. Zheng et al. (2018) introduced the seminal NOTEARS [72] method, which enforces an analytic acyclicity constraint to learn DAGS via standard gradient-base solvers. This breakthrough spurred many extensions to capture nonlinear relations and improve scalability. For example, GraN-DAG 351 extended NOTEARS to nonlinear settings by using feedforward networks model causal mechanisms, and DAG-GNN [69 embedded a variational autoencoder with a graph neural network to learn DAG structures from data. Subsequent works proposed alternative differentiable formulations of the acyclicity constraint for stability and speed. ENCO (38 optimizes independent edge likelihoods without explicit acyclicity constraints, while Nazaret al. introduce a spectral constraint a Stable DCD 44 approach improve numerical stability and scale causal discovery to thousands of variables.
Reinforcement learning (RL) has also been applied to structure learning - e.g., Zhu and Chen 77 use an RL. agent to sequentially build the DAG, avoiding expensive combinatorial searches. Beyond point estimation of single graph, recent variational methods estimate a distribution over plausible DAGs. For instance, BCD-Nets [16] perform Bayesian causal discovery via variational inference, yielding a posterior over DAG structures and quantifying uncertainty. Similarly, Differentiable Causal Discovery with Interventions (DCDI) 12 integrates normalizing flows to model complex causal mechanisms under interventions.
Mediation Analysis. Deep learning has also been applied to mediation analysis to estimate path- specific effects (natural direct and indirect effects) in complex settings. DeepMed (65] is a recent example that uses deep neural networks to flexibly estimate the necessary nuisance functions in the efficient influence function for mediation, achieving semiparametric-optinaa estimation of direct and indirect effects. By cross-fitting DNNs for propensity and outcome models, DeepMed debiases mediator and outcome predictions and attains the efficiency bound without restrictive model assump- tions. Other works leverage deep variable models to handle mediators and hidden confounding.
35

---

CMA-VAE [13] introduced variational autoencoder for causal mediation that accounts for latent con- founders, enabling estimation of natural direct and indirect effects even with unobserved confounders. Building on this, Disentangled Mediation Analysis VAE (DMAVAE) [66] separates latent factors into mediating, confounding, and irrelevant components. This disentanglement allows DMAVAE to consistently estimate direct and indirect effects under weaker assumptions than traditional sequential ignorability, significantly improving mediation effect estimation in the presence of multiple types of latent confounders.
In parallel, deep models have been used for path-specific effect estimation in nigh-dimensional settings like text and faimess. CausalBERT [63] fine-tunes transformer language models to learn low-dimensional text embeddings that retain information needed for causal adjustment. This approach enables estimation of causal effects from text (e.g., the effect of review's content on outcomes) by adjusting for textual confounders within the embedding space. Similarly, algorithmic faimess, researchers have examined path-specific effects via deep nets to detect mediated discrimination. For example, Chiappa [1] studies counterfactaal fairness by removing information about protected attributes along certain causal paths, and Nabi & Shpitser [43) develop methods to identify path- specific effects under mediation and confounding. These efforts underscore a growing theme: integrating deep representation learning with causal mediation analysis to uncover nuanced causal pathways (e.g., via stability analysis or regularization to select stable mediators) in complex, high- dimensional data.
Deep Treatment Effect Estimation. Estimating individualized treatment effects (ITEs) or hetero- geneous causal effects has benefited greatly from deep learning, particularly through representation learning and efficient influence function (EIF) techniques. A landmark work by Johansson et al. [3 introduced learning balanced representations for treatment and control groups to improve ITE gener- alization. This approach, extended by Shalit et al. [53] with integral probability metric regularization, underlies methods like TARNeI/CFR that reduce selection bias by aligning latent representations across treatment groups Subsequent methods added adversarial objectives to enforce balance: Yao et al. [67] propose a SITE approach that preserves local similarity while using adversarial training to minimize covariate distribution differences, and GANITE [68] directly employs GANs to generate counterfactual outcomes for ITE estimation.
Another line of work incorporates semiparametric theory to improve estimation efficiency and robustness. DragonNet (54) a neural architecture that adapts the network training to target the efficient influence function of the average treatment effect. By jointly learning the outcome and propensity predictor with a targeted regularization term, DragonNet yields doubly-robust estimates that empirically outperform earlier architectures. In a related vein, CEVAE (39 uses deep variational inference to handle unobserved confounders, learning a latent variable model that estimates treatment effects under a causal graphical model assumption. Extensions of CEVAE have further improved disentanglement; for instance, TEDVAE [70] factorizes latent confounders into multiple parts to better separate causal effects.
To inject inductive biases about potential outcomes, FlexTENet 17 adaptively learns what represen- tations to share between the treated and control outcome functions, reflecting the assumption that many predictive factors are common between potential outcomes. Beyond standard observational studies, deep models have tackled more complex causal inference tasks: DeepIV (5 uses deep nets for instrumental variable regression to estimate counterfactual predictions in the presence of hidden confounding, and other works address time-series and panel data with recurrent networks for sequential treatment effects (5)
Causal Inference under Structural Uncertainty A distinguishing challenge for methods like SIGMA is causal inference under structural uncertainty, where the true causal graph is unknown or ambiguous. Recent research has begun to explicitly account for this uncertainty using Bayesian deep learning, model averaging, and ensembles of causal models. One approach is to perform Bayesian model averaging (BMA) over DAGs: rather than committing to a single estimated graph. one can weight effect estimates by each graph's 's posterior probability. Cundy et al. [16] implement this via their variational BCD-Nets, which output a distribution over possible DAG structures. By sampling DAGs from the learned posterior, BCD-Nets can compute expectation of causal effects across many graph realizations, inherently providing credible intervals for effects. DECI (2), a deep end-to-end causal inference framework, uses normalizing flows within variational Bayesian scheme.
36

---

DECI simultaneously learns a DAG posterior and estimates causal quantities, essentially performing causal discovery and inference in one unified model.
Furthermore, techniques like stability selection have been adapted to causal structure learning to identify consensus edges under data perturbations, thereby informing model averaging strategies (e.g., selecting only edges that consistently appear across high-scoring DAGs) [44]. Such ideas align with the stability analysis component of SIGMA, which seeks robust mediation pathways by exploring variations in the learned structure. Overall, the literature suggests that accounting for structural uncertainty - through Bayesian posterior averaging. DAG ensembles, or stability- -bseed model selection can substantially improve the reliability of causal effect estimates in complex models. SIGMA builds on these insights by integrating BMA of mediation pathways with EIF-based effect estimation, ensuring that uncertainty in the causal graph is properly propagated into uncertainty in estimated direct and indirect effects.
D Human Phenotype Project Dataset
D.1 Description of HPP Cohort
The Human Phenotype Project (HPP) is deeply phenotyped, longitudinal cohort study conducted in Israel, designed to investigate the interplay between various physiological systems and environmental exposures. The present study analyzed data collected from 6,748 participants aged 40-75 years between January 2019 and December 2022. All participants provided written informed consent, and the study protocol was approved by the institutional review board of the Weizmann Institute of Science (IRB no. 1719-1). Participants underwent comprehensive profiling, including clinical, physiological, behavioral, and multi-omics assessments, which were categorized into body systems [44 one of which focused on sleep. Sleep characteristics were derived from home sleep apnea testing (HSAT) using the WatchPAT 300 device (ZOLL Itamar), with each participant completing up to two monitoring series, each comprising three nights of testing over two-week period. A total of 16,812 nights of valid HSAT recordings from 6,366 individuals (47.8% male, 52.2% female; mean age 52.4 + 7.7 years; mean BMI 26.1 + 4.1 kg/m²) were included in the final analysis.
The cohort primarily consisted of healthy, educated individuals of European (Ashkenazi) Jewish ancestry residing in a relatively homogenous environmental context. Exclusion criteria included severe medical conditions at enrollment. Participants were scheduled for biennial follow-up assess- ments over a 25-y period. This extensive data resource enables high-resolution investigation of associations between sleep and systemic health.
D.2 Body System- Features
We utilize data from 11 body systems, comprising over 500 raw clinical and physiological attributes. To facilitate tractable and clinically meaningful mediation analysis, we derived set of 62 aggregated features by iamaiin-nformmed grouping of related variables across systems. These features are con- structed to capture coherent physiological constructs suitable for pathway-level causal interpretation. Please see Table for the detailed description of our aggregated features.
Sleep Characteristic Sleep characteristics were captured using the WatchPAT 300 home sleep mon- itoring device across up to three nights per participant. A total of 448 features were extracted, compris- ing two main subgroups: 100 sleep test measurements and 348 pulse rate variability (PRV) features. Sleep test features included metrics related to respiratory events (e.g., peripheral apnea-hy -ymnaee index [pAHI]), oxygen desaturation (e.g., mean SpO2 nadir, time below 90% saturation), snoring intensity, sleep stage distribution (light, deep, REM), sleep efficiency, and body position during sleep. PRV features were derived from the peripheral arterial tonometry signal using the NeuroKit2 li- brary (2). spanning time- and frequency-domat indices, entropy, and nonlinear dynamics, computed across REM, NREM, wake, and full-night segments.
Metabolic Pathways Metabolic pathway features were derived from metagenomic profiling of the gut microbiome using the HUMAnN3 4 functional annotation pipeline. Relative abundances of microbial metabolic pathways were quantified based on high-throughput sequencing data, followed by taxonomic and functional mapping against curated reference database. Only pathways present in
37

---

Table 3: Summary of physiological systems, features and medical explanations
System
Explanation
bodyy for pencninge rato
high valtrr sm licked to abouty and metabelic risk
Body Compesition
Refects android gyand disribution, mic
sugests
candiousealur
aound organe: resecissed Measures linc music mass a key andicato for scroopenia and Etaiyy
Quannittes
informmation.
apperdicola loan mass
severity ssucing scvaiity
Iodictes obetnattvve spaca (OSA), levels suggest Serute OSA. Key indicator of OSA severity bacd in staning loadses and frequency.
OSN Reiazed Phenorypes
enuring pilo
Prepertion of snoring time dering tost sleep time; linked to CSA inck.
slorp
disoddrs. Rato actual slep time to time in belt reffect steep mecnnenance quality. Meusurzs Eragneenttion lower continitty indicttes more fisgrnnted sleep.
skopp laggh:
secscbeall
Sleep Stuuture Quality
sleep effeiency
sloep continuity
slepp
rotteaise dhap siagaa
clags nelleats slotp quality
night mece HH
Anarage sightime hane mit: reflass astonamic balance.
Sleep HRV
fex Index
HRV indet during steep: lower HRY licked stress and metabolic
stage vatiation fev wake deta
vuriabllay scross steep stgee, reflects sleep structure quility.
HRV diffrrence berwees sleep and wake, indiores recovery capacity.
Nusturall Hypusii Bunden
paygen desauuation uegge hanken
Measurss oeyges severity
slaeps indicaics hypoccia level. Comelative oxygen henden from all eonne aeflocs physiiaggical liod.
Iying disaolic prosem::
Distriiz blod
mesomrd while
Blood lying
lying paise
Palse messured birg.
Iying aystoli pressure
Sysolic booo presure nessured while Iying.
stting diastlic pressure
Diatolic blood presuuee seasured while
Blood Pressure Sitting
vising pulse mac sittig yyteeli presoue
Pole nie meaeum
diting
Sycolie food presurr meesured while satingg
stnddnn Imin diasatiic prosure starding Lais pulse rue standng lmis ealic pezsare susding Jain dissallc presuse steningg Imin pube stndiing Imin gyoole pressure
Diedelic presare mesound ator sonndnng for
Pulie nue nessured after standing for I ninue
Blood Provor: Standing
Diztoli blood presuuee
Syenlic bood prosure mexsered alter standing foe I minate.
sunling for
Pole náe meaemm sfee sonding minnto.
Syeolie food prestrr meesured ofher standing for 3 minates.
reiig sycollc prosore resting distclic presuze resing pule
Syeolic possrre
Blood Provamm Rexting
Distelic blood pressure rest.
Pulic nue res.
rrlitstatic SBP drup rethvetatic hRP change Iiin octhontrie SUP Imin
Drop systelic blood prosuee minute afler standing.
10
Blood Pressure Orthostane
Chongs in diseolic hirdd proczr: I minst: aheer standing
systelic bloed presure minates standing.
change
Change disonlic blond prosere
after stunding Incease pubce I maute afte standing,
11
Blood Presours Puke
pulie rute increse Lais pale
pulce miatce slir standing.
ABl min
Askle-bnchial indo; valuss -019 sugerst peripheral antery discase. Poto velocity indicatss anenial sithnee and cartinasoute rak
12
Viacular Healds
PWV mean
mux
Maximon systeli blood pressure indicaee peck uterial load.
int left
Intimo nedi diickness (IMT) of left caroid attry.
int
Intino dicckees (ITT) right
ini windor keft int windon walls right
Window with of mosumement on the
Window with of IMT moseeeneni on the right side.
15
Caroid Ulioonnd INT
int left
IMU value en the left ede..
right
value ou right side:
mas cint
Average cardid IMT acnos mocsred segions.
cint
Masinum czetid IMT otscrved.
shnormd
Indicais preceees
cacerling chinccl theodalds.
14
Caroid Ulsawwnnd Plague
plaque presence by tit plopee prosence thickness
Indicnos
Indiates presence plaque based on IMT Ir. criteria
plaque haod IMT micknes threchodd
buyasz produces abandance popimaic prodacen abondlance
Abundaace butynue produciu bacisii, lisked to asti inflamnstrry effects and gut banier ingarit. repporraic-peesccing battrna invlval in appeiic reguaiion and mutsholiom
peciotic
ahundasee
Probiatic feraaciing
produce lactic
ofter beneficial metabciies.
viamin prodhcere ahurdbnne
Basterio eapobic of syrhhoiingg vioeins each viseins
15
Gut Micnoiocte Abundance
Mikemasiia abandance
Abundance of "Akttmmasiia, associated with netabolic healt.
produces shandance
Busteria produciu sacgellsaccaaaiie (LPS), may promcte infirmattonn.
retucess abuudiace
Sulfite reducing basteri potertall hannfll
TMA produces abundanee
Busteria producing TMA gprecane of TTAA linked to candoosacaee dissac.
pestedyyic hocteiia shondanee
Roctori that degnake provine: inclhed produccion of lock mussholiza.
shannon divesity index
Shaznon index, messues dicunss everness of gut nicoobes a key makker gat bealt).
at least 5% of participants were retained for analysis, and abundance values were orgio-transformed to stabilize variance. These pathways represent a functional summary of microbial activity, covering key domains such as amino acid biosynthesis, vitamin and cofactor metabolism, carbohydrate degradation, and chort-chain fatty acid production. They were used to probe host-microbiome interactions relevant to systemic metabolic health and potential associations with sleep-disordered breathing phenotypes.
Cardiovascular System Cardiovascular features encompassed a comprehensive set of hemody- namic and structural measurements obtained through non-invasive techniques. Peripheral and central
38

---

blood pressures were measured in sitting, supine, and standing positions. Arterial stiffness was evalu- ated using pulse wave velocity (PWV) recorded with the Falcon device (Viasonix), and endothelial function was indirectly assessed through carotid intima-media thickness using high-resolution ultra- sonography (SuperSonic Aixplorer MACH 30). Retinal microvaseular parameters were extracted via fundus imaging (iCare DRSplus) and processed with the Muth (74] Python package. Cardiac electrical activity was characterized by 12-lead resting ECG (PC-ECG 1200, NORAV), from which standard and derived indices (e.g., heart rate, QT intervals) were calculated. The cardiovascular dataset thus integrates structural, functional, and clectrophysiological dimensions of cardiovascular health, enabling multiscale association testing with sleep traits.
Inflammation The inflammation system features were derived from standard clinical laboratory tests obtained from participants' health maintenance organizations (HMOs) as part of their routine medical care. Participants were encouraged to upload results from various tests, including those indicative of systemic inflammation. Key biomarkers relevant to the inflammation system included C-reactive protein (CRP), white blood cell (WBC) count, and differential counts such as neutrophils, lymphocytes, monocytes, eosinophils, and basophils. These markers provide insights into both acute and chronic inflammatory states. CRP, an acute-phast protein synthesized by the liver, serves as a sensitive indicator of systemic inflammation and is commonly elevated in response to infection or tissue injury. The WBC count and its differentials offer information on the immune system's cellular components, with variations potentially reflecting underlying inflammatory immune responses. All laboratory values were standardized and stored in a structured format to facilitate downstream analyses.
Baseline Characteristics (BMI) Body mass index (BMI) was calculated as weight in kilograms divided by height in meters squared (kg/m²) and was treated baseline covariate alongside age and sex. Though BMI is not classified as a standalone body system in this analysis, serves as a critical covariate due to its known influence both sleep-disordered breathing and cardiometabolic risk. It was excluded from the body composition system to avoid redundancy in statistical modeling.
E
Synthetic Datasets Details
E.1
Generation Process
We construct synthetic datasets using a modular simulation framework that mimics realistic causal sys- tems with configurable complexity. The generation process is parameterized to control graph structure, variable types, functional complexity. effect magnitudes, missing data patterns, and more, making the synthetic data suitable for benchmarking causal discovery and mediation analysis algorithms. The entire dataset generation procedure is outlined in Algorithm 2
Causal Graph Construction Each dataset is generated from DAG G over p variables and contains n observational samples. The graph is sampled from an Erdös-Rén (ER) model with a specified expected average in-degree followed by random topological sort to ensure acyclicity. The DAG defines the full set of causal relationships among the p variables, with no predefined roles assigned to individual variables. To support the construction of realistic and diverse causal structures, we explicitly control two key properties of the graph: the chain path length E, which determines the depth of the longest directed chain in the DAG, and the number of structurally parallel paths k, which specifies how many disjoint or partially overlapping causal routes exist across the graph. This parameterization allows us to generate graphs that reflect varying levels of structural complexity, enabling robust evaluation of causal discovery and inference methods under different topological regimes.
Structural Equation Modeling Let v = (u1,.... Ppp) be a collection of random variables. Each node is randomly assigned a variable type according to a multinomial distribution over three categories: continuous with probability n., binary with ma and categorical (with up to Leat levels) with probability Teat where we + my + Tent = 1. Once the types are assigned, we generate variables according to structural equation modeling (SEM) [46] framework based on a given DAG. Specifically, each variable in is assigned structural equation of the form:
1/ F: (Vpa(1.))) 21),
(43)
39

---

where is an exogenous noise term independent of all other variables, pa(i; G) denotes the set of parent nodes of the graph G, and F; is the deterministic function governing the influence of parent variables on u4 In this work, we focus special class of SEMs known as additive noise models (ANMs), where the structural equations take the simplified form:
U2 fa (M)()))))) 2is in vector form: v == ==f(v +
(44)
The functional form f; may be linear nonlinear, including polynomial (c.g., v²), sinusoidal (e.g.. sin) or exponential (e.g., exp(-lul)) transformations. To control the proportion of nonlinear mechanisms, we define nonlinearity ratio parameter Pnonlin E [0., 1], which specifies the expected fraction nodes governed by nonlinear functions. For nodes not selected as nonlinear, fr defaults to linear combination of parent variables. This stochastic selection enables systematic benchmarking across datasets with different levels of structural complexity. Pairwise interactions between parent variables also be incorporated increase structural complexity. The exogenous noise z; is sampled from a mixture of distributions Gaussian, exponential, or Student-t- wwtt user-defined mixture weights simulate heteroskedasticity and non-Gaussian variability, For discrete variables, logistic or softmax functions are used to generate probabilistic outcomes conditioned on parent variables.
Causal Effeet Scaling To enable controlled experimentation with mediation analysis, the generator includes an effect-scaling mechanism that adjusts the magnitude of causal effects along specific paths. Structural coefficients on the direct path from treatment to outcome, and on the indirect paths through mediators, are scaled such that the resulting absolute values of the natural direct effect (NDE) and natural indirect effect (NIE) fall within user-specified intervals, denoted by ('DDE 'NDE and [rmin max), respectively. This ensures that generated datasets reflect wide spectrum of causal regimes, from dominantly direct to heavily mediated effects.
mAX]
Missingness Injection and Ground Truth Estimation To further enhance realism, the framework optionally injects missingness into the dataset under both missing completely at random (MCAR) and missing at random (MAR) mechanisms. The overall missingness rate Tmiss determines the proportion of missing entries. In addition, the framework allows configuration of the proportion of variables and the conditional dependencies under which missingness occurs. Ground-truth causal effects, including total effect (TE), NDE, and NIE, are estimated via Monte Carlo simulation using large number of samples (e.g., 10,000 or more) by evaluating counterfactual outcomes under interventions on treatment and mediator variables.
Algorithm 2 Synthetic Dataset Generation Procedure
1: Require: Number of variables p, samples n, average in-degree k
2:
Target effect ranges TNDE, TNIE
3:
Number of parallel paths in DAG:
4:
Length of the chain path in DAG:
5:
Variable type proportions: Tes Tb. Teat
6:
Nonlinearity proportion Pronlin
7:
Noise distribution mix Duuise missingness rate Tiiss
8: Ensure: Synthetic dataset D with DAG G, SEM specs, and ground-trutt effects
9: Generate DAG G using Erdös-Rényí model with p nodes and average in-degree k.
10: Verify that G contains & parallel directed paths of length <
11: Assign variable types each 1y E
12: Continuous with Te; Binary with Th: Categorical (up to Lea)) with Teat-
for each node E v do
14: Define: + fr (Vpu(i:G)) + ž;i
Choose fa as linear or nonlinear (e.g., poly2, sin, exp, interaction) with probability Pronlin:
16: Sample noise 2; ~ Duoise
end for
18: Generate samples by topological traversal of G
19: Inject missing values using Timiss and conditional rules (e.g., MAR/MCAR).
20: Estimate ground-truth effects (TE, NDE, NIE) via Monte Carlo simulation.
21: return /
40

---

E.2 Mediation Pathway Structures
To evaluate mediation analysis methods under diverse structural assumptions, we incorporate range of causal graphs with varying mediation configurations into the synthetic data generator. Figure several representative DAG structures that are automatically instantiated within our framework. These examples cover key topologies relevant for causal mediation, including both simple and complex pathways, as well as potential nonlinear effects.
X1
X1
Mi
M
Y
Y
X2
X2
M2
(a) Two node
(b) Single Mediator with Confounding
(c) Dual Independent Mediators
M1
M1
X
M3
M2
M2
(d) Nonlinear Dual Mediators
(e) Hybrid Chain-and Mhaan-and-Paraallll Mediation
Figure 4: DAG graphs. Unless otherwise stated, we take X as the treatment, Y as the outcome, and M as mediator variable.
1. (a) Direct-only structure ("Two-Node Path"): This minimal structure includes only direct effect from treatment to outcome Y, with no mediators. It serves as a baseline for evaluating spurious mediation detection.
2. (b) Single mediator with confounding inputs ("Single Mediator with Confounding "): A single mediator lies between X and but additional variables (e.g., X2) influence the mediator, simulating confounded or high-dimensional mediator function.
3. (c) Parallel mediation ("Dual Independent Mediators"). Two mediators Mi and M2 are both affected by X and influence independently, forming parallel mediation paths.
4. (d) Nonlinear parallel mediation ("Nonlinear Dual Mediators")) Similar to (c), but the structural functions between nodes involve nonlinearities such as interactions or sine/exponential transforms.
5. e) Mixed mediation ("Hybrid Chain-and and-Parallel Mediation more complex topology where two mediators M, and M2 form parallel channels, and a third node M3 forms a chain extending toward Y. This structure tests the ability of algorithms disentangle indirect effects across multiple pathways.
E.3 Dataset Benchmarks
Table summarizes the configurations used to generate representative datasets in our experiments. These settings allow us to benchmark the performance of causal inference methods across regimes with varying dimensionality, mediation complexity, and data quality. Specifically, we hold constant the overall variable type proportions (л, 0.6, = 0.3, Tear 0.1). category levels (Leat 5), missingness (rmiss 0. 1), average in-degree (k = 3) and the target standard deviation ranges of the causal effects (TNDE [0.2, 0.4), TNIE E (0.15, 0.3]. This design allows for controlled benchmarking across varying causal and statistical challenges.
To further validate the quality of mediation structures embedded in the synthetic datasets, we visualize the variable relationships of sampled mediation pathways. Specifically, we examine the empirical dependencies among treatment, mediator, and outcome variables to see whether the intended causal effects are faithfully reflected in the generated data. We analyze four types of associations: (i) Treatment + Mediator, (ii) Mediator Outcome, (iii) Treatment - Outcome, and (iv) Mediator
41

---

Dataset ID
P
Proalia
Configuration Focus
LowDim-L
20
6000
1
Baseline (Low-dim, linear, simple graph)
LowDim-N
20
6000
0.5
1
Baseline (Low-dim, nonlinear, simple graph) Complex path (Low-dim, parallel, nonlinear) Complex path (Low-dim, long chain, nonlinear) Mid-dim effect (long chain. nonlinear, simple graph) Mid-dim effect (long chain, nonlinear, complex graph) High-dim effect (monlinear, simple grapgh)
LowDim-P
20
6000
0.5
2
LowDim-D
20
6000
0.5
MidDim-S
100
6000
0.5
MidDim-D
100
6000
0.5
HigDim-S
200
HigDim-D
200
6000
High-dim effect. (long chain, nonlinear, upper bound)
MidDim-P
100
6000
0.5
2
Model misspecification test (nonlinear, parallel)
MidDim-C
50
6000
0.5
Mid-dim effect (nonlinear, long chain, smoothed comparison)
Table 4: Parameter configurations for synthetic datasets used in our experiments.
+ Outcome under different treatment. These visualizations provide intuitive verification that the synthetic data preserves the expected mediation mechanisms and highlight the nonlinear or conditional dependencies that may arise under the configured data-generating process.
LowDim-L: LowDim-L serves as the baseline configuration in the benchmark suite, featuring low dimensionality (p = 20), fully linear structural specification (Phonlia = 0.0), and a simple graph topology with single short causal path (x I, l 3). The graph structure is minimal and interpretable, offering a clean testbed for evaluating algorithmic correctness in effect estimation under idealized, statistically well-behaved conditions. The causal graph consists variable V8, continuous mediator V2, categorical mediator V18, and a binary outcome V13. This setup enables clean interpretability of causal effects and serves as reference for assessing performance under idealized conditions. LowDim-L provides controlled environment for validating correctness of pathway identification and effect estimation algorithms under structurally simple and statistically well-behaved conditions.
binary treatment
vB.- 12
VB V1a
Medistor V2 Dutcome V15
Mediate's V18 Outcome V13
V18-0.0
Vid=1.1
V12-2.0
V18-3.0
g
13 LB
2012
V2 V9 V15
V2 V13 uB
V18
V/3 v
VB+0.0
VB-0.0
VB=10
VS-10
Mediation Pathway
Treatment: V8:
Mediator: V2, V18;
Outcome: V13;
Figure 5: Visualization of variable relationships sampled mediation pathway in LowDim-L. The pathway comprises binary treatment variable V8, a continuous mediator V2, categorical mediator V18, and binary outcome V13.
42

---

LowDim-N: LowDim-N mirrors the low-dimensional structure of LowDin-L (p 20, n 6000).
but introduces nonlinear functional mechanisms (Pnonlin 0.5) while maintaining structurally simple graph with one short directed path (x 1. {= = 3). This setting incorporates nonlinear mappings within the graph's 's structural equations, allowing us evaluate whether causal inference methods remain robust to nonlinearity. The variable types are consistent: treatment V8 is binary, mediators V2 and V18 are continuous and categorical respectively, and outcome V13 is binary. In contrast to LowDim-L, the functional dependencies here include non-additive effects and interactions, enabling the examination empirical recoverability under nonlinear data-generating processes. In particular, the direction of the -utiite relationship is reversed: V2 now has positive rather than negative effect on V13, consistent with the design of its generating function. Conditional plots further highlight interactions modulated by treatment levels, providing visual confirmation of nonlinear mediation behavior.
Medisty V2 V15
Mediaeer ViB Outoome
V18-2.0
V1d+3.0
Probobiity (-1)
VB 122 V13 VB
A-13
V2 V8 - VIS
V18
VTB VI3 vB
VB-0.0
Va-00
VS-10
VB-10
Mediation Pathway
XM
Treatment:
Mediator: V2, V18;
byVV10 1/000-0 Probate Ryy 10-8-11
Outcome: V13;
V15
Figure 6: Visualization of variable relationships of sampled mediation pathway in LowDim-N. The pathway comprises binary treatment variable V8, a continuous mediator V2, a categorical mediator V18, and binary outcome V13.
43

---

LowDim-P: LowDin-P is designed to introduce structural complexity through multiple causal pathways, while maintaining low dimensionality (p 20, n 6000) and moderate nonlinearity (Pnonlin = 0.5). The graph contains 2 parallel paths of moderate depth (E = 3), representing con- current channels influence. The variable configuration includes binary treatment V8, continuous mediator V2, categorical mediator V18, binary mediator V7, and binary outcome V13. This setup is designed to evaluate the interaction and combined influence multiple concurrent mediators under nonlinear conditions.
VidH0.3
Vb+0.0
ProtasilyV/1311)
V1d+1
VBer
V18-2.0
VO
Via Via
V7 V13
VBI
VI
V/16 V1ai Vi
V/H0.
Pr (jsee ())
ViHH0
V/+10
Protonilt(Vaarrpeee ()
Vile 1.0
V15
V13
12
V15
We+10
Vi+10
Mediation Pathway
Treatment: V8;
Mediator: V2, V7, V18;
Outcome: V13;
Figure 7: Visualization of variable relationships a sampled mediation pathway in LowDim-P.
The pathway comprises binary treatment variable V8, continuous mediator V2, categorical mediator V18, a binary mediator V7 and a binary outcome V13.
44

---

LowDim-D: LowDim- D represents a more structurally complex synthetic setting with p 20, 6000, moderate nonlinearity (Pmonlin 0.5), and a long-chain causal structure characterized by single deep path (N e = 6). It includes binary treatment variable V8, continuous mediators V22 V16, and V19, a categorical mediator V18, binary mediator V11, and a binary outcome V13. This configuration captures extended sequential dependencies, where indirect effects propagate through multiple intermediate variables. It provides a testbed for assessing robustness in effect estimation under deeper, multi-step causal pathways.
18
VB.-
V1G+BL
V10+1
V18-2.0
V1d+30
VS V18 VT3
1B:03
VS
VS
P5 - VIT
V2 - V13
Vid - VIS
VD+0
VB-1.0
Probade R61V155-11
Probal Ry1V13H11)
12
V18 V2 - V131 VB
v15 V13
VY V1
V95 - Vi31 V9
ViT+000
Vit=10
Protab 1y
Probul 9B110-1-0y111-1----0-000
Photesh 1y
Photwo spttt
spotes
Ve
V18 V13
V'9 V13) VS
V11 VB
V0+0.0
VB+0.0
v8-0.0 VB-1
VG-1133
V8=1.0
Pestoedd by Vieco
Probote Ry V8:100
Figure 8: Visualization of variable relationships sampled mediation pathway in LowDim-D. The pathway comprises binary treatment variable V8, continuous mediators V2, V and V 19, categorical mediator V /188, a binary mediator V V11 and a binary outcome V 13.
45

---

MidDim-S: MidDin-S increases both dimensionality and structural complexity, with = 100, 6000, and a moderate level of nonlinearity (Proalin 0.5), while preserving simple causal topology with single short directed path (k = 1 3). The treatment variable V35 is binary, mediators include a categorical variable V2, a binary variable V13, and continuous variable V14, with V17 serving as the continuous outcome. This setup is designed to assess algorithm performance in mid-dimensional regimes where the causal structure remains sparse and shallow, isolating the effect of dimensional scaling under otherwise controlled conditions.
V25 12
V35 V13
V35 - Vid
V2-V17
V2-04
V25-0.0
12-10
Vitbe
V2-2.
Pora 100
V2+30
V25
VIS
V26
V2
V13
V35+0
Mediation Pathway
Treatment: V35;
Mediator: V2,
VSS-00 O
Outcome: V17;
Vitel
VS6
Figure 9: Visualization of variable relationships of sampled mediation pathway in MidDim-S. The pathway comprises binary treatment variable V35, continuous mediator 14, categorical mediator V2, a binary mediator V13 and continuous outcome 17.
46

---

MidDim-D: D: MidDin-D extends the structural complexity of the mid -iddesmnnni setting with p 100, n 6000, a nonlinearity ratio of Proalin 0.5, and a long-chain causal configuration characterized by single extended path (N = 1, l 6). The DAG includes multiple layers of intermediate variables connected through deep sequential structure, challenging algorithms to recover indirect effects that accumulate across heterogeneous transformations. The treatment variable V35 is binary, and the mediation pathway involves heterogeneous mediators: categorical variable V2, binary variables /3, V15 and V89, and continuous variables V 14 and V40, with continuous outcome V17. This dataset well-suited evaluating a method's ability to capture layered, nonlinear, and high-dimensional mediation dynamics.
V25 32
V25 140
V89
V15 - VIS
V2-00
V35-0
V05+0
V2-10
V95-1
VSS-1.0
V2-2.
V2+3.
VI5
VIB
V20
VIS
V2
V40
V35-0.0
V85-0.0
V171
VS5
V17
Mediation Pathway
v35e11
Treatment: V35;
Mediator: V2, V13, V14, V15, V89;
V25HO.0
Outcome: V17:
V35+110
RHQ
Figure 10: Visualization of variable relationships sampled mediation pathway MidDim-D The pathway comprises a binary treatment variable V/3,, continuous mediators V14 and V40, categorical mediator V2, binary mediators V13, V and V889, and a continuous outcome 17.
47

---

HigDim-S: HighDin-S explores high-ding digh--iiensional setting with p 200, n n 6000, moderate nonlinearity (Pnontin 0.5), and sparse causal structure characterized by single short path é t= 3). The causal pathway consists of binary treatment variable V71,, continuous mediator V35, a binary mediator V28, and a binary outcome V83. Despite the limited path depth, the nigh-dimensional background introduces substantial noise and potential confounding, making this configuration well-suited for testing robustness in pathway recovery under low-signal regimes. is designed to assess performance when precise signal extraction is required in the presence of many irrelevant or weakly related variables.
V20H0
V28:110
V80-1
-2
V11
VES
V71:000
Mediation Pathway
Treatment: V71;
Mediator: V28, V35;
Probabilty 11
Outcome: V83;
Figure 11: Visualization of variable relationships of a sampled mediation pathway in HigDim-S. The pathway comprises binary treatment variable V71, continuous mediator V35, binary mediator V28, and a binary outcome V88.
48

---

HigDim-D: HighDim-D represents the most structurally complex setting in the benchmark, with high dimensionality (p 200), moderate nonlinearity (Pronlin 0.5), and deep causal pathway characterized by single extended chain (xs = 1, 6). The underlying DAG contains multiple layers of variables with diverse data types, connected through long directed path. The causal structure includes a binary treatment V71, continuous mediators V55 and V176, categorical mediators V5 and V57, binary mediators V28, V15, and V89, and binary outcome V83. This design models rich causal hierarchies and nonlinear dependencies, mimicking the structure of intricate biological or behavioral systems.
V170
VS5
V21
VS-0.
VSR-0.
Viel
VS:2.
VS7-2.0
V7
Vtr=1.0
P1200-11
Probabilty
1y1V59:1/77110
Probubely448:11171-000 Phibssetyllt
VS7
V178
V25
171:1100
V57
V22
V23
Figure 12: Visualization of variable relationships of sampled mediation pathway in HigDim-D. The pathway comprises a binary treatment variable V71, continuous mediators V35 and V176, categorical mediators V5 and 577, binary mediators V28, V15 and V89, and a binary outcome V88.
49

---

MidDim-P:
MidDin-P is designed to evaluate robustness under potential model misspecification.
featuring moderate dimensionality (p 100), nonlinearity (Pronlin 0.5), and parallel causal structure with overlapping pathways (K = 2, f 3). The causal graph includes binary treatment V35, continuous mediators V14, V68, and V72, categorical mediator V2, and a continuous outcome V17. Unlike prior configurations emphasizing either deep chains or sparse connections, this setting introduces partially redundant, nonlinear paths that may obscure structural signals and challenge estimation accuracy when assumptions are violated.
V05
VAS
V2-2.
R000 1A
V2-0
60:33
69-12
V12
12- VIT
VIT
-VT
V12
V2
vea
V14
VITI
V12 V36 V11
V53
V14
Modiation Pathway
Treatment:
Mediator: V2, V14, V68, V72:
Outcome: V17;
Figure 13: Visualization of variable relationships of a sampled mediation pathway in MidDin-P. The pathway comprises binary treatment variable V/3,, continuous mediators V14, V68 and V72, a categorical mediator V2, and a continuous outcome V17.
50

---

MidDim-C: MidDim-C explores mid-dimensional, high-complexity causal dynamics under nonlin- ear transformations, with p 50, n 6000, and deep chain structure defined by a single extended path (k = 1, l 6). The treatment variable V35 is binary, the mediators-Vi V14, V33, V42, and V48- are all continuous, and the binary outcome V45. This configuration emphasizes smooth, continuous interactions across multiple intermediate nodes, where causal influence accumulates gradually along the chain. It is designed to test methods' ability to detect and quantify multi-step effects under structurally rich yet wwwpaarllellsm conditions.
VSS V35
V05 56
182
V35
VB5 148 - y45
V85
V95
V25 - V14
VIS Vi5
VS - ViB
VAC
VAC
Probuc Probuo 10
Probets Ry1145411
Probub-iy1145-11
VI5
V03
114 - V45
Vic - V45
V48 V451 V/55
V20 V451 VIS
1465
0
-01 -54
-06
18
V48 129
V42 Wid V25
V/54 V45 V99
V29 V45
V35-0
v35-1
Pistetiityyttt 10225-0000000 ())
Proteds 12220110 (1225-0
Reoobillyy45-11 10-36-1-111000
Proceti 21-44451133-5111111
ProbadilyV45 1113--1----
Figure 14: Visualization of variable relationships a sampled mediation pathway in MidDim-C. The pathway comprises a binary treatment variable V35, continuous mediators V6, V14, V33, V42 and V 48, and binary outcome V /5..
51

---

F
Experiment
E1
Experiment Setup
In the structure discovery phase, we configure the Flow-SEM model with a two-layer MLP (hid- don_dim multiplier=2)) as the conditional distribution model. The model is trained using the Adam optimizer with a learning rate of 0,001, an initial DAG constraint penalty coefficient À = 0.1, which increases by a factor of 10 when the acyclicity constraint is violated (h threshold=10 up to maximum of 10 Sparsity is controlled through L1 regularization (a = 0.01), with maximum gra- dient norm limit of 1.0. The training iterates for 1000 epochs, with early stopping if no improvement occurs for 100 consecutive epochs.
In the posterior sampling phase, we extract 1000 DAG samples from the learned structure. During sampling, each node retains only the top 3 highest probability incoming edges to control sparseness, and sigmoid function converts the weight matrix to edge existence probabilities. Graph structures are randomly sampled according to Bernoulli distribution, validated for acyclicity through topological sorting, and cyclic graphs are discarded.
For path identification, we set the frequency threshold to 5%, meaning paths appearing with frequency above this threshold in the DAG sample set identified stable paths. The maximum path length is limited to 10, and multi-processor parallel computing is employed to accelerate the identification process.
The mechanism modeling phase implements 5-fold cross-validation. Each model uses uniform parameters: 50 iterations, batch size of 128, learning rate of 10-4, and hidden layer dimension of 128. For continuous treatment variables, intervention levels are set at the mean (4) and the mean plus one standard deviation (μ + o), while binary variables use intervention levels {0.1).
Effect estimation employs a mixed approach: paths of length 3 use Efficient Influence Function estimation, while longer paths use plugin estimation. Multiple DAG effects are aggregated through Bayesian Model Averaging, considering both within-DAG estimation uncertainty (average variance) and between-DAG structural uncertainty (effect variation), constructing 95% confidence intervals.
We test SIGMA 10 synthetic datasets with feature dimensions p E (20, 50, 100, 200), nonlinearity degree E {0, 0.5), containing parallel (K E {2, 3, 6}) and chain (f {1,2)) structures.
F.2 Computation of Ground Truth Mediation Effects
To validate SIGMA's estimation accuracy, we compute the true NDE and NIE values for the two mediation pathways through counterfactual simulation. For a mediation pathway X - M + Y,
we define three counterfactual data sets: Dref [(X; Freff Ma( (1 Me( })}}
D Deomp {(X; = Tcomps Mi(Deomp)) Yi(Feomp: Mc(Feomp)) H= and Def Frompy Ma ) Y.Ccompp Mi(FFref ))) In 2= where M1xx) denotes the potential mediator value when treatment X; = x, and V,(x,, m. represents the potential outcome when X, = I and M, = m.
{(x;
Based on these counterfactual sets, the true mediation effects are calculated as: NDE EE[Y (Iref; (Xyer)) and NIE = E[Y (Fcomps M (Icomp)) E[Y(Acoopp: M (Iref))) In Monte Carlo simulation, these expectations are estimated via sam- 1 A (x(coomp M4 (Zve:)) - - Yi(xre,, ((rer))) and NIE 1 ET Ye(Coopp: M, (Tcomp) E Y( (arrer))
BI (Tcomp: M (Xres)))
ple means: NDE
We compute these effects using structural equation model specifications with 5000 simulation samples, 115-+V19-V1666 at Øre/ = 0.012 and Zcomp = 1.253, and
setting intervention levels for pathway for pathway (V12 (112+V10~+V6) at Tref 0.243 and Zeomp 0.926. To generate Def, we extract M; (Zref) values from Dref. then set X, = I comp and fix M; = M3 (x ). employing topological sorting ensure proper variable generation order in the DAG.
52

---

F33
Causal Discovery Evaluation
Table 5: Causal discovery on benchmark datasets.
Method
Metrios
Datavet
Datastt 2
Datrect
Datnet 4
Datuset
Datrset 6
Daiset 7
Dotestt 8
Dainset 9
Dataset 10
Adjacency Fl
0.5217
0.4427
0.3498
0.3475
0.5333
Orientation
0.2526
0.2246
B:2564
0.1143
0.2369
0.2922
0.1771
0.3605
0.2823
0.3551
Cunsal Accuracy
0.577
0.5735
0.5849
0.5286
0.6168
0.6243
0.5883
0.5806
0.6247
0.6478
Adjacency F1
0.3728
0.4061
0.471
0.5727
0.4018
0.441
0.32
0.3194
0.4639
0.4829
Noteirs ti
Orientation Fl
0.2158
0.2065
02213
0.1056
0.2636
0.2695
0.1601
0. 0.1488
0.2576
0.3261
Cansal Accuracy
0,5291
0,5301
0,5407
0.4828
0.5745
0.5692
0.5388
0.5244
0.577
0,594
Adjacney Fl
0.3922
0.416
0.5075
0.3935
0/4166
0.416
0.3296
0.3283
0.4186
0.5225
Grandag 3V
Otienation
0.2247
0.2172
0.2485
0.1119
0.2738
0.2359
0.1701
0.1571
022633
0.3364
Cuasal Accuracy
0.5485
0,5056
0.5915
0.598
0.6026
Fl
0377
04168
0.4897
0.364
0.4171
0.413
0.3278
0.3276
0,4102
0.4967
ICANLIGAM
Orientation
0.2432
0,1062
0.2676
0.1641
0,1514
0.2628
0.3323
Cansal Accuracy
0.5476
0.5571
0.5401
0.4257
0.5792
0.583
0.5456
0.5417
0.5879
0.5998
Adjueency Fl
0.431
0.4453
0.5188
0,4019
0/4512
0.4477
0.3326
0
0.4869
0,5077
DeC 22
Otiemation Pl
0,2282
0.2256
0.25%
0.1141
0.2858
0.2991
0.1298
111551
02255
0.355
Guaaa Accuracy
0.5833
0.5661
0.623
0.5366
0.6129
0.5955
0. 0.5872
0.5649
0.6776
0.6625
Adjacency Fl
0,5263
0.4969
0,5602
0.4512
0,4939
0.4998
0,3998
0L3741
0.4754
0.5918
SIGMA (Ours)
Orientation Fl
0,4474
0.3601
IL2784
0.1266
1,3185
0.3163
0.1925
IL1782
0.3162
0.386
Cansal Accurary
0.7695
0.643
0,6719
0,5712
0,6808
0.6747
0.6437
0.6561
0.6914
0.7259
F4 Real Data HPP Validation
Table 6: Estimated NDE, NIE and TE with the corresponding 95% confidence intervals (CI) across all the aub-pathways within the validated mediation pathways shown in Figure [5]b)) Estimates are obtained using the proposed SIGMA framework with structural uncertainty propagation.
TB
NIE
Pathways
Estimate
CI Lower
CI Upper
Fstimate
CI I Lower
CI Upper
Estimate
CI Lower
CI Upper
4-+ 1 1
0,1880
0.1292
0.2432
0.1422
0,0902
0.1995
0.0458
0.0112
0.0823
4-11+12 6
0.1880
0.1306
0.2111
0.1636
0.0818
0,2006
0.0243
0.0036
0.0597
4-11-2-6
0,0462
0.0072
0.0737
-0.0776
-0.0614
0,0502
0.0538
0.0127
0.0921
4-1--2-11--
0.0554
.000566
-0.0310
0.0244
0.0006
0.0604
0.0310
.0.0034
0.0790
4-11-134-6
0.0554
-0.0255
-0.0313
0.0391
+00077
0.0758
0.0163
+0.0149
0.0580
3-1-122-6
0.0965
0.0638
0.1384
0,0821
-0.0024
0.1100
0.0144
.0,0075
0.0543
3-1-2-+12-4
0.0965
0.0472
0,1198
0.0192
-0.0033
0.0575
0.0773
0.0502
0.1023
2++1+6
0,1584
0.0853
0.2267
-0.0252
-0.0884
0.0522
0.1836
0.0777
0,2409
2-1-12-46
0.0323
-0.0580
0.0811
0,1314
0,0816
0.1402
-0.0991
-0.1346
-0.0124
2-1-2-
-0.0397
-0.0653
-0.0036
-0,0462
-0.0809
-0.0093
0.0065
0.0132
0.0502
2-1-2-6
0.1584
0.0834
0.1839
0.1410
0,0023
0.2262
0.0174
0.0000
0.0408
Table Estimated NDE, NIE and TE with the corresponding 95% confidence intervals (CI) across
all the sub-pathways within the validated mediation pathways shown in Figure b).. Estimates are obtained using classical regression-
mediation analysis without structural uncertainty modeling.
fe
NDE
NIE
Pathways
Estimate
CI Lower
CI Upper
Estimate
CI Lower
CI Upper
Estimate
CI Lower
CI Upper
4+ 12 + 6
0.1880
0.1624
0.2133
0.1875
0,1618
0.2126
0.0005
-0.0063
0.0017
4- 12
0,1880
0.1623
0,2124
0.1684
0.1414
0.1950
0.0196
0.0115
0.0285
4-2 12 + 6
-0.0554
-0.0871
-0.0302
-0.0502
-0.0767
-0.0249
-0.0052
-0.097
-0.0035
4-42 12 13
0.0554
0.0800
-00073
0.0626
0.0866
+0.0381
0.0072
0.0031
0.0118
0.0462
0.0210
0.0713
0.0340
0.0074
0.0614
0.0123
0.0042
0.0203
3-11-12 12
0.0964
0.0736
0.1240
0.0482
0,0204
0.0752
0.0483
0.0334
0.0638
3-11+6
0.0858
0.0624
0.1
0,0853
0.0574
0.1123
.0,0005
000147
0.0157
3-1
-0.0730
-0.1087
-0.0545
-0.0821
-0.1091
-0.0553
0.0000
-0.0002
0.0002
3 1-1-22-1--6
0.0964
0.0745
0.1188
0.0956
0.0612
0.1297
0.0004
-0.0011
0.0040
3-1-2-6
0.1853
0.1537
0.2165
0.15
0,1519
0.2150
0.0019
-0.0551
0.0078
2-1-6
0.0323
0.0100
0.0706
0.0183
0.0578
0.0198
0.0506
0.0195
0.0863
2-1-124--13
-0,0397
-00655
-0.0124
-0.065
-0.0459
-0.0285
0.0077
0.0029
0.0131
2 2-1412-6
0.1584
0.1338
0.1828
0.1839
0.1430
0.2250
0.0008
0.0004
0.0023
53

---

Table 8: Estimated NDÉ, NIE and TE with the corresponding 95% confidence intervals (CI) across all the sub- pathways within the novely discovered mediation pathways shown in Figure c). Estimates are obtained using the proposed SIGMA framework with structural uncertainty propagation.
TB
NDE
NIB
Patlways
Estimate
C1 Lower
CI Upper
Estimate
CI Lower
CI Upper
Estimate
Lower
CI Upper
9
-0.720
-0.0997
-0.0499
0.0534
0.0014
0.0815
-0.1254
-0.1774
-0.9917
4- 9 - -13
0.0858
0.1084
0.0624
0.0756
0.0315
0.1028
0.1614
0.2148
00.1172
-
-0.5544
0.1036
10 10-+-13
-0.0359
-0.0528
-0.0128
0.0306
0.0008
0.0526
-0,0665
-0.0904
-0.0206
-12
-0.0270
-0.0442
-0.0072
0.0425
0.0127
0.0568
-0.0695
-0.0917
-0.0326
10 11-112--14--11
000441
000652
-0.0209
0.0387
0,0095
0.0652
-0,0828
-0.1152
-0.0411
- --11
-00333
0..0088
.0.0126
0.0315
0.0011
0.0452
-0.0648
-0.0842
-0.0221
8+111 -12---13
-0.0550
-0.0770
-0.0291
0.0506
0.0018
0.0814
-0.1056
-0.1515
-0.0633
8 +-15
0.0763
0.0501
0.1011
-0.0516
-0.1002
-0.0173
0.1279
0.0814
0.1919
Table 9: Estimated NDE, NIE and TE with the corresponding 95% confidence intervals (CI) across all the sub-pathway within the novely discovered mediation pathways shown in Figure [5 Estimates are obtained using classical regression- -bseed mediation analysis without structural uncertainty modeling.
NIE
Pathways
Estimate
CT Lower
CI Upper
Estimate
CI Lower
CI Upper
Estimate
CI Lower
CI Upper
3 9 + -13
-0.0723
-0,0955
-0.0500
-0.0764
-0.0987
-0.0533
0.0030
0.0014
0.0054
4 -13
0.0858
0.0614
0.1094
0.0983
0.0734
0.1221
-0.0096
0.0002
00,088
4- 9-9--111 -2
0.0554
0,0800
0.0306
0.0576
00.841
000326
0.0022
0.0007
0.0041
10 -13
0.0462
0.0212
0.0712
0.0554
0,0312
0.0809
-0,0083
0.0010
-0.0633
10 10-+-13
0.0360
0.0584
0.0080
0.0531
0.0053
0.0010
10 00-+13-2 -12
-0.0270
000150
00,0442
-0.0271
-0.0944
-0.0043
0.0001
0,0010
0.0005
13
-0.0441
-0.0681
-0.0205
-0.0521
-0.0748
-0.0295
0.0077
0.0054
0.0116
2-+1
-0.0167
-0.0391
0.0049
-0.0254
-0.0464
-0.0036
0.0090
0.0025
0.0156
2 2-11- +115
-0.0325
-0.0594
-0.0085
-0.0113
-0.0330
0.0104
-0.0212
-0,0393
-0.0051
44-2411--0
-0,0325
-0.0694
0.0028
-0.0087
-0.0315
0.0140
-0.0239
-0.0510
0.0032
8 8--13-4 +15
000325
-.0666
0.0015
0.0014
-0.0225
0.0246
-0.0340
0.07702
0.0019
0.0066
115
0.0763
00,0011
0.0515
-0.0765
-0.1014
-0,0477
0.0002
-0,0013
0.0019
Plot Mediation Effects
2-51-56
Jub 13
3-x11->12
4---12-2-
2-011-> 6
3-5 6
4=2:100-12--:1
1-
1=>1-=2-1---11
12
4-3
3-01
Total Effect
Direct Fffect
12
-0.10
-0.05
0.00
0.05
0.10
0.15
0.20
Effect Size
Figure 15: Visualization of the effect estimation on HPP data.
54


Reviews:
[{"review_id": "FTkMlZw2Ne", "reviewer": "Reviewer_6ZQR", "summary": "This paper presents SIGMA, a unified framework for causal mediation analysis. The authors identify that traditional causal mediation suffers from the issue of trying to balance *principled mechanisms for quantifying structural uncertainty* (interpretability/prior-informed) with *efficient Pathway Identification* (accuracy/performance). They build SIGMA by integrating three key components/functionalities together: structure discovery, pathway identification, and effect estimation. \n\nSIGMA uses flow-based structural equation models to learn uncertain causal structures, samples an ensemble of DAGs, and then introduces a Path Stability Score to identify high-confidence mediation pathways across those DAGs—avoiding the need for manual path specification. For each pathway, the authors also provide a Bayesian sampling inspired approach for quantifying structural uncertainty in mediation effects. The authors provide an interesting simulation study as well as a real-world application to HPP data, revealing novel mediation pathways linking sleep quality to cardiovascular health.", "questions": "I have four main concerns, which - if addressed - can substantially improve the paper, especially from an applied perspective: \n\n**Simplistic SEMs in Simulation Studies:**\nThe additive noise structural equation models (SEMs) used in the simulations (Appendix E1) appear to have relatively simplistic functional forms. This raises the concern that SIGMA's strong performance may be driven in part by the particular choice of SEMs, rather than its general robustness. It would be helpful if the authors could test SIGMA under more complex or nonlinear SEMs, or at least provide a stronger justification for why the current SEM choices reflect realistic data-generating processes.\n\n**Unclear Robustness to Noise:**\nMethods of this kind often degrade with increasing observational noise, yet it's not clear how sensitive SIGMA is to increasing noise levels. Including such a robustness analysis would help clarify how noise-tolerant SIGMA is, and where its limitations lie. \n\n**Scalability to High-Dimensional Settings:**\nWhile the simulation studies do explore dimensionalities up to 200 variables, many real-world systems—such as gene regulatory networks—involve thousands of variables. It's unclear whether SIGMA remains computationally feasible or accurate at such scales. I encourage the authors to evaluate its performance on higher-dimensional data (e.g., 500, 1000, or 5000 variables), as this would greatly help readers assess its applicability across domains.\n\n**Real-World Application:**\nFinally, while I appreciate the depth of the theoretical development, I think the paper would benefit from a second, more large-scale real-world application. For instance, applying SIGMA to a large-scale gene regulatory dataset could be a strong demonstration of its practical utility. In particular, evaluating whether SIGMA can recover known mediation pathways in well-studied biological circuits would offer a more tangible result for potential users in biology.", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 4, "strengths_and_weaknesses": "Strengths: \n- The authors present a strong theoretical foundation (sections 3 and 4) for the pathway, mediation, and DAG estimation\n- Interesting design of simulation experiments\n- Meaningful comparisons to competing methodologies\n\nWeaknesses (elaborated below): \n-  The additive noise SEMs in simulation studies have a simplistic form (Appendix E1) which may be unrealistic \n-  The robustness of the method to noise is not clear\n-  The scalability of the model is questionable\n-  Need more real world evidence", "was_revised": true, "final_justification": "The authors have provided a comprehensive and thoughtful rebuttal that addresses nearly all of the concerns I previously raised. In particular, the new experiments on robustness to systematic complexity and scalability are much appreciated and significantly strengthen the empirical foundation of the work. These additions demonstrate methodological rigor and a strong commitment to validation. Given these improvements and clarifications, I believe the paper is stronger and have accordingly raised my score.", "rebuttal": "## Q1 and Q2\n\n### **Clarification of Existing Simulation Design Complexity**\n\nWe may not have explicitly emphasized the degree of complexity already implemented. However, our simulation design incorporates substantial complexity to closely reflect real-world data-generating processes:\n\n**Functional Complexity (Appendix E.1):**\n\n* Nonlinear transformations: polynomial (v²), sinusoidal (sin(v)), and exponential (exp(-|v|)) functions (lines 1541-1542)\n* Nonlinearity ratio control: ρnonlin ∈ [0,1] systematically modulates nonlinearity degree (lines 1543-1544)\n* Parent variable interactions: pairwise interaction terms added to increase structural complexity (lines 1546-1547)\n* Non-Gaussian noise: exogenous noise zi sampled from Gaussian, exponential, and Student-t mixture distributions, simulating heteroscedasticity and non-Gaussianity (lines 1547-1549)\n\n**Data Realism (Appendix E.1):**\n\n* Missing data patterns: MCAR and MAR mechanisms simulate non-random measurement gaps common in clinical studies (lines 1560-1561)\n* Heterogeneous variable types: mixed continuous, binary, and categorical variables (πc=0.6, πb=0.3, πcat=0.1) reflecting real biomedical data characteristics (lines 1532-1534, Table 3)\n* Complex mediation structures: parallel, chain, and hybrid mediation pathways systematically embedded in data (lines 1113-1114)\n\n**Structural Diversity (Appendix E.2, E.3):**\n\n* Graph topology variation: Erdős-Rényi model configured with different chain lengths ℓ∈{1,2,6} and parallel paths κ∈{1,2,3} (lines 1521-1527, Table 3)\n* Dimensional scaling: evaluation across p∈{20,50,100,200} variables (lines 1712-1713)\n* Effect magnitude control: systematic scaling ensures identifiable direct and indirect effects (lines 1552-1558)\n\n### **Systematic Complexity and Robustness Verification**\n\nNevertheless, we completely agree with the reviewers' perspective—more systematic complexity verification will further enhance the method's credibility. Therefore, we conducted supplementary experiments:\n\n**Experimental Configuration:**\n\n* **Dataset scale**: n=6000 samples, p=50 variables\n* **Graph parameters**: ER edge density 0.3; SF attachment m∈{2,3}; Hierarchical 3-layer structure\n* **Mediation structure**: κ=2 parallel paths, ℓ=2 chain length\n* **Missing data**: 15% MCAR + 10% MAR\n* **Effect sizes**: Direct effects η ∈ [0.2, 0.6], Indirect effects θ ∈ [0.1, 0.4]\n\n**Table R1: SIGMA Robustness under Complex SEMs, Graph Structures, and Noise Levels**\n\n| Graph Type | SEM Level | SNR | Adj F1 (Mean) | Adj F1 (SE) | Orient F1 (Mean) | Orient F1 (SE) | Path F1 (Mean) | Path F1 (SE) |\n|--|---|---|----|----|------|---|---|---|\n| **ER** | Current | 4.0 | 0.672 | 0.031 | 0.618 | 0.037 | 0.651 | 0.034 |\n| **ER** | Current | 2.0 | 0.658 | 0.034 | 0.601 | 0.041 | 0.634 | 0.038 |\n| **ER** | Current | 1.0 | 0.631 | 0.042 | 0.572 | 0.049 | 0.605 | 0.045 |\n| **ER** | Current | 0.5 | 0.582 | 0.055 | 0.518 | 0.063 | 0.549 | 0.059 |\n| **ER** | High | 4.0 | 0.643 | 0.036 | 0.587 | 0.042 | 0.621 | 0.039 |\n| **ER** | High | 2.0 | 0.628 | 0.041 | 0.571 | 0.047 | 0.604 | 0.044 |\n| **ER** | High | 1.0 | 0.594 | 0.051 | 0.535 | 0.058 | 0.567 | 0.055 |\n| **ER** | High | 0.5 | 0.531 | 0.068 | 0.471 | 0.077 | 0.498 | 0.073 |\n| **SF (m=2)** | Current | 2.0 | 0.634 | 0.038 | 0.577 | 0.044 | 0.611 | 0.041 |\n| **SF (m=2)** | Current | 0.5 | 0.558 | 0.061 | 0.497 | 0.069 | 0.524 | 0.065 |\n| **SF (m=2)** | High | 2.0 | 0.615 | 0.043 | 0.556 | 0.050 | 0.588 | 0.047 |\n| **SF (m=2)** | High | 0.5 | 0.507 | 0.075 | 0.443 | 0.084 | 0.471 | 0.080 |\n| **SF (m=3)** | Current | 2.0 | 0.621 | 0.041 | 0.563 | 0.048 | 0.597 | 0.045 |\n| **SF (m=3)** | High | 1.0 | 0.574 | 0.054 | 0.513 | 0.062 | 0.547 | 0.058 |\n| **Hierarchical** | Current | 2.0 | 0.649 | 0.035 | 0.593 | 0.041 | 0.626 | 0.038 |\n| **Hierarchical** | High | 1.0 | 0.601 | 0.048 | 0.541 | 0.055 | 0.574 | 0.052 |\n\n**Extended Testing Dimensions:**\n\n* **SEM complexity progression**: Current level (original design) vs. High level (heteroscedastic noise + composite nonlinear functions)\n* **Graph structure diversity**: ER networks vs. Scale-Free networks (m∈{2,3}) vs. Hierarchical network structures\n* **Noise robustness analysis**: SNR ∈ {4.0, 2.0, 1.0, 0.5} systematic testing of observational noise impact\n\n**Key Component Mathematical Definitions:**\n\n* **SNR definition**: `SNR = Var(fi(PaG(Vi))) / Var(σi(PaG(Vi); θi) · Ui)`\n* **Current SEM**: `Vi = fi(PaG(Vi)) + σ·εi`, where `fi(x) = Σ(αjxj + βjxj² + γj sin(xj) + δj exp(-|xj|))`, `εi ~ 0.6·N(0,1) + 0.3·Exp(1) + 0.1·t(3)`\n* **High SEM**: `Vi = fi(PaG(Vi)) + σi(PaG(Vi))·εi`, where `fi(x) = tanh(Σ wjxj) × exp(Σ ujxj) + log(1 + |Σ vjxj|)`, `σi²(x) = exp(Σ γjxj)`, `εi ~ Cauchy(0,1)`\n* **Graph structures**: ER graphs `P(eij = 1) = 0.3`; SF graphs `P(ki) ∝ ki^(-γ)`; Hierarchical graphs `G = (L₁ ∪ L₂ ∪ L₃, E)` where `Li → Li+1`\n\n## Q3\n\n### Experimental Setup\n\nWe extended testing of SIGMA on 3 high-dimensional synthetic datasets with feature dimensions p ∈ {500, 1000, 2000}, nonlinearity degree ρ = 0.5, containing parallel (κ = 2) and chain (ℓ = 2) structures. Each configuration replicated 10 independent runs with sample size n = 6000.\n\n**Table R2: SIGMA High-Dimensional Scalability Analysis**\n\n| Dimension (p) | Adj F1 (Mean) | Adj F1 (SE) | Orient F1 (Mean) | Orient F1 (SE) | Path F1 (Mean) | Path F1 (SE) |\n|--|--|--|--|--|--|--|\n| 500| 0.603| 0.044| 0.541| 0.051| 0.587| 0.047|\n| 1000| 0.548| 0.058| 0.479| 0.066| 0.521| 0.062|\n| 2000| 0.467| 0.074| 0.401| 0.083| 0.429| 0.079|\n\n### Computational Complexity Analysis\n\nSIGMA's computational complexity comprises four main components:\n\n1. **Flow-SEM training phase**: O(T × p² × N), where T is the number of training iterations, p is the number of variables, and N is the sample size.\n2. **DAG posterior sampling**: O(N_DAG × p²), requiring cycle detection and removal for each sampled structure.\n3. **Path enumeration**: O(p!) in the worst case, but controllable to O(p^L) through path length constraint L.\n4. **Bayesian model averaging**: O(N_path × N × K), where N_path is the number of stable paths and K is the number of cross-validation folds.\n\n### Overall Complexity\n\nWhen implementing path length constraint L, SIGMA's total time complexity is:\n\n**O(T × p² × N + N_DAG × p² + p^L + N_path × N × K)**\n\n## Q4\n\nWe sincerely appreciate the reviewer’s thoughtful suggestion. It is important to highlight that the core design philosophy of SIGMA is precisely targeted at the identification and quantification of **long-chain complex mediation pathways**, which fundamentally distinguishes it from traditional single-step or simple mediation analyses.\n\n**Long-chain Mediation Pathways: The Core Advantage of SIGMA**\n**Design Comparison:**\n\n- **Traditional Mediation Analysis:**\n  Treatment → Mediator → Outcome (L=3)\n- **SIGMA’s Design Goal: Supporting Complex Structures:**\n  \n  - **Long-Chain Sequences:** Treatment → M1 → M2 → M3 → M4 → Outcome (L≥4)\n  - **Parallel Paths:**\nTreatment → M1 → Outcome\n\nTreatment → M2 → Outcome\n  - **Hybrid Networks:**\n\nTreatment → M1 → M3 → Outcome\n\nTreatment → M2 ↗\n  - **Branching and Convergence:**\n\nTreatment → M1 ↘ M4 → Outcome\n\nTreatment → M2 → M3 ↗\n\n**Methodological Advantages:**\n\n- **Path Stability Score:** Specifically designed to evaluate the stability of various complex pathway structures under uncertainty\n- **Bayesian Model Averaging:** Effectively addresses accumulated uncertainty in long-chain and network mediation processes\n- **Cross-level Identification:** Capable of spanning multiple biological layers for complex pathway identification\n\n**Fundamental Limitations of Existing Genomic Datasets**\n\nCurrent gene regulatory network datasets face several critical limitations that hinder their applicability for long-chain mediation analysis:\n\n**Structural Constraints:**\n\n- **Cross-sectional Limitations:** Typically cover only a single omics layer (e.g., gene expression), lacking cross-omics integration\n- **Single-Level Hierarchy:** Fail to provide a complete biological cascade from gene → protein → metabolite → phenotype\n- **Path Length Constraints:** Mostly limited to simple mediation with L=3, insufficient for supporting long-chain analysis\n\n**Practical Comparison:**\n\n- **Typical Genomic Study Capability:**\n  Gene A → Gene B → Phenotype (L=3, single omics)\n- **SIGMA’s Design Requirement:**\n  Gene → Protein → Metabolite → Biomarker → Phenotype (L=5+, cross-omics)\n\n**HPP Long-Chain Mediation Experiment Plan and Technical Support**\n\nWe are currently extending SIGMA to genomic data. The **HPP dataset** provides a multi-omics data foundation, including genomic SNPs and expression profiles, proteomic biomarkers, metabolomic product profiles, and clinical phenotypes. This multi-omics integration creates the conditions for identifying **gene-to-phenotype long-chain pathways**.\n\n**Effectiveness of Cross-sectional Data in Long-chain Mediation:**\n\nSIGMA’s theoretical foundation is built on causal graph structure rather than temporal sequence. Statistical dependencies in cross-sectional data are sufficient to reveal underlying causal pathways. Long-chain mediation paths reflect biological hierarchy and functional relationships. Our Flow-SEM can recover latent causal dependencies from cross-sectional correlational structure, while Path Stability Score ensures the statistical robustness of identified pathways.\n\n**Theoretical and Algorithmic Suitability:**\n\nSIGMA’s core algorithms are suitable for multi-omics cross-sectional data: 1) Flow-SEM handles heterogeneous data types;2) DAG posterior sampling quantifies causal relationships across omics layers;3) Path Stability Score evaluates cross-biological-level pathway robustness;4) Bayesian Model Averaging integrates effect uncertainty across omics.\nAs shown in the scalability analysis (Q3), the algorithm remains computationally feasible at p=2000, supporting integrative multi-omics analysis."}, {"review_id": "n595YSLPoU", "reviewer": "Reviewer_rfqt", "summary": "The paper proposes a new method for mediation analysis called SIGMA, which has three main steps: causal structure learning, mediation path extraction, and estimation of the natural direct and indirect effects. SIGMA is end-to-end, and the estimated effects account for uncertainty about the causal structure. SIGMA is compared to other algorithms on synthetic data and also used to analyze a real health dataset.", "questions": "1. What does \"identify\" (and all the variations on this root) mean in this paper? Within the causal inference and discovery literature, I'm used to seeing \"identify\" in the formal sense of expressing a causal estimand in terms of observed distributions. This paper seems to mostly (but not always) mean something more like \"estimate\". While identifiability is required for unbiased estimation, these are distinct concepts, and conflating them detracts from the rigor and clarity of the paper. (A similar point can be made for \"bias\" versus \"error\", e.g., on line 237 and in Table 1.)\n\nAs a follow-up, is the paper presenting/making use of new identifiability results for mediation paths beyond those implied by the usual causal discovery setup? I don't see an explicit formal statement of this anywhere, but it seems to be implied (e.g. top of Page 2).\n\n2. Can the authors give some sense of the complexity or runtime of the proposed method? This is important to understand how practical it actually is. Ideally, there would be explicit statements/proofs of the complexity along with runtime comparisons to other algorithms.\n\n3. How do the SIGMA results on the real dataset compare to other methods? The current results suggest SIGMA is reasonable but they don't demonstrate its significance compared to existing methods.\n\n*Minor comments*:\n1. L. 111: should cite the original Erdős and Rényi paper rather than an arbitrary recent causal discovery paper.\n2. Section 3: would be helpful to relate these to the Markov and faithfulness assumptions---for example, isn't conditional ignorability satisfied when conditioning on the parents of A as described in Section 4.4.1?\n3. L. 158: cofounders -> covariates?\n4. L. 322--323: \"successfully identified multiple physiologically meaningful mediation pathways\" seems a bit strong since there's no ground truth; something along the lines of \"learned multiple physiologically plausible pathways\" seems to be all that's justfied based on the current results\n5. I don't really get the distinction between DECI and the others in lines 349--351. For example, GES uses the MLE (via the BIC score) for structure learning, so isn't it also jointly inferring causal structure and structural equations?", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "*Quality*:\n- (+) the paper is comprehensive, touching upon many ideas from causal discovery and inference as well as its health application\n- (-) however, technical depth and rigor is lacking, with some important concepts left ambigous, undefined, or conflated (see Question 1)\n- (-) state-of-the-art causal discovery methods (like BOSS and order MCMC; see [here](https://benchpressdocs.readthedocs.io/en/latest/examples.html#study-4-fpr-tpr-pattern-png) for more details) aren't used in comparisons; for the effect estimation task, there's no comparison to an obvious (causal discovery-based) baseline like using the MLE for the output of LiNGAM\n\n*Clarity*: (+) the paper is generally well-written ([-] with the important exception of some technical imprecision)\n\n*Significance*: (+) there's a lot of potential value in an end-to-end mediation analysis that accounts for structural uncertainty, but (-) I don't see much theoretical significance here, and (-) the practical significance remains unclear based on the current experiments\n\n*Originality*: (+) the approach for explicitly accounting for structural uncertainty in the causal effect estimate is novel (and seems to sensibly builds off/connects to some existing work in an interesting way); but (-) it's not always so clear exactly what is novel versus what is existing work, e.g. Section 4 contains no references at all", "was_revised": true, "final_justification": "Through the rebuttal process, I understand the submission better and think it's probably free of any serious mathematical mistakes (so I increase my rating to 3)---but the current presentation in the submission does not make this sufficiently clear, so I do not increase my rating more unless/until the authors address the rest of my remaining concerns.\n\n**Update after another iteration:** The authors have committed to a number of changes that will significantly improve the presentation clarity and evaluation of the work/comparison to existing methods, so I increase my rating to 5 (assuming the authors will successfully implement the changes as promised).\n\n**another update:** The authors have provided some of the requested results, making their promises to implement the rest more trustworthy.", "rebuttal": "## **Q1**\n\nWe thank the reviewer but must point out a disciplinary understanding discrepancy. \n\n### Our use of \"identify\" follows standard practice in causal analysis and mediation analysis literature:\n\n**Statistical identification level**\n\n- Pearl (2001): \"identification of causal effects transmitted through intermediate variables\"\n- Robins & Greenland (1992): \"identifiability of direct and indirect effects\"\n- VanderWeele (2015): \"identification of mediation effects\"\n- Our line 139: \"Identifying NDE and NIE\"\n\nOther standard usage: Imai et al. (2010), Tchetgen & Shpitser (2012), Baron & Kenny (1986), Sobel (1982), Holland (1988), etc.\n\n**Pathway identification level**\n\n- Daniel et al. (2015): \"identification of path-specific effects\"\n- Avin et al. (2005): \"identifiability of path-specific effects\"\n- VanderWeele & Vansteelandt (2014): \"identifying multiple mediation pathways\"\n- Our line 206: \"identifies mediation pathways\"\n\nOther standard usage: Miles et al. (2017), Steen et al. (2017), Nguyen et al. (2021), Bind et al. (2016), etc.\n\n**Bias vs Error terminology**\n\n- Chernozhukov et al. (2018): explicitly distinguish \"bias\" (systematic deviation) and \"error\" (total error)\n- Kennedy (2020): \"bias-variance decomposition\"\n- Our Table 1: Bias refers to E[θ̂]-θ₀, SE refers to standard error, CI refers to confidence intervals\n\nOther standard practices: Bickel et al. (1993), Tchetgen & Shpitser (2012), van der Laan & Rose (2011), etc.\n\n### Question about new identifiability theory:\n\nSIGMA has never claimed new identifiability theory. Our contribution is algorithmic/methodological innovation, addressing practical limitations of existing methods in real applications: systematic identification of multi-step pathways, pathway selection under structural uncertainty, complex pathway effect estimation.\n\nSIGMA identified previously unspecified multi-step pathways connecting sleep quality to cardiovascular health, which traditional methods cannot systematically discover.\n\n## **Q2:**\n\nSIGMA's computational complexity comprises four main components: (1) Flow-SEM training phase: O(T × p² × N), where T is the number of training iterations, p is the number of variables, and N is the sample size; (2) DAG posterior sampling: O(NDAG × p²), requiring cycle detection and removal for each sampled structure; (3) Path enumeration: O(p!) in the worst case, but controllable to O(p^L) through path length constraint L; (4) Bayesian model averaging: O(Npath × N × K), where Npath is the number of stable paths and K is the number of cross-validation folds.\n\n**Overall Complexity:** When implementing path length constraint L, SIGMA's total time complexity is:\n**O(T × p² × N + NDAG × p² + p^L + Npath × N × K)**\n\n## **Q3:**\n\n### Existing methods are technically unable to handle our problem setting:\n\n**1. Fundamental differences in method capabilities:**\n\n**Traditional mediation analysis methods** (DeepMed)  requires pre-specified pathways, computationally infeasible in high-dimensional data\n\n**Causal discovery methods** (PC, NOTEARS, ICALiNGAM) can only learn structures but cannot quantify mediation effects, unable to handle pathway selection under structural uncertainty\n\n**2. Real data lacks ground truth making quantitative comparison impossible. SIGMA provides end-to-end analysis while existing methods handle individual components.**\n\n### Our validation strategy follows field standards:\n\n**Synthetic data quantitative comparison + real data domain validation is established practice in causal analysis:**\n\nSuch as Chernozhukov et al. (2018), Kennedy (2020), VanderWeele (2015), Tchetgen & Shpitser (2012) all adopt this strategy\n\n### We adopted standard biomedical validation:\n\n**For effect estimation validation, we adopted the gold standard approach used in top journals like Nature Medicine and Nature Communications:** Using traditional regression mediation analysis (linear regression with bootstrapping) to independently validate SIGMA-identified pathways, by separately establishing X→M and M→Y regression models, calculating indirect effect ab and direct effect c', and using bootstrap confidence intervals for statistical testing. This method is the established standard in medical journals, aligned with biological research requirements. **Real data validation focused on biological rigor:** (1) **Known mechanism reproduction**: validating established medical pathways; (2) **Biological plausibility**: new pathways conform to physiological mechanisms; (3) **Independent statistical validation**: traditional regression mediation analysis confirms effect consistency. Our validation shows SIGMA results are highly consistent with traditional methods in effect direction and magnitude, proving the method's reliability. **In medical data lacking ground truth, domain expertise validation is the standard scientifically valid approach.**\n\n### **MC1:**\n\nWe thank the reviewer for this correction. We will cite in the revision.\n\n### **MC2:**\n\nThis requires conceptual clarification. **The Markov and Faithfulness assumptions mentioned by the reviewer and the conditional ignorability assumptions we discuss in Section 3 operate at different levels of causal inference: the former are used in the structure discovery stage, while the latter are used for effect identification given a structure.** Regarding the key question of \"whether conditioning on the parents of A automatically satisfies conditional ignorability,\" we address this in Appendix A.3.1, and **the answer is no.** Simply constructing adjustment sets based on learned graph structures cannot automatically guarantee satisfaction of ignorability assumptions; misspecified structures violate key identification assumptions. **A core contribution of SIGMA is precisely handling the impact of such structural uncertainty on effect estimation.** We propagate uncertainty through structural posterior sampling and Bayesian model averaging, rather than assuming correct structure. This theoretical treatment ensures our method's inference validity under real-world structural uncertainty.\n\n### **MC3:**\n\n**\"Confounders\" is the theoretically necessary choice.** In the mediation analysis framework of Pearl (2009) and VanderWeele (2015), \"confounders\" has a clear operational definition: variables that simultaneously affect relationships between multiple variables in the causal chain, which would lead to systematic bias if not properly adjusted. **This differs from 'covariates,' which refers to arbitrary predictor variables.** In the context of line 158, we discuss \"observational data over treatment, outcome, mediators, and confounders\", a classic causal inference setting where the core issue is identification and control of causal confounding, representing the direct manifestation of three key no-confounding assumptions in mediation analysis. **A core theoretical contribution of SIGMA is precisely handling these causal confounding problems under structural uncertainty conditions.** Using the statistical term \"covariates\" would **obscure the causal inference nature of our method and weaken the clear expression of its theoretical contributions.**\n\n### **MC4:**\n\n**We use \"successfully identified\" for mechanisms already reported in medical literature.** We explicitly state: \"the pathway (sleep HRV index → pulse wave velocity → carotid IMT) verified a known mechanism whereby sleep quality affects cardiovascular health via vascular function,\" and cite relevant medical studies (Bourdillon et al. 2021; Kadoya & Koyama 2019; Saz-Lara et al. 2022) as validation evidence. **For these pathways with existing ground truth literature support, using \"successfully identified\" is well-justified.**\n\n**For newly discovered mediation pathways, we adopted more cautious wording while employing standard medical validation methods.** As shown in line 328, we use \"identified\" rather than \"successfully identified,\" and as stated in line 339: \"Classical medical regression mediation analysis further validated these pathways, with an overall TE error of 0.2734 and consistent NIE and NDE estimation directions.\"\n\n### **MC5:**\n\n**While GES uses MLE and BIC scoring, it differs fundamentally from DECI in \"joint inference.\"** GES employs conditional inference: for each candidate structure G, it first fixes G to estimate parameters θ̂_G = argmax P(D|θ,G), then uses θ̂_G to calculate BIC scores for structure evaluation—this is a typical P(θ|G,D) conditional estimation process. **In contrast, DECI and other deep learning methods directly optimize the joint posterior P(G,θ|D), simultaneously learning structural distributions and functional parameters through variational inference, achieving true end-to-end optimization.** Additionally, GES relies on preset parametric models (typically linear Gaussian) and searches in discrete DAG spaces, while DECI uses neural networks to learn arbitrary nonlinear relationships and optimizes in continuous parameter spaces.\n\n### **W1:**\n\n**BOSS and order MCMC methods are incompatible with mixed data types.** BOSS is designed for linear Gaussian continuous data (Andrews et al., 2022), Order MCMC's BiDAG documentation states: \"acceptable data matrices are homogeneous with all variables of the same type.\" Our datasets contain heterogeneous variable types (continuous, binary, categorical), making these methods' core assumptions invalid. Recent mixed data causal discovery work consistently excludes these methods: Raghu et al. (2018) ACM SIGKDD, NOTEARS-M (2024), DECI (NeurIPS 2022) all systematically exclude BOSS/order MCMC, reflecting established practice.\n\n### **W2:**\n\nWe cited SIGMA's core components in earlier sections: Flow-SEM [69] cited on line 59, EIF [6] cited on lines 47 and 64, BMA [53] cited on line 67. Section 6 discussed existing research (NOTEARS, GraN-DAG, DECI, etc). Section 4 focuses on theoretical innovation and mathematical proofs based on these components according to AI conference paper standards."}, {"review_id": "iz2UANMZXK", "reviewer": "Reviewer_MyHV", "summary": "This paper addresses the challenge of handling structural uncertainty in Causal Mediation Analysis (CMA). Traditional CMA methods either assume the causal structure is known, making them unable to handle unknown structures, or they fail to provide consistent inference that properly accounts for structural uncertainty, even when the structure is estimated. To overcome these limitations, the authors propose a novel framework called SIGMA (Structure-Informed Guided Mediation Analysis). SIGMA leverages a Flow-Structural Equation Model (Flow-SEM) to model the posterior probabilities of edges in possible Directed Acyclic Graphs (DAGs). It introduces a new metric, Path Stability Score (PSS), to identify high-confidence mediation pathways. Effect estimation is then performed based on these stable pathways.\n\nExperiments on synthetic datasets demonstrate that SIGMA achieves more accurate structure recovery than conventional causal discovery methods and more reliable mediation path identification compared to existing CMA approaches. Furthermore, application to real-world data related to sleep and cardiovascular health shows that SIGMA can successfully reconstruct known mediation mechanisms and uncover previously unreported pathways, highlighting its potential for new scientific discovery.", "questions": "1. Would it be possible to compare the proposed method against conventional structure learning methods that are enhanced with bootstrap-style techniques to account for structural uncertainty? If the proposed approach can be shown to outperform such alternatives, it would significantly strengthen the overall evaluation of the paper.\n\n2. For the novel mediation pathway shown in Figure 3(c), is it possible to cite any prior studies, however minimal, that suggest the plausibility of the entire pathway as a coherent mechanism? While the paper acknowledges that individual links have been previously studied, it also states that the complete pathway has not been reported as such. If even partial evidence supporting the full path exists, the method could be seen as successfully identifying an under-recognized yet plausible hypothesis. This would highlight the proposed approach as a powerful tool for supporting scientific discovery and could substantially elevate the contribution of the paper.", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\nThe paper presents a highly robust and reliable framework for Causal Mediation Analysis (CMA) by incorporating structural uncertainty into effect estimation through the combination of Flow-SEM and Path Stability Score (PSS), effectively leveraging a relaxed acyclicity constraint. This integration enables a novel end-to-end treatment of uncertainty that surpasses prior approaches. The effectiveness of the proposed method is quantitatively demonstrated through comprehensive experiments on synthetic data. Furthermore, its applicability and potential for scientific discovery are substantiated via real-world evaluation using sleep-related physiological data.\n\nWeaknesses:\nWhile the proposed treatment of structural uncertainty is conceptually compelling and methodologically sound, the paper does not include a comparison against baseline methods that incorporate structural uncertainty via bootstrap-style approaches, which are feasible even for conventional causal discovery algorithms. Although it is intuitively and theoretically plausible that modeling structural uncertainty directly within the learning framework, as done in this paper, yields superior results compared to such post-hoc strategies, empirical validation of this advantage would substantially strengthen the paper's claims. Including such a comparison in the synthetic experiments would significantly enhance the persuasiveness and completeness of the evaluation.", "was_revised": true, "final_justification": "The comparison experiments with the bootstrap baselines have fully addressed my concerns regarding the novelty of the proposed method, and I am now convinced of its innovative nature. Furthermore, I believe you have provided clear and sufficient biological evidence to support your claims. I am confident that incorporating these points into the manuscript will further enhance its persuasiveness. In light of these considerations, I have decided to raise my rating by one point.", "rebuttal": "## Q1\n\nWe are very grateful to the reviewer for this valuable suggestion. The comparison with the bootstrap-based method for handling structural uncertainty that you pointed out is indeed an important and valuable evaluation dimension.\n\n### Methodological differences\n\nSIGMA and the bootstrap method have fundamental methodological differences in handling structural uncertainty:\n\n**SIGMA's Methodology:**\n- **End-to-end posterior learning of structures:** Directly learns the posterior distribution of structures `P(G|D)` through Flow-SEM to generate a DAG ensemble `{G_s}`.\n- **Native uncertainty modeling:** Structural uncertainty is directly encoded in the model's parameterization.\n- **Integrated propagation:** Uncertainty is propagated throughout the entire process, from structure discovery to path identification to effect estimation.\n\n**Limitations of the Bootstrap Method:**\n- **Two-stage process:** First obtains a point estimate of the structure, then estimates uncertainty through data resampling.\n- **Indirect uncertainty quantification:** Relies on data perturbation to indirectly infer structural changes.\n- **Potential underestimation:** May systematically underestimate the true structural uncertainty, especially in high-dimensional sparse scenarios.\n\n### Theoretical Framework for Bayesian Quantification of Structural Uncertainty\n\nBased on our theoretical analysis in Appendix B.1, we can argue for the superiority of the SIGMA method from a Bayesian inference perspective.\n\n#### Direct Modeling of Posterior Probability\n\nAccording to the theoretical results in Appendix B.1.1, our Path Stability Score (PSS) is an unbiased and consistent estimator of the true posterior inclusion probability:\n\n**PSS formula:** `PSS(π) = (1/N_valid) Σ_{s=1}^{N_valid} I(π ⊆ G_s)`\n\nwhere Theorem 1 proves its **unbiasedness**: `E[PSS(π)] = P(π ⊆ G|D)`, and **consistency**: `PSS(π) →^p P(π ⊆ G|D)` as `N_valid → ∞`.\n\n#### Theoretical Advantages of Bayesian Model Averaging\n\nOur BMA estimator is directly based on the learned structural posterior:\n\n**BMA estimator:** `θ̂_BMA(π) = (1/N_path) Σ_{s=1}^{N_path} θ̂_s(π)`\n\nAccording to the convergence analysis in Appendix B.1.2, as the number of DAG samples increases, our variance estimation:\n\n**Variance decomposition:** `Var(θ̂_BMA(π)) = (1/N_path) Σ_{s=1}^{N_path} Var(θ̂_s(π)) + (1/N_path) Σ_{s=1}^{N_path} (θ̂_s(π) - θ̂_BMA(π))²`\n\ncan simultaneously capture both **estimation uncertainty** (the first term) and **structural uncertainty** (the second term).\n\n### Bootstrap Baseline Methods\n\n For each of the `B=100` bootstrap iterations, we first generate a data subset `D_b` by sampling with replacement from the original dataset `D`. Subsequently, the corresponding base learner is used to learn a causal graph `G_b = StructureLearner(D_b)`, from which all potential mediation paths `Π_b = PathIdentifier(G_b)` are identified. For each path `π ∈ Π_b`, its effect is then estimated as `θ̂_b(π) = EffectEstimator(π, G_b, D_b)`. \n\nAfter completing all `B` iterations, the results are aggregated:\n- **Path stability:** `PSS_boot(π) = (1/B) Σ_{b=1}^B I(π ∈ Π_b)`\n- **Path effect estimate:** `θ̂_boot(π) = (1/B) Σ_{b=1}^B θ̂_b(π)`\n- **Variance estimation:** `Var_boot(π) = (1/(B-1)) Σ_{b=1}^B (θ̂_b(π) - θ̂_boot(π))²`\n\nRegarding the experimental configuration, we selected four representative synthetic datasets: LowDim-L, LowDim-N, MidDim-S, and MidDim-D. The sample size was set to `n = 6000` (consistent with the original experiments), with data dimensions of `p ∈ {20, 100}` and a non-linearity proportion of `ρ_nonlin ∈ {0, 0.5}`. Bootstrap parameters were kept consistent with SIGMA, with the number of resamples `B = 100` and a path stability threshold of `τ = 0.05`. Each bootstrap iteration used a sample size equal to that of the original data to ensure comparability. \n\nFor the base method parameters, the PC algorithm used a significance level of `α = 0.05` with a Gaussian conditional independence test. The NOTEARS method used a sparsity parameter `λ_1 = 0.01` and a DAG constraint parameter `λ_2 = 0.01`, with a maximum of 100 iterations. The GraN-DAG method was configured with a temperature parameter `τ_g = 0.2` and was trained for 1000 epochs.\n\n**Table r1: Bootstrap Baseline Comparison**\n\n|Dataset|Method|Adjacency F1|Orientation F1|Path Recovery Rate|False Discovery Rate|Avg PSS (True)|Avg PSS (False)|\n|---|---|---|---|---|---|---|---|\n|**LowDim-L**|PC+Bootstrap|0.4082|0.2326|0.347|0.312|0.425|0.089|\n||NOTEARS+Bootstrap|0.3728|0.2158|0.321|0.298|0.398|0.076|\n||GraN-DAG+Bootstrap|0.3922|0.2247|0.334|0.287|0.412|0.081|\n||**SIGMA(Ours)**|**0.5263**|**0.4474**|**0.521**|**0.124**|**0.645**|**0.052**|\n|**LowDim-N**|PC+Bootstrap|0.4287|0.2286|0.298|0.367|0.384|0.105|\n||NOTEARS+Bootstrap|0.4061|0.2131|0.276|0.345|0.356|0.098|\n||GraN-DAG+Bootstrap|0.4161|0.2172|0.287|0.321|0.371|0.091|\n||**SIGMA(Ours)**|**0.4909**|**0.2601**|**0.456**|**0.143**|**0.578**|**0.063**|\n|**MidDim-S**|PC+Bootstrap|0.4427|0.2869|0.234|0.445|0.312|0.134|\n||NOTEARS+Bootstrap|0.4018|0.2636|0.218|0.398|0.289|0.118|\n||GraN-DAG+Bootstrap|0.4216|0.2738|0.231|0.376|0.305|0.112|\n||**SIGMA(Ours)**|**0.4939**|**0.3185**|**0.387**|**0.187**|**0.498**|**0.074**|\n|**MidDim-D**|PC+Bootstrap|0.4479|0.2929|0.189|0.498|0.267|0.156|\n||NOTEARS+Bootstrap|0.4141|0.2695|0.174|0.452|0.245|0.142|\n||GraN-DAG+Bootstrap|0.4360|0.2859|0.185|0.424|0.258|0.135|\n||**SIGMA(Ours)**|**0.4998**|**0.3163**|**0.312**|**0.234**|**0.423**|**0.089**|\n\n**Table r2: Comparison of Effect Estimation (NDE / NIE)**\n\n|Dataset|Method|Bias(NDE/NIE)|SE(NDE/NIE)|Coverage Probability(NDE/NIE)|Interval Width(NDE/NIE)|\n|---|---|---|---|---|---|\n|**LowDim-L**|PC+Bootstrap|0.0847/0.0693|0.1523/0.1298|0.891/0.902|0.421/0.387|\n||NOTEARS+Bootstrap|0.0672/0.0548|0.1389/0.1176|0.908/0.921|0.398/0.364|\n||GraN-DAG+Bootstrap|0.0634/0.0512|0.1301/0.1134|0.915/0.928|0.387/0.351|\n||**SIGMA(Ours)**|**0.0421/0.0347**|**0.1087/0.0923**|**0.947/0.954**|**0.356/0.318**|\n|**LowDim-N**|PC+Bootstrap|0.1156/0.0934|0.1876/0.1567|0.847/0.863|0.467/0.428|\n||NOTEARS+Bootstrap|0.0923/0.0756|0.1654/0.1389|0.873/0.887|0.445/0.401|\n||GraN-DAG+Bootstrap|0.0851/0.0687|0.1598/0.1323|0.886/0.899|0.431/0.389|\n||**SIGMA(Ours)**|**0.0583/0.0468**|**0.1245/0.1045**|**0.932/0.941**|**0.389/0.342**|\n|**MidDim-S**|PC+Bootstrap|0.1453/0.1278|0.2234/0.1945|0.824/0.836|0.523/0.476|\n||NOTEARS+Bootstrap|0.1187/0.1021|0.1987/0.1723|0.851/0.864|0.498/0.451|\n||GraN-DAG+Bootstrap|0.1089/0.0934|0.1876/0.1634|0.867/0.878|0.482/0.434|\n||**SIGMA(Ours)**|**0.0724/0.0621**|**0.1467/0.1267**|**0.918/0.926**|**0.421/0.378**|\n|**MidDim-D**|PC+Bootstrap|0.1789/0.1567|0.2678/0.2345|0.789/0.803|0.587/0.534|\n||NOTEARS+Bootstrap|0.1456/0.1289|0.2341/0.2067|0.816/0.829|0.556/0.503|\n||GraN-DAG+Bootstrap|0.1324/0.1156|0.2198/0.1934|0.834/0.847|0.534/0.485|\n||**SIGMA(Ours)**|**0.0892/0.0785**|**0.1723/0.1456**|**0.893/0.904**|**0.467/0.421**|\n\n## Q2\n\nWe sincerely thank the reviewer for this excellent suggestion, which significantly strengthens our contribution. Recent studies published in 2025 and substantial evidence from recent years potentially support the biological plausibility of this pathway.\n\n### Recent Clinical Evidence\n\nThe HPP study published in Nature Medicine in January 2025[1], based on sleep monitoring data from 6,366 participants across 16,812 nights, found that visceral adipose tissue was the body characteristic most strongly correlated with the peripheral apnea-hypopnea index, and sleep characteristics could predict over 15% of physiological features across 15 body systems. Simultaneously, the HPP deep phenotyping study from July 2025[2] further revealed complex associations among sleep, metabolic, cardiovascular, and gut microbiome systems.\n\n### Sleep Apnea-Gut Microbiome-Cardiovascular Pathway\n\nZhang et al.[3] demonstrated that OSA leads to gut dysbiosis through intermittent hypoxia, causing cardiovascular complications via inflammatory cascades. Recent research directly proved gut microbiota mediate OSA-induced atherosclerosis[4]. Liu et al.[5] confirmed causal relationships between gut microbiome and OSA through Mendelian randomization.\nDurgan et al.[6] showed OSA-induced gut dysbiosis promotes hypertension by reducing butyrate-producing bacteria. Visceral adipocytes secrete pro-inflammatory factors (TNF-α, IL-6, IL-1β) that participate in atherosclerosis development, with visceral adipose tissue predicting subclinical carotid atherosclerosis more accurately than waist circumference[7,8].\n\nReferences:\n\n[1] Phenome-wide associations of sleep characteristics in the Human Phenotype Project\n\n[2] Deep phenotyping of health–disease continuum in the Human Phenotype Project\n\n[3] Recent progresses in gut microbiome mediates obstructive sleep apnea‐induced cardiovascular diseases\n\n[4] Gut microbiota and derived metabolites mediate obstructive sleep apnea induced atherosclerosis\n\n[5] Causal relationships between gut microbiome and obstructive sleep apnea: a bi-directional Mendelian randomization\n\n[6] Role of the Gut Microbiome in Obstructive Sleep Apnea-Induced Hypertension\n\n[7] Visceral Adipose Tissue Is a Better Predictor of Subclinical Carotid Atherosclerosis Compared with Waist Circumference\n\n[8] Visceral adipose tissue as a source of inflammation and promoter of atherosclerosis"}, {"review_id": "KpQcpljrpn", "reviewer": "Reviewer_APS7", "summary": "The paper proposes an end-to-end framework called SIGMA (Structure-Informed Guided Mediation Analysis), which aims to address the disconnection between structure discovery and the understanding of mediation mechanisms in existing causal mediation analysis. SIGMA integrates three key steps: probabilistic causal structure discovery, automated mediation pathway identification, and effect estimation, while systematically quantifying and propagating structural uncertainty. Specifically, the framework first utilizes differentiable Flow-Structural Equation Models (Flow-SEM) to learn the structural posterior distribution, generating diverse Directed Acyclic Graphs (DAGs) to quantify structural uncertainty. Based on these DAGs, the paper introduces the Path Stability Score (PSS) to evaluate the marginal probability of pathways, thereby identifying high-confidence mediation paths. For the identified pathways, SIGMA combines Efficient Influence Functions (EIF) and Bayesian model averaging (BMA) to estimate pathway effects, propagating both within-structure estimation uncertainty and between-structure effect variation to the final effect estimates.", "questions": "1.  As mentioned in the weaknesses, some experiments could be supplemented.\n2.  In Appendix F1, it is stated that the frequency threshold is set to 5%. Is this \"frequency threshold\" the same as the τ mentioned earlier in the text? If so, why was τ valued at 0.15 previously and 0.05 here?\n3.  I have difficulty understanding the descriptions related to Section 5.3 and Figure 3. For example, I did not see the path \"sleep HRV index → pulse wave velocity\" 【line:324】described in the text explicitly labeled in Figure 3. Could the description be made clearer?\n4.  \"structures uncertainty\" 【line:22】should be \"structural uncertainty\".\n5. The author is requested to carefully check and confirm the usage of specialized terms within the manuscript. For example, a given specialized term should be written in its full form only at its first mention (as is the case with BMA and PSS in the text). Additionally, the capitalization (of initial letters) for these terms should be consistent (e.g., Path Stability Score, NOTEARS), and so on.", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths:**\n\n1.  This paper proposes a unified framework (SIGMA) that integrates causal structure discovery, mediation pathway identification, and effect estimation, effectively addressing the problems caused by the separation of these stages in traditional methods.\n2.  The SIGMA framework combines multiple advanced statistical and machine learning techniques, such as Flow-SEM, EIF, BMA, and a Variational Autoencoder (VAE) imputation module, providing a relatively solid theoretical foundation. The authors have also provided quite detailed proofs for the core theoretical basis of this paper.\n3.  This paper conducts experiments on both real-world and synthetic datasets. The experiments are presented from relatively sufficient perspectives, including linear/non-linear scenarios and different data scales, which can effectively demonstrate the feasibility of its method.\n\n**Weaknesses:**\n\n1.  The SIGMA framework may contain multiple computationally intensive steps, making an analysis of its time complexity necessary. The authors should at least provide running times on datasets of different scales to illustrate the relationship between time consumption and data size.\n2.  The synthetic datasets mentioned in the paper are all constructed based on the ER model; however, real-world data often do not entirely conform to the distribution of ER models. The authors should attempt experiments on different models (e.g., SF).", "was_revised": true, "final_justification": "The authors have convincingly addressed my questions. I have raised my score.", "rebuttal": "## Q1-W1\n\n**We thank the reviewer for this important concern regarding time complexity.**\n\n### **Theoretical Time Complexity Analysis**\n\nSIGMA's computational complexity comprises four main components: (1) Flow-SEM training phase: O(T × p² × N), where T is the number of training iterations, p is the number of variables, and N is the sample size; (2) DAG posterior sampling: O(NDAG × p²), where NDAG is the number of sampled DAGs; (3) Path enumeration: controllable to O(p^L) through path length constraint L, where L is the maximum path length; (4) Bayesian model averaging: O(Npath × N × K), where Npath is the number of stable paths and K is the number of cross-validation folds.\n\n**Overall Complexity:** O(T × p² × N + NDAG × p² + p^L + Npath × N × K)\n\n### **Experimental Setup and Runtime Verification**\n\nWe conducted systematic runtime testing on Intel Xeon E5-2690 v4 processor (28 physical cores, 2.6GHz) with 128GB DDR4 memory. The experiments used 5 synthetic datasets with feature dimensions p∈{50,100,200,500,1000}, nonlinearity degree ρ=0.5, containing parallel (κ=2) and chain (ℓ=2) structures, sample size n=6000, variable types with continuous variable proportion πc=0.6, binary variable πb=0.3, categorical variable πcat=0.1.NDAG=100 DAG samples, T=1000 training iterations, path length constraint L=3, and 5-fold cross-validation.\n\n### **Parallel Computing Implementation and Optimization**\n\nSIGMA's four main components support different degrees of parallelization. In the Flow-SEM training phase, sample batches can be processed in parallel across cores, while different blocks of the weight matrix can compute gradients in parallel. The DAG posterior sampling phase exhibits parallel characteristics, as NDAG DAG samples are mutually independent, enabling parallel processing with each core handling NDAG/cores samples. Cycle detection and sample validation algorithms can be executed independently in parallel.\n\nPath enumeration and effect estimation phases demonstrate the algorithm's parallel potential. Path enumeration for different DAGs can be processed simultaneously, and Path Stability Score calculations can be executed in parallel across multiple candidate paths. In effect estimation, different folds of K-fold cross-validation can run independently, effect estimation for different stable pathways can be computed in parallel, and bootstrap resampling for uncertainty quantification also supports parallel processing.\n\n### **Performance Evaluation Results**\n\n**Table R1: SIGMA Runtime Analysis (hours)**\n\n| Dimension (p)  | Single-core | Multi-core (8 cores) | Multi-core (16 cores) |\n| -- |-- |--|--|\n| 50| 1.02| 0.39| 0.26|\n|100 | 4.15| 1.59| 1.14|\n|200 | 16.12| 6.18| 4.83 |\n|500 | 100.84| 38.46 | 23.11|\n|1000| 399.72|154.08 | 109.51|\n\n### **GPU Acceleration Optimization**\n\nMatrix multiplication operations in Flow-SEM training can be migrated to GPU, utilizing CUDA cores to parallelize forward and backward propagation computations of weight matrix Wω. Gradient computations can be accelerated through GPU tensor operations, particularly when computing ∂L/∂Wω for large-scale matrix derivatives. Probability computations in the DAG sampling phase can leverage GPU's parallel architecture to simultaneously process multiple Bernoulli sampling operations. Bootstrap resampling in effect estimation can be parallelized through GPU random number generation, while influence function calculations can utilize GPU's vector processing capabilities to handle large batches of samples.\n\n## Q1-W2\n\n**We appreciate the reviewer's important observation regarding graph structure diversity.**\n\n### **Scale-Free Network Experiments**\n\nTo address this limitation, we conducted supplementary experiments testing SIGMA on Scale-Free networks using the Barabási-Albert model with different attachment parameters. As detailed in our response to Reviewer 6ZQR\n\n**For the Experimental Configuration, Extended Testing Dimensions, and Key Component Mathematical Definitions, please see our rebuttal for reviewer 6ZQR due to word count restrictions.**\n\n**Table R2: SIGMA Robustness under Complex SEMs, Graph Structures, and Noise Levels**\n\n| Graph Type | SEM Level | SNR | Adj F1 (Mean) | Adj F1 (SE) | Orient F1 (Mean) | Orient F1 (SE) | Path F1 (Mean) | Path F1 (SE) |\n|--|---|---|----|----|------|---|---|---|\n| **ER** | Current | 4.0 | 0.672 | 0.031 | 0.618 | 0.037 | 0.651 | 0.034 |\n| **ER** | Current | 2.0 | 0.658 | 0.034 | 0.601 | 0.041 | 0.634 | 0.038 |\n| **ER** | Current | 1.0 | 0.631 | 0.042 | 0.572 | 0.049 | 0.605 | 0.045 |\n| **ER** | Current | 0.5 | 0.582 | 0.055 | 0.518 | 0.063 | 0.549 | 0.059 |\n| **ER** | High | 4.0 | 0.643 | 0.036 | 0.587 | 0.042 | 0.621 | 0.039 |\n| **ER** | High | 2.0 | 0.628 | 0.041 | 0.571 | 0.047 | 0.604 | 0.044 |\n| **ER** | High | 1.0 | 0.594 | 0.051 | 0.535 | 0.058 | 0.567 | 0.055 |\n| **ER** | High | 0.5 | 0.531 | 0.068 | 0.471 | 0.077 | 0.498 | 0.073 |\n| **SF (m=2)** | Current | 2.0 | 0.634 | 0.038 | 0.577 | 0.044 | 0.611 | 0.041 |\n| **SF (m=2)** | Current | 0.5 | 0.558 | 0.061 | 0.497 | 0.069 | 0.524 | 0.065 |\n| **SF (m=2)** | High | 2.0 | 0.615 | 0.043 | 0.556 | 0.050 | 0.588 | 0.047 |\n| **SF (m=2)** | High | 0.5 | 0.507 | 0.075 | 0.443 | 0.084 | 0.471 | 0.080 |\n| **SF (m=3)** | Current | 2.0 | 0.621 | 0.041 | 0.563 | 0.048 | 0.597 | 0.045 |\n| **SF (m=3)** | High | 1.0 | 0.574 | 0.054 | 0.513 | 0.062 | 0.547 | 0.058 |\n| **Hierarchical** | Current | 2.0 | 0.649 | 0.035 | 0.593 | 0.041 | 0.626 | 0.038 |\n| **Hierarchical** | High | 1.0 | 0.601 | 0.048 | 0.541 | 0.055 | 0.574 | 0.052 |\n\n## Q2\n\nThank you for this methodological question.\n\n**Our Adaptive Threshold Selection Strategy**\n\nAs you noted, we indeed mentioned in Appendix B.2 (lines 1344-1348):\n\n> \"Therefore, in practice, the recommended strategy is to select a moderate intermediate value (e.g., τ = 0.15), and then make fine adjustments based on the specific characteristics of the data and prior knowledge in the relevant domain. We explicitly adopted this strategy in the numerical experiments...\"\n\n **τ = 0.15 serves as the theoretical starting recommendation**​, which is then adaptively adjusted based on specific data characteristics.\n\n**Long-Chain Design Validation in Synthetic Datasets**\n\nWe considered this in our dataset design. In lines 1524-1528, we described:\n\n> \"the chain path length ℓ, which determines the depth of the longest directed chain in the DAG... This parameterization allows us to generate graphs that reflect varying levels of structural complexity\"\n\nIn Appendix E.3 Table 3, we designed configurations with different chain lengths (ℓ):\n\n* **Short-chain datasets** (ℓ = 3): LowDim-L, LowDim-N, MidDim-S, HigDim-S\n* **Long-chain datasets** (ℓ = 6): LowDim-D, MidDim-D, HigDim-D, MidDim-C\n\n**Frequency Characteristics Analysis of Long vs. Short-Chain Mediation Pathways**\n\nLong-chain pathways naturally have lower occurrence frequencies in DAG ensembles because they require multiple consecutive edges to exist simultaneously. This design enables systematically validate optimal τ selection under different structural complexities.\n\n**Special Properties of the HPP Dataset**\n\nThe HPP dataset, as we described in lines 92-104 and Section 5.3, contains complex multi-step biomedical pathways. For such real-world biomedical data with long-chain characteristics, τ = 0.05 better captures important pathways that have clinical significance but relatively lower statistical frequency.\n\n**Medical Validation Support**\n\nWe validated these pathways using standard medical mediation analysis methods in Section 5.3 (lines 322-340) and Appendix F.4 (Tables 5-8), confirming their clinical validity.\n\n**Our Data-Driven Guidelines**\n\nBased on our research, we recommend:\n\n* ​**Short-chain and gene network data**​: τ = 0.15 provides balance between sensitivity and specificity\n* ​**Complex biomedical data with long-chain pathways**​: τ = 0.05 is more suitable for capturing meaningful mechanisms\n\n\n## Q3\n\nThank you for raising this important clarification question. \n\n**Specific Locations of Variable Mapping and Mediation Pathway Effects:**\n\nThe specific variable mapping and effect estimates for the pathway mentioned in line 324 have been provided in detail in the appendices:\n\n1. **Variable Mapping**: **Table 2** (Appendix D.2, line 1479, page 37) provides detailed descriptions of all body system variables\n2. **Specific Mediation Pathway Effects**: **Table 5** (Appendix F.4, line 1732, page 52) and **Table 6** (page 52) provide complete pathway effect estimates\n\n**Specific Correspondence of Variables Mentioned in Section 5.3:**\n\nAccording to the detailed definitions in Table 2 (Appendix D.2, line 1479, page 37) , the clinical terms mentioned in Section 5.3 correspond as follows:\n\n- **Sleep HRV index** → **hrv index** (Sleep HRV system): HRV index during sleep; lower HRV linked to stress and metabolic issues\n- **Pulse wave velocity** → **PWV mean** (Vascular Health system): Pulse wave velocity; indicates arterial stiffness and cardiovascular risk\n- **Carotid IMT** → **mean cimt** (Carotid Ultrasound IMT system): Average carotid IMT across measured regions\n\n  The complete mapping of all clinical terms to their corresponding body system variables can be found in Table 2.\n\n**Improvements in the Revision:**\n\n**We will add more specific appendix references in the main text**, ensuring that each pathway mentioned in Section 5.3 clearly references:\n\n1. The corresponding Table 2 variable definitions and their associated body systems\n2. The corresponding Table 5/6 effect estimation results\n3. Clear Figure 3 variable position identification\n\n## Q4&5\n\nWe sincerely thank the reviewer for identifying these editorial issues and will correct all in the revised manuscript. We commit to conducting a comprehensive review of the entire 53-page manuscript to ensure consistent acronym usage, capitalization, and technical term formatting throughout all sections."}]

Author rebuttal:

Meta-review: The authors study the problem of causal mediation analysis. In particular, they propose SIGMA, which combines causal discovery, mediation pathway identification, and effect estimation. All reviewers agree that this is an interesting and relevant problem, and that the proposed approach is novel and exciting. The rebuttal and active discussion that followed addressed the main concerns of the reviewers, leading them all to recommend acceptance. I'm happy to follow their recommendation, on the expectation that the authors include the clarifications and additional results in the camera-ready copy.
Decision: Accept (poster)