1 Introduction↩︎

Retinal imaging provides a unique, non-invasive window into systemic cardiovascular health, motivating its widespread use in modern AI-driven diagnostic systems [1]. Recent deep learning approaches have demonstrated strong predictive performance for cardiovascular risk estimation directly from retinal images [2]. However, these models are predominantly associative and often operate as black-box predictors, offering limited insight into the underlying biological mechanisms or how upstream physiological variation may be reflected in retinal microvascular structure and related to downstream vascular outcomes [3]. In particular, they provide limited explanation of how upstream physiological variation may be reflected in retinal microvascular structure and related to downstream vascular outcomes.

In this work, the retinal block is not interpreted as a direct biological cause of downstream cardiovascular or vascular outcomes. Instead, retinal traits are treated as an accessible phenotype layer that may statistically reflect upstream systemic vascular, metabolic, and inflammatory processes. The proposed mediator structure is therefore used as a hypothesis-testing strategy: it evaluates whether retinal microvascular features carry measurable mediator-like signal between systemic exposures and vascular outcomes under explicit modelling assumptions. This framing is consistent with the view of retinal microvascular signs as a non-invasive window onto systemic vascular dysfunction, rather than as isolated causal drivers of cardiovascular disease [4]–[6]. Beyond prediction, a central scientific question is whether retinal traits behave as passive correlates, downstream biomarkers, or phenotype-level indicators that may carry mediator-like statistical signal between systemic exposures and vascular outcomes.

This question is difficult to address in real-world biomedical settings because multimodal data are often fragmented. Clinical variables, molecular measurements, and retinal imaging are typically collected in independent cohorts and are rarely available jointly at the participant level [7]. For example, NHANES provides extensive clinical and physiological information [8], while retinal imaging datasets such as APTOS 2019 are collected separately [9]. As a result, causal modelling must operate under limited-resource conditions, requiring harmonisation of partially observed modalities rather than assuming access to a complete participant-linked multimodal dataset.

Classical causal discovery methods, including PC [10], NOTEARS [11], and LiNGAM [12], provide principled frameworks for learning directed relationships from observational data. However, in high-dimensional and noisy biomedical settings, these approaches can produce unstable or biologically implausible structures and may generalise poorly across heterogeneous data regimes [13]. This motivates causal modelling frameworks that incorporate domain structure while retaining interpretable pathway-level analysis.

To address this gap, we propose RetiSEM, a domain-constrained structural equation modelling framework for fragmented multimodal biomedical data. RetiSEM enforces an ordered variable structure across genetics, molecular traits, retinal features, and vascular outcomes, while restricting implausible edges through a forbidden-edge masking mechanism. The resulting pathway prior is illustrated in Fig. 1.

RetiSEM is designed not only for structure learning but also for pathway-level interpretation. By incorporating causal mediation analysis, the framework decomposes estimated effects into total, direct, and indirect components, enabling systematic evaluation of whether retinal variables behave as passive indicators, reflective biomarkers, or mediator-like statistical components under explicit modelling assumptions. The overall workflow is shown in Fig. 2.

We validate the framework in both synthetic and real-world fragmented settings. The synthetic benchmark suite contains ten scenarios with varying dimensionality, nonlinearity, causal depth, and parallel pathway structure, allowing assessment of whether domain constraints improve graph recovery as the data-generating process becomes more complex. In the real-world analysis, NHANES clinical variables are combined with externally derived retinal representations and proxy population-structure variables, reflecting a low-resource setting in which biomedical modalities are not fully participant-linked.

The main contribution of this work is a domain-constrained causal modelling framework for disjoint multimodal biomedical data. RetiSEM combines biologically ordered variable blocks, forbidden-edge masking, structural equation modelling, and causal mediation analysis to recover interpretable pathway structures under limited-resource conditions. Collectively, the experiments evaluate whether domain constraints improve causal graph recovery and whether retinal phenotypes carry passive, reflective, or mediator-like pathway signal in vascular outcome modelling.

2 Method↩︎

3pt

We propose a domain-constrained SEM framework for causal graph recovery and pathway-level interpretation in data-scarce multimodal biomedical settings. The central objective is to recover biologically consistent causal structure while characterising retinal variables as passive indicators, reflective biomarkers, or mediator-like statistical components within cardiovascular pathways. The overall causal hypothesis is illustrated in Fig. 1, where upstream genetic and molecular factors propagate through retinal features before influencing vascular outcomes.

2.1 Problem Setup and Domain Structure↩︎

Let \(X \in \mathbb{R}^{n \times d}\) denote the observed data matrix, where variables are partitioned into biologically meaningful blocks: \[(G, Z) \prec L \prec R \prec V,\] with \(G\) denoting genetic or ancestry-proxy features, \(Z\) covariates, \(L\) molecular traits, \(R\) retinal microvascular features, and \(V\) vascular outcomes. This ordering defines a domain-informed statistical pathway prior in which retinal variables form an interpretable phenotype layer between upstream biological context and vascular outcome modelling. Within RetiSEM, this structure is used to evaluate whether retinal microvascular features carry mediator-like statistical signal, reflective biomarker information, or only weak associative signal under the specified assumptions. The imposed ordering constrains the DAG search space, reduces biologically implausible reverse-direction edges, and improves interpretability in limited-resource biomedical settings. Compared with unconstrained DAG learning, this design uses prior biomedical knowledge to guide structure recovery while retaining RetiSEM as a structured hypothesis-testing framework.

2.2 Domain-Constrained Graph Recovery↩︎

We aim to estimate a directed acyclic graph (DAG) represented by a weighted adjacency matrix \(W \in \mathbb{R}^{d \times d}\). To prevent invalid relationships, we define a forbidden-edge mask \(F \in \{0,1\}^{d \times d}\) such that:

\[\hat{W} = W \odot (1 - F),\]

where \(F_{ij} = 1\) indicates disallowed edges (e.g., \(V \rightarrow L\) or \(R \rightarrow G\)). This enforces strict adherence to the domain ordering \((G, Z) \prec L \prec R \prec V\).

This constraint effectively reduces the hypothesis space, improving stability under limited data and preventing structurally invalid graph recovery. In particular, it ensures acyclicity by construction and avoids reverse-direction edges that frequently arise in unconstrained methods. The mask does not force every variable to participate in a single fixed path. Instead, it defines a permissible pathway space: edges consistent with the RetiSEM prior can be estimated when supported by the data, whereas reverse or biologically implausible directions are excluded. In this sense, the framework can follow the prior pathway organisation illustrated in Fig. 1, while still allowing the learned graph to select a sparse subset of supported directed relationships.

2.3 Structural Equation Modelling and Estimation↩︎

The pathway structure is modelled using SEM over the ordered domain blocks: \[L = f_L(G, Z) + \epsilon_L,\] \[R = f_R(L, G, Z) + \epsilon_R,\] \[V = f_V(R, L, G, Z) + \epsilon_V.\] Here, \(L\) denotes molecular or systemic exposure variables, \(R\) denotes retinal microvascular phenotype variables, \(V\) denotes vascular outcomes, and \((G,Z)\) denote genetic or ancestry-proxy variables and covariates. This formulation allows upstream biological context to influence molecular traits, retinal phenotypes, and vascular outcomes while preserving the domain-informed ordering introduced above. In practice, we adopt a linear SEM approximation as an interpretable and statistically stable estimator for causal graph recovery and pathway-level effect decomposition. Although biomedical relationships may involve nonlinear and interaction-driven mechanisms, the linear approximation provides a tractable first-order model for estimating total, direct, and indirect components across high-dimensional and noisy biomedical datasets. This choice also enables consistent comparison across simulated and real-data settings, while keeping the estimated pathways transparent and reproducible.

Under this approximation, each structural component is estimated using regression-based models subject to the domain-constrained adjacency mask. The learned coefficients are then used to quantify pathway-specific effects through the mediation formulation described below.

2.4 Retina-Hub Mediation Modeling↩︎

The model uses causal mediation analysis to quantify how much pathway-level signal passes through the retinal phenotype block. For an upstream exposure \(L\), retinal phenotype layer \(R\), vascular outcome \(V\), and adjustment variables \((G,Z)\), let \(R(l)\) denote the retinal phenotype value under exposure level \(l\), and let \(V(l,r)\) denote the vascular outcome under exposure level \(l\) and retinal state \(r\). We write \(V(l)=V(l,R(l))\) for the potential outcome when both the exposure and retinal phenotype take their natural values under \(l\). Comparing exposure levels \(l\) and \(l^{\ast}\), the total effect (TE), natural indirect effect (NIE), and natural direct effect (NDE) are defined as: \[\mathrm{TE}(l,l^{\ast}) = \mathbb{E}\{V(l)-V(l^{\ast})\},\] \[\mathrm{NIE}(l,l^{\ast}) = \mathbb{E}\{V(l,R(l))-V(l,R(l^{\ast}))\},\] \[\mathrm{NDE}(l,l^{\ast}) = \mathbb{E}\{V(l,R(l^{\ast}))-V(l^{\ast},R(l^{\ast}))\}.\]

Under consistency, positivity, and conditional exchangeability assumptions, these quantities can be identified as pathway-level summaries of direct and indirect association components [14], [15]. Let \(C=(G,Z)\) denote the adjustment set, and define \[\mu_V(l,r,c)=\mathbb{E}[V \mid L=l, R=r, C=c].\] The nested counterfactual mean is identified by the mediation functional: \[\begin{align} \mathbb{E}\{V(l,R(l^{\ast}))\} &= \int \int \mu_V(l,r,c)\, p(r \mid L=l^{\ast}, C=c)\, p(c)\,dr\,dc . \end{align}\] This gives the standard decomposition: \[\mathrm{TE}(l,l^{\ast}) = \mathrm{NDE}(l,l^{\ast}) + \mathrm{NIE}(l,l^{\ast}).\]

In the linear SEM implementation, these counterfactual quantities are operationalised through a product-of-coefficients estimator: \[R = \alpha L + \gamma^\top [G,Z] + \epsilon_R,\] \[V = c'L + \beta R + \delta^\top [G,Z] + \epsilon_V.\] Here, \(\alpha\) captures the exposure–retina association, \(\beta\) captures the retina–outcome association conditional on \(L\), \(G\), and \(Z\), and \(c'\) captures the direct exposure–outcome component. For a single retinal mediator, the pathway components are estimated as \[\mathrm{NIE}=\alpha\beta, \qquad \mathrm{NDE}=c', \qquad \mathrm{TE}=c' + \alpha\beta.\] For a retinal block with multiple features \(R_1,\ldots,R_m\), the indirect component is estimated by summing feature-level pathways: \[\mathrm{NIE}=\sum_{j=1}^{m}\alpha_j\beta_j, \qquad \mathrm{TE}=c' + \sum_{j=1}^{m}\alpha_j\beta_j.\] This decomposition allows RetiSEM to distinguish weak associative signals from more stable mediator-like pathway contributions based on the magnitudes and robustness of the estimated direct and indirect components.

2.5 Fragmented Multimodal Evaluation↩︎

To reflect real-world biomedical constraints, RetiSEM is evaluated under both controlled synthetic regimes and a fragmented applied setting. In the applied analysis, clinical and physiological variables are obtained from NHANES, retinal representations are derived from external fundus-image resources, and proxy genetic or ancestry-related variables are used for population-structure adjustment where individual genotype data are unavailable. Retinal representations are extracted using the AutoMorph [18] retinal image-analysis pipeline, producing quantitative microvascular descriptors such as vessel calibre, tortuosity, vessel density, and branching-related morphology, which define the retinal phenotype block \(R\).

The ten synthetic regimes are generated using the block-structured mechanism summarised in Table 1. Each regime follows the ordered pathway prior \((G,Z) \prec L \prec R \prec V\), while varying dimensionality, nonlinearity, causal depth, and parallel pathway structure to test the framework under progressively more complex generative conditions. Fig. 3 illustrates one representative LowDim-N pathway structure, showing how lipidomic, retinal, and vascular feature blocks are organised within the constrained mediation design.

This evaluation design tests whether the imposed pathway prior remains informative when moving from known synthetic structures to realistic data fragmentation. Although estimation uses a linear SEM approximation for interpretability and identifiability, the experiments assess whether the constrained pathway design remains useful across nonlinear synthetic regimes and partially observed multimodal data.

3 Experiments↩︎

We evaluated the proposed framework on both synthetic and real-world multimodal datasets to assess its ability to recover causal structure and characterize retinal-mediated pathways. Our experimental design focuses on three key aspects: (i) structural recovery accuracy, (ii) robustness under diverse generative conditions, and (iii) interpretation of retinal variables within cardiovascular pathways.

3.1 Synthetic Benchmark Evaluation↩︎

To systematically evaluate generalization, we construct a suite of ten synthetic benchmark scenarios (D1–D10) that vary along multiple axes, including dimensionality (low, medium, high), causal depth (shallow vs.deep), and functional form (linear vs.nonlinear). Each scenario follows the biologically motivated ordering \((G,Z) \prec L \prec R \prec V\), ensuring consistency with the domain assumptions of RetiSEM.

2pt

3.1.0.1 Synthetic Scenario Design.

The synthetic benchmarks are designed to systematically stress-test the proposed framework across increasing levels of complexity. As summarized in Table 1, the scenarios vary along three primary axes: (i) functional complexity through nonlinear transformations (\(\rho_{\text{nonlin}}\)), (ii) structural complexity via parallel paths (\(\kappa\)) and causal depth (\(\ell\)), and (iii) scalability through increasing dimensionality (\(p\)). This structured design ensures that the evaluation does not rely on idealized assumptions. In particular, the inclusion of both linear (D1) and nonlinear scenarios (D2–D10) allows us to assess whether RetiSEM remains effective when the underlying data-generating mechanisms deviate from linearity. These benchmarks are designed to test whether the proposed framework can reliably recover graph structure and pathway-level relationships under increasingly challenging conditions. Table 2 reports the full quantitative comparison across all ten synthetic benchmark scenarios. Across D1-D10, RetiSEM achieves the lowest SHD and the highest causal accuracy in every setting, indicating more accurate recovery of the underlying graph structure under both simple and challenging regimes. Some baseline methods achieve strong adjacency F1 in selected scenarios, but their performance is less stable across metrics and regimes. Collectively, RetiSEM provides the most consistent balance between low structural error, reliable orientation recovery, and high causal accuracy across the benchmark suite.

Beyond graph recovery, we examine whether the recovered structures preserve meaningful pathway-level effects. Figure 3 provides an illustrative example of the mediation pathway structure in the LowDim-N synthetic setting. Because the full mediation outputs across all ten scenarios are difficult to visualise in a single compact figure, this representative low-dimensional case shows how RetiSEM organises lipidomic, retinal, and vascular indicator blocks into a constrained pathway model. The displayed coefficients indicate loading and pathway strengths within this example, while the corresponding TE, NDE, and NIE decompositions across all benchmark scenarios are reported in Table 3.

3.2 Real-World Multimodal Analysis↩︎

To evaluate RetiSEM under realistic fragmentation, we combine NHANES clinical and physiological variables with externally derived retinal representations and proxy population-structure variables. This setting does not assume complete participant-level multimodal alignment; instead, it tests whether a domain-constrained causal framework can operate when biomedical modalities are partially observed and harmonised through feature definitions and covariate structure.

Across 96 exposure–retinal mediator–vascular outcome combinations, TE and NDE dominate most relationships, whereas NIE values are generally smaller but non-zero. Figure 4 ranks the strongest NHANES-based pathways by \(|\mathrm{NIE}|\) to identify where the retinal block contributes the largest mediator-like statistical signal. Each pathway follows the adjusted complete-case specification \(L \rightarrow R \rightarrow V \mid Z,G\) (CC), where \(L\) denotes an upstream exposure, \(R\) denotes a retinal mediator, and \(V\) denotes a vascular outcome.

For example, LBDSTRSI (refrigerated serum triglycerides, mmol/L) \(\rightarrow\) Vein_Distance_tortuosity \(\rightarrow\) BPXPLS can be interpreted as a lipid-related systemic exposure associated with a vascular examination outcome through a retinal tortuosity feature, after adjustment for \(Z\) and \(G\). This example illustrates how the model reports both the direct pathway component and the retinal indirect component within the same adjusted specification. The overall pattern in Fig. 4 shows that even among the highest-ranked indirect pathways, NIE remains smaller than TE and NDE in most cases. This supports the interpretation that retinal traits mainly behave as downstream vascular-reflective biomarkers in this cohort, while retaining limited but detectable mediator-like statistical contributions in selected pathways.

3.3 Analysis of Retinal Roles↩︎

The mediation results support the interpretation of retinal variables along a continuum rather than as fixed biological categories. Features with negligible indirect effects are best understood as weak associative indicators. Features with stable non-zero NIE values provide a mediator-like statistical signal under the specified pathway assumptions. Features with stronger direct associations with vascular outcomes may serve as reflective downstream biomarkers of systemic vascular state. Across both synthetic and real-data analyses, most retinal variables fall closer to the biomarker or weak-mediator end of this continuum, reinforcing the view that retinal morphology reflects systemic vascular variation more strongly than it acts as a dominant causal driver of downstream vascular outcomes. Overall, the experiments show that RetiSEM improves structural recovery in synthetic settings with increasing dimensionality, nonlinearity, and causal depth, while preserving interpretable TE/NDE/NIE decomposition. The real-data analysis further demonstrates how the framework can be applied under fragmented multimodal resources, where retinal traits provide mainly biomarker-like information with smaller but measurable indirect effects. These findings support the use of domain-constrained causal priors as a practical strategy for interpretable pathway modelling in limited-resource biomedical AI.

4 Discussion↩︎

The results show that domain-constrained structure learning can improve causal graph recovery in fragmented multimodal biomedical settings. Table 1 defines a benchmark suite that varies dimensionality, nonlinearity, causal depth, and parallel pathway structure, providing a stress test beyond a single idealised simulation. The quantitative results in Table 2 show that RetiSEM achieves the lowest SHD and highest causal accuracy across D1-D10. This suggests that biologically informed masks can reduce unstable or implausible edge orientations when data are limited, noisy, or generated under heterogeneous regimes. The mediation results further show that RetiSEM provides more than edge-level graph recovery. Fig. 3 illustrates how the LowDim-N setting is represented as a constrained pathway model, while Table 3 reports TE, NDE, and NIE decompositions across all synthetic scenarios. Together, these results support the methodological claim that domain priors and SEM-based decomposition can be combined to evaluate whether a retinal phenotype layer carries pathway-level statistical signal, even when the underlying generative process is not perfectly linear. In the real-world analysis, Figure 4 indicates that indirect effects through retinal variables are generally smaller than total and direct effects. This pattern prevents over-interpretation of the retina as a dominant causal mediator. Instead, the findings suggest that retinal microvascular features behave primarily as downstream biomarker-like indicators of systemic vascular variation, with selected pathways showing detectable mediator-like statistical signal. This interpretation is consistent with oculomics work that views retinal morphology as a non-invasive window into systemic circulation and cardiometabolic status [6].

Several limitations should guide interpretation. The real-data setting remains an observational, low-data regime, and the retinal representations are externally derived rather than fully participant-linked to all NHANES variables. The linear SEM estimator provides interpretability and stability but cannot capture all nonlinear or interaction-driven biomedical mechanisms. Estimated indirect effects should therefore be viewed as pathway-level statistical summaries under structural assumptions, not as definitive biological proof. Future work should evaluate RetiSEM in participant-linked multimodal cohorts, extend the estimator to nonlinear structural models, and test the stability of inferred pathways across independent biomedical datasets.

5 Conclusion↩︎

In this work, we introduced RetiSEM, a domain-constrained structural equation modelling framework for causal graph recovery and mediation analysis under limited biomedical resources. By combining ordered biomedical variable blocks, forbidden-edge masking, and TE/NDE/NIE effect decomposition, RetiSEM supports interpretable pathway modelling when complete multimodal data are unavailable. Synthetic experiments across ten benchmark scenarios demonstrate improved structural recovery and interpretable mediation estimates under varying dimensionality, nonlinearity, and causal depth. The real-world analysis further suggests that retinal features behave mainly as downstream biomarker-like indicators, with a smaller but detectable mediator-like pathway signal. Overall, RetiSEM provides a practical and reproducible framework for testing structured causal hypotheses with limited resources in multimodal health data.

References↩︎