Organizations: Kyoto University Graduate School of Medicine, Kyoto, Japan · SANKEN, The University of Osaka, Osaka, Japan · Faculty of Data Science, Shiga University, Shiga, Japan · AIP, RIKEN, Tokyo, Japan · Faculty of Medicine, University of Miyazaki, Miyazaki, Japan · Hiroshima University Graduate School of Biomedical and Health Sciences, Hiroshima, Japan
Causal discovery has achieved substantial theoretical progress, yet its deployment in large-scale longitudinal systems remains limited. A key obstacle is that operational data are generated under institutional workflows whose induced partial orders are rarely formalized, enlarging the admissible graph space in ways inconsistent with the recording process. We characterize a workflow-induced constraint class for longitudinal causal discovery that restricts the admissible directed acyclic graph space through protocol-derived structural masks and timeline-aligned indexing. Rather than introducing a new optimization algorithm, we show that explicitly encoding workflow-consistent partial orders reduces structural ambiguity, especially in mixed discrete--continuous panels where within-time orientation is weakly identified. The framework combines workflow-derived admissible-edge constraints, measurement-aligned time indexing and block structure, bootstrap-based uncertainty quantification for lagged total effects, and a dynamic representation supporting intervention queries. In a nationwide annual health screening cohort in Japan with 107,261 individuals and 429,044 person-years, workflow-constrained longitudinal LiNGAM yields temporally consistent within-time substructures and interpretable lagged total effects with explicit uncertainty. Sensitivity analyses using alternative exposure and body-composition definitions preserve the main qualitative patterns. We argue that formalizing workflow-derived constraint classes improves structural interpretability without relying on domain-specific edge specification, providing a reproducible bridge between operational workflows and longitudinal causal discovery under standard identifiability assumptions.
Figures & tables
Figure 1: Bootstrap distributions of lagged total effects from health guidance in 2020 to outcomes measured in 2021–2023 (lags 0–2), under the workflow-derived constraints ( B=1000 ). Panels are arranged in a 3×5 grid: rows correspond to lag and columns to outcomes (BMI, SBP, DBP, HbA1c, LDL). Dashed vertical lines indicate the 95% bootstrap percentile interval, and solid vertical lines mark zero.
Figure 2: Compact recurring subgraph (motif) summarizing within-time relations among the five continuous health screening outcomes. Directed edges indicate directions that are consistent across time points 1–3 under the workflow-derived constraints. The undirected SBP–DBP link indicates a recurring adjacency whose direction varies across time points. The motif is a descriptive summary of recurrent within-time structure, not the full longitudinal graph.
Appendix figures & tables14 assets
Supplementary material from the paper’s appendix.
Appendix
Figure A.1: Variable set and time-point alignment used in this study. Each time point contains 15 variables: health guidance, five continuous health screening outcomes (BMI, SBP, DBP, HbA1c, LDL), three medication indicators, three lifestyle indicators, demographics (Age, Sex), and Check_num. (attendance history in the three years preceding the measurement year). Cells marked with † are shown only to illustrate the fixed tensor layout across time points and are excluded from model fitting.
Figure A.2: Schematic of the health guidance assignment logic used to construct the assignment indicator for sensitivity analysis. The selection proceeds from waist circumference screening to assessment of risk factor status, with medication-treated individuals excluded from eligibility (details are simplified for exposition).
Figure B.1: Conceptual diagram of the workflow-derived prior knowledge. The prior knowledge encodes (i) time ordering across annual visits, (ii) within-time block ordering consistent with the recording protocol, and (iii) admissible cross-time links restricted to one-year lag ( t−1→t ). The prior knowledge is derived from generic invariances and recording-protocol constraints, not from medical causal assumptions.
Figure D.1: Sensitivity analysis (assignment indicator): bootstrap distributions of lag-specific total effects anchored at 2020. Rows correspond to outcome measurement years 2021–2023 (lags 0–2), and columns correspond to outcomes (BMI, SBP, DBP, HbA1c, LDL). Dashed vertical lines indicate the 95% bootstrap percentile interval; the solid vertical lines mark zero when it lies within the plotted range (otherwise zero is indicated by an arrow). Distributions are computed by subject-level bootstrap resampling ( B=1000 ) under the workflow-derived prior-knowledge constraints.
Figure D.2: Sensitivity analysis (waist circumference): bootstrap distributions of lag-specific total effects anchored at 2020 when replacing BMI with waist circumference. Rows correspond to outcome measurement years 2021–2023 (lags 0–2), and columns correspond to outcomes (Waist, SBP, DBP, HbA1c, LDL). Dashed vertical lines indicate the 95% bootstrap percentile interval, and solid vertical lines mark zero. Distributions are computed by subject-level bootstrap resampling ( B=1000 ) under the workflow-derived prior-knowledge constraints.
Figure D.3: Sensitivity analysis (body weight): bootstrap distributions of lag-specific total effects anchored at 2020 when replacing BMI with body weight. Rows correspond to outcome measurement years 2021–2023 (lags 0–2), and columns correspond to outcomes (Weight, SBP, DBP, HbA1c, LDL). Dashed vertical lines indicate the 95% bootstrap percentile interval, and solid vertical lines mark zero. Distributions are computed by subject-level bootstrap resampling ( B=1000 ) under the workflow-derived prior-knowledge constraints.
Figure E.1: Full learned longitudinal directed graph (point estimate) under the primary specification. The graph is estimated using workflow-constrained longitudinal causal discovery with the prior-knowledge mask described in Appendix B . Time points follow measurement-year indexing, and cross-time links are restricted to one-year lag ( t−1→t ). The figure shows the full structure from which the compact motif summary in the main text is extracted. Edges with ∣β^∣<0.01 are omitted by the package’s default visualization.
Figure E.2: Zoomed view (crop) of Figure E.1 for the first block ( Health-guidance (2020) and outcomes measured in 2021).
Figure E.3: Zoomed view (crop) of Figure E.1 for the second block ( Health-guidance (2021) and outcomes measured in 2022).
Figure E.4: Zoomed view (crop) of Figure E.1 for the third block ( Health-guidance (2022) and outcomes measured in 2023).
Figure F.1: Overview of the practitioner-facing simulator (user interface in Japanese; English callouts added for the paper).
Figure F.2: Forward prediction (What-if): changing a baseline variable X and predicting outcome Y after k years under the fitted longitudinal system.
Figure F.3: Goal-seeking (Target setting): specifying a desired outcome value and solving for the required baseline value under the fitted system.
Figure F.4: Decision messaging (guardrail): suppressing recommendations when the bootstrap 95% CI includes 0.
Causal Bayesian Networks (CBNs) are a powerful tool for reasoning under uncertainty about complex real-world problems. Such problems evolve over time, responding to external shocks as they occur. To support decision-making, CBNs require a cause-and-effect map of the variables under consideration, known as the network's structure. Learning the graphical structure of a causal model from data remains challenging; learning it from time-series data is even harder because dependencies may arise at different time lags. Existing time-series causal discovery methods often assume a fixed lag window and do not explicitly optimise edge-specific lags. We propose a Tabu-based structure learning algorithm that searches for a time-ordered directed structure (i.e., where every edge respects time) while allowing edge-specific lags up to a specified maximum lag. The approach uses a decomposable BIC-based score with node-specific effective sample sizes and an explicit lag-length penalty encouraging parsimonious delay assignments while preserving efficient local score updates. We provide theoretical guarantees of validity and local optimality, and we also describe a parallel implementation for improved scalability. In simulations, the method recovered graph structure competitively and estimated lags accurately when true adjacencies were recovered. On a real-world UK COVID-19 policy dataset, the learnt structure was dominated by short delays while retaining a substantial minority of longer-lag dependencies, consistent with delayed behavioural and epidemiological effects.
Bruno Petrungaro, Anthony C. Constantinou
Bayesian Artificial Intelligence research lab, MInDS research group, Queen Mary University of London, Mile End Road, London, E1 4NS, UK.
Unobserved common causes are pervasive in real-world time series and can induce spurious associations that causal discovery methods mistake for direct edges. We propose LUCID (Learning Under Confounding for Inference and Discovery, a regime-adaptive deconfounding layer that first estimates the confounding regime from data using a Marčenko--Pastur spectral router, then applies a deconfounding strategy matched to that regime. When the spectrum indicates pervasive factor confounding, LUCID attenuates factor-dominated variation and recovers contemporaneous (lag-0) structure from the resulting innovations, with edge selection calibrated against a data-driven edge-free null. Rather than being tied to a particular discovery algorithm, it can wrap existing discovery engines; we demonstrate consistent improvements across three such methods. On a diverse synthetic out-of-distribution benchmark spanning changes in confounder strength and sparsity, loading density, lag structure, volatility dynamics, edge heterogeneity, persistence, intermittency, and tail behavior, LUCID achieves the best family-weighted directed, lag-resolved graph F1 (0.60), improving over the strongest baseline by 0.19 absolute (≈46% relative). Its advantage widens relative to looser lag-collapsed scoring, and remains robust under intermittent and heavy-tailed confounding. Code reproducing the method, the benchmark generators, and every reported experiment is available at https://github.com/bloomberg/causal-ts.
Causal discovery with instantaneous effects in multivariate time series is challenging, as the instantaneous structure must be acyclic. Prior methods enforce this by either separating instantaneous and lagged estimation into multi-stage pipelines or imposing algebraic acyclicity constraints via complex augmented Lagrangian optimization, both of which incur high computational cost. In this work, we propose a different approach: we learn a differentiable permutation of variables using the Gumbel--Sinkhorn operator and triangularize the instantaneous coefficient matrix of a Structural Vector Autoregressive (SVAR) model in the learned order. This converts acyclicity from a hard constraint into a parameterization and keeps it valid throughout optimization. In doing so, our method enables unified, continuous optimization with gradient-based learning, leading to improved efficiency in time--series causal discovery. Across three real-world benchmarks, our method achieves the best overall performance compared with 12 baselines in both discovery accuracy and efficiency. On the large-scale benchmark, it further demonstrates strong scalability, achieving more than a 6x speedup over competing methods.