arXiv CSopen access

Joint Treatment Effect Estimation from Incomplete Healthcare Data: Temporal Causal Normalizing Flows with LLM-driven Evolutionary MNAR Imputation

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

artificial intelligence, reasoning, knowledge representation

arXiv:2605.05125v1 [cs.LG] 6 May 2026

Joint Treatment Effect Estimation from Incomplete Healthcare Data: Temporal Causal Normalizing Flows with LLM-driven Evolutionary MNAR Imputation Olivia Jullian Parra∗ Department of Mathematical Modeling and Machine Learning University of Zürich Winterthurerstrasse 190, Zürich 8057 [email protected]

Sara Zoccheddu∗ Department of Mathematical Modeling and Machine Learning University of Zürich Winterthurerstrasse 190, Zürich 8057 [email protected]

David Catalan Cerezo Department of Mathematics ETH Zürich Rämistrasse 101, Zürich 8092 [email protected]

Tom Forzy Epidemiology, Biostatistics & Prevention Institute (EBPI), University of Zürich Hirschengraben 84, Zürich 8001 [email protected]

Franziska Ulrich Department of Mathematics ETH Zürich Rämistrasse 101, Zürich 8092 [email protected]

William Sutcliffe Department of Physics University of Zürich Winterthurerstrasse 190, Zürich 8057 [email protected]

Jakob Martin Burgstaller Institute of Primary Care University of Zürich Sonneggstrasse 6, Zürich 8091 [email protected]

Oliver Senn Institute of Primary Care University of Zürich Sonneggstrasse 6, Zürich 8091 [email protected]

Patrick Owen Department of Physics University of Zürich Winterthurerstrasse 190, Zürich 8057 [email protected]

Nicola Serra Department of Physics & Department of Mathematical Modeling and Machine Learning University of Zürich Winterthurerstrasse 190, Zürich 8057 [email protected]

Abstract Target trial emulation (TTE) provides a framework for answering causal questions using observational data when randomized controlled trials (RCTs) are infeasible. However, standard methods for treatment effect estimation have been developed in isolation, failing to jointly address the compounding challenges inherent to analyses of observational data such as electronic health records (EHRs). In particular, these challenges include time-varying confounding and missing-not-at-random (MNAR) missingness reaching 50%–80% for critical biomarkers. To address this gap, we propose a two-stage pipeline that jointly handles MNAR missingness and causal ∗ Equal contribution.

Preprint.

structure across time. CausalFlow-T, a Directed Acyclic Graph (DAG)-constrained normalizing flow with Long Short-Term Memory (LSTM)-encoded patient history, performs exact invertible counterfactual inference, eliminating the approximation errors and confounding biases where existing variational and adversarial methods fail silently. Ablations on four synthetic and one semi-synthetic dataset with known counterfactuals confirm that its two core design choices address strictly non-overlapping failure modes (DAG constraints for confounding separation, exact inference for structural propagation) with neither compensating for the absence of the other. To handle the incomplete data CausalFlow-T receives as input, we propose an LLM-driven evolutionary imputer and evaluate it with three LLM backends, including two open-source models. Across 30%–80% MNAR missingness, the imputer achieves the best pooled rank across biomarker and causal metrics, leading on point-wise accuracy and temporal extrapolation while maintaining average treatment effect (ATE) recovery where statistical methods progressively degrade. Applied to a cohort of adults with type 2 diabetes in Swiss primary care initiating a GLP-1 receptor agonist or SGLT-2 inhibitor, the pipeline recovers a per-protocol weight-loss difference of −0.98 kg [95% CI −1.01, −0.96] favoring GLP-1 receptor agonists, consistent with RCT evidence and estimated directly from realistically incomplete real-world EHR data.

Introduction

Type 2 diabetes (T2D) is one of the most common chronic conditions affecting over half a billion adults worldwide [International Diabetes Federation, 2025]. To optimize glycemic control and reduce the risk of complications (e.g. cardiovascular disease), effective management requires a multifactorial approach [American Diabetes Association, 2025]. Recent therapeutic advances, particularly two drug classes, GLP-1 receptor agonists (GLP-1RAs; including semaglutide, marketed as Ozempic) and sodium-glucose co-transporter-2 inhibitors (SGLT-2is), have transformed T2D management by demonstrating cardiovascular, renal, and weight loss benefits beyond glycemic control in randomized controlled trials (RCTs) [American Diabetes Association, 2025]. However, RCTs are limited by high costs, restricted sample sizes, and strict eligibility criteria, making them insufficient to address the complexity of T2D treatment, where therapies are numerous, combined, and evolve over time [Schneeweiss and Patorno, 2021]. Observational data (e.g., electronic health records [EHRs]) can complement RCTs by providing real-world evidence across diverse patient populations and care settings [Schneeweiss and Patorno, 2021]. Yet, causal inference from EHRs remains challenging due to confounding, selection bias, and other structural biases arising from non-randomized treatment assignment. Target trial emulation (TTE) addresses this by explicitly framing observational analyses to mimic a hypothetical RCT, thereby improving causal validity [Hernán and Robins, 2016]. Despite this framework, existing methods for treatment effect estimation have been developed in isolation, and none jointly addresses the structural challenges of real EHR data [Shalit et al., 2017, Louizos et al., 2017, Bica et al., 2020, D’Amour et al., 2021]. First, treatment decisions are driven by the patient’s evolving clinical profile, creating time-varying confounding that existing methods handle only under parametric assumptions or without explicit causal structure [Shalit et al., 2017, Louizos et al., 2017, Bica et al., 2020, D’Amour et al., 2021]. Second, EHR data exhibit substantial missingness (often exceeding > 50%) that arises under missing-not-at-random (MNAR) mechanisms (clinical decisions drive data collection) while imputation strategies have been typically validated at far lower missingness rates [Curnow et al., 2024, He et al., 2025, Mangussi et al., 2026]. These challenges interact: MNAR missingness distorts biomarker (e.g., blood pressure)–outcome relationships, compounding bias in downstream causal estimation. We address this gap with a two-stage pipeline that jointly handles time-varying confounding and MNAR missingness within the TTE framework. We apply this pipeline to estimate the per-protocol effect of GLP-1RAs versus SGLT-2is on 1-year body weight in a Swiss primary care cohort of adults with T2D, extending existing RCT evidence from specific agents in selected populations to the full drug classes used in real-world clinical practice. Our contributions are as follows: CausalFlow-T: A directed acyclic graph (DAG)-constrained normalizing flow (NF) conditioning a causal masked autoregressive flow (CausalMAF) on a LSTM-encoded patient history to model the joint distribution of covariates, treatment, and outcomes over time. Experiments on four synthetic 2

datasets and one semi-synthetic evaluation dataset with known counterfactuals reveal two distinct failure modes: (i) DAG constraints are necessary for confounding separation (unconstrained models maintain systematic bias in low-effect subgroups even despite low factual prediction error); (ii) exact inference is required for structural propagation (variational methods fail to recover the true treatment effect under complex mediation paths even with correctly specified causal graphs). LLM-driven evolutionary imputation: A large language model (LLM)-driven evolutionary pipeline benchmarked against LOCF [Fitzmaurice et al., 2012], MissForest [Stekhoven and Bühlmann, 2012], and a DAG-aware flow matching (CFM) baseline (inspired by the m-graph framework [Mohan and Pearl, 2021]) across 30–80% MNAR missingness using three LLM backends (GPT-5.4 and GPTOSS-120b from OpenAI, Qwen3.5-Plus from Alibaba Cloud). Evaluated via reconstruction, causal correlation structure preservation, and downstream ATE recovery, LLM-based methods dominate, with GPT-5.4 achieving the best overall performance and open-source models confirming robustness and reconstruction–causal trade-offs. Real-world application: When applied to a Swiss primary care T2D cohort initiating GLP-1RAs or SGLT-2is, our two-stage pipeline produces robust per-protocol estimates of 1-year body weight change (kg) from incomplete EHR data.

Related Work

Treatment effect estimation and time-varying confounding: Marginal structural models (MSMs) with inverse probability weighting (IPW) [Hernán and Robins, 2020] address time-varying confounding but rely on parametric assumptions that break down in high-dimensional, nonlinear settings [D’Amour et al., 2021]. Deep learning (DL) approaches such as R-MSN [Lim et al., 2018], CRN [Bica et al., 2020], and the Causal Transformer [Melnychuk et al., 2022] improve flexibility but rely on approximate inference or lack explicit causal structure, making them prone to spurious associations under interventions [Javaloy et al., 2023]. Generative models such as CEVAE [Louizos et al., 2017] and GANITE [Yoon et al., 2018] introduce additional approximation gaps via variational or adversarial objectives. All have been primarily evaluated on synthetic datasets with limited missingness and large sample sizes [Shalit et al., 2017], conditions rarely met in real-world clinical data [Sun et al., 2024, Barrett et al., 2020]. Normalizing flows for causal inference: [Khemakhem et al., 2021] establish identifiability via nonlinear independent component analysis for autoregressive flows; [Javaloy et al., 2023] extend this to general triangular mappings with a full do-operator formulation. [Chao et al., 2023] use diffusion models at the node level, but requiring a separate network per variable limits scalability for highdimensional EHR data. All existing methods assume static data. We extend DAG-constrained flows to temporal settings by conditioning on LSTM-encoded patient histories, capturing distributional shifts over time in a unified model for longitudinal treatment effect estimation. Missing data handling in causal inference: Missing data in EHRs can bias treatment effect estimation even after imputation [Zhou et al., 2023]. Effective imputation must leverage temporal disease trajectories and account for the underlying missingness mechanism [Zhou et al., 2023]. Standard methods such as Multiple Imputation by Chained Equations (MICE) [Sterne et al., 2009] and missing forest (MissForest) [Stekhoven and Bühlmann, 2012] assume missing-at-random (MAR) and can amplify bias under MNAR, while LOCF distorts temporal correlation structure by construction. Recent LLM-based approaches show improved reconstruction, but have been validated mainly at regimes of 5–40% missingness [He et al., 2025, Mangussi et al., 2026], below the high-rate MNAR missingness often encountered in real-world EHR data; in parallel, LLMs have been used as mutation operators in evolutionary program search for mathematical and algorithmic discovery [RomeraParedes et al., 2024, Liu et al., 2024, Novikov et al., 2025, Lange et al., 2025].

Methodology

We propose a two-stage pipeline for treatment effect estimation from incomplete longitudinal observational data. CausalFlow-T (Section 3.2) addresses time-varying confounding and exact counterfactual inference, while the LLM-driven evolutionary imputer (Section 3.3) handles MNAR missingness and feeds completed data into CausalFlow-T. Validating the pipeline sequentially isolates the contribution 3

of each stage, ensuring that downstream performance differences reflect imputation quality rather than estimator failure. 3.1

Problem formulation and assumptions

Consider N patients over T discrete time steps. At each time t, patient i has covariates xi,t ∈ Rd (indexed by j ∈ {1, . . . , d}), treatment ai,t ∈ {0, 1}, outcome yi,t ∈ R, and missingness mask mi,t ∈ {0, 1}d . Causal relationships between (x, a, y) are encoded in an expert-specified DAG G = (V, E). The incomplete longitudinal dataset is e = {(x̃i,t , ai,t , yi,t , mi,t )}N,T D i=1,t=1

where entry (i, t, j) is observed if and only if mi,t,j = 0. Missingness is MNAR: the probability of a value being unobserved depends on the unobserved value itself, P (mi,t,j = 1 | xi,t,j , xi,t,−j , ai,t , yi,t ) = fj xi,t,j , xi,t,−j , ai,t , yi,t , (1) where fj depends on xi,t,j itself (the MNAR component), as well as on other observed covariates xi,t,−j and on the outcome state yi,t , mirroring the clinical pattern where tests ordered on suspicion produce structured, not random, absence. Our goal is to estimate the ATE = E yi,t (1) − yi,t (0) , where yi,t (a′ ) denotes the potential outcome under do(A = a′ ) [Pearl, 2009], as well as individual treatment effects τi,t = yi,t (1) − yi,t (0). Assumptions. We require (i) sequential ignorability At ⊥⊥ Yt (a′ ) | ht , G, where ht is the LSTM-encoded patient history serving as the adjustment set; and (ii) correct DAG specification: G correctly encodes the causal structure of (x, a, y); misspecification propagates systematic bias through the abduction-action-prediction (AAP) procedure regardless of estimator quality (sensitivity analysis in Appendix I). Two-stage pipeline. We first validate CausalFlow-T on complete-data settings D with known counterfactuals, establishing its causal estimation properties independently of any imputation choices. Once validated, the LLM-driven evolutionary imputer produces D̂, enabling CausalFlow-T to operate on realistic MNAR data with following two-stage pipeline (Figure 1): ⋆

g CausalFlow-T e− D → D̂ −−−−−−−−−→ {ŷi,t (a′ ), τ̂i,t }

3.2

CausalFlow-T for Exact Counterfactual Inference

CausalFlow-T estimates ŷi,t (a′ ) and τ̂i,t via a DAG-constrained normalizing flow conditioned on longitudinal patient history, validated first on D̂ where ground-truth counterfactuals are known. DAG-Constrained Causal MAF. The core generative model is a Causal Masked Autoregressive Flow (CausalMAF) [Javaloy et al., 2023]. For variable vector vt = [xt , at , yt ] ∈ Rd the flow learns an invertible mapping fθ : vt 7→ zt with zt ∼ N (0, I), with autoregressive structure constrained to a topological sort of G: vj = fj−1 (zj ; vpa(j) , ht ), j = 1, . . . , d, (2) where pa(j) denotes the causal parents of node j in G. Intervening on at therefore propagates exclusively through causal descendants while non-descendants are held fixed, implementing P (Y | do(a′ ), X) directly. Without this constraint, an unconstrained flow can propagate interventions through arbitrary learned dependencies (including anti-causal directions) producing systematically biased counterfactuals even when the factual distribution is correctly fitted (Section 4.2). θ The flow maximizes the exact log-likelihood log p(vt ) = log pz (fθ (vt )) + log det ∂f ∂vt eliminating the ELBO approximation gap that accumulates under variational inference. CVAE and GNN-CVAE collapse to a hazard ratio (HR) ≈ 1.0 on CVD Risk despite a true 18% hazard reduction, consistent with posterior collapse under the Kullback-Leibler(KL) regularization (Section 4.2, Appendix K.1).

Temporal encoding via LSTM. The CausalMAF is conditioned on an LSTM encoder: vt ∼ CausalMAF( · | ht ),

ht = LSTM(xt , ht−1 ), H

(3)

The hidden state ht ∈ R summarizes the full observed covariate history up to time t, capturing disease progression, prior treatment exposure, and time-varying confounding. The DAG G governs the contemporaneous factorization P (Xt , At , Yt | ht ), enabling valid do-calculus [Pearl, 2009]. 4

Counterfactual inference via AAP. Given the trained flow, counterfactuals are computed via Pearl’s AAP procedure [Pearl, 2009]. Abduction: invert the flow exactly to recover the exogenous noise zt = fθ (vt ; ht ), capturing all individual-specific factors not explained by observed covariates. Action: replace at ← a′ ; because the autoregressive ordering respects G, the intervention propagates only to causal descendants of at ; zt is held fixed across both arms, implementing the surgical intervention of the do-operator. Prediction: propagate through fθ−1 via the DAG ordering to obtain vtcf and τ̂i,t = ŷi,t (1) − ŷi,t (0). Because flow inversion is exact, the same zt represents the individual-specific exogenous noise under both treatment arms, implementing the twin-network assumption and enabling individual-level counterfactuals, a guarantee CVAE-based approaches cannot provide (Appendix K.1). 3.3

LLM-driven evolutionary imputation under MNAR

e by evolutionary search over imputation operators. Crucially, the This stage constructs D̂ from D LLM does not fill missing entries directly; it proposes candidate Python imputers {g (k) }K k=0 , each a e 7→ D̂(k) (full algorithm in Algorithm 1, Appendix L). complete module g (k) : D Self-supervised proxy score. Let Ωobs = {(i, t, j) : mi,t,j = 0}. We draw a proxy holdout Ωp ⊂ Ωobs by masking each observed cell independently with probability ρ ∈ (0, 1) (contiguous-run e− . The composite proxy score penalizes both pointwise masking in real EHR; Appendix L), forming D reconstruction error and distortion of covariate–treatment and covariate–outcome correlations: v u u 1 s(g) = t |Ωp | |

2 e− ) − x̃ g(D +λY ∆Y (g) + λT ∆T (g), i,t,j i,t,j

(4)

(i,t,j)∈Ωp

}

RMSE(g)

where ∆Y,j (g) = |Corr(x̂j , y) − Corr(x̃j , y)| and ∆T,j (g) = |Corr(x̂j , a) − Corr(x̃j , a)| measure how much imputation distorts biomarker–outcome and biomarker–treatment correlations, with column averages ∆Y , ∆T (weights λY and λT reported in Table 6). Penalizing ∆Y directly targets the correlation structure that CausalFlow-T’s do-calculus depends on; penalizing ∆T targets the propensity-related correlations that drive confounding adjustment. Evolutionary loop. Starting from a deterministic seed imputer g (0) , at each iteration k the LLM ⋆ ⋆ receives the current-best candidate gk−1 , its score s(gk−1 ), and a compact search history summary, (k) then proposes a variation g . The update rule is (k) ⋆ g if s(g (k) ) < s(gk−1 ), ⋆ gk = (5) ⋆ gk−1 otherwise. ⋆ e to produce D̂. We use a single-parent evolutionary scheme, After K iterations, g ⋆ = gK is applied to D maintaining a single current-best imputer; implementation details are reported in Appendix L.1.

Experiments

We evaluate the pipeline in two steps: (i) CausalFlow-T on complete data with known counterfactuals, and (ii) four imputation strategies with CausalFlow-T fixed under 30%, 50%, and 80% MNAR missingness, including three LLM-based variants. Three findings emerge: (1) CausalFlow-T is the only model satisfying all reliability criteria, revealing two distinct failure modes (Section 4.2); (2) LLM-driven imputation achieves the best pooled performance across biomarker and causal metrics, with GPT-5.4 strongest overall and open-source LLMs demonstrating robustness and reconstruction–causal tradeoffs (Section 4.3); (3) the full pipeline recovers a robust per-protocol treatment effects from incomplete real-world EHR data (Section 5). Experiments were run on a single HPC node (Appendix L.3). 4.1

Experimental setup and evaluation protocol

Datasets. We evaluate CausalFlow-T on four synthetic datasets with ground-truth counterfactuals of increasing structural complexity: (1) Simple 3-Node (N =10k, T =5; linear, no confounding; 5

Stage 1: LLM-driven Evolutionary Imputation

Self-supervised score s(g (k) )

iterate k = 1, . . . , K

LLM proposes imputer g (k)

Stage 2: CausalFlow-T Expert-Based DAG G = (V, E)

Update gk⋆

Incomplete dataset e D

Final imputer g⋆

e g ⋆ (D)

Completed Patient Data D̂ (x1:T , a1:T , y1:T )

x̂1:T

Temporal Encoder (LSTM)

ht / τ t

Causal MAF (DAG-constrained)

AAP

Counterfactual ŷ(a′ )

(vt )

Figure 1: Stage 1 (LLM-driven Evolutionary Imputation): an LLM iteratively proposes candidate imputers g (k) , scores them via a self-supervised proxy s(g (k) ), and updates the running best gk⋆ ; after K rounds, g ⋆ produces D̂. Stage 2 (CausalFlow-T): a temporal encoder conditions a DAGconstrained Causal MAF on patient history; counterfactual outcomes ŷ(a′ ) are obtained via AAP.

positive control); (2) LDL Toy (N =10k, T =5; heterogeneous saturating effect, autoregressive AR(1) dynamics, selection-on-gain confounding; primary stress-test for subgroup calibration); (3) Cox Survival (N =30k, T =11; time-dependent hazard, bimodal covariate heterogeneity); and (4) CVD Risk Toy (N =50k, T =10; 17-node DAG, fully mediated treatment effect via systolic blood pressure (SBP) lags, HR ≈ 0.82; hardest benchmark). For imputation, we use a semi-synthetic dataset built on a real-world EHR backbone [Chmiel, 2011] with ten synthetic longitudinal biomarkers, a known benchmark ATE (−3.484), and MNAR missingness introduced at 30%, 50%, and 80%. All datasets are described in Appendix A. Causal inference baselines. We compare against CVAE, GNN-CVAE, and TARNet, all sharing the same LSTM temporal encoder; discriminative sequence models (Causal Transformer, R-MSN, CRN) are excluded by construction as they fail the three jointly necessary criteria for distributional counterfactual evaluation (Appendix B.1). Imputation baselines. Six imputation strategies are evaluated with CausalFlow-T held fixed: LOCF (carry-forward heuristic), MissForest (non-parametric random-forest), CausalCFM (a DAG-aware flow matching baseline, Appendix C), and LLM-driven evolutionary imputation with three backends (GPT-5.4, Qwen3.5-Plus, and GPT-OSS-120b). Evaluation protocols. CausalFlow-T is validated on four metrics exposing failure modes invisible to predictive error: subgroup calibration |BiasQ1 |/MAEQ1 , where MAE is the mean absolute error (ratio ≈ 1 signals systematic confounding failure), arm reconstruction error Erra , tail variance ratio VRQ4 (collapse < 1, overdispersion > 1), and hazard ratio recovery |HRtrue −HRpred |. Imputation is evaluated on two complementary layers: biomarker quality (MAE, root mean squared error RMSE, autocorrelation error |∆AC|, consecutive-step error CSE, and correlation inflation |infY |) d − ATE∗ |), with CausalFlow-T held fixed ¯ a , and |ATE and downstream causal quality (Q1 MAE, Err so differences reflect imputation alone (formal definitions in Appendix D). 4.2

CausalFlow-T: Two non-overlapping failure modes

Table 1 summarizes per-metric ranks and binary reliability criteria across all four synthetic datasets; full numerical results per dataset are in Appendix G (Tables 10–11). All models perform comparably on the positive control (Simple 3-Node), confirming the DAG constraint does not degrade performance when confounding is absent; two non-overlapping failure modes emerge on harder benchmarks. Failure mode 1: Confounding separation requires causal structure. On LDL Toy, GNN-CVAE, CVAE, and TARNet all reach |Bias/MAE|Q1 ≈ 1.0 (every low-effect error is systematic), while NF (no DAG) reduces this to 0.195 via exact inference alone and CausalFlow-T to 0.270, both clearing the 0.5 threshold; CausalFlow-T further achieves VRQ4 = 1.049 ± 0.033 (rank 1), the only model preserving tail variance without collapse or inflation. TARNet’s lowest absolute MAE (0.392 ± 0.081) masks a bias ratio of ≈ 1.0 across all quartiles and VRQ4 = 0.639, illustrating why MAE rank alone is an insufficient criterion. 6

Table 1: Per-metric ranks (1 = best, 5 = worst) across four synthetic datasets of increasing structural complexity, and binary reliability criteria evaluated across all datasets simultaneously. BiasQ1 : ratio |BiasQ1 |/MAEQ1 . VR: |VRQ4 − 1|. Arm: mean arm error. HR: |HRtrue − HRpred |. Reliability criteria (✓/✗/–): Bias <0.5; Tail = |VRQ4 − 1|<0.5 on all benchmarks; HR = correct direction; Arm = best on ≥2 benchmarks; Stable = no explosion or inversion [Austin and Stuart, 2015]. Simple Model CausalFlow-T NF (no DAG) GNN-CVAE CVAE TARNet

LDL

Cox

CVD Risk

MAE VR MAE Bias VR Arm MAE Arm HR MAE VR HR 3 2 1 3 4

1 2 3 5 4

2 3 5 4 1

2 1 4 3 5

1 2 3 4 5

3 4 5 2 1

3 2 4 1 5

1 2 5 4 3

2 4 3 1 5

1 2 3 3 5

2 5 3 4 1

1 2 3 4 5

Mean rank 1.83 2.58 3.58 3.17 3.67

Reliability Bias Tail HR Arm Stable ✓ ✓ ✗ ✗ ✗

✓ ✗ ✗ ✗ ✓

✓ ✓ ✗ ✗ ✗

✓ – ✗ ✗ –

✓ ✗ ✗ ✗ ✗

Failure mode 2: Structural propagation requires exact inference. CVD Risk stress-tests DAG constraints because the treatment effect is entirely mediated, so any model that cannot respect causal ordering will collapse or diverge. GNN-CVAE and CVAE encode the graph yet rely on ELBO approximation, and the resulting inference gap causes posterior collapse to HR ≈ 1.005 (null effect); NF (no DAG) performs exact inference but without structural constraint, routing the treatment signal through arbitrary pathways and yielding HR = 0.834 ± 0.324 (clinically meaningless); TARNet, lacking both, inverts the effect entirely (HR = 1.133 ± 0.148, predicting harm where there is benefit). Only CausalFlow-T, combining DAG factorization with exact normalizing-flow inference, forces the treatment signal through the correct mediation pathways and recovers HR = 0.786 ± 0.051 vs. true 0.831, with arm errors an order of magnitude below all competitors. On the other hand, the dequantization of binary survival outcomes (Cox benchmark) introduces a calibration cost that explains CVAE’s MAE advantage on that benchmark; on structurally meaningful metrics CausalFlowT leads with best arm-1 error (0.014 ± 0.002) and closest HR recovery (0.866 ± 0.011 vs. true 0.887); discrete flow extensions are a natural direction for future work (Section 6). Semi-synthetic validation. CausalFlow-T and GNN-CVAE are the only models passing the bias threshold (0.316 and 0.351); CausalFlow-T is further the only model with near-perfect variance preservation (VRQ4 = 1.006 ± 0.009) while GNN-CVAE collapses (VR = 0.616), confirming that low systematic error alone is insufficient for individual-level reliability. TARNet achieves the best arm reconstruction (0.081 ± 0.012 / 0.038 ± 0.003) but fails the bias threshold (0.708) and collapses variance (VR = 0.608), reinforcing that factual accuracy does not imply counterfactual reliability. Removing the DAG constraint degrades the bias ratio to 0.883 with the largest ATE deviation (−3.382 ± 0.251), underscoring that structural supervision is necessary under realistic clinical distributions. Full results can be found in Appendix F. CausalFlow-T for causal inference. CausalFlow-T uniquely resolves the bias–variance–calibration tradeoff across all four synthetic datasets (Table 1; mean rank 1.83), satisfying all five reliability criteria. Similar results are found when stress-testing in the real-world EHR semi-synthetic dataset. Competing models fail along at least one axis (trading bias for variance collapse, capturing trajectories with subgroup error, or achieving low MAE while reversing treatment direction) and none meets more than two criteria. We therefore fix CausalFlow-T in subsequent imputation experiments so that differences in ATE recovery reflect imputation quality alone. 4.3

Imputation under MNAR missingness

Full numerical results are in Appendix H (Tables 8–9) while Table 2 summarizes per-metric ranks across both evaluation layers. Biomarker reconstruction. GPT-5.4 achieves the lowest pointwise MAE and RMSE at every missingness level with only modest degradation from 30% to 80% (MAE 1.914→2.059), and the narrowest normalized error distribution across all missing positions (nMAE 0.55, 0.56, 0.64; Appendix H). Qwen3.5-Plus best preserves lag-1 autocorrelation at 30%, while MissForest is closest at 50% and 80%; CausalCFM best preserves biomarker (outcome correlation at 50% and 80% [|infY | ≈ 0.01]). GPT-OSS-120b tracks GPT-5.4 closely on pointwise reconstruction at 50% and 80%, supporting robustness of the LLM-driven strategy beyond a single model. 7

Table 2: Per-metric ranks (1 = best, 6 = worst) across the three missingness levels. Biomarker block: pointwise MAE, pointwise RMSE, lag-1 autocorrelation error |∆AC|, consecutive-step error (CSE), and absolute biomarker–outcome correlation inflation |inf Y |. Causal block: Q1 MAE, mean arm ¯ a , and absolute ATE residual |ATE [ − ATE∗ | (ATE∗ = −3.484). reconstruction error Err Biomarker quality Miss. Method

Downstream causal

¯ a |∆ATE| MAE RMSE |∆AC| CSE |inf Y | Q1 MAE Err

Mean rank

LOCF MissForest CausalCFM 30% GPT-5.4 Qwen3.5-Plus GPT-OSS-120b

6 4 5 1 3 2

6 2 5 4 1 3

6 4 5 1 3 2

6 5 2 4 1 3

5 6 1 4 3 2

4 6 1 3 5 2

4 6 2 1 5 3

5.38 4.62 3.25 2.38 3.00 2.38

LOCF MissForest CausalCFM 50% GPT-5.4 Qwen3.5-Plus GPT-OSS-120b

6 3 5 1 4 2

6 1 3 2 5 4

6 3 5 1 4 2

6 5 1 3 2 4

4 6 5 3 1 2

3 6 4 2 1 5

1 5 3 4 2 6

4.75 4.00 3.88 2.12 2.88 3.38

LOCF MissForest CausalCFM 80% GPT-5.4 Qwen3.5-Plus GPT-OSS-120b

6 3 5 1 4 2

6 1 2 5 4 3

5 3 6 1 4 2

6 2 1 5 3 4

3 1 6 4 5 2

1 4 6 3 5 2

2 4 6 1 5 3

4.38 2.62 4.62 2.62 4.25 2.50

Method

Mean rank

LOCF MissForest CausalCFM GPT-5.4 Qwen3.5-Plus GPT-OSS-120b

4.83 3.75 3.92 2.38 3.38 2.75

Downstream causal recovery and the reconstruction: causal tradeoff. Pointwise reconstruction quality and causal correlation preservation are not aligned objectives, and Table 2 exposes this tradeoff directly. At 30%, GPT-5.4 gives the closest ATE recovery (|∆ATE| = 0.024) while CausalCFM leads on subgroup and arm reconstruction; at 50%, Qwen3.5-Plus gives the strongest causal reconstruction ¯ a = 0.069) with near-best ATE residual (0.013 vs. 0.007 for LOCF); at 80%, GPT-5.4 again gives (Err the closest ATE recovery (|∆ATE| = 0.013). Above 50% missingness (the regime of most real EHR data [Sun et al., 2024]) LLM imputation is the only strategy preserving causal correlation structure and ATE recovery where statistical methods progressively degrade. LLM-driven evolutionary imputation. LLM-driven imputation achieves the strongest pooled performance across 30–80% MNAR (mean rank: GPT-5.4 2.38, GPT-OSS-120b 2.75, Qwen3.5-Plus 3.38; MissForest 3.75, CausalCFM 3.92, LOCF 4.83). GPT-5.4 is selected as the default; CausalCFM remains a useful comparator for causal-fidelity diagnostics at lower missingness. Appendix Figure 7 shows GPT-5.4 search progress, and Appendix Tables 8–9 compare the final imputer with the first executable GPT-5.4 proposal; bootstrap sensitivity is shown in Appendix Table 13).

Real-world application: Causal inference from healthcare data

TTE in incomplete EHR data: Semi-synthetic benchmarks provide controlled settings with known ATEs, enabling direct evaluation of estimator bias against ground-truth counterfactuals. While our pipeline recovers the ATE under MNAR missingness in this setting, real-world EHR data introduces additional challenges: unlike RCTs, it is collected for clinical and administrative purposes rather than to answer specific causal questions; treatment is not randomized, follow-up is irregular, treatment deviations are informative, and missingness is rarely MAR. To bridge this gap, we apply the two-stage pipeline (LLM imputation followed by CausalFlow-T) to the FIRE database, the most comprehensive primary care EHR database in Switzerland [Chmiel, 2011] (Appendix A.6). Using the TTE framework [Hernán and Robins, 2016], we estimate the per-protocol effect among initiators of GLP-1RAs (n = 2,392) versus SGLT-2is (n = 3,722) on body weight loss (kg) at 1-year in adults with T2D (details in Appendix A.6.2). Causal assumptions are encoded in an expert-based DAG and missing measurements are imputed via our validated LLM-based strategy (Appendix A.6.2). Due to privacy and regulatory constraints, the dataset used in this study cannot be publicly released. 8

Per-protocol effect estimation and external consistency: Figure 2 shows both treatment arms exhibit monotonic weight loss from initiation to 1-year, estimated using the proposed pipeline d GLP-1RA − 4.02 kg [95% CI −4.15, −3.88] and (Appendix A.6.2). At 1-year, the model recovers ∆W d SGLT-2i − 3.04 kg [95% CI −3.15, −2.91], yielding a per-protocol effect of −0.98 kg [95% CI ∆W −1.01, −0.96], favoring GLP-1RAs (bootstrap uncertainty analysis in Appendix Table 13). This estimate is consistent with the SUSTAIN 8 trial (an RCT comparing a GLP-1RA [semaglutide] vs (an SGLT-2i [canagliflozin]) in T2D, reporting 1-year body weight difference of −1.06 kg (95% CI −1.76, −0.36) [Lingvay et al., 2019]. Our study extends RCT evidence from specific agents in selected trial participants to the full GLP-1RA and SGLT-2i classes used in routine Swiss primary care, where mixed agents and > 50% MNAR missingness are expected to attenuate the contrast compared to the SUSTAIN 8 trial estimate [Lingvay et al., 2019]. While ground-truth counterfactuals are unavailable, consistency with RCT evidence supports the plausibility of the per-protocol effects estimated by the proposed pipeline from incomplete real-world EHR data.

Figure 2: Weight change over 1-year for adults with type 2 diabetes initiating GLP-1 receptor agonists vs SGLT-2 inhibitors, estimated from real-world healthcare data using the two-stage pipeline.

Discussion

This work addresses a key gap in causal inference from observational data: jointly handling timevarying confounding and high-rate MNAR missingness, two compounding challenges that may substantially distort treatment effect estimation in observational data. We propose a two-stage pipeline in which CausalFlow-T, a DAG-constrained normalizing flow with exact counterfactual inference, is validated on complete data before an LLM-driven evolutionary imputer enables deployment under realistic missingness, isolating each component’s performance impact. A central finding is not that CausalFlow-T outperforms competing models, but that DAG constraints and exact inference address strictly non-overlapping failure modes, neither compensates for the absence of the other. Existing methods evaluated on simpler benchmarks may be systematically overestimating their robustness: models that perform well on factual prediction can still fail on structurally meaningful criteria such as subgroup calibration, tail variance, and mediation-consistent effect propagation. Our five-criterion reliability framework provides a more stringent evaluation template than any single metric and may serve as a useful benchmark for future work in longitudinal causal inference. The imputation results reveal a complementary insight: pointwise reconstruction and causal correlation preservation are not aligned objectives. LLM-based methods, especially GPT-5.4, dominate on reconstruction and temporal extrapolation, while CausalCFM performs better on selected causalfidelity diagnostics, especially biomarker–outcome correlation preservation. This tradeoff is invisible to standard imputation benchmarks, motivating joint biomarker–causal evaluation as a necessary standard for causal pipelines under missingness. Future methods could explicitly optimize this tradeoff instead of separating reconstruction from causal fidelity. Applied to 6,114 adults with T2D in Swiss primary care, the pipeline recovers a per-protocol effect of −0.98 kg [95% CI −1.01, −0.96] favoring GLP-1RAs over SGLT-2is, consistent with RCT evidence [Lingvay et al., 2019]. Although ground-truth counterfactuals are unavailable, this suggests that the proposed pipeline recovers clinically plausible treatment effect estimates from incomplete EHR data, extending causal evidence beyond selected participants and conditions typical for RCTs. Several limitations suggest directions for future work. First, CausalFlow-T uses dequantization for binary outcomes, incurring calibration costs on survival endpoints that explain CVAE’s lower 9

Cox Survival MAE despite poorer hazard-ratio recovery and arm reconstruction (Appendix K); discrete or hybrid flow-based architectures for survival outcomes are a natural extension. Second, the DAG is assumed fixed and expert-specified. Clinically motivated single-edge removals reveal a bias–instability tradeoff: deleting unique causal pathways primarily increases bias, whereas deleting redundant pathways increases estimation instability, suggesting that cross-seed instability may provide a differentiable signal for data-driven DAG refinement without ground-truth counterfactuals (Appendix I). However, edge-direction errors and latent confounding remain important limitations for future work. Third, although we evaluate three LLM backends, all use the same evolutionary search protocol, and performance may therefore depend on the search design. Future work should explore alternative evolutionary strategies, including recombination and multi-objective acceptance rules. Finally, like all observational studies, our real-world data analysis assumes conditional exchangeability given the expert-specified adjustment set, a standard but untestable assumption. Future work could incorporate sensitivity analyses for unmeasured confounding.

Conclusion

We introduced a two-stage pipeline for treatment effect estimation from incomplete longitudinal EHR data, combining CausalFlow-T (a DAG-constrained normalizing flow for exact counterfactual inference) with an LLM-driven evolutionary imputer for MNAR missingness. Controlled ablations show that DAG constraints and exact inference address distinct failure modes in longitudinal causal inference, with CausalFlow-T the only model satisfying all five reliability criteria simultaneously. Joint biomarker–causal evaluation further reveals causal implications of imputation choices that reconstruction-only metrics overlook. Applied to real-world EHR data, the pipeline recovers a perprotocol effect consistent with RCT evidence despite substantial MNAR missingness. Our findings demonstrate that jointly addressing causal structure and MNAR missingness is both necessary and feasible, supporting ML-based treatment effect estimation from observational data as a practical complement to RCTs and extending causal inference to real-world populations and care settings where trials are infeasible or insufficient and real-world evidence is needed to inform clinical and policy decision-making.

References American Diabetes Association. 10. Cardiovascular Disease and Risk Management: Standards of Care in Diabetes—2025. Diabetes Care, 48(Supplement_1):S207–S238, 1 2025. ISSN 0149-5992. doi: 10.2337/DC25-S010. URL https://dx.doi.org/10.2337/dc25-S010. Peter C Austin and Elizabeth A Stuart. Moving towards best practice when using inverse probability of treatment weighting (iptw) using the propensity score to estimate causal treatment effects in observational studies. Statistics in medicine, 34(28):3661–3679, 2015. Donika Balaj, Jakob M Burgstaller, Audrey Wallnöfer, Katja Weiss, Oliver Senn, Thomas Rosemann, Thomas Grischott, Stefan Markun, FIRE research group, et al. Leveraging free-text diagnoses to identify patients with diabetes mellitus, obesity or dyslipidaemia–a cross-sectional study in a large swiss primary care database. Swiss Medical Weekly, 155(2):3360–3360, 2025. James E Barrett, Aylin Cakiroglu, Catey Bunce, Anoop Shah, and Spiros Denaxas. Selective recruitment designs for improving observational studies using electronic health records. Statistics in Medicine, 39(19):2556–2567, 2020. Ioana Bica, Ahmed M. Alaa, James Jordon, and Mihaela van der Schaar. Estimating counterfactual treatment outcomes over time through adversarially balanced representations. In International Conference on Learning Representations (ICLR), 2020. Patrick Chao, Patrick Blöbaum, Sapan Patel, and Shiva Prasad Kasiviswanathan. Modeling causal mechanisms with diffusion models for interventional and counterfactual queries. arXiv preprint arXiv:2302.00860, 2023. Moshinsky Chmiel. The fire project. Swiss medical weekly, 141(0304):w13142–w13142, 2011. Elinor Curnow, Rosie P Cornish, Jon E Heron, James R Carpenter, and Kate Tilling. Multiple imputation using auxiliary imputation variables that only predict missingness can increase bias due to data missing not at random. BMC medical research methodology, 24(1):231, 2024. Alexander D’Amour, Peng Ding, Avi Feller, Lihua Lei, and Jasjeet Sekhon. Overlap in observational studies with high-dimensional covariates. Journal of Econometrics, 221(2):644–654, 2021. Garrett M Fitzmaurice, Nan M Laird, and James H Ware. Applied longitudinal analysis. John Wiley & Sons, 2012. Xinrui He, Yikun Ban, Jiaru Zou, Tianxin Wei, Curtiss Cook, and Jingrui He. Llm-forest: Ensemble learning of llms with graph-augmented prompts for data imputation. In Findings of the Association for Computational Linguistics: ACL 2025, pages 6921–6936, 2025. Miguel A Hernán and James M Robins. Practice of Epidemiology Using Big Data to Emulate a Target Trial When a Randomized Trial Is Not Available. American Journal of Epidemiology, 183 (8):758–764, April 2016. doi: 10.1093/aje/kwv254. URL https://academic.oup.com/aje/ article/183/8/758/1739860. Miguel A. Hernán and James M. Robins. Causal inference: What if. Boca Raton: Chapman & Hall/CRC, 2020. Miguel A Hernán, Wei Wang, and David E Leaf. Target trial emulation: a framework for causal inference from observational data. Jama, 328(24):2446–2447, 2022. International Diabetes Federation. IDF Diabetes Atlas, 11th edition, 2025. URL https://diabetesatlas.org/media/uploads/sites/3/2025/04/IDF_Atlas_11th_ Edition_2025.pdf. Accessed 9 September 2025. Adrián Javaloy, Pablo Sánchez-Martín, and Isabel Valera. Causal normalizing flows: from theory to practice. In Advances in Neural Information Processing Systems (NeurIPS), 2023. Ilyes Khemakhem, Ricardo Pio Monti, Robert Leech, and Aapo Hyvärinen. Causal autoregressive flows. In International Conference on Artificial Intelligence and Statistics (AISTATS), pages 3520–3528, 2021. 11

Robert Tjarko Lange, Yuki Imajuku, and Edoardo Cetin. ShinkaEvolve: Towards open-ended and sample-efficient program evolution. arXiv preprint arXiv:2509.19349, 2025. Bryan Lim, Ahmed M. Alaa, and Mihaela van der Schaar. Forecasting treatment responses over time using recurrent marginal structural networks. In Advances in Neural Information Processing Systems (NeurIPS), volume 31, 2018. Ildiko Lingvay, Andrei-Mircea Catarig, Juan P Frias, Harish Kumar, Nanna L Lausvig, Carel W le Roux, Desirée Thielke, Adie Viljoen, and Rory J McCrimmon. Efficacy and safety of once-weekly semaglutide versus once-daily canagliflozin as add-on to metformin in patients with type 2 diabetes (SUSTAIN 8): a double-blind, phase 3b, randomised controlled trial. The Lancet Diabetes & Endocrinology, 7(11):834–844, 2019. doi: 10.1016/S2213-8587(19)30311-0. Fei Liu, Xialiang Tong, Mingxuan Yuan, Xi Lin, Fu Luo, Zhenkun Wang, Zhichao Lu, and Qingfu Zhang. Evolution of heuristics: Towards efficient automatic algorithm design using large language model. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 32201–32223. PMLR, 2024. Christos Louizos, Uri Shalit, Joris M. Mooij, David Sontag, Richard Zemel, and Max Welling. Causal effect inference with deep latent-variable models. In Advances in Neural Information Processing Systems (NeurIPS), volume 30, 2017. Arthur Dantas Mangussi, Ricardo Cardoso Pereira, Ana Carolina Lorena, and Pedro Henriques Abreu. Large language models for missing data imputation: Understanding behavior, hallucination effects, and control mechanisms. arXiv preprint arXiv:2603.22332, 2026. Valentyn Melnychuk, Dennis Frauen, and Stefan Feuerriegel. Causal transformer for estimating counterfactual outcomes. In International Conference on Machine Learning (ICML), pages 15293–15329, 2022. Karthika Mohan and Judea Pearl. Graphical models for processing missing data. Journal of the American Statistical Association, 116(534):1023–1037, 2021. Alexander Novikov, Ngân Vũ, Marvin Eisenberger, Emilien Dupont, Po-Sen Huang, Adam Zsolt Wagner, Sergey Shirobokov, Borislav Kozlovskii, Francisco J. R. Ruiz, Abbas Mehrabian, M. Pawan Kumar, Abigail See, Swarat Chaudhuri, George Holland, Alex Davies, Sebastian Nowozin, Pushmeet Kohli, and Matej Balog. Alphaevolve: A coding agent for scientific and algorithmic discovery. arXiv preprint arXiv:2506.13131, 2025. Judea Pearl. Causality. Cambridge university press, 2009. Bernardino Romera-Paredes, Mohammadamin Barekatain, Alexander Novikov, Matej Balog, M. Pawan Kumar, Emilien Dupont, Francisco J. R. Ruiz, Jordan S. Ellenberg, Pengming Wang, Omar Fawzi, Pushmeet Kohli, and Alhussein Fawzi. Mathematical discoveries from program search with large language models. Nature, 625(7995):468–475, 2024. doi: 10.1038/s41586-023-06924-6. Sebastian Schneeweiss and Elisabetta Patorno. Conducting Real-world Evidence Studies on the Clinical Outcomes of Diabetes Treatments. Endocrine reviews, 42(5):658–690, 10 2021. ISSN 1945-7189. doi: 10.1210/ENDREV/BNAB007. URL https://pubmed.ncbi.nlm.nih.gov/ 33710268/. Uri Shalit, Fredrik D Johansson, and David Sontag. Estimating individual treatment effect: Generalization bounds and algorithms. In International Conference on Machine Learning (ICML), 2017. Daniel J Stekhoven and Peter Bühlmann. Missforest—non-parametric missing value imputation for mixed-type data. Bioinformatics, 28(1):112–118, 2012. Jonathan AC Sterne, Ian R White, John B Carlin, Michael Spratt, Patrick Royston, Michael G Kenward, Angela M Wood, and James R Carpenter. Multiple imputation for missing data in epidemiological and clinical research: potential and pitfalls. Bmj, 338, 2009. 12

Minghui Sun, Matthew M Engelhard, Armando D Bedoya, and Benjamin A Goldstein. Incorporating informatively collected laboratory data from ehr in clinical prediction models. BMC Medical Informatics and Decision Making, 24(1):206, 2024. 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. Yizhao Zhou, Jiasheng Shi, Ronen Stein, Xiaokang Liu, Robert N Baldassano, Christopher B Forrest, Yong Chen, and Jing Huang. Missing data matter: an empirical evaluation of the impacts of missing ehr data in comparative effectiveness research. Journal of the American Medical Informatics Association, 30(7):1246–1256, 2023.

Dataset Descriptions

All four benchmarks are fully synthetic, so both ground-truth outcomes Y (0) and Y (1) are available for every test patient, enabling exact evaluation of counterfactual predictions. Train/test splits of 80/20 are used throughout. A.1

Simple 3-Node

The dataset contains N patients observed over T timesteps (default N =10,000, T =5). Three variables are simulated: (t)

• Confounder X1 ∼ N (t/4, 1): a single continuous covariate whose mean drifts linearly with time, inducing non-stationarity. • Treatment T ∈ {0, 1}: a binary treatment assigned once at baseline and kept constant. The assignment probability is a nonlinear function of X1 : 0.2 if X̄1 > 2.5, P (T = 1 | X1 ) = 0.8 otherwise, (t)

where X̄1 = Et [X1 ] and the intermediate variable X12 − sin(X1 ) + ε (ε ∼ N (0, 0.25)) mediates the propensity. This creates nonlinear confounding with a treatment probability that flips between 20% and 80%. • Outcome Y (t) : a continuous outcome generated as ( (t) 3X1 + 0.25 (t/T ) + ε0 (t) Y = (t) 3X1 − 0.50 (t/T ) + ε1

T = 0, T = 1,

where ε0 , ε1 ∼ N (0, 1). Treatment reduces the time trend by 0.75 (t/T ) relative to control, inducing a modest, time-growing ITE. Causal graph. The DAG is X1 → T → Y , X1 → Y , with self-loops at each timestep, encoded as a 3×3 adjacency matrix. Purpose. A positive control with a low confounded structure (linear dependencies mostly). Because the causal ordering matches the variable index order and confounding is low, all five models are expected to perform comparably. The dataset verifies that the framework performs well even in the simplest case. A.2

LDL Toy

Generated dataset with N =10,000 patients, T =5 timesteps. The outcome Y represents a continuous biomarker (analogous to LDL cholesterol) subject to a heterogeneous drug response, log-saturation pharmacokinetics, autoregressive dynamics, and selection-on-gain confounding. Specifically, the response magnitude at baseline is δi ∼ Uniform(0.20, 0.60), giving each patient an individual treatment-effect scale (20%–60% reduction). The dose-response saturates over time according to log(1 + λt) sat(t) = , λ = 0.4, log(1 + λTmax ) so early timesteps carry most of the treatment signal. Temporal dependence is introduced via a first-order AR(1) process with coefficient ρ = 0.75, such that each observation exhibits strong dependence on its immediate predecessor. Confounding is selection on gain: patients who respond more strongly (δi above median) are preferentially assigned to treatment, so the treated group has a higher true effect than the control group. Purpose. Tests calibration of heterogeneous treatment effects. The combination of log-saturation (non-constant effect over time), AR(1) dynamics (temporal correlation), and selection-on-gain confounding makes this the primary benchmark for the structural-failure and variance-collapse diagnostics. 14

A.3

Cox Survival Toy

The dataset contains N =30,000 patients, T =11 timesteps (including t=0), p=5 covariates. 1. Covariates. Five covariates are drawn from bimodal Gaussian mixtures: ( (j) N (µ1 , 0.01) Zij = 1, Xij | Zij = (j) N (µ2 , 0.04) Zij = 0, (1)

(1)

with Zij ∼ Bernoulli(0.5) and (µ1 , . . .) = (3, 1, 2, 5, 0), (µ2 , . . .) = (0, 2, 4, 5, 5). 2. Hazard function. Per-patient log hazard ratios are log HRi = β ⊤ Xi ,

β = (0.2, −0.4, −0.3, 0.3, −0.5),

giving HRi = exp(β ⊤ Xi ). The baseline hazard is h0 = 0.5. 3. Treatment assignment (confounding). Patients with HRi ≥ 0.5 (i.e., higher risk) are treated with probability 0.70; the remainder with probability 0.30. 4. Treatment effect. Treatment reduces the individual hazard ratio by 30%: HRtreated = i 0.70 HRi . 5. Event times. Survival times are drawn from an exponential distribution, E(1/[h0 HRi ]), and converted to a binary event indicator at each of the 11 timesteps. Censoring is applied only at the end of the study (all weight is on the final timestep), so within-study censoring is minimal. 6. Counterfactuals. The smooth individual cumulative incidence function Fi (t) = 1 − exp(−h0 HRi t) is computed for both the always-treated and never-treated hazards and stored as ground-truth potential outcomes for evaluation. Purpose. Tests survival-model calibration under continuous confounding. The bimodal covariate structure creates two patient clusters. The correct arm reconstruction requires the model to disentangle hazard ratio heterogeneity from the selection-on-risk confounding. A.4

CVD Risk Toy

The CVD Risk Toy is the most complex benchmark, designed to emulate a real-world antihypertensive target-trial setting. It is semi-synthetic: the causal structure, covariate relationships, treatment dynamics, and outcome model are grounded in established cardiovascular epidemiology and RCT evidence, while outcomes are re-generated under a fully specified data-generating process to provide a known ground truth (including exact counterfactual trajectories and a precise hazard ratio (HR ≈ 0.82)). The dataset comprises N =50,000 patients over T =10 observation years, with 15 time-varying covariates and a 17-node DAG encoding the full mediation structure. 1. Baseline demographics. Age ∼ N (55, 100); BMI from age with Gaussian noise; total cholesterol (TC) from a log-normal with age/BMI dependence; high-density lipoprotein (HDL) inversely related to age and BMI; smoking status ∼ Bernoulli(0.6); diabetes probability proportional to BMI; race and sex drawn from fixed population proportions. 2. Time-varying dynamics. At each year t, age increments by 1; BMI, TC, and HDL evolve with Gaussian noise; hypertension (htn) transitions to 1 irreversibly with a probability driven by age, TC, and smoking; diabetes is absorbing once acquired. Systolic blood pressure (SBP) at time t is: SBP(t) = 70

SBP(t−1) + 0.5 age(t) + 0.15 TC(t) + 10 smoker + ε, ε ∼ N (0, 400), SBP0

clipped to [80, 200]. 3. Treatment assignment (absorbing). At t=0, patients with hypertension are treated with probability 0.70; non-hypertensive patients with probability 0.30. Treatment is absorbing: once assigned it persists. At subsequent time steps, untreated hypertensive patients initiate treatment with probability 0.01 per year. 15

4. Mediated treatment effect. The antihypertensive reduces SBP via lagged mediation: (t) SBPfinal = SBP(t) 1 − [0.20 T (t−1) − 0.05 T (t−2) − 0.02 T (t−3) ] , with lags 1-3. This is the only pathway through which treatment affects the outcome; the complete causal chain is T → SBPfinal → Y . 5. CVD outcome. Binary CVD event at each timestep from a logistic model with 12 risk-factor coefficients: p log = −10 + 0.005 age + 0.15 sex + 0.03 BMI + 0.015 SBPfinal 1−p − 0.01 HDL + 0.01 TC + 0.25 smoker + 0.30 diabetes + 0.20 fam_hx + 0.10 ⊮[race = Black] − 0.05 ⊮[race = Asian]. The outcome is absorbing: rows after the first event are excluded. Both always-treated and never-treated probability trajectories are stored as counterfactual targets. 17-node DAG. The causal graph encodes the full mediation structure: T → SBPfinal → Y , plus direct paths from all 15 covariates to Y and temporal self-edges. Purpose. Strong confounding (hypertensive patients preferentially treated), a fully mediated treatment pathway (SBP lags 1-3), 15 correlated time-varying covariates, and an absorbing binary outcome combine to make this the hardest benchmark. Correctly recovering the true HR ≈ 0.82 requires disentangling the mediation chain from the confounded observational distribution (precisely the setting where the expert DAG prior is expected to matter most). A.5

Summary and main purpose

The four benchmarks vary deliberately in complexity. Simple 3-Node is a positive control (minimal confounding, three variables, known closed-form ITE) and is used to verify that all the models perform well in a simple, unrealistic case. Cox Survival Toy introduces a survival endpoint with bimodal covariate heterogeneity and selectionon-risk confounding, testing whether models correctly recover per-arm cumulative incidence functions under censoring. LDL Toy and CVD Risk Toy are the two hardest benchmarks and carry the most diagnostic weight in our evaluation. LDL Toy combines log-saturation pharmacokinetics, AR(1) temporal dynamics, and selection-on-gain confounding, making it the primary stress-test for subgroup calibration: a model that merely minimizes factual mean squared error (MSE) will exhibit structural bias (bias/MAE ≈ 1) and variance collapse (VRQ4 < 1) under this design. On the other hand, CVD Risk Toy is the most demanding benchmark overall: 15 time-varying covariates, an absorbing binary outcome, a purely mediated treatment effect through SBP (lags 1–3), and a 17-node expert DAG create a setting where correctly recovering the interventional distribution p(Y | do(T )) requires explicit encoding of the causal graph. It is the only benchmark where direction d < 0 or > 1) is observed for models lacking structural constraints, making it the inversion (HR definitive test of causal identifiability. Together, this benchmark aims to cover the three main failure modes of longitudinal causal models (confounding, variance miscalibration, and mediation) while providing exact ground-truth potential outcomes for all evaluation metrics. A.6

Real-World EHR Database Description and Use

The Family Medicine Research using Electronic Medical Records (FIRE) database (managed by the Institute of Primary Care at the University of Zurich, at the University Hospital Zurich, Switzerland) is the most comprehensive Swiss primary-care database of longitudinal (time-stamped) EHR from general practitioners across Switzerland. [Balaj et al., 2025] The FIRE database contains individuallevel information on demographics, consultations, medication prescriptions, laboratory values, vital signs, and diagnoses [Balaj et al., 2025, Chmiel, 2011]. FIRE is used in two ways throughout this work: as a clinical backbone for generating the semi-synthetic evaluation dataset (Section A.6.1), and to create the study population for the real-world application (Section A.6.2). 16

A.6.1

Semi-Synthetic MNAR Dataset

The semi-synthetic evalution dataset is built on a patient-month panel structure from FIRE, keeping all observed real-world information on demographics, comorbidities, comedication, and calendar time. Synthetic covariates: Ten continuous longitudinal variables are generated recursively on top of this backbone, indexed by patient and by month. The variables share latent factors and cross-variable dependencies, but in order to be realistic, the variables are designed to be non-smooth across time: they rely on near-memoryless shared latent drivers, external shocks, and noise. The generation proceeds as follows: 1. A patient-month panel is constructed from selected demographic and clinical covariates. 2. Within-patient time indices are created. 3. Three patient-level latent traits are drawn once per patient. 4. Two weakly autocorrelated shared state variables drive autocorrelation across series. 5. External shocks introduce abrupt level changes. 6. Ten variable-specific noise terms complete the stochastic structure. 7. For each patient, a recursive simulation initializes each variable from covariates, latent traits, and noise, then iterates forward introducing autocorrelation, cross-variable dependence, treatment indicators, latent states, and shocks. Synthetic outcome: The continuous outcome is a nonlinear function of the synthetic covariates, clinical covariates, treatment assignment, time gap, and lagged terms, with added individual noise and occasional outliers. Ground-truth counterfactuals and benchmark average treatment effect (ATE): Potential outcomes under both treatment arms are computed directly from the structural outcome function, providing exact individual-level counterfactuals and a known benchmark ATE. MNAR missingness mechanism: Missingness is introduced into the ten synthetic covariates at 30%, 50%, and 80% following the MNAR mechanism specified in Equation 1. This mirrors the clinical pattern where measurements are triggered by patient condition rather than scheduled, meaning the probability of observing a value depends on the underlying patient state. This controlled setup allows us to directly evaluate imputation quality and the subsequent robustness of causal estimators against a known ground truth. Directed acyclic graph (DAG) and variable distributions: Figure 3 illustrates the DAG of the generative process. We also show the resulting distribution of the synthetic variables in Figure 4. Lagged covariates are included for realism: lag(Synthetic Variable) at time t is derived from the same variable at time t − 1, so the causal direction is past value → lagged feature at the current step. Although the naming may seem counterintuitive, the causal direction is straightforward: the past value of a variable drives its own lagged feature at the current time step. Figure 4 shows the resulting variable distributions. A.6.2

Target trial emulation using Real-World EHR data

Study design and patient cohort construction: Using the FIRE database, we design an activecomparator, new-user cohort study (emulating a target trial; see Table 3) to estimate the effect of initiating a GLP-1 receptor agonist (GLP-1RA) or an SGLT-2 inhibitor (SGLT-2i) on body weight change from treatment initiation over 1-year. Cohort entry (time-zero) is defined as the date of the first-ever prescription for a GLP-1RA or an SGLT-2i as first second-line glucose-lowering therapy between 01-Jan-2015 and 07-Dec-2024. To be eligible, individuals are aged 18-years or older at cohort entry and have a diagnosis of type 2 diabetes (T2D) before or at cohort entry (see detailed eligibility criteria in Table 3). Treatment strategies and follow-up: We use a per-protocol analysis in which eligible individuals are followed from cohort entry until the earliest of: deviation from the assigned treatment strategy 17

Figure 3: Directed-Acyclic-Graph of the generation of the synthetic variables and the outcome. (defined as a gap between consecutive prescriptions exceeding 365-days, with a 30-day grace period); initiation of the comparator drug class; death; or end of the pre-specified follow-up period (390 days after cohort entry). Outcome: The outcome is the absolute change in bodyweight (kg) from baseline over 1-year. The outcome is extracted in 30-day intervals throughout follow-up. Covariates: Time-fixed baseline covariates measured at cohort entry include demographic characteristics (age, sex assigned at birth). Time-varying covariates updated from baseline every 30-day interval throughout follow-up include current comedications (antihypertensives, statins), comorbidities (cardiovascular disease, chronic kidney disease, hypertension, dyslipidemia), and laboratory measurements (Hemoglobin A1c , eGFR, LDL cholesterol, systolic blood pressure). Confounders and DAG: Using the covariates defined above, causal assumptions are encoded in an expert-specified DAG developed with pharmacoepidemiology experts (Figure 5), following standard practice in target trial emulation [Hernán et al., 2022]. The DAG identifies the minimal sufficient adjustment set required for conditional exchangeability between treatment arms. Real-world application of the proposed ML pipeline: The final cohort comprised 6,114 adults with T2D initiating a GLP-1RA (n = 2,392) or an SGLT-2i (n = 3,722) as their first second-line glucose-lowering therapy. Following our proposed pipeline (section 3), missing covariates or outcome measurements are imputed using the validated LLM-based strategy across follow-up. CausalFlow-T was then fitted on the completed longitudinal panel.

Baseline Architectures for Causal Inference Estimation

B.1

Baseline Selection Criteria

We restrict causal inference baselines to generative distributional models satisfying three jointly necessary criteria: 18

Table 3: Specification and emulation of the target trial emulation protocol aiming to estimate the perprotocol effect of GLP-1 receptor agonist (GLP-1RA) versus SGLT-2 inhibitor (SGLT-2i) initiation on body weight at one year in adults with type 2 diabetes (T2D). Protocol component

Target trial specification

Target trial emulation

Eligibility criteria

Inclusion criteria: Adults (≥18 yrs) initiating a GLP-1RA or an SGLT-2i; T2D diagnosis. Exclusion criteria: Prescription for the comparator drug class before treatment initiation; type 1 diabetes diagnosis; last recorded eGFR ≤30ml/min/1.73m2.

Inclusion criteria: Adults (≥18 yrs) initiating a GLP-1RA or an SGLT-2i between 01-Jan-2015 and 07-Dec-2024 (cohort entry = first prescription date); T2D diagnosis before or at cohort entry; at least one weight measurement recorded within the 365 days before cohort entry. Exclusion criteria: Prescription for the comparator drug class before or at cohort entry; type 1 diabetes diagnosis before or at cohort entry; eGFR ≤30 before or at cohort entry. Missingness in variables need to assess eligibility are imputed using the LLM-based strategy prior to eligibility assessment.

Treatment strategies

(1) Initiate and sustain a GLP-1RA receptor agonist; (2) initiate and sustain an SGLT-2i inhibitor; as add-on to first-line glucose-lowering therapy over 12 months.

Same strategies; patients classified by observed prescription at cohort entry.

Treatment assignment

Random assignment at baseline (open-label); participants not blinded to assigned strategy.

Non-random; conditional exchangeability emulated by adjusting for the minimal sufficient adjustment set identified from the expert-specified DAG (Figure 5).

Outcome

Change in body weight (kg) from baseline to 1-year (12 months), measured in 30-day intervals throughout follow-up.

Same. Missing outcome measurements are imputed using the LLM-based strategy.

Follow-up

From baseline until the earliest of: per-protocol deviation from assigned strategy; initiation of comparator drug class; death; or 390 days after cohort entry.

Same. Per-protocol deviation defined as a gap between consecutive prescriptions exceeding 365 days (and adding a 30-day grace period) or initiation of the comparator drug class.

Causal contrast

Per-protocol effect: treatment effect for all randomized individuals if all adhered to treatment and did not initiate the an agent of the comparator drug class.

Observational analog of the per-protocol effect, estimated under the assumption of conditional exchangeability given the adjustment set identified from the DAG(Figure 5).

Statistical analysis

ATE estimation via marginal structural model with inverse probability of censoring weighting to account for informative censoring due to non-adherence.

Following our proposed pipeline (section 3), missing covariates or outcome measurements are imputed using the validated LLM-based strategy across follow-up. CausalFlow-T was then fitted on the completed longitudinal panel.

Figure 4: Distribution of the synthetic variables and outcome.

Figure 5: Expert-specified directed acyclic graph for the active-comparator, new-user cohort study assessing the effect of GLP-1 receptor agonist versus SGLT-2 inhibitor initiation on weight loss. The DAG is a simplified version with no differentiation between baseline and time varying-confounders. 1. Generative mechanism. The model must define a joint distribution p(Y (0) , Y (1) | X) rather than only conditional means E[Y | X, A], as the latter is insufficient for individual-level counterfactual inference via the AAP procedure [Pearl, 2009]. 2. Invertible abduction. The model must support exact or approximate recovery of patientspecific exogenous noise z from observations, implementing the twin-network assumption required for individual-level counterfactuals; without this, ŷ (0) and ŷ (1) are population-level contrasts rather than structural counterfactuals for the same individual. 3. Distributional evaluability. The model must produce outputs compatible with our five reliability criteria (subgroup calibration, tail variance ratio, arm reconstruction error, HR 20

recovery, and training stability) all of which require access to the full potential outcome distributions p(ŷ (a) ) rather than point predictions. Table 4 summarizes how each candidate model maps onto these criteria and the resulting inclusion decision. Table 4: Baseline selection against the three jointly necessary criteria. A model is included only if all three are satisfied. Model

Generative mechanism

Invertible abduction

Distributional evaluability

Included

✓ ✓ ✓ ✗ ✗ ✗ ✗

✓ ✗ ✗ ✗ ✗ ✗ ✗

✓ ✓ ✓ ✗ ✗ ✗ ✗

✓ ✓ ✓ ✓ ✗ ✗ ✗

CausalFlow-T (ours) CVAE [Louizos et al., 2017] GNN-CVAE TARNet [Shalit et al., 2017] Causal Transformer [Melnychuk et al., 2022] R-MSN [Lim et al., 2018] CRN [Bica et al., 2020]

Discriminative sequence models (including the Causal Transformer [Melnychuk et al., 2022], RMSN [Lim et al., 2018], and CRN [Bica et al., 2020]) fail all three criteria and are therefore excluded by construction, independently of their factual prediction performance. We select CVAE [Louizos et al., 2017], GNN-CVAE, and TARNet [Shalit et al., 2017] to span the key design axes underlying the proposed criteria. CVAE provides a variational inference baseline without causal structure; GNN-CVAE augments this with DAG-structured encoding under the same variational objective, isolating the contribution of causal structure within approximate inference. TARNet, in contrast, serves as a purely discriminative reference point, failing all three criteria. Together, these baselines expose the two failure modes: confounding separation and structural propagation. B.2

Shared Components

All models share: an outcome-reader LSTM (hYt , used only in the encoder to prevent outcome ∗ leakage), a treatment-encoder LSTM (hA t , replaced from t onwards at counterfactual generation cov time), and a covariate aggregator (linear projection to ht ). B.3

Conditional CVAE

Encoder. qϕ (z | x1:T , y1:T ) = N (µϕ , diag(σϕ2 )), cov Y MLPϕ ([flatten(h1:T ), h̄ ]). Latent dimension L = 64.

where

[µϕ , log σϕ2 ]

Decoder. ŷ1:T = MLPθ ([z, flatten(hA 1:T )]). Objective. LCVAE = −Eqϕ [log pθ (y1:T | z, hA 1:T )] + βt KL[qϕ ∥N (0, I)], with linear βt warm-up (0→1 over first 20 epochs). Key limitation. The variational posterior is an approximation; the abducted noise is not the exact exogenous noise of the patient, and this approximation error propagates into counterfactual predictions. B.4

GNN-CVAE

Replaces the flat LSTM covariate encoder with a DAG-structured recurrent encoder performing message passing in topological order: (j) (j),in (j) (j) ht = LayerNorm GRU ht + mt , ht−1 , (6) 21

Transformer Encoder xns 1:T → ht

(L=2, H=4)

Conditioning [xst , xns t , Tt , Yt , mt , τt , tflow , ht ] f

MLP Vector Field uθ outputs velocity v̂ v̂

Euler ODE Integration x0 ∼ N (0, I) → x̂smiss

Figure 6: CFM imputer. The Transformer Encoder produces context ht from always-observed xns 1:T . The MLP uθ is conditioned on ht and DAG evidence (Tt , Yt ), trained with masked MSE on v ∗ =x1 −x0 . At inference, missing values are recovered by ODE integration from Gaussian noise. P (j) (k) (j) where mt = |pa(j)|−1 k∈pa(j) MLPk→j ([ht , ht−1 ]). The causal structure is enforced in the encoder only. However the ELBO approximation gap in the abduction step remains, setting a ceiling that only exact inference can overcome.

Baseline model for MNAR imputation: Conditional DAG-aware Flow Matching

We implement a new Conditional DAG-aware Flow Matching (CFM) imputer as a baseline for MNAR imputation of the synthetic biomarker variables (SynthVars). CFM learns the conditional transport pnoise → p(xsmiss | xsobs , xns , T, Y, τ ) via flow matching, with a loss evaluated only at missing positions. Architecture. A Transformer encoder (L=2, H=4) processes the always-observed non-SynthVar sequence xns 1:T with a key-padding mask on unobserved timesteps, producing a leakage-free context vector ht that interpolates from truly-observed neighbors on both sides of each gap. A three-layer MLP vector field uθ then maps the concatenation s s xtf ∥ xns t ∥ Tt ∥ Yt ∥ mt ∥ xt−1 ∥ mt−1 ∥ ht ∥ τt ∥ tflow to a velocity v̂ ∈ R|S| . Conditioning on Yt and Tt ensures that downstream DAG evidence propagates into every imputed value. The architecture is illustrated in Figure 6. Training. On fully-observed rows a random binary mask mmiss is applied and a linear interpolant xtf = (1 − tf )x0 + tf x1 , x0 ∼ N (0, I) is constructed, with target velocity v ∗ = x1 − x0 . The objective is a masked mean-squared error at missing positions only, requiring no Jacobian: L(θ) =

∥mmiss ∥1

2 mmiss v̂i − vi∗ . i

Inference. The encoder runs once over the full dataset to produce h1:T . Missing values are then recovered by Euler integration of the ordinary differential equation (ODE) dx/dtf = uθ (xtf , ct , tf ) from Gaussian noise to the imputed target, with observed dimensions held fixed throughout integration. 22

Evaluation Protocol

D.1

CausalFlow-T Evaluation Protocol

(i) Subgroup calibration. MAEq =

1X [ ATEq (t) − ATEq (t) , T t (1)

Biasq =

1X [ ATEq (t) − ATEq (t) , T t

(0)

[ q (t) = Ei∈q [ŷi (t)] − Ei∈q [ŷi (t)] and analogously for the ground truth. Quartiles are where ATE P (1) (0) defined by the per-patient mean true ITE, τ̄i = T1 t [yi (t) − yi (t)], ensuring identical patient subgroups across all models. (ii) Arm reconstruction error. Erra =

1X E[ŷ (a) (t)] − E[y (a) (t)] , T t

a ∈ {0, 1},

measured in original outcome units, matching the scale of the mean trajectory plots. (a)

(a)

(iii) Variance calibration. VRQ4 = Var(ŷQ4 )/Var(yQ4 ), averaged over both arms and all time steps within Q4. Values below 1 indicate distributional collapse (the model predicts near-constant outputs for high-effect patients); values above 1 indicate overdispersion. D.2

Imputer Evaluation Protocol

Let M = {(i, t, j) : mi,t,j = 1} denote the set of MNAR-missing biomarker cells, x̂i,t,j the imputed value, and x∗i,t,j the ground truth (available by construction on the FIRE oracle, Section A.6). Section 4.3 reports two layers: biomarker quality (i)–(v) on the reconstruction itself, and downstream causal quality (vi)–(ix) on the fixed CausalFlow-T estimator (Section 4.2). Throughout, d=10 is the number of synthetic biomarkers, N the cohort size, and T the number of timesteps. D.2.1

Biomarker quality

(i) Pointwise MAE. MAE =

1 |M|

x̂i,t,j − x∗i,t,j ,

(i,t,j)∈M

in original biomarker units, restricted to masked positions. (ii) Pointwise RMSE. v u u 1 RMSE = t |M|

2 x̂i,t,j − x∗i,t,j .

(i,t,j)∈M

Like MAE, RMSE is computed in original biomarker units and restricted to masked positions. (iii) Lag-1 autocorrelation error. For biomarker j and patient i, let ri,j (z) = Corrt zi,t,j , zi,t+1,j denote the Pearson correlation between consecutive timesteps of a series z. The mean lag-1 autocorrelations on the imputed and ground-truth datasets are, AC =

1 X ri,j (x̂), N d i,j

AC∗ =

1 X ri,j (x∗ ), N d i,j

and we report |AC − AC∗ |, with AC∗ = 0.366. LOCF inflates AC toward 1; shrinkage-to-mean imputers deflate it toward 0. 23

(iv) Consecutive-step error. 0, mi,t+1,j = 1},

On observed-to-missing transitions T = {(i, t, j) : mi,t,j =

CSE =

1 |T |

x̂i,t+1,j − x∗i,t+1,j ,

(i,t,j)∈T

isolating one-step extrapolation from an observed clinical anchor. (v) Biomarker–outcome correlation inflation. For each biomarker j, let κj (z) = Corr z·,·,j , y·,· be the Pearson correlation between the (imputed or true) biomarker series and the outcome, pooled over all (i, t) pairs. We define inf Y =

d 1X κj (x̂) − κj (x∗ ) , d j=1

with inf ∗Y = 0. Positive values inflate the biomarker–outcome association, while negative values indicate deflation, in which a real signal is washed out. D.2.2

Downstream causal quality

CausalFlow-T is held fixed (Section 4.2) so reported differences reflect imputation alone. (vi) Subgroup calibration MAE. Following Appendix D (i), for each true-ITE quartile q ∈ {Q1 , Q2 , Q3 , Q4 } we compute T

MAEq =

1X [ ATEq (t) − ATEq (t) . T t=1

We focus on the two tail quartiles, which expose complementary failure modes invisible in aggregate metrics: MAEQ1 on low-effect patients, where confounding separation is hardest and systematic bias is most exposed; and MAEQ4 on high-effect patients, where variance collapse manifests as flattening of predicted treatment effects. Table 2 ranks methods by MAEQ1 ; Appendix Table 9 reports both tail quartiles. ¯ (vii) Mean per-arm reconstruction error Erra . As in Appendix D (ii), averaged across the two ¯ a = 1 Erra=0 + Erra=1 , summarizing population-level trajectory accuracy under both arms, Err 2 treatment assignments. (viii) Variance calibration.

VRQ4 as defined in Appendix D (iii).

(ix) Absolute ATE residual. The ATE is the population-level expected difference between the potential outcome under treatment and under control, averaged across patients and timesteps: X (1) 1 X (1) (0) (0) [= 1 ATE ŷi (t) − ŷi (t) , ATE∗ = y (t) − yi (t) , N T i,t N T i,t i [ − ATE∗ |. with ATE∗ = −3.484. We report |ATE

Training and Hyperparameter Details

Tables 5 and 6 define the hyperparameter details for all the experimental parts.

FIRE Semi-Synthetic Oracle Results

Table 7 reports full results on the FIRE semi-synthetic oracle (no missingness), used as a controlled bridge between synthetic benchmarks and the imputation experiments (Section 4.3). True ATE ≈ −3.48; n=10 seeds throughout. 24

Table 5: Hyperparameters across datasets Component

Shared training

Hyperparameter

Synthetic

FIRE semi

FIRE real

Optimizer Learning rate Batch size Data split (train(train/val)/test)

Adam 10−3 512 80(80/20)/20

Adam 10−3 256 80(80/20)/20

Adam 10−4 64 80(80/20)/20

LSTM encoder

Hidden dimension H Dropout

128 -

256 -

256 0.3

CausalMAF

Coupling layers Hidden dimension

1 128

TARNet

Hidden dimension MLP depth

128 2

128 1

Not Used Not Used

GNN-CVAE

Encoder layers Latent dimension

2 64

Not Used Not Used

CVAE

Architecture

No graph encoder

Not Used

Normalization

Covariates & treatments Outcomes

Joint standardization Per-timestep

Per-timestep Per-timestep

Table 6: Main hyperparameters of the LLM-driven evolutionary imputation search. Hyperparameter

Value

LLM model Iteration budget K Proxy holdout fraction ρ Holdout base seed σ Composite-score weight λY Composite-score weight λT Selection rule History window W Per-dataset timeout τrun Memory cap Mrun

gpt-5.4, qwen3.5-plus, gpt-oss-120b 20 0.10 fixed across candidate evaluations 2.0 0.5 minimize s = RMSE + λY ∆Y + λT ∆T 3 180 s 6144 MB

CausalFlow-T (0.316) and GNN-CVAE (0.351) are the only two models passing the <0.5 reliability threshold, indicating that low-responder errors are predominantly random rather than systematic; NF (no DAG) (0.883), TARNet (0.708), and CVAE (0.938) all fail, confirming that explicit causal structure is necessary for confounding separation under realistic clinical covariate distributions. CausalFlow-T is further distinguished as the only model with near-perfect variance preservation (VR = 1.006±0.009), while GNN-CVAE collapses (VR = 0.616) and CVAE collapses more severely (VR = 0.482), confirming that passing the bias threshold alone is insufficient for individual-level reliability. TARNet achieves the best arm-1 reconstruction (0.081) and closest ATE point estimate (−3.415) but fails the ratio threshold (0.708) and collapses variance (VR = 0.608), reinforcing that factual accuracy does not imply counterfactual reliability. Removing the DAG constraint (NF no DAG) degrades the bias ratio to 0.883 and produces the largest ATE deviation (−3.382, ∆ = 0.102) with substantial cross-seed instability (±0.251), underscoring that structural supervision is necessary for stable estimation under clinical covariate distributions.

Full Numerical Results

Tables 10 and 11 report complete per-model, per-benchmark results for the four synthetic benchmarks; they are the primary evidence base for Table 1 in the main body. FIRE semi-synthetic oracle results are in Table 7 (Appendix F). • Simple 3-Node. On Simple 3-Node (linear, no confounding), all models show |Biasq | ≈ MAEq with near-zero standard deviations throughout. This is expected: in a confoundingfree linear setting, residual errors are purely systematic rather than variance-driven. Seed 25

Table 7: FIRE semi-synthetic oracle: systematic-error ratio |BiasQ1 |/MAEQ1 (↓; <0.5 = predominantly random errors), arm reconstruction, variance ratio, and ATE recovery (n=10). VRQ4 → 1. True ATE ≈ −3.48. Bold: best per metric. |BiasQ1 |/MAEQ1 ↓

Model CausalFlow-T (ours) NF (no DAG) TARNet GNN-CVAE (DAG) CVAE

0.316 0.883 0.708 0.351 0.938

Erra=1

Erra=0

VRQ4

ATE

0.088±0.032 0.063±0.022 1.006±0.009 −3.568±0.112 0.161±0.073 0.084±0.044 1.016±0.010 −3.382±0.251 0.081±0.012 0.078±0.003 0.608±0.003 −3.415±0.015 0.216±0.048 0.202±0.038 0.616±0.021 −3.395±0.113 0.134±0.045 0.127±0.047 0.482±0.055 −3.451±0.075

consistency confirms training stability across all architectures, so differences on harder benchmarks reflect genuine architectural sensitivity. GNN-CVAE achieves the lowest Q1 MAE (0.425) in this regime, illustrating that explicit graph encoding can reduce absolute error when confounding is absent (but this advantage disappears under the more complex generating processes below). • LDL Toy. TARNet achieves the lowest absolute MAE across all quartiles (0.392–2.107), but maintains a |Biasq |/MAEq ratio at or near 1.0 throughout (every patient in a given quartile is wrong in the same direction), indicating systematic confounding failure rather than calibrated uncertainty. CausalFlow-T achieves VRQ4 = 1.049 ± 0.033 (closest to 1 across all models), confirming that tail variance is faithfully preserved rather than collapsed or inflated, while maintaining sign-consistent positive bias (+0.547–+0.963 across quartiles) indicating mild but coherent overestimation rather than directionally inconsistent errors. NF (no DAG) systematically underestimates across all quartiles (−0.402 to −0.604), confirming that causal factorization is necessary for directionally correct inference. GNN-CVAE and CVAE show variance collapse (VRQ4 ≤ 0.85) and substantially larger arm errors, revealing distributional narrowing in the high-effect tail. • CVD Risk Ground-truth HR = 0.831 (protective treatment). CausalFlow-T recovers HR = 0.786 ± 0.051 (∆ = 0.045), the only model within clinically meaningful range with arm errors an order of magnitude below all competitors. GNN-CVAE and CVAE collapse to HR ≈ 1.005 ± 0.030 across all seeds (predicting a null effect despite a true 17% hazard reduction) consistent with posterior collapse under the ELBO objective. TARNet predicts HR = 1.133 ± 0.148, inverting the effect direction on average. NF (no DAG) yields HR = 0.834 ± 0.324: the large standard deviation (related to the explosion in the variance inf Q4) reveals that without the DAG constraint the flow is fundamentally unstable under strong mediation (sometimes approximately correct, sometimes catastrophically wrong), making it unreliable for clinical deployment regardless of average performance. • Cox Survival. CVAE achieves lower Q1–Q4 MAE on Cox (0.023 ± 0.001–0.020 ± 0.006) because variational models handle binary outcomes without continuous relaxation. However, CausalFlow-T achieves the best arm-1 reconstruction error (0.014 ± 0.002) and closest HR recovery (0.866 ± 0.011 vs. true 0.887), while CVAE and GNN-CVAE show arm errors 4–6× higher. TARNet achieves the best arm-0 error (0.005 ± 0.001) but overshoots the HR (0.966 ± 0.005), indicating factual accuracy without causal calibration. The Cox MAE advantage of CVAE does not generalize to more complex DAGs, as the CVD Risk results confirm. CVAE vs. GNN-CVAE. Across benchmarks, CVAE consistently outperforms GNN-CVAE despite the latter incorporating explicit graph structure in the encoder. This contrasts with the normalizing flow setting, where the DAG constraint yields consistent improvement. Both CVAE variants rely on variational inference: the ELBO approximation introduces an irreducible inference gap that dominates performance, limiting the benefit of improved structural inductive bias. The GNN encoder does not translate into better counterfactual estimates and introduces additional optimization complexity, evidenced by GNN-CVAE’s substantially higher LDL arm-1 error (4.854 ± 0.451) relative to plain CVAE (0.535 ± 0.245). Normalizing flows perform exact likelihood-based inference, so architectural improvements such as DAG factorization directly improve abduction quality. This explains the pattern CVAE > GNN-CVAE, whereas CausalFlow-T > NF (no DAG). 26

Full Imputation Results

Tables 8 and 9 report complete per-method, per-missingness-level results underlying Table 2. Figure 7 summarizes the normalized proxy-score trajectories during the selected GPT-5.4 evolutionary search, showing rapid early improvement followed by smaller accepted refinements. To separate the effect of the evolutionary loop from a single executable LLM proposal, the tables also report “GPT-5.4 first-valid”, defined as the first GPT-5.4 candidate that passed static and runtime checks and was evaluated. This row is diagnostic only and is not included in the pooled ranks used for model selection. The first-valid rows show that executable LLM proposals are not sufficient on their own, especially at 50%–80% missingness, where the final selected imputer substantially improves biomarker reconstruction and ATE recovery. Biomarker-level metrics are computed once over all imputed cells across the ten semi-synthetic biomarkers and are therefore reported without seed variability; downstream causal metrics are averaged over n = 10 training seeds.

Figure 7: Evolutionary search progress for the LLM-driven imputer across 30%, 50%, and 80% MNAR missingness, shown as normalized proxy scores over the 20-iteration budget.

DAG Sensitivity Analysis

We evaluate CausalFlow-T under two partial-graph misspecifications on the FIRE semi-synthetic dataset (true ATE = −3.484; n = 10 seeds). Three baselines serve as reference points: NF (no DAG) shares CausalFlow-T’s exact normalizing flow inference but uses no causal structure, isolating the contribution of graph specification within the exact-inference family; CVAE and TARNet use neither exact inference nor graph structure and are fully DAG-invariant, producing identical results across all graph variants by construction. All covariate-to-outcome edges are preserved across variants, isolating confounding path misspecification from outcome model misspecification. We evaluate two misspecification types that expose complementary robustness properties: (i) a unique-pathway removal (Asthma), where the removed covariate has no correlated substitute in the adjustment set, producing maximal bias ratio deterioration; and (ii) a redundant-pathway removal (Hypertension), where correlated covariates partially absorb the misspecification, producing instability rather than systematic bias as the primary signal. Results can be found in Table 12 The two misspecification types evaluated here are structurally motivated to expose complementary failure modes. Peripheral misspecification (unique pathway removal) represents the worst case for bias: the removed covariate has no correlated substitute, so confounding separation degrades maximally. Central misspecification (redundant pathway removal) represents the worst case for stability: correlated covariates partially absorb the missing path, producing cross-seed instability rather than systematic bias as the primary signal. Together these two cases bound the space of singleedge removals along the bias-instability tradeoff axis. While multi-edge removals, edge direction errors, and latent confounders represent additional misspecification regimes not evaluated here, the 27

Table 8: Biomarker-level imputation quality across missingness levels, averaged over the semisynthetic biomarkers. Pointwise MAE and RMSE are computed at missing positions only in original biomarker units. Lag-1 autocorrelation AC is compared against the ground-truth value AC∗ = 0.375; CSE denotes consecutive-step error; and inf Y is biomarker–outcome correlation inflation (ground truth: 0).“GPT-5.4 first-valid” denotes the first executable GPT-5.4 proposal evaluated during the evolutionary search.

Missing

Method

MAE ↓

RMSE ↓

AC (→ 0.375)

30%

LOCF MissForest CausalCFM GPT-5.4 Qwen3.5-Plus GPT-OSS-120b

3.314 2.090 2.755 1.914 2.055 1.972

4.399 2.674 3.523 2.461 2.659 2.548

0.473 0.384 0.350 0.392 0.375 0.386

2.781 2.073 2.708 1.893 1.956 1.914

−0.113 +0.036 −0.007 +0.034 +0.003 +0.026

GPT-5.4 first-valid

2.026

2.598

0.389

2.009

+0.036

LOCF MissForest CausalCFM GPT-5.4 Qwen3.5-Plus GPT-OSS-120b

3.378 2.219 2.833 2.052 2.461 2.059

4.463 2.839 3.623 2.628 3.163 2.646

0.554 0.390 0.334 0.413 0.314 0.423

2.773 2.185 2.758 2.003 2.349 2.011

−0.197 +0.059 −0.009 +0.052 −0.013 +0.054

GPT-5.4 first-valid

2.805

3.741

0.619

2.343

−0.138

LOCF MissForest CausalCFM GPT-5.4 Qwen3.5-Plus GPT-OSS-120b

3.407 2.504 2.875 2.059 2.788 2.147

4.509 3.187 3.679 2.634 3.645 2.745

0.713 0.380 0.348 0.596 0.208 0.514

2.760 2.468 2.798 2.035 2.612 2.109

−0.347 +0.072 −0.011 +0.159 −0.080 +0.129

GPT-5.4 first-valid

3.146

4.193

0.729

2.442

−0.318

50%

80%

CSE ↓ inf Y (→ 0)

instability diagnostic identified in the central case (detectable through cross-seed variance without access to ground-truth counterfactuals) generalizes as a practical signal across misspecification types, since any structural misspecification that disrupts the adjustment set will manifest as estimation instability across random initializations Interpretation. A partially correct DAG is better than no DAG. Comparing CausalFlow-T against NF (no DAG) within the exact-inference family isolates the contribution of causal structure from inference quality. Under central misspecification, CausalFlow-T’s bias ratio (0.517) still improves on NF (no DAG)’s DAG-invariant 0.883, demonstrating that even an incomplete causal graph provides meaningful confounding separation beyond what exact inference alone achieves. Distributional reliability is preserved under misspecification. CausalFlow-T maintains near-perfect variance calibration across both removals (VRQ4 = 0.996 and 1.005), while CVAE (0.482) and TARNet (0.608) collapse regardless of graph quality. NF (no DAG) also preserves VRQ4 (1.016) but at the cost of a consistently high bias ratio (0.883), confirming that variance calibration and confounding separation are non-redundant properties requiring both exact inference and causal structure jointly. Misspecification is detectable through instability. The central removal produces a 4.5× increase in CausalFlow-T’s cross-seed standard deviation relative to the full DAG (0.504 vs. 0.112). DAGinvariant models show no such signal (TARNet’s seed SD remains 0.015 and NF (no DAG)’s 0.251 across all conditions) meaning they fail silently under the same misspecification that CausalFlow-T 28

Table 9: Downstream causal inference quality. All methods use the same fixed CausalFlow-T estimator (Section 4.2); only the imputation strategy varies. Q1 MAE and Q4 MAE (lowest- and ¯ a = 1 (Erra=0 + Erra=1 ), tail variance highest-ITE quartiles), mean per-arm reconstruction error Err 2 ratio VRQ4 (target value 1), and absolute ATE residual against ATE∗ = −3.484. Means ± standard ¯ a= deviations over n = 10 seeds. † Oracle reference: CausalFlow-T on fully observed data, Err 0.058 (Table 7). “GPT-5.4 first-valid” denotes the first executable GPT-5.4 proposal evaluated during the evolutionary search. Missing

30%

50%

80%

Q1 MAE ↓

Q4 MAE ↓

¯ a↓ Err

VRQ4

|∆ATE| ↓

LOCF MissForest CausalCFM GPT-5.4 Qwen3.5-Plus GPT-OSS-120b

0.224 ± 0.123 0.416 ± 0.137 0.161 ± 0.073† 0.223 ± 0.120 0.201 ± 0.107 0.180 ± 0.040

0.153 ± 0.022 0.401 ± 0.144 0.118 ± 0.015 0.159 ± 0.036 0.220 ± 0.054 0.181 ± 0.034

0.079 ± 0.023 0.177 ± 0.093 0.059 ± 0.006 0.071 ± 0.016 0.086 ± 0.035 0.065 ± 0.010

1.005 ± 0.006 0.998 ± 0.003 1.008 ± 0.005 1.007 ± 0.011 1.004 ± 0.004 1.005 ± 0.007

0.090 ± 0.113 0.331 ± 0.223 0.027 ± 0.066 0.024 ± 0.120 0.145 ± 0.092 0.075 ± 0.053

GPT-5.4 first-valid

0.363 ± 0.239

0.556 ± 0.159

0.233 ± 0.073

1.008 ± 0.009

0.464 ± 0.199

LOCF MissForest CausalCFM GPT-5.4 Qwen3.5-Plus GPT-OSS-120b

0.224 ± 0.074 0.305 ± 0.143 0.232 ± 0.105 0.219 ± 0.067 0.187 ± 0.099 0.194 ± 0.095

0.146 ± 0.049 0.342 ± 0.192 0.206 ± 0.105 0.199 ± 0.043 0.132 ± 0.034 0.283 ± 0.120

0.084 ± 0.031 0.155 ± 0.084 0.096 ± 0.048 0.083 ± 0.029 0.069 ± 0.022 0.116 ± 0.056

1.008 ± 0.008 0.999 ± 0.006 1.007 ± 0.006 1.012 ± 0.007 1.008 ± 0.006 1.010 ± 0.006

0.007 ± 0.174 0.184 ± 0.306 0.048 ± 0.198 0.131 ± 0.085 0.013 ± 0.114 0.199 ± 0.149

GPT-5.4 first-valid

0.263 ± 0.146

0.285 ± 0.191

0.129 ± 0.063

1.005 ± 0.008

0.216 ± 0.201

LOCF MissForest CausalCFM GPT-5.4 Qwen3.5-Plus GPT-OSS-120b

0.235 ± 0.045 0.175 ± 0.044 0.548 ± 0.358 0.311 ± 0.133 0.370 ± 0.220 0.222 ± 0.105

0.176 ± 0.050 0.253 ± 0.169 0.440 ± 0.244 0.183 ± 0.086 0.219 ± 0.080 0.211 ± 0.153

0.088 ± 0.031 0.102 ± 0.070 0.226 ± 0.132 0.093 ± 0.050 0.133 ± 0.063 0.088 ± 0.044

1.009 ± 0.005 1.001 ± 0.002 1.037 ± 0.024 1.016 ± 0.005 1.032 ± 0.007 1.019 ± 0.008

0.052 ± 0.163 0.145 ± 0.188 0.408 ± 0.332 0.013 ± 0.191 0.190 ± 0.185 0.072 ± 0.172

GPT-5.4 first-valid

0.283 ± 0.150

0.327 ± 0.182

0.140 ± 0.064

1.010 ± 0.006

0.201 ± 0.205

Method

flags through detectable instability. This diagnostic property is practically valuable: a practitioner monitoring cross-seed variance can detect central misspecification without access to ground-truth counterfactuals. Bias is directionally consistent and bounded. Both misspecifications produce attenuation towards zero (peripheral: −3.852; central: −3.896; true: −3.489), consistent with standard epidemiological reasoning where uncontrolled confounding attenuates protective effects. The direction of bias is therefore predictable from domain knowledge even when its magnitude is not, providing a practical prior for sensitivity reasoning in real-world deployments where some DAG uncertainty is unavoidable.

Sensitivity Analysis for the LLM Imputer

Table 13 reports reports bootstrap sensitivity intervals for the primary endpoint estimates obtained after applying the selected LLM imputer. After imputation, patients were resampled with replacement. For each bootstrap replicate, the weighting models were re-estimated and the ATE was recomputed. We used 500 bootstrap replicates and percentile-based 95% confidence intervals. For computational reasons, we applied this bootstrap to the completed imputed datasets rather than repeating the full imputation procedure inside each bootstrap replicate. Therefore, these intervals quantify sampling and downstream weighting uncertainty conditional on the completed imputed dataset; they do not capture full re-imputation uncertainty. 29

−0.005±0.005 0.007±0.007 0.011±0.000 0.012±0.000 0.021±0.012

−0.007±0.008 0.021±0.011 0.027±0.001 0.027±0.000 0.080±0.011

CausalFlow-T 0.010±0.005 0.007±0.004 0.004±0.002 0.003±0.001 NF (no DAG) 0.023±0.010 0.009±0.004 0.006±0.003 0.557±1.750 GNN-CVAE 0.027±0.001 0.011±0.000 0.008±0.000 0.004±0.000 CVAE 0.027±0.000 0.012±0.000 0.008±0.000 0.004±0.000 TARNet 0.080±0.011 0.022±0.011 0.011±0.007 0.013±0.010 † GNN-CVAE and CVAE evaluated on the full LDL / full Cox simulation variant.

CVD Risk Toy

0.025±0.003 0.024±0.006 0.020±0.007 0.014±0.007 0.047±0.002

0.027±0.003 0.028±0.009 0.015±0.005 0.014±0.007 0.042±0.004

2.393±0.449 2.752±0.442 7.563±0.658 1.477±0.418 1.003±0.065

3.564±0.626 3.769±0.453 1.729±0.565 4.216±0.499 2.107±0.114

−0.003±0.003 0.003±0.007 0.008±0.000 0.008±0.000 0.001±0.013

−0.011±0.007 −0.016±0.015 −0.002±0.007 −0.002±0.008 0.042±0.004

−0.010±0.008 −0.012±0.013 −0.006±0.012 0.002±0.011 0.047±0.002

−0.005±0.018 −0.011±0.011 −0.014±0.017 0.011±0.018 0.047±0.003

0.025±0.008 0.022±0.005 0.026±0.012 0.023±0.001 0.047±0.003

2.326±0.289 2.449±0.427 13.111±0.833 1.540±0.242 0.842±0.071

0.000±0.002 0.555±1.751 0.004±0.000 0.004±0.000 −0.007±0.014

−0.021±0.013 −0.029±0.012 −0.000±0.003 −0.002±0.006 0.029±0.003

0.963±1.020 −1.004±0.650 1.677±0.568 −4.212±0.494 2.107±0.114

0.034±0.008 0.038±0.009 0.018±0.002 0.020±0.006 0.029±0.003

CausalFlow-T NF (no DAG) GNN-CVAE† CVAE† TARNet

Cox Survival

2.029±0.201 2.057±0.320 23.056±1.133 7.005±0.554 0.392±0.081

0.781±0.657 −0.773±0.641 7.516±0.653 −1.477±0.418 0.885±0.068

Q1 0.433±0.002 0.499±0.000 0.425±0.000 0.514±0.000 0.549±0.000 0.677±0.569 −0.713±0.734 13.071±0.854 1.296±0.311 0.751±0.074

Q4 0.553±0.001 0.586±0.000 0.659±0.000 0.566±0.000 0.546±0.000 0.547±0.639 −0.402±0.356 23.035±1.170 6.809±0.543 −0.546±0.126

CausalFlow-T NF (no DAG) GNN-CVAE† CVAE† TARNet

LDL

Q3 0.151±0.005 0.182±0.000 0.248±0.000 0.156±0.000 0.133±0.000

Q2 0.127±0.001 0.093±0.000 0.110±0.000 0.121±0.000 0.147±0.000

−0.353±0.001 −0.586±0.000 −0.659±0.000 −0.566±0.000 −0.546±0.000

0.533±0.002 0.499±0.000 0.425±0.000 0.514±0.000 0.549±0.000

Quartile Bias −0.148±0.001 −0.182±0.000 −0.248±0.000 −0.156±0.000 −0.133±0.000

CausalFlow-T NF (no DAG) GNN-CVAE CVAE TARNet

Simple 3-Node

Quartile MAE 0.110±0.001 0.093±0.000 0.028±0.000 0.121±0.000 0.147±0.000

Model

Dataset

1.036±0.001 expl. expl. expl. 1.032±0.051

1.528±0.016 1.522±0.009 0.916±0.035 0.849±0.038 0.855±0.058

1.049±0.033 1.098±0.025 0.846±0.033 0.698±0.026 0.639±0.003

0.991±0.001 0.991±0.000 0.852±0.000 0.763±0.000 0.810±0.000

VR Q4

Table 10: Subgroup calibration and tail variance ratio across all benchmarks. GNN-CVAE and CVAE on LDL and Cox use the full simulation variant. Bold: best value(s) per dataset per metric. “expl.” = numerical explosion.

Table 11: Arm reconstruction errors and ATE/HR recovery. Bold: best per dataset per metric. † GNN-CVAE and CVAE on LDL/Cox use the full simulation variant. Erra=1

Erra=0

ATE or HR (true / pred)

CausalFlow-T NF (no DAG) GNN-CVAE CVAE TARNet

0.025±0.004 0.028±0.000 0.203±0.000 0.200±0.000 0.081±0.000

0.019±0.003 0.022±0.000 0.239±0.000 0.204±0.000 0.064±0.000

−0.940 / −0.950±0.001 −0.940 / −0.984±0.000 −0.942 / −1.055±0.000 −0.942 / −0.964±0.000 −0.942 / −0.938±0.000

LDL

CausalFlow-T NF (no DAG) GNN-CVAE† CVAE† TARNet

0.963±0.119 1.035±0.154 4.854±0.451 0.535±0.245 0.192±0.076

1.549±0.240 1.697±0.218 6.671±0.773 0.786±0.372 0.757±0.075

−28.351 / −27.609±0.579 −28.351 / −28.874±0.459 −28.351 / −17.027±0.789 −28.351 / −27.747±0.337 −28.351 / −27.527±0.066

Cox Survival

CausalFlow-T NF (no DAG) GNN-CVAE† CVAE† TARNet

0.014±0.002 0.015±0.005 0.076±0.016 0.072±0.017 0.040±0.002

0.013±0.002 0.012±0.002 0.071±0.009 0.074±0.007 0.015±0.001

0.887 / 0.866±0.011 0.887 / 0.857±0.021 0.887 / 0.865±0.020 0.887 / 0.879±0.022 0.887 / 0.966±0.005

CVD Risk Toy

CausalFlow-T NF (no DAG) GNN-CVAE CVAE TARNet

0.003±0.001 0.006±0.002 0.053±0.000 0.052±0.000 0.040±0.005

0.002±0.001 0.143±0.438 0.065±0.000 0.065±0.000 0.017±0.008

0.831 / 0.786±0.051 0.831 / 0.834±0.324 0.831 / 1.006±0.030 0.831 / 1.005±0.007 0.831 / 1.133±0.148

Dataset

Model

Simple 3-Node

Table 12: Sensitivity to DAG misspecification on FIRE semi-synthetic (true ATE = −3.489, n = 10 seeds). NF (no DAG), CVAE, and TARNet are DAG-invariant and produce identical results across all graph variants by construction († ). VRQ4 → 1; |Bias|/MAE < 0.5 indicates predominantly random errors.

|Bias|/MAEQ1 ↓

VRQ4 → 1

ATE

Seed SD

Full DAG

CausalFlow-T NF (no DAG)† CVAE† TARNet†

0.316 0.883 0.938 0.708

1.006 1.016 0.482 0.608

−3.568 −3.382 −3.451 −3.415

0.112 0.251 0.075 0.015

Peripheral (Asthma removed)

CausalFlow-T NF (no DAG)† CVAE† TARNet†

0.955 0.883 0.938 0.708

0.996 1.016 0.482 0.608

−3.852 −3.382 −3.451 −3.415

0.145 0.251 0.075 0.015

Central (HTN removed)

CausalFlow-T NF (no DAG)† CVAE† TARNet†

0.517 0.883 0.938 0.708

1.005 1.016 0.482 0.608

−3.896 −3.382 −3.451 −3.415

0.504 0.251 0.075 0.015

DAG variant

Model

Theoretical Aspects of CausalFlowT

K.1

Exact Likelihood vs. ELBO: Why the Gap Matters

CVAEs maximize LELBO = Eqϕ [log pθ (v | z)] − KL[qϕ (z | v)∥p(z)] ≤ log pθ (v). The gap log pθ (v) − LELBO ≥ 0 is the variational approximation error. Normalizing flows compute the exact log-likelihood via the change-of-variables formula: log p(vt ) = log pz (fθ (vt )) + log | det ∂fθ /∂vt |. For MAF-type flows the Jacobian is triangular and computable in O(D). The AAP procedure is therefore: (1) Abduction: zt = fθ (vt ) (exact and deterministic), (2) Action: at ← a′ , (3) Prediction: decode all descendants of A in G, holding non-descendants fixed. Because inversion is exact, the same zt is used under both arms, implementing the twin-network assumption underlying individual-level counterfactual validity [Pearl, 2009]. 31

Table 13: Bootstrap sensitivity intervals for selected LLM-imputed primary endpoint estimates. IPTW denotes inverse probability of treatment weighting; IPCW denotes inverse probability of censoring weighting.

K.2

Miss.

Estimator

30% 50% 80% Real

IPTW IPTW IPTW IPTW–IPCW

Point ATE

95% CI

-3.230 -3.228 -3.172 -0.936

[−3.354, −3.101] [−3.354, −3.095] [−3.323, −3.011] [−1.335, −0.553]

The Role of the DAG in Counterfactual Propagation

For an unconstrained autoregressive flow with ordering (A, X, Y ), replacing A ← a′ causes the decoder to adjust X along the learned autoregressive path, traversing the A → X direction which does not exist in the causal graph. This anti-causal propagation produces systematically biased counterfactuals even when the factual distribution is perfectly fitted. CausalFlow-T fixes the autoregressive ordering to a topological sort of G, so the intervention propagates only through causal descendants, implementing do-calculus exactly. K.3

Limitations of Continuous Relaxation for Binary Outcomes

Normalizing flows require dequantization of binary targets, effectively relaxing the discrete problem into a continuous approximation. This enables exact likelihood-based inference but introduces a modeling approximation that can affect calibration for survival outcomes. On the Cox benchmark, this explains why variational approaches can match or exceed flows in MAE and VRQ4 , while CausalFlow-T maintains advantage in hazard-ratio recovery and arm errors (the structural benefit persists at a calibration cost). Future work on discrete or hybrid flow-based models for survival outcomes will address this limitation.

Theoretical Aspects of LLM-driven Evolutionary Imputation

L.1

Search Algorithm

The search maintains a single current-best candidate and accepts a new proposal only if it strictly improves the composite score s(g) defined in Eq. (4). Starting from a deterministic seed imputer g (0) , Algorithm 1 repeats for K iterations. The same proxy holdout Ωp is used at every iteration so that all candidates are scored on identical hidden cells. Real-world proxy holdout. For the semi-synthetic missingness benchmarks, the proxy holdout Ωp is sampled independently at the observed-cell level. In the real FIRE application, the imputation target set includes all incomplete analysis variables, including the outcome. Because raw clinical measurements are expanded onto a 30-day longitudinal target grid, the same observed measurement can appear across adjacent rows; a cellwise proxy holdout would therefore allow trivial reconstruction from neighboring rows via forward/backward copying, artificially favoring LOCF-like candidates. We therefore sample Ωp at the level of contiguous within-patient observed-value runs: for each patient and target variable, consecutive rows with the same observed value are treated as one segment, and selected segments are masked in full with probability ρ. The selected segment-level holdout is fixed across candidate evaluations. Single-parent scheme. We deliberately avoid population-based search. The evaluation cost is e− , which is cheap relative to the cost of an LLM call. A singledominated by the proxy fit on D parent scheme keeps the prompt focused on a single, well-characterized current best, which we found empirically to converge faster than larger populations under a fixed call budget. The strictimprovement rule is equivalent to a deterministic acceptance step in a single-parent evolution strategy without recombination. 32

Algorithm 1 LLM-driven evolutionary imputation search (high-level view). ed }d∈{30,50,80} ; seed imputer g (0) ; holdout fraction ρ; correlation Require: Incomplete datasets {D weights λY , λT ≥ 0; budget K; history length W . Ensure: One imputer gd⋆ per missing-rate dataset. 1: for d ∈ {30, 50, 80} do ▷ three independent runs (d) ed independently with 2: Sample proxy holdout Ωp by masking each observed target cell of D probability ρ (fixed seed) 3: gd⋆ ← g (0) ; s⋆d ← S CORE(g (0) , d); H ← ∅ 4: for k = 1, . . . , K do 5: g (k) ← LLM P ROMPT gd⋆ , s⋆d , H[−W :] 6: if g (k) fails static check or runtime guard then 7: continue ▷ log failed 8: end if 9: sk ← S CORE(g (k) , d) 10: if sk < s⋆d then ▷ minimize scalar s 11: gd⋆ ← g (k) ; s⋆d ← sk ▷ accept 12: end if 13: H ← H ∪ {(k, g (k) , sk )} 14: end for 15: end for 16: return {gd⋆ }d∈{30,50,80} 17: function S CORE(g, d) (d)

18: Compute RMSE(g), ∆Y (g), ∆T (g) on Ωp 19: return s ← RMSE(g) + λY ∆Y (g) + λT ∆T (g) 20: end function

L.2

▷ plus MAE(g) as diagnostic for H

Prompting

At iteration k, the prompt pk assembled in line 4 of Algorithm 1 contains three blocks, in order: ⋆ • Current-best source code. The full Python module implementing gk−1 , including any helper functions. The LLM is asked to return a complete replacement module rather than local edits; this avoids whole classes of patch-application errors and makes acceptance a binary check on the new module. ⋆ • Proxy summary of gk−1 . A short text block stating the selection rule (minimize s = RMSE + λY ∆Y + λT ∆T ) and the current-best values of RMSE, MAE, ∆Y , and ∆T on the proxy holdout. Each missingness regime is searched independently, so the summary is regime-specific.

• Search history H. A short log of the last W proposals from the same run, summarizing for each one whether it was accepted, rejected without improvement, or failed. Accepted and rejected proposals contribute their proxy values RMSE, ∆Y , and ∆T , with MAE included only as diagnostic feedback; failed proposals contribute a short reason (e.g. banned import, timeout, malformed output). This lets the LLM see which directions have already been tried and which kinds of mistakes to avoid in the next proposal. L.3

Constraints and Failure Handling

Each candidate g (k) is screened in two stages before its proxy score is computed. Static checks. Before execution, the candidate source is parsed and rejected if it tries to access the file system, the network, or any form of dynamic code execution (e.g. eval, exec, disk I/O via pandas). Runtime guard. Candidates that pass the static check run in an isolated child process with a per-dataset time and memory budget (Table 6). Any crash, timeout, or attempt to overwrite observed 33

values or non-target columns aborts the candidate; the failure is recorded in H and the search moves on.

Compute Resources All experiments were run on a single node of a HPC cluster, equipped with 4× NVIDIA L40S GPUs (46 GB VRAM each, 350 W TDP), CUDA 12.8, and driver version 570.172.08. In practice, all experiments ran on a single GPU, with peak memory usage well below 2 GB, confirming the pipeline is lightweight relative to available hardware. Synthetic and semi-synthetic benchmark experiments (Sections 4.2–4.3) require on the order of a few minutes per run; the LLM evolutionary imputation search (K = 20 iterations) takes approximately 00:57 hours per missingness level; and the real-world TTE experiment (Section 5) requires approximately 00:58 hours end-to-end. Each fully imputed dataset with GPT-5.4 incurred an estimated API cost of approximately $2–3. All models are implemented in PyTorch. Total compute across all experiments, including n = 10 seed repetitions through the full pipeline, amounts to approximately 10 GPU-hours on a single L40S.

Record · ID 158543 · SHA-256 0eefe9e8a0f2990c
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.