arXiv CSopen access

Differentiable latent structure discovery for interpretable forecasting in clinical time series

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

machine learning, deep learning, neural networks

Differentiable latent structure discovery for interpretable forecasting in clinical time series Ivan Lerner1,2*, Jean Feydy1 , Alexandre Kalimouttou2 , Anita Burgun2 , Francis Bach3,4,5

arXiv:2604.27967v1 [cs.LG] 30 Apr 2026

Université Paris Cité and Université Sorbonne Paris Nord, INRIA, Inserm, HEKA, Paris, 75015, France. 2 Department of Medical Informatics, Assistance Publique Hôpitaux de Paris (AP-HP), Georges Pompidou European Hospital, Paris, France. 3 SIERRA, Inria Paris, Paris, 75015, France. 4 Département d’Informatique de l’École Normale Supérieure, France. 5 PSL Research University, Paris, France.

*Corresponding author(s). E-mail(s): [email protected]; Abstract ∗

Background: Timely, uncertainty-aware forecasting from irregular electronic health records (EHR) can support critical-care decisions, yet most approaches either impute to a grid or sacrifice interpretability. We introduce StructGP, a continuous-time multi-task Gaussian process that couples process convolutions with differentiable structure learning to uncover a sparse, ordered directed acyclic graph (DAG) of inter-variable dependencies while preserving principled uncertainty. We further propose LP-StructGP, which augments StructGP with latent pathways—shared, temporally shifted trajectories inferred via subject-specific coupling filters and a softmax gating mechanism—to capture cross-patient progression patterns. Both models are trained under sparsity and acyclicity constraints (augmented Lagrangian, Adam) using scalable low-rank updates: StructGP via exact marginal likelihood maximization and LP-StructGP via an online conditional pseudo-marginal likelihood. Results: In simulations, the approach reliably recovers ground-truth graphs (Structural Hamming Distance approaching 0 as cohorts grow) and pathway assignments (high Adjusted Rand Index). On a MIMIC-IV septic shock cohort (n=1,008; norepinephrine, creatinine, mean arterial pressure), StructGP ∗ This manuscript is under review at BioData Mining.

improves short-horizon (6 h) forecasting over independent-task baselines (average RMSE 0.68 [95%CI: 0.63–0.74] vs. 0.88 [0.83–0.94]) and, with 15 additional inputs, markedly outperforms unstructured kernels (0.63 [0.58–0.69] vs. 3.02 [2.85–3.18]) with superior calibration (coverage 0.96 vs. 0.84). For long horizons (up to 6 days), LP-StructGP further reduces error for creatinine and norepinephrine (e.g., creatinine RMSE 0.95 [0.87–1.05] vs. 1.17 [1.08–1.25]) and improves overall coverage (0.94 vs. 0.91). On the PhysioNet Challenge (12k patients, 41 variables), StructGP attains competitive accuracy (MAE 3.72×10−2 ) relative to a state-of-the-art graph neural model while maintaining calibrated uncertainty. Conclusion: These results show that structured process convolutions with latent pathways deliver interpretable, scalable, and well-calibrated forecasting for irregular clinical time series. Keywords: Clinical time series , Multi-task Gaussian Process , Differentiable structure learning

1 Introduction Clinical decision-making relies not only on assessing a patient’s current state but also on anticipating its future course. Care plans, for example, are adjusted when renal function declines, and more aggressive interventions may be required if cardiovascular stability cannot be restored. Predictive models that provide model-based, data-driven forecasts can therefore play a crucial role in supporting clinicians. Their potential impact is considerable: up to 11% of deaths in critically ill patients are attributable to delayed care [1], underscoring the importance of timely and accurate predictions. Electronic health records (EHRs) capture a wide range of information during hospital admissions, including sensor measurements, biomarkers, and treatments. Excluding unstructured text and imaging data, these records can be viewed as collections of multivariate time series, one per patient’s stay. Unlike controlled experimental settings, clinical observations are collected opportunistically rather than on a regular grid, making them difficult to model with classical frameworks (e.g., vector autoregressive models). The resulting patterns are irregular and heterogeneous (Supplementary Figure S2), with variables measured at different times, frequencies, and subsets of patients. A common strategy to address irregular sampling is to resample the data onto a fixed grid and impute the missing values, after which standard time-series methods can be applied [2, 3]. While simple, this approach is computationally intensive and introduces potential biases through resampling and imputation. A more principled alternative is to work directly in continuous time, thereby avoiding these intermediate steps. Continuous-time modeling can be approached through a variety of frameworks: multi-task Gaussian processes (MTGPs) [4–11]; neural ODEs and related stochastic differential equation models [12–15]; adaptations of deep learning architectures, such as recurrent [16] or attention-based networks [17–20]; and, more recently, graph neural networks (GNNs) [21–23]. Although GNN-based methods offer some interpretability

and, in the case of GraFITi [23], state-of-the-art predictive performance, they generally lack principled uncertainty quantification, a key strength of MTGPs. In this work, we revisit the classical framework of process convolutions [24–26], which define the covariance structure of a multi-task Gaussian process (MTGP) through kernel convolutions. We combine this framework with modern differentiable structure learning [27, 28], enabling the discovery of interpretable dependency structures between time series in continuous time. A preliminary version of this approach was presented at the Bayesian Decision-making and Uncertainty Workshop (NeurIPS 2024) [29]. The present work extends it by introducing internal standardization, latent pathways, an online learning algorithm, and an extensive biomedical evaluation. In our formulation, the sparsity pattern of the cross-task convolution filters is represented as a directed acyclic graph (DAG), encoding ordered conditional independence relations and causal effects among time series. The corresponding adjacency matrix and variable ordering can be inferred directly from data by solving a continuous optimization problem [27, 28]. Compared to classical MTGPs, our approach provides both flexible probabilistic modeling and interpretable, actionable dependency structures between time series. In addition, we address the question of modeling dependencies across patients by introducing the concept of a latent pathway: a partially shared trajectory that groups of patients may follow. This assumption is broad enough to capture diverse phenomena, ranging from different patients experiencing the same underlying pathological process (e.g., urinary tract– versus lung infection–induced sepsis) to patients undergoing similar treatment strategies (e.g., progressive catecholamine withdrawal). At test time, our model infers both the latent pathway along which a patient is evolving and their position within that pathway. In contrast to previous approaches that tackled patient-level dependencies [8–10], we integrate latent pathways directly into the covariance design. This yields calibrated uncertainty estimates while preserving interpretable dependency structures across tasks. With this study, our objectives were to design, implement, and evaluate a process convolution model, StructGP, along with an extension that incorporates latent pathways, LP-StructGP, together with their associated learning algorithms. We assessed their properties using both realistic simulations and real-world data, specifically a MIMIC-IV septic shock cohort from a previous study [30, 31]. For comparison with state-of-the-art deep learning methods, we further evaluated our models on the PhysioNet Challenge [32], using the preprocessing pipeline described in [18]. Our contributions are thus twofold: (i) StructGP, which unifies process convolutions with differentiable structure learning, and (ii) LP-StructGP, which extends this framework to capture patient-level latent pathways.

2 Methods 2.1 Process convolution models We begin by recalling the process convolution formulation of Gaussian processes [33], which provides a flexible framework for constructing structured covariance functions. Let y (x) = {y (x) : x ∈ X } denote a stochastic process, i.e., a collection of random 3

variables indexed by x ∈ X . Process convolution models [24, 25] adopt a signalprocessing perspective, where y (x) is viewed as the output of a linear time-invariant (LTI) system driven by a latent noise process. An LTI system is entirely characterized by its impulse response function h(x), which represents the system’s response to a Dirac input. If the input is a Gaussian white noise process w(ξ ) with covariance E[w(ξ )w(ξ ′ )] = δ (ξ − ξ ′ ), then by the closure of Gaussian processes under linear operations, y (x) is also a Gaussian process provided that h has compact support or is square-integrable:

Z ∞ y (x) =

h(x − ξ )w(ξ ) dξ = (h ∗ w)(x),

(1)

h(x − ξ )h(x′ − ξ ) dξ = (h ∗ h)(x − x′ ),

(2)

−∞

k (x, x′ ) =

Z ∞ −∞

where k (x, x′ ) denotes the covariance function of y (x). By appropriately designing these filters, one can derive analytical covariance functions that naturally capture the smoothness and cross-dependence of observed clinical time series, while remaining well-defined under irregular sampling.

2.2 Structured process convolution models 2.2.1 StructGP: an inter-task graphical structure We now consider StructGP, a k -dimensional Gaussian process Y(t) obtained by convolving a multivariate white noise process w(t) with a matrix-valued filter H(t):

Y(t) = (H ∗ w)(t), ′

⊤

(3) ′

Cov[Y(t), Y(t )] = (H ∗ H )(t − t ),

(4)

where w(t) is a k -dimensional Gaussian white noise process and H(t) is a sparse, lowertriangular, matrix-valued impulse response function. This formulation generalizes the univariate process convolution to a structured multi-task setting. The notion of conditional independence between stationary Gaussian process components can be characterized in the frequency domain through their spectral density matrix [34–37]. In our previous work [29], we showed how this concept enables learning ordered conditional relations between time series directly from data. Accordingly, we parameterize the filter matrix as

H(t) = (I + S) ◦ L(t),

(5)

where ◦ denotes the Hadamard product, S is a sparse lower-triangular matrix (up to permutation) encoding inter-task dependencies, and L(t) is a positive matrix of filter coefficients with elements t2 Lvu (t) = exp − , (6) ℓvu 4

with ℓvu > 0 denoting learnable length scales. The nonzero entries of S define the adjacency matrix of a directed acyclic graph (DAG) G , which encodes ordered conditional dependencies among tasks:

Svu ̸= 0

⇔

Yu → Yv

⇔

Yu ̸⊥⊥ Yv | Y1:u−1 .

(7)

Each filter element Hvu (t) thus represents the continuous-time influence of process Yu on process Yv , yielding an interpretable, structured covariance that captures directed inter-task relations (see Additional file 1, Section S1.1). Finally, in order to detach the signal structure from the scale of each variable, we follow [38] and introduce an internal standardization step such that the marginal variance of each process is always equal to one (see Additional file 1, Section S1.2).

2.2.2 LP-StructGP: an inter-subject latent pathway structure

Fig. 1 Sketch of the latent pathway model for one task

To capture shared temporal dynamics across subjects, we extend StructGP with an additional latent layer representing group-level trajectories, referred to as latent pathways (see Figure 1). Each latent pathway corresponds to a multivariate Gaussian process

Qu (t) = (H ∗ wu )(t),

(8)

where wu (t) ∈ Rk is a latent multivariate white noise process and H(t) is the shared inter-task filter. Individual trajectories are modeled as mixtures of these latent 5

pathways, together with a subject-specific component:

Yi (t) =

p X

(Giu ∗ Qu )(t) + (H ∗ wi )(t),

(9)

u=1

where Giu (t) is a subject-specific coupling filter. To account for differences in disease progression, Giu (t) includes a temporal shift parameter, while subject membership to each pathway is controlled through a softmax-based gating mechanism [39, 40]. This formulation allows subjects to be probabilistically assigned to shared latent trajectories while preserving subject-specific variability. Full mathematical details and the corresponding covariance expressions are provided in Additional file 1, Section S1.3.

2.3 Learning and inference As the time series are irregularly sampled across patients, we adopt a set-based indexing scheme that treats all observations as elements of a common collection rather than as regularly spaced sequences. Owing to the marginalization properties of multivariate Gaussians, the covariance structure defined above remains valid under this irregular sampling. This allows us to directly apply standard Gaussian process learning and inference procedures: (i) maximizing the process marginal likelihood to learn the convolution filter parameters, and (ii) computing posterior predictive distributions at test points (see Additional file 1, Section S1.4).

2.3.1 Differentiable structure learning To uncover interpretable dependencies between tasks, we learn the adjacency matrix S under the constraint that it forms a weighted directed acyclic graph (DAG) G . We adopt the differentiable structure learning framework of NOTEARS [27], which enforces acyclicity through a smooth matrix function based on the trace of the matrix exponential (see Additional file 1, Section S1.5.1). The parameters θ = {S, L} are estimated by minimizing the process negative marginal likelihood with sparsity and acyclicity constraints:

h i θ∗ = argmin − log L(y, X, θ) + Pλ (S) s.t. h(S) = 0,

(10)

where Pλ is a sparsity penalty and h(S) enforces the DAG constraint. Optimization is carried out with the augmented Lagrangian method [41], using Adam [42] as the gradient-based optimizer. A smooth approximation of the ℓ1 norm ensures differentiability everywhere, and the sparsity weight λ is selected by minimizing an AIC-like or validation criterion. Finally, a hard threshold is applied to S to obtain an acyclic structure. Full algorithmic details are provided in Additional file 1, Section S1.5.

2.3.2 Low-rank plus block-sparse scheme for scalability A key limitation of Gaussian process regression is its O(n3 ) computational complexity, with n the number of observations, due to the covariance matrix inversion. To scale to 6

large cohorts, we exploit independence across subjects and implement all models with batched Cholesky solvers in PyTorch, leveraging GPU-based parallel computation. The resulting complexity of matrix inversion scales with the number of observations per subject, making this approach efficient in regimes with many subjects but fewer than ∼1000 observations per subject. Independence across subjects further allows scaling through mini-batch gradient steps (e.g., using the Adam optimizer). In StructGP, subjects are independent by design, and mini-batches are constructed by padding subjects with different sequence lengths. By contrast, LPStructGP couples subjects through shared latent pathways. In the online implementation, we replace the full-data negative marginal log likelihood (NMLL) by an online conditional pseudo-marginal likelihood (see Additional file 1, Section S1.8), which conditions on the current empirical estimate of the latent pathways. This surrogate objective can be interpreted as an approximation to the true marginal likelihood in which the posterior over the shared latent process is replaced by its current empirical mean. Conditioning on this fixed estimate makes subjects effectively independent, yielding a tractable optimization that becomes increasingly accurate as the posterior over the latent pathways concentrates with more data. In practice, this is achieved by applying a low-rank approximation to the intersubject covariance using the Hilbert Space Gaussian Process (HSGP) method [43, 44] (see Additional file 1, Section S1.6). The resulting kernel follows a low-rank plus block-sparse structure, enabling efficient inference via the Woodbury matrix identity. Matrix inversion and log-determinant computations required for marginal likelihood evaluation and posterior inference are performed using a batch online algorithm with geometrically discounted accumulators (see Additional file 1, Section S1.7 and S1.8). The main accumulator, which acts as a memory of the latent pathways accumulated across subjects, has dimension mk × mk , where m is the number of basis functions in the HSGP approximation and k is the number of tasks. Its inclusion in the model therefore scales computationnaly as O((mk )3 ).

2.4 Simulation study We evaluated the ability of StructGP and LP-StructGP to recover both the inter-task graph and the latent pathway structure from simulated data. For each simulation, a ground-truth graph, covariance parameters, and observations were sampled from the corresponding generative model (see 2.2). Models were then fitted using the learning procedure described in 2.3, and the inferred structures were compared to the ground truth. Graph recovery accuracy was assessed using the Structural Hamming Distance (SHD), which counts the number of edge insertions, deletions, or reversals needed to match the true graph. Latent pathway assignments were evaluated with Adjusted Rand Index (ARI) and Normalized Mutual Information (NMI). Each setting was repeated 30 times, and we report median and interquartile ranges across runs. Full simulation settings, parameter sampling schemes, and additional results are provided in Additional file 1, Section S1.9 and in Lerner et al. (2024) [29].

2.5 Real-life evaluation 2.5.1 MIMIC-IV septic shock cohort We evaluated StructGP and LP-StructGP on a septic shock cohort extracted from the MIMIC-IV database [30, 31]. MIMIC-IV is a publicly available collection of deidentified ICU records from the Beth Israel Deaconess Medical Center (Boston, US) between 2008 and 2019. Septic shock was defined according to Sepsis-3 criteria [45]. Only the first ICU stay per patient was retained. Detailed inclusion criteria and preprocessing steps are provided in Additional file 1, Section S1.10. We aimed to jointly forecast creatinine (Cr), norepinephrine (NE) dose, and mean arterial pressure (MAP) over the first seven ICU days. Two forecasting regimes were considered: (i ) short-term, where within each 24-hour window the model predicts the last 6 hours from the preceding 18 hours of observations, and (ii ) long-term, predicting the remainder of the stay from the first 24 hours. All variables were standardized using a quantile transform fitted on the training set. Data were split by subject into training, validation, and test sets (70%/15%/15%). We compared kernel architectures (StructGP, LP-StructGP, and an independenttask baseline) under identical hyperparameter settings. Performance was assessed using the Root Mean Square Error (RMSE) and the empirical 95% coverage of the predictive distribution. Additional preprocessing and feature-engineering details are given in Additional file 1, Section S1.10. To assess scalability and multi-task generalization, we conducted additional large-scale experiments including 18 physiological and treatment variables spanning cardiovascular, respiratory, renal, and metabolic systems (see Additional file 1, Section S1.10 for the complete list).

2.5.2 PhysioNet Challenge dataset We further evaluated StructGP on the PhysioNet Challenge dataset [32], following the preprocessing pipeline of Zhang et al. [18]. The dataset includes approximately 12,000 ICU patients, each with 41 irregularly sampled clinical variables recorded during the first 48 hours after admission. For each patient, the first 24 hours were used as observations to predict values over the following 24 hours. Data were split into training, validation, and test sets (60%/20%/20%). For comparison, we also refitted the GraFITi model [23], which currently achieves state-of-the-art performance on this benchmark.

3 Results 3.1 Simulation study

Fig. 2 Accurate structure recovery with increasing number of subjects.

Simulation parameters are detailed in Additional file 1, Section S1.9. F1: harmonic mean of precision and recall for edge recovery. SHD: Structural Hamming Distance. ARI: Adjusted Rand Index. NMI: Normalized Mutual Information.

Our simulations (Figure 2) show that the learning algorithm accurately recovers the latent structures of StructGP and LP-StructGP under realistic, albeit low-noise, parameter regimes. Performance improves with the number of subjects, achieving perfect graph recovery (SHD = 0) and consistent identification of latent pathways (high ARI and NMI). Supplementary Figure S3 shows how the predictive mean initially combine subject-specific and group-level components, before converging to the group-level forecast at longer horizons in a simulated cohort.

3.2 MIMIC-IV septic shock cohort We evaluated short-term forecasting on a MIMIC-IV septic shock cohort comprising 1,008 ICU stays, split into training (705), validation (151), and test (152) subjects. Each record included three target variables—norepinephrine administration rate (34,411 observations), serum creatinine (8,959), and mean arterial pressure (109,635)—measured over the first seven ICU days (time window 0–7 days). Measurements were sampled with a minimum sampling interval of one hour. When including 15 additional physiological and treatment variables, the dataset totaled 511,186 observations. All time series were standardized after preprocessing as described in Section 2.5.1.

3.2.1 Short-term forecast Models Average

Output Tasks Cr

MAP

RMSE [95% CI] Input Tasks: NE, Cr, MAP No Structure Independent Tasks StructGP StructGP + 2h lag features StructGP + task-spe noise LP-StructGP Input Tasks: + 15 tasks No Structure StructGP StructGP + 2h lag features StructGP + task-spe noise

0.65 [0.60, 0.71] 0.88 [0.83, 0.94] 0.68 [0.63, 0.74] 0.65 [0.60, 0.71] 0.63 [0.57, 0.68] 0.66 [0.60, 0.71]

0.67 [0.61, 0.73] 0.76 [0.70, 0.82] 0.65 [0.60, 0.71] 0.66 [0.61, 0.72] 0.58 [0.53, 0.64] 0.63 [0.58, 0.69]

0.43 [0.26, 0.58] 0.94 [0.79, 1.11] 0.44 [0.30, 0.58] 0.46 [0.32, 0.63] 0.43 [0.26, 0.60] 0.44 [0.30, 0.60]

0.87 [0.85, 0.90] 0.95 [0.93, 0.98] 0.96 [0.93, 0.98] 0.84 [0.82, 0.87] 0.86 [0.84, 0.88] 0.90 [0.88, 0.92]

3.02 [2.85, 3.18] 0.63 [0.58, 0.69] 0.66 [0.60, 0.72] 0.62 [0.57, 0.67]

2.84 [2.67, 3.03] 0.64 [0.58, 0.70] 0.65 [0.60, 0.71] 0.61 [0.55, 0.67]

3.25 [2.80, 3.77] 0.42 [0.28, 0.59] 0.42 [0.27, 0.59] 0.38 [0.26, 0.54]

2.97 [2.85, 3.09] 0.84 [0.81, 0.86] 0.91 [0.88, 0.93] 0.85 [0.83, 0.87]

Coverage [95% CI] Input Tasks: NE, Cr, MAP No Structure Independent Tasks StructGP StructGP + 2h lag features StructGP + task-spe noise LP-StructGP Input Tasks: + 15 tasks No Structure StructGP StructGP + 2h lag features StructGP + task-spe noise

0.96 [0.95, 0.97] 0.95 [0.93, 0.96] 0.97 [0.96, 0.98] 0.96 [0.95, 0.97] 0.96 [0.95, 0.97] 0.97 [0.96, 0.97]

0.95 [0.94, 0.96] 0.97 [0.96, 0.98] 0.96 [0.95, 0.97] 0.96 [0.94, 0.97] 0.95 [0.93, 0.96] 0.96 [0.95, 0.97]

0.97 [0.95, 0.99] 0.95 [0.92, 0.98] 0.97 [0.94, 0.99] 0.96 [0.93, 0.98] 0.97 [0.95, 0.99] 0.97 [0.94, 0.99]

0.96 [0.95, 0.96] 0.92 [0.91, 0.93] 0.98 [0.98, 0.99] 0.96 [0.95, 0.97] 0.96 [0.95, 0.96] 0.97 [0.97, 0.98]

0.84 [0.82, 0.86] 0.96 [0.95, 0.97] 0.95 [0.94, 0.97] 0.96 [0.95, 0.97]

0.86 [0.84, 0.88] 0.96 [0.95, 0.97] 0.95 [0.94, 0.96] 0.93 [0.92, 0.95]

0.83 [0.77, 0.88] 0.97 [0.95, 0.99] 0.96 [0.92, 0.98] 0.98 [0.97, 1.00]

0.84 [0.82, 0.85] 0.96 [0.95, 0.96] 0.96 [0.95, 0.97] 0.96 [0.95, 0.97]

Table 1 Short-term (6 hour) horizon forecasting Root mean squared error (RMSE) and predictive coverage with bootstrap 95% confidence intervals. Metrics are reported over a six-hour forecasting horizon for Norepinephrine (NE), Creatinine (Cr), and mean arterial pressure (MAP), using 18 hours of observations within overlapping 24-hour ICU stay windows, with and without 15 additional tasks.

In the short-term forecasting task, StructGP outperformed the independent-task baseline, achieving an average RMSE of 0.68 [95% CI: 0.63–0.74] compared to 0.88

[0.83–0.94] (Table 1). With three input tasks (NE, Cr, MAP), the explicit graphical structure provided no benefit over an unconstrained model (average RMSE 0.68 [0.63–0.74] vs. 0.65 [0.60–0.71]). However, when 15 additional input tasks were included, the structured formulation became essential, yielding a markedly lower average RMSE (0.63 [0.58–0.69]) compared to the unstructured variant (3.02 [2.85–3.18]), and substantially better calibration (0.96 [0.95–0.97] vs. 0.84 [0.82–0.86]). Finally, relaxing the assumption of shared observation noise and learning task-specific noise variances consistently improved performance across all settings.

Fig. 3 Graphical representation of the peak total effect Total effect is the sum of direct and indirect effect. Graph representing the models reported in Table 1. A: StructGP. B: StructGP + 2h lag features. C: StructGP with 15 additional tasks (see Additional file 1, Section S1.10). Black and red edges represent positive and negative influences, respectively, and edge thickness indicates the strength of the effect.

In Figure 3, we illustrate the learned dependency structures for three StructGP variants reported in Table 1. In panel A, a red link between MAP and NE, indicate that an increase in MAP is associated with a decrease in NE, whereas a black link between Cr and NE indicates that an increase in Cr is associated with a increase in NE. In panel B, a model with additional feature tasks, where “1” denotes a constant (intercept) pseudo-task, and “NE–2h” represents the 2-hour lagged version of NE. Here, NE–2h is associated with increases in both Cr and NE.

Models Average

Output Tasks Cr

MAP

RMSE [95% CI] Input Tasks: NE, Cr, MAP StructGP LP-(fixed)StructGP LP-StructGP

1.04 [1.00, 1.07] 1.04 [1.01, 1.07] 0.97 [0.94, 1.00]

0.97 [0.94, 1.00] 0.97 [0.94, 1.00] 0.93 [0.90, 0.96]

1.17 [1.08, 1.25] 1.16 [1.08, 1.25] 0.95 [0.87, 1.05]

0.97 [0.96, 0.98] 0.97 [0.96, 0.98] 1.01 [1.00, 1.02]

Coverage [95% CI] Input Tasks: NE, Cr, MAP StructGP LP-(fixed)StructGP LP-StructGP

0.91 [0.91, 0.92] 0.87 [0.86, 0.87] 0.94 [0.94, 0.95]

0.89 [0.88, 0.90] 0.83 [0.82, 0.85] 0.93 [0.92, 0.94]

0.90 [0.89, 0.92] 0.86 [0.84, 0.87] 0.94 [0.92, 0.95]

0.95 [0.94, 0.95] 0.91 [0.90, 0.91] 0.96 [0.96, 0.97]

Table 2 Long-term (six-day) horizon forecasting Root mean squared error (RMSE) and predictive coverage with bootstrap 95% confidence intervals. Metrics are reported for up to a six-day forecasting horizon for Norepinephrine (NE), Creatinine (Cr), and mean arterial pressure (MAP), after the first 24 hours of the ICU stay.

3.2.2 Long-term forecast Table 2 compares StructGP with its latent-pathway extensions, LP-(fixed)StructGP and LP-StructGP, using three latent pathways. In LP-(fixed)StructGP, the convolution filter parameters H were kept fixed to those learned by StructGP for computational efficiency. LP-StructGP achieved the best long-term performance, improving forecasts for creatinine (RMSE 0.95 [0.87–1.05] vs. 1.17 [1.08–1.25]) and norepinephrine (0.93 [0.90–0.96] vs. 0.97 [0.94–1.00]). Predictive calibration also improved, with coverage increasing from 0.91 [0.91–0.92] to 0.94 [0.94–0.95]. Visualization of the learned latent pathways (Figure 4) revealed two broad temporal patterns, primarily differing in their creatinine trajectories. Qualitatively, these latent pathways appear to improve long-term forecasts for certain patients, as illustrated for two representative cases in Figure 5.

3.3 PhysioNet Challenge dataset

GraFITi (SOTA) StructGP StructGP + 3lat + time feat.

MAE ×10−2 [95% CI]

MSE ×10−3 [95% CI]

3.50 [3.39, 3.62] 4.06 [3.92, 4.21] 3.72 [3.59, 3.86]

4.79 [4.27, 5.58] 7.55 [6.67, 8.48] 6.74 [6.05, 7.63]

Table 3 Comparison with SOTA on the PhysioNet challenge Mean absolute error (MAE) and Mean squared error (MSE) with bootstrap 95% confidence intervals. Metrics are reported for up to a 1-day forecasting horizon macro-average for 41 tasks, after the first 24 hours of the ICU stay.

The dataset comprises 41 time series from a total of 12,000 subjects (7,200 for training, 2,400 for development, and 2,400 for testing), with approximately 5.26 million observations in total (3,149,079 for training, 1,054,084 for development, and

Fig. 4 Latent patient pathways visualization Learned from the LP-StructGP model reported in Table 2. Each trajectory represents one inferred latent pathway across the 7-day ICU stay.

Fig. 5 Latent patient pathways six-day horizon forecasts Forecasts produced by the LP-StructGP model from Table 2, illustrating subject-level predictions along the inferred latent pathways.

1,056,460 for testing). Table 3 compares StructGP with the state-of-the-art model on this dataset, GraFITi [23]. The best-performing StructGP variant—incorporating time-feature pseudo-tasks and additional flexibility via three latent white-noise processes per task with a nested graphical structure—achieved an average MAE of 3.72 [3.59–3.86] × 10−2 , compared to 3.50 [3.39–3.62] × 10−2 for GraFITi. However, the gap was larger in MSE, with StructGP reaching 6.74 [6.05–7.63] × 10−3 versus 4.79 [4.27–5.58] × 10−3 . A detailed breakdown of per-variable performance is provided in Supplementary Table S1, listing the 18 tasks where StructGP underperformed, and Supplementary Table S2, which reports the 19 tasks where StructGP outperformed GraFITi. For context, all other models evaluated by Yalavarthi et al. [23]—except Zhang et al. [18] (MAE = 3.72 × 10−2 )—performed worse than both GraFITi and StructGP. This includes NeuralFlow [15] (MAE = 4.67 × 10−2 ), Latent-ODE [12] (4.23 × 10−2 ), and CrossFormer [17] (4.81 × 10−2 ) 1 . Moreover, when restricting evaluation to the top tercile of predictions with highest confidence (lowest posterior variance), StructGP slightly outperformed GraFITi, with an MAE of 2.46 [2.34–2.60] × 10−2 versus 2.53 [2.44–2.64] × 10−2 (Supplementary Table S3).

4 Discussion We developed StructGP, a process-convolution model with an inter-task differentiable graphical structure, and LP-StructGP, an extension that embeds inter-subject latent pathways directly within the covariance. StructGP enables learning a shared, ordered dependency structure across time series in a differentiable manner, enhancing forecasting performance as the number of time series increases. LP-StructGP identifies latent trajectory motifs shared across patient subgroups, yielding more accurate longhorizon forecasts. In realistic simulations, both the inter-task graph and pathway assignments were reliably identifiable from data. On a MIMIC-IV septic shock cohort, incorporating structure became increasingly valuable with higher task dimensionality, improving both error and calibration; over longer horizons, latent pathways further enhanced creatinine and norepinephrine forecasts and coverage. On the PhysioNet Challenge, StructGP achieved accuracy close to the current state of the art while maintaining calibrated uncertainty. Together, these contributions support interpretable and transparent transfer across tasks and patients, advancing an uncertainty-aware, self-supervised structure learning framework for medium-sized, noisy clinical datasets. As shown in Figure 2, the proposed models accurately recovered the underlying graph structure—even in the presence of latent pathways—with as few as approximately 1,000 subjects. Although the simulated graphs were relatively small (10 nodes) and sparse (mean degree of 2, i.e., on average two incoming or outgoing edges per node), this represents a notable result given that observations were irregularly sampled, and that the models were internally standardized—preventing them from exploiting scale information to recover variable ordering. In addition, the latent pathway (cluster) assignments were jointly recovered with high reliability. 1

The complete benchmark results are available at https://github.com/yalavarthivk/GraFITi/tree/main.

As shown in Table 1, incorporating the structural constraint improved forecasting performance in the setting with 18 tasks. This gain likely arises from the regularizing effects of the acyclicity and sparsity constraints. It may also reflect the invariance properties conferred by identifying a latent structure that approximates the true causal graph [46]. In Figure 3, the model identified an association in which an increase in mean arterial pressure (MAP) was linked to a decrease in norepinephrine (NE), whereas physiologically NE increases MAP. This apparent reversal is in fact consistent with the difference in timescales: NE’s hemodynamic effect occurs over minutes [47], whereas our dataset is sampled hourly, capturing instead the clinical adjustment of NE dosage in response to changes in MAP. Together, these findings highlight the importance of allowing models to represent different directional effects at different temporal resolutions or time shifts. A simple strategy is to include lagged “feature tasks,” as illustrated in Panel B of Figure 3, where a two-hour lag of NE (NE–2h) leads to an increase in creatinine—consistent with the known renal side effects of NE—while a contemporaneous link shows that increases in creatinine are associated with decreases in NE, potentially reflecting real corrective clinical interventions. For long-term forecasts, LP-StructGP outperformed StructGP on the creatinine (Cr) and norepinephrine (NE) tasks (Table 2). It also revealed latent pathways characterized by increasing mean arterial pressure (MAP), concurrent decreases in NE, and variable Cr trajectories—either gradually decreasing or decreasing followed by a late increase (Figure 4), consistent with different patterns of septic shock recovery. Interestingly, the task dependency structure was shared across latent pathways, preserving interpretability within this more flexible model.

4.1 Related Works Compared with leading graph neural network models such as GraFITi [23], the best StructGP variant underperformed on 18 tasks and outperformed on 19 tasks based on the MAE metric. In addition, its predictions provide calibrated uncertainty estimates and remain interpretable, confirming convolution process models as a promising direction for continuous-time modeling of clinical time series. From a scalability perspective, the current implementation is tailored for large cohorts with many patients but relatively few observations per patient. Under this regime, StructGP successfully fit 5.26 million observations from the PhysioNet dataset (41 tasks) in approximately 45 minutes on a single NVIDIA A40 GPU, including 100 grid-search iterations and repeated minimization of the NMLL to satisfy the acyclicity constraint at each step (Additional file 1, Section S1.5). Future work could extend scalability to settings with longer individual trajectories or larger numbers of tasks while remaining within an exact GP framework, by leveraging dedicated block-sparse routines for Gaussian process solvers [48]. From a structure learning standpoint, StructGP can be viewed as a continuous-time extension of DYNOTEARS [28], a necessary step for modeling clinical time series. Dallakyan et al. [49] proposed a continuous-time structure learning approach based on a Gaussian linear additive SCM in the Fourier domain. Unlike their method, which estimates separate weight matrices at each frequency, StructGP assumes a single invariant structure across frequencies. 15

Gaussian processes (GPs) have long been applied to clinical time series. Early work by Durichen et al. [4] introduced process convolution kernels with time-shift parameters, and later studies addressed population- versus patient-level modeling [8– 11]. These approaches typically separate inter-task and inter-patient dependencies or rely on fixed hierarchical structures. In contrast, StructGP jointly learns both withinand between-patient relationships through a unified covariance model optimized via (approximate) marginal likelihood.

5 Limitations and Perspectives The main limitation of the current implementation is that, since this study focused on developing the model’s structural assumptions, we restricted our experiments to Gaussian likelihoods. Although the data were standardized to roughly satisfy this assumption, more flexible likelihood formulations [50–53] would be necessary to capture the heterogeneity of clinical time series and to jointly model discrete patient states, such as intubation status, missingness indicators, or mortality. From a structure learning perspective, however, the fully Gaussian case remains the most challenging, as it requires strong assumptions on variances for identifiability [54, page 138]. Our learned graphs represent total-effect structures, where nonzero amplitude parameters in the transfer function matrix correspond to the combined direct and indirect effects between tasks. Because different length-scale parameters are allowed for each task, improving flexibility, the model no longer distinguishes direct from indirect effects. Nonetheless, this design has the desirable property that, even in large graphs, relationships do not become nearly deterministic, as occurs in structural causal models [55]. From a design standpoint, LP-StructGP currently faces scalability limitations. In practice, one must choose either a small number of latent pathways or a small number of tasks, since the feature dimension of the inter-subject component scales as m × k × p, where m is the number of basis functions in the HSGP approximation, k the number of tasks, and p the number of latent pathways. A simple remedy is to assign latent pathway components only to a subset of nodes (e.g., root nodes). An even more appealing direction would be to restrict latent pathways to latent unobserved tasks, although inference and identifiability in that setting would be more challenging. The group-level component of the model provides greater flexibility than the individual component in extrapolation settings, as it interpolates time series across training patients to identify shared patterns and trends. However, it remains limited to patients following a single latent pathway within the observed time period. Introducing additional patient-specific parameters at selected time points could enable transitions between different dynamical motifs. Further experiments are needed to fully characterize the model’s properties, including its causal properties, robustness to missingness patterns, transfer learning capabilities, and comparative performance across heterogeneous datasets. Nevertheless, this study introduces, for the first time, two continuous-time models that leverage latent structure in clinical data to improve forecasting, paving the way toward algorithmic approaches that bridge the gap from real-life data to knowledge.

6 Conclusion We developed StructGP, a process-convolution model that learns dependency structures among irregularly sampled time series, and its extension LP-StructGP, which captures dependencies across subjects through shared latent pathways. We showed that incorporating these structures in continuous-time models improves forecasting of real-world clinical time series while maintaining interpretability and calibrated uncertainty estimates. Future work should further characterize the conditions under which a causal interpretation is valid, paving the way toward a self-supervised framework of structured random functions that algorithmically bridges the gap from natural data streams to knowledge.

Acknowledgements. We would like to thank INRIA, INSERM, and AP-HP for their administrative, technical, and material support.

Declarations Ethics Permission to use the data was obtained for MIMIC-IV (No. 45398938). Because of the de-identified nature of the data, informed consent was waived. Consent was obtained for the original data collection.

Data availability This study was performed with the data from the MIMIC-IV version 1.0. Even though datasets are de-identified, restrictions have been imposed on data sharing since they contain sensitive information. Conventions are signed by researchers before any access to the data. For data access, interested researchers must fulfill all of the following requirements: be a credentialed user of https://physionet.org/, finish required training, and sign the data use agreement for the project. Reproducing the PhysioNet Challenge experiment with GraFiti model can be made easily at https://github.com/yalavarthivk/GraFITi.

Declaration of generative AI and AI-assisted technologies in the manuscript preparation process. Statement: During the preparation of this work the authors used GPT-5 (OpenAI) in order to proofread and reformulate sentences, paragraphs, or math formulas in latex. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the published article.

Contributions IL had full access to all of the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis. Concept and design: IL, FB. Acquisition, analysis, or interpretation of data: AK, IL, AB, JF, FB. Drafting of the manuscript: IL. Critical revision of the manuscript for important intellectual content: 17

IL, FB, JF, AB. Administrative, technical, or material support: AB. Supervision: FB, AB.

Funding This work was supported in part by the French government under the management of ANR as part of the “Investissements d’avenir” program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute).

Declaration of interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Additional file Additional file 1: Supplementary methods, supplementary results, supplementary figures, and supplementary tables. PDF file containing detailed mathematical derivations, algorithmic details, additional experimental details, Supplementary Figures S1–S3, and Supplementary Tables S1–S3.

References [1] Thomson R, Britain) NPSAG, Observatory PS. Safer Care for the Acutely Ill Patient: Learning from Serious Incidents. National Patient Safety Agency; 2007. Available from: https://books.google.fr/books?id=sSxyHwAACAAJ. [2] Lu F, Zheng Y, Cleveland H, Burton C, Madigan D. Bayesian hierarchical vector autoregressive models for patient-level predictive modeling. PLOS ONE. 2018 12;13(12):1–27. https://doi.org/10.1371/journal.pone.0208082. [3] Purushotham S, Meng C, Che Z, Liu Y. Benchmarking deep learning models on large healthcare datasets. Journal of Biomedical Informatics. 2018;83:112–134. https://doi.org/https://doi.org/10.1016/j.jbi.2018.04.007. [4] Dürichen R, Pimentel MA, Clifton L, Schweikard A, Clifton DA. Multitask Gaussian processes for multivariate physiological time-series analysis. IEEE Transactions on Biomedical Engineering. 2014;62(1):314–322. [5] Ghassemi M, Pimentel M, Naumann T, Brennan T, Clifton D, Szolovits P, et al. A multivariate timeseries modeling approach to severity of illness assessment and forecasting in ICU with sparse, heterogeneous clinical data. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 29; 2015. . [6] Zhang K, Karanth S, Patel B, Murphy R, Jiang X. A multi-task Gaussian process self-attention neural network for real-time prediction of the need for mechanical ventilators in COVID-19 patients. Journal of Biomedical Informatics. 2022;130:104079. https://doi.org/https://doi.org/10.1016/j.jbi.2022.104079. [7] Liu Z, Wu L, Hauskrecht M. Modeling clinical time series using gaussian process sequences. In: Proceedings of the 2013 SIAM International Conference on Data Mining; 2013. p. 623–631. [8] Liu Z, Hauskrecht M. Learning adaptive forecasting models from irregularly sampled multivariate clinical data. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 30; 2016. . [9] Cheng LF, Dumitrascu B, Darnell G, Chivers C, Draugelis M, Li K, et al. Sparse multi-output Gaussian processes for online medical time series prediction. BMC Medical Informatics and Decision Making. 2020;20(1):1–23. [10] Cui S, Yoo EC, Li D, Laudanski K, Engelhardt BE. Hierarchical Gaussian processes and mixtures of experts to model Covid-19 patient trajectories. In: Pacific Symposium on Biocomputing 2022. World Scientific; 2021. p. 266–277. [11] Karch JD, Brandmaier AM, Voelkle MC. Gaussian process panel modeling—machine learning inspired analysis of longitudinal panel data. Frontiers in Psychology. 2020;11:351.

[12] Rubanova Y, Chen RT, Duvenaud DK. Latent ordinary differential equations for irregularly-sampled time series. Advances in Neural Information Processing Systems. 2019;32. [13] De Brouwer E, Simm J, Arany A, Moreau Y. GRU-ODE-Bayes: Continuous modeling of sporadically-observed time series. Advances in neural information processing systems. 2019;32. [14] Scholz R, Born S, Duong-Trung N, Cruz-Bournazou MN, Schmidt-Thieme L.: Latent Linear ODEs with Neural Kalman Filtering for Irregular Time Series Forecasting. OpenReview preprint. [15] Biloš M, Sommer J, Rangapuram SS, Januschowski T, Günnemann S. Neural flows: Efficient alternative to neural ODEs. Advances in neural information processing systems. 2021;34:21325–21337. [16] Che Z, Purushotham S, Cho K, Sontag D, Liu Y. Recurrent neural networks for multivariate time series with missing values. Scientific reports. 2018;8(1):6085. [17] Zhang Y, Yan J. Crossformer: Transformer utilizing cross-dimension dependency for multivariate time series forecasting. In: The eleventh international conference on learning representations; 2023. . [18] Zhang W, Yin C, Liu H, Zhou X, Xiong H. Irregular multivariate time series forecasting: A transformable patching graph neural networks approach. In: Fortyfirst International Conference on Machine Learning; 2024. . [19] Liu Y, Hu T, Zhang H, Wu H, Wang S, Ma L, et al. itransformer: Inverted transformers are effective for time series forecasting. arXiv preprint arXiv:231006625. 2023;. [20] Zhang J, Zheng S, Cao W, Bian J, Li J. Warpformer: A multi-scale modeling approach for irregular clinical time series. In: Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining; 2023. p. 3273– 3285. [21] Wu Z, Pan S, Long G, Jiang J, Chang X, Zhang C. Connecting the dots: Multivariate time series forecasting with graph neural networks. In: Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining; 2020. p. 753–763. [22] Cao D, Wang Y, Duan J, Zhang C, Zhu X, Huang C, et al. Spectral temporal graph neural network for multivariate time-series forecasting. Advances in neural information processing systems. 2020;33:17766–17778.

[23] Yalavarthi VK, Madhusudhanan K, Scholz R, Ahmed N, Burchert J, Jawed S, et al. Grafiti: Graphs for forecasting irregularly sampled time series. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 38; 2024. p. 16255–16263. [24] Boyle P, Frean M. Dependent Gaussian processes. Information Processing Systems. 2005;17:217–224.

Advances in Neural

[25] Higdon D. Space and space-time modeling using process convolutions. In: Quantitative Methods for Current Environmental Issues. Springer; 2002. p. 37–56. [26] Alvarez MA, Luengo D, Lawrence ND. Linear latent force models using Gaussian processes. IEEE transactions on pattern analysis and machine intelligence. 2013;35(11):2693–2705. [27] Zheng X, Aragam B, Ravikumar P, Xing EP. Dags with no tears: Continuous optimization for structure learning. arXiv preprint arXiv:180301422. 2018;. [28] Pamfil R, Sriwattanaworachai N, Desai S, Pilgerstorfer P, Georgatzis K, Beaumont P, et al. Dynotears: Structure learning from time-series data. In: International Conference on Artificial Intelligence and Statistics. PMLR; 2020. p. 1595–1605. [29] Lerner I, Burgun A, Bach F. Spectral structure learning for clinical time series. arXiv preprint arXiv:250211680. 2025;. [30] Kalimouttou A, Lerner I, Cheurfa C, Jannot AS, Pirracchio R. Machine-learningderived sepsis bundle of care. Intensive Care Medicine. 2023;49(1):26–36. [31] Johnson A, Bulgarelli L, Pollard T, Horng S, Celi LA, Mark R. Mimic-iv. PhysioNet Available online at: https://physionet org/content/mimiciv/10/(accessed August 23, 2021). 2020;. [32] Silva I, Moody G, Scott DJ, Celi LA, Mark RG. Predicting in-hospital mortality of icu patients: The physionet/computing in cardiology challenge 2012. In: 2012 computing in cardiology. IEEE; 2012. p. 245–248. [33] Rasmussen CE. Gaussian processes in machine learning. In: Summer school on machine learning. Springer; 2003. p. 63–71. [34] Eichler M. Fitting graphical interaction models to multivariate time series. arXiv preprint arXiv:12066839. 2012;. [35] Brillinger DR. Remarks concerning graphical models for time series and point processes. Brazilian Review of Econometrics. 1996;16(1):1–23. [36] Dahlhaus R. Graphical interaction models for multivariate time series. Metrika. 2000;51:157–172. 21

[37] Bach FR, Jordan MI. Learning graphical models for stationary time series. IEEE transactions on signal processing. 2004;52(8):2189–2199. [38] Ormaniec W, Sussex S, Lorch L, Schölkopf B, Krause A. Standardizing structural causal models. arXiv preprint arXiv:240611601. 2024;. [39] Jacobs RA, Jordan MI, Nowlan SJ, Hinton GE. Adaptive mixtures of local experts. Neural computation. 1991;3(1):79–87. [40] Shazeer N, Mirhoseini A, Maziarz K, Davis A, Le Q, Hinton G, et al. Outrageously large neural networks: The sparsely-gated mixture-of-experts layer. arXiv preprint arXiv:170106538. 2017;. [41] Nemirovsky A. Optimization II. Numerical methods for nonlinear continuous optimization. Technion–Israel Institute of Technology; 1999. [42] Kingma DP, Ba J. Adam: A method for stochastic optimization. arXiv preprint arXiv:14126980. 2014;. [43] Solin A, Särkkä S. Hilbert space methods for reduced-rank Gaussian process regression. Statistics and Computing. 2020;30(2):419–446. [44] Riutort-Mayol G, Bürkner PC, Andersen MR, Solin A, Vehtari A. Practical Hilbert space approximate Bayesian Gaussian processes for probabilistic programming. Statistics and Computing. 2023;33(1):17. [45] Singer M, Deutschman CS, Seymour CW, Shankar-Hari M, Annane D, Bauer M, et al. The third international consensus definitions for sepsis and septic shock (Sepsis-3). Jama. 2016;315(8):801–810. [46] Bühlmann P. Invariance, causality and robustness. 2020;35(3):404–426.

Statistical Science.

[47] Wiggins BS, Sanoski CA. Emergency cardiovascular pharmacotherapy: a pointof-care guide. vol. 4. ASHP; 2012. [48] Charlier B, Feydy J, Glaunes JA, Collin FD, Durif G. Kernel operations on the GPU, with autodiff, without memory overflows. The Journal of Machine Learning Research. 2021;22(1):3457–3462. [49] Dallakyan A. On Learning Time Series Summary DAGs: A Frequency Domain Approach. arXiv preprint arXiv:230408482. 2023;. [50] Chan AB, Dong D. Generalized Gaussian process models. In: CVPR; 2011. p. 2681–2688. [51] Bonilla EV, Krauth K, Dezfouli A. Generic inference in latent Gaussian process models. arXiv preprint arXiv:160900577. 2016;. 22

[52] Murray S, Kjellström H. Mixed likelihood gaussian process latent variable model. arXiv preprint arXiv:181107627. 2018;. [53] Ramchandran S, Koskinen M, Lähdesmäki H. Latent Gaussian process with composite likelihoods and numerical quadrature. In: International Conference on Artificial Intelligence and Statistics. PMLR; 2021. p. 3718–3726. [54] Peters J, Janzing D, Schölkopf B. Elements of Causal Inference: Foundations and Learning Algorithms. The MIT Press; 2017. [55] Reisach A, Tami M, Seiler C, Chambaz A, Weichwald S. A scale-invariant sorting criterion to find a causal order in additive noise models. Advances in Neural Information Processing Systems. 2024;36.

S1 Supplementary methods S1.1 Detailed StructGP covariance The StructGP model encodes directed task dependencies through the matrix-valued filter H(t) = (I + S) ◦ L(t). For notational simplicity, we redefine S := I + S in what follows, so that

H(t) = S ◦ L(t). This convention absorbs the identity self-connections into S and keeps the notation lighter without affecting the subsequent derivations. Note that acyclicity constraints or sparsity penalties, when applied, concern only the off-diagonal part of S, corresponding to inter-task connections. Each task Yv (t) then results from the convolution of independent latent white noise processes {wu (t)}ku=1 :

Yv (t) =

k X

(Hvu ∗ wu )(t),

(1)

u=1

so that the covariance between tasks v and v ′ is given by

kvv′ (t, t′ ) =

k X

(Hvu ∗ Hv′ u )(t − t′ ).

(2)

u=1

Since the white noise processes are mutually independent, only shared sources u contribute to the covariance. We assume here that all filters Hvu (t) are Gaussian functions of the form t2 Hvu (t) = Svu exp − , (3) ℓvu where Svu encodes the connectivity strength from task u to task v , and ℓvu > 0 controls its temporal smoothness. Under this assumption, each convolution term has a closed-form expression: (Hvu ∗ H

v′ u

′

)(t − t ) = Svu S

v′ u

ℓvu ℓv′ u (t − t′ )2 π exp − . ℓvu + ℓv′ u ℓvu + ℓv′ u

(4)

The overall covariance between tasks v and v ′ is thus obtained by summing these terms over their common parents in the graph, that is,

kvv′ (t, t′ ) =

(Hvu ∗ Hv′ u )(t − t′ ),

u∈Pa(v)∩Pa(v ′ )

(5)

where Pa(v ) denotes the set of parent tasks of v in the DAG G . The learnable kernel parameters are grouped into θ = {S, L}, with S ∈ Rk×k encoding the graph structure and L ∈ Rk×k the temporal scales. + We also include a (potentially task-specific) regularization term σj δx=x′ to the kernel definition which models observation-level noise.

S1.2 Internal standardization The concept of internal standardization was originally introduced by Ormaniec et al. [1] in the context of linear additive structural causal models (SCMs). In their formulation, the structural coefficients are rescaled such that each variable has unit variance, thereby disentangling causal dependencies from the arbitrary measurement scales of the variables. We extend this principle to the Gaussian process setting: instead of linear instantaneous effects, dependencies between tasks are mediated through continuoustime convolutional filters. This normalization ensures that the learned amplitudes Svu represent only the relative strength of dependencies between tasks, independent of their individual measurement scales. For Gaussian filters of the form t2 Hvu (t) = Svu exp − , (6) ℓvu the marginal variance of process Yv (t) is obtained by integrating the squared filters over time: Var[Yv (t)] =

k Z ∞ X u=1

2 Hvu (t) dt =

−∞

k X

r 2 Svu

u=1

πℓvu . 2

(7)

We define the normalization factor for each task v as

sv =

k X

r 2 Svu

u=1

πℓvu 2

!1/2 .

(8)

The standardized amplitudes are then given by

S̃vu =

Svu , sv

(9)

and the corresponding normalized filters

t2 H̃vu (t) = S̃vu exp − ℓvu

(10)

satisfy Var[Yv (t)] = 1 for all v . This internal standardization stabilizes learning and isolates the structural information carried by the amplitude parameters from the raw scale of the observed variables. 2

S1.3 Detailed LP-StructGP formulation To capture shared temporal dynamics across subjects, we extend StructGP with an additional latent layer representing group-level trajectories, referred to as latent pathways. We introduce p latent multivariate white noise processes {wu (t)}pu=1 , with wu (t) ∈ Rk , each giving rise to a multivariate latent process Qu (t) ∈ Rk through the shared inter-task filter H(t):

Qu (t) = (H ∗ wu )(t).

(11)

All latent pathways therefore share the same inter-task dependency structure, consistent with a biological invariance hypothesis that assumes common task dependencies mechanisms across subjects. Each subject i expresses an individualized trajectory Yi (t) as a subject-specific mixture of these latent pathways, combined with a subject-specific residual component:

Yi (t) =

p X

(Giu ∗ Qu )(t) + (H ∗ wi )(t),

(12)

u=1

where wi (t) is a subject-specific white noise process. The resulting covariance between subjects i and i′ factorizes into inter-task and inter-subject components: ! p X ′ ′ ⊤ ′ k [(i, t), (i , t )] = (H ∗ H )(t − t ) ∗ δii′ + Giu ∗ Gi′ u (t − t′ ). (13) u=1

Each subject-level filter Giu (t) is parameterized as a shifted Gaussian function: (sub) Giu (t) = Siu exp

−

(t − τ iu )2

(sub)

(14)

ℓiu

(sub)

where Siu controls the coupling strength between subject i and latent pathway u, (sub) ℓiu > 0 determines the temporal smoothness, and τ iu represents a subject-specific temporal shift capturing progression along the shared latent trajectory. To encourage each subject to follow predominantly one latent pathway, we apply a softmax transformation to the coupling amplitudes: (sub)

π i = softmax(Si

(15)

yielding normalized weights πiu ∈ [0, 1] that sum to one. This mechanism acts as a soft gating function [2, 3], assigning each subject probabilistically to the most relevant latent pathway. The complete parameter set of LP-StructGP is thus

θ = {S, L, S(sub) , L(sub) , τ (sub) }, 3

(16)

where (S, L) define the shared inter-task structure, and (S(sub) , L(sub) , τ (sub) ) govern the subject–pathway interactions. In practice, the inter-task and inter-subject components are combined through a fixed gating coefficient, controlled by a hyperparameter that balances their relative contributions. For clarity, this mixing term is omitted from the above equations.

S1.4 Gaussian process learning and inference We consider a collection of n scalar observations y (x), drawn from r individuals and k tasks: {y (x) : x = (i, j, t), i ∈ {1, . . . , r}, j ∈ {1, . . . , k}, t ∈ R+ }. Each observation is indexed by the triplet x = (i, j, t), corresponding to the subject index, task index, and observation time, respectively. Without loss of generality, we assume a zero-mean prior, E[y (x)] = 0, since the posterior mean will be data-dependent. The process is then fully characterized by its covariance function:

y (x) ∼ GP (0, k (x, x′ )),

k (x, x′ ) = E[y (x)y (x′ )].

Gaussian process regression consists of two main steps. First, we place a prior distribution on y (x), with modeling assumptions encoded in the kernel k . The kernel hyperparameters θ (which include the convolution filter parameters and the noise variance σ 2 ) are learned by maximizing the log marginal likelihood of the training observations y given the inputs X [4]: 1 n 1 log p(y|X, θ) = − y⊤ (K + σ 2 I)−1 y − log |K + σ 2 I| − log(2π ), 2 2 2

(17)

where K is the n × n covariance matrix with entries Kuv = k (xu , xv ) for xu , xv ∈ X. This objective balances data fit (the quadratic term) and model complexity (the logdeterminant term), corresponding geometrically to finding the smallest covariance ellipsoid that adequately explains the data. In the second step, we define the joint prior distribution between n training observations y and m test observations y∗ , indexed by X ∈ (N, N, R+ )n and X∗ ∈ (N, N, R+ )m , respectively:

K + σ2 I K∗ y ∼ N 0, , y∗ K⊤ K∗∗ + σ 2 I ∗

(18)

where K, K∗ , and K∗∗ denote the covariance matrices over training, cross, and test inputs, respectively. Conditioning on y yields the posterior predictive distribution:

y∗ | y, X, X∗ ∼ N (µ∗ , Σ∗ ), where

2 −1 µ∗ = K⊤ y, ∗ (K + σ I)

(19) (20)

2 −1 Σ∗ = K∗∗ − K⊤ K∗ . ∗ (K + σ I)

(21)

In practice, the kernel kθ (x, x′ ) encodes the structured covariance induced by either StructGP or LP-StructGP, with shared parameters across tasks and patients.

S1.5 Differentiable structure learning: algorithmic details In StructGP, we impose that S encodes the weighted adjacency matrix of a sparse directed acyclic graph (DAG) G , i.e., that S is sparse and lower-triangular up to a permutation of the task order. To jointly learn the graph structure, the task ordering, and the sparsity pattern, we adapt the NOTEARS algorithm [5]. NOTEARS introduces a differentiable characterization of acyclicity based on the trace of the matrix exponential of the adjacency matrix (see S1.5.1). Our goal is therefore to solve the following constrained optimization problem:

h i θ∗ = argmin − log L(y, X, θ) + Pλ (S) θ s.t. h(S) := Tr exp(S ◦ S) − k = 0,

(22)

where Pλ is a sparsity-inducing penalty with strength λ, L is the process marginal likelihood and θ the set of learnable parameters. This problem is solved using the augmented Lagrangian method [6], which converts the constrained problem into a sequence of unconstrained subproblems (the primal ) and updates a dual variable by gradient ascent (see S1.5.2). The primal objective is the penalized likelihood augmented with a quadratic penalty term enforcing the acyclicity constraint. To make Pλ differentiable everywhere, we use the SmoothL1-SC approximation to the ℓ1 norm [7]:

Pλ (S) = λ

X 1 h i log 1 + eβL1 Si + log 1 + e−βL1 Si , βL1 i

where a large βL1 yields a close approximation to the standard ℓ1 penalty. We use Adam [8] as the gradient-based optimizer. The sparsity weight λ∗ is selected by grid search to minimize either an AIC-like criterion, AIC = 2∥S∥0 − 2 log L(y, X, θ), or a validation loss. The search is conducted on a logarithmic scale from λmax to λmin with warm starts. Because solving the augmented Lagrangian problem to high precision is computationally expensive and may cause numerical instability when the penalty coefficient ρ becomes large, we only solve it approximately, with a loose tolerance for the acyclicity constraint (typically ϵ = 0.01, see S1.5.2). To guarantee acyclicity in the final graph, we apply a hard-thresholding step to S, removing all entries below the smallest value that ensures h(S) = 0. 5

All implementation details and code are available at gitlab.

S1.5.1 The matrix exponential acyclicity constraint

b12

Y0 b2 0

 0 b01 0 B =  0 0 b12  b 0 0   20 0 0 b01 b12 0 0  B2 = b12 b20 0 b01 b20 0   b01 b12 b20 0 0 3  0 b01 b12 b20 0 B = 0 0 b01 b12 b20 

Fig. S1 Counting cycles in adjacency matrices using the matrix exponential.

The NOTEARS formulation expresses the acyclicity of a directed graph as a differentiable function of its weighted adjacency matrix [5]. Let B ∈ Rk×k denote the weighted adjacency matrix of a graph G with nodes {Y1 , . . . , Yk }. An edge Yu → Yv exists if and only if buv ̸= 0. The graph G (B) is acyclic if and only if

h(B) := Tr exp(B ◦ B) − k = 0, where ◦ denotes the Hadamard (elementwise) product. Expanding the matrix exponential, 1 Tr(exp S) = Tr(I ) + Tr(S) + 2! Tr(S2 ) + · · · , shows that Tr(Sm ) counts weighted cycles of length m (see Figure S1). Thus h(B) = 0 ensures that all such cycles vanish, yielding a differentiable characterization of acyclicity. For a detailed discussion, see [9].

S1.5.2 The augmented Lagrangian optimization algorithm Given:

• Objective function: f (θ) = − log L(y, X, θ) + Pλ (S) • Constraint: g (θ) = Tr(exp(S ◦ S)) − k = 0 • Convergence criteria: constraint tolerance ϵ and ρmax • Primal solver: Solver Define:

• Lagrange multiplier α(k) at step k • Augmented Lagrangian: L(θ, α(k) , ρ) = f (θ) + α(k) g (θ) + ρ2 g (θ)2 Do: 1. Choose initial guess θ(0) , Lagrange multipliers α(0) = 0 and ρ = 1 2. For k = 0, 1, 2, . . . (a) While ρ < ρmax update θ(k+1) . • Minimize the augmented Lagrangian function:

θ(k+1) = Solver L(θ, α(k) , ρ) • Break if g (.) sufficiently decreases: g (θ(k+1) ) < 0.25 g (θ(k) ) • Else augment ρ: ρ = 10 ρ (b) Update Lagrange multipliers:

α(k+1) = α(k) + ρ · g (θ(k+1) ) (c) Check convergence criteria. If satisfied, stop.

g (θ(k+1) ) < ϵ or ρ ≥ ρmax

S1.6 Hilbert-space Gaussian process approximation We build on the Hilbert-space Gaussian process (HSGP) approximation introduced by Solin and Särkkä [10], which provides a computationally efficient low-rank representation of stationary Gaussian process (GP) covariance functions. The key idea is to project the GP onto the eigenbasis of the negative Laplacian operator defined on a bounded domain Ω = [−L, L] with Dirichlet boundary conditions. Let {ϕj , λj }j≥1 denote the eigenfunctions and eigenvalues of the Laplacian, satisfying −∇2 ϕj (x) = λj ϕj (x), x ∈ [−L, L], ϕj (±L) = 0. (23)

In one dimension, these eigenpairs admit closed-form expressions:

λj =

πj 2L

r ,

ϕj (x) =

1 πj (x + L) sin , L 2L

j = 1, 2, . . .

(24)

forming an orthonormal basis of L2 ([−L, L]). For a stationary kernel k (x, x′ ) with spectral density S (ω ) (via Bochner’s theorem), the covariance function can be approximated as

k (x, x′ ) ≈

m p X S λj ϕj (x) ϕj (x′ ).

(25)

j=1

Defining feature functions: Φj (x) =

r S

p λj ϕj (x),

(26)

the kernel matrix admits a low-rank factorization

K ≈ ΦΦ⊤ ,

(27)

where Φ ∈ Rn×m is the feature matrix evaluated at the training inputs.

S1.6.1 Squared-exponential kernel in one dimension For the squared-exponential (SE) kernel

(t − t′ )2 kSE (t, t′ ) = α exp − , 2ℓ2

(28)

the corresponding spectral density is Gaussian:

√ S (ω ) = α 2π ℓ exp − 21 ℓ2 ω 2 . Substituting ωj = features:

(29)

λj = πj/(2L) into (25) yields the one-dimensional HSGP q √ q1 πj(t+L) Φj (t) = α 2π ℓ exp − 14 ℓ2 λj sin . (30) L 2L

Truncating the basis to m terms provides a controlled low-rank approximation whose accuracy increases with both m and the ratio L/ℓ. This approximation is applied to each component of the structured covariance function (see S1.1)

S1.7 Online Woodbury Updates with Geometrically Discounted Accumulators We consider the problem of computing x = K −1 v , where

K = ΦΦ⊤ + M, and M is block-diagonal. New blocks of data arrive sequentially. For each new minibatch of k blocks, {(vi , Φi , Mi )}ki=1 , we update both the inverse action and the log-determinant of K using low-rank Woodbury corrections and geometrically discounted accumulators. For each block i: Ai = Mi−1 vi , Bi = Mi−1 Φi , and the corresponding local contributions:

di = Φ⊤ i Ai ,

ci = Φ⊤ i Bi .

The global accumulators (D, C ) are then updated using a fading-memory rule with decay parameter β ∈ [0, 1):

D ← βD +

C ← βC +

di ,

ci .

Finally, the correction term and log-determinant are updated as:

E = (I + C )−1 D,

log det K = log det M + log det(I + C ).

The inverse action for the latest mini-batch is then obtained as

xi = Ai − Bi E,

for i = 1, . . . , k.

This provides the blockwise components of K −1 v corresponding to the most recent data. This update implements a geometric fading memory : past contributions are scaled by β t , so β controls how quickly old information decays. Setting β = 0 uses only the current batch (no memory), whereas β → 1 approximates full memory. This allows us to emulate continuous learning across epochs during the optimization of the kernel parameters.

S1.8 Online conditional pseudo-marginal likelihood For LP-StructGP, the joint covariance of all subjects or tasks is

Kθ = Φθ Φ⊤ θ + Mθ ,

where Mθ is block-diagonal and encodes subject-specific noise or short-range dependencies, while Φθ collects the low-rank HSGP basis features that couple all subjects through shared latent pathways. The negative log marginal likelihood (NMLL) of all observations y ∈ RN under this model is 1 N 1 log(2π ), L(θ) = y ⊤ Kθ−1 y + log det Kθ + 2 2 2 where θ denotes all kernel and noise parameters. When the data are processed sequentially, this global objective can be decomposed exactly into a sum of conditional contributions: log p(y1:T | θ) =

T X

log p(yt | y1:t−1 , θ),

t=1

so that the full NMLL can be written as

L(θ) =

T X

Lt (θ),

Lt (θ) =

t=1

Nt 1 ⊤ −1 1 y K yt + log det Kt − log det Kt−1 + log(2π ), 2 t t 2 2

where Kt denotes the covariance restricted to all data observed up to batch t. Each term Lt (θ) is therefore the conditional negative log marginal likelihood of the new block yt given all previous blocks y1:t−1 , and their sum exactly recovers the full-data NMLL. This conditional factorization forms the basis for our online updates. By maintaining sufficient statistics of the accumulated data (Ct , Dt ), we can evaluate each incremental term Lt (θ) without revisiting earlier batches.

Batch structure The dataset is organized into T mini-batches, each containing kt subject-level blocks: t Bt = {(yi , Φi , Mi )}ki=1 ,

t = 1, . . . , T.

Each block corresponds to a subject (or subject–task group) with local covariance Mi and basis design matrix Φi . Although these blocks are not independent, their crosscovariances are represented implicitly through the shared low-rank term Φθ Φ⊤ θ . During sequential processing, we maintain global low-rank accumulators

−1 Φ⊤ i M i Φi ,

−1 Φ⊤ i Mi yi ,

which summarize the contributions of all observed subjects and encode the current estimate of the shared latent structure. These accumulators are updated after each new batch and enter the Woodbury updates described below and detailed in S1.7, where in practice we apply a geometric fading memory to smooth the loss across epochs.

Batchwise conditional NMLL The quadratic fit term of the full NMLL can be written as T 1 X X ⊤ −1 1 ⊤ −1 y Kθ y = yi (Kθ y )i . 2 2 t=1 i∈Bt

For each batch t, the online Woodbury solver provides the inverse actions (t)

xi = Kt−1 yi = Ai − Bi Et ,

Ai = Mi−1 yi ,

Bi = Mi−1 Φi ,

Et = (I + Ct )−1 Dt .

The batch contribution to the fit term is then Fitt =

1 X ⊤ (t) yi xi . 2 i∈Bt

The log-determinant of the covariance decomposes as −1 Cθ = Φ⊤ θ M θ Φθ =

log det Kθ = log det Mθ +log det(I +Cθ ),

−1 ci , ci = Φ⊤ i Mi Φi .

When data are processed sequentially, we track the incremental log-determinant contribution of each batch as X ∆ log det Kt = log det Mi + log det(I + Ct ) − log det(I + Ct−1 ) . i∈Bt

Combining both components, the conditional negative log marginal likelihood for batch t reads

Lt (θ | Ct , Dt ) =

Nt 1 X ⊤ (t) 1 yi xi + ∆ log det Kt + log(2π ), 2 2 2 i∈Bt

where Nt =

i∈Bt dim(yi ) is the total number of observations in the batch.

Interpretation for LP-StructGP In LP-StructGP, subjects are coupled through shared latent pathways encoded by the low-rank term Φθ Φ⊤ θ . As a result, the full marginal likelihood cannot be written as a sum over subjects: X L(θ) ̸= Li (θ), i

and its gradient does not decompose into independent per-subject contributions. During mini-batch training, we therefore optimize a conditional pseudo-marginal likelihood Lt (θ | Ct , Dt ) = − log p(yt | Ct , Dt , θ),

which conditions on the current empirical estimate of the latent process summarized by the accumulators (Ct , Dt ). In the large-sample limit, this conditional objective can be viewed as an exact approximation of the full-data NMLL: Z p(yt | y1:t−1 , θ) = p(yt | f, θ) p(f | y1:t−1 , θ) df ≈ p(yt | fˆt , θ), where the posterior distribution over the latent process p(f | y1:t−1 , θ) is replaced by its mode (or mean) fˆt , whose empirical estimate is encoded by (Ct , Dt ). As more data are processed, the accumulators provide increasingly accurate estimates of fˆt , making this approximation effectively exact and yielding a scalable surrogate objective that asymptotically matches the true marginal likelihood.

S1.9 Simulation details We summarize here the data generation procedures and parameter settings used in the simulation study.

Graph and kernel parameters For each simulation, the support of S was sampled from a random Erdős–Rényi graph with mean degree 2. Nonzero edge weights of S were uniformly drawn from [−1.5, −0.5] ∪ [0.5, 1.5]. Log-lengthscale parameters log ℓu were uniformly sampled in [0, 1]. Observation times t were drawn uniformly in [0, 10]. The noise variance was fixed to σ 2 = 10−2 and provided as an oracle during model fitting. Latent pathway parameters In LP-StructGP simulations, each subject trajectory was generated as a mixture of one shared latent pathway and an individual component, Yi (t) = 0.3 (Gi ∗ Qui )(t) + 0.7 (H ∗ wi )(t), where ui denotes the latent pathway assigned to subject i, and the coefficients 0.3 and 0.7 control the relative contributions of the latent and subject-specific components, respectively.

Design Each simulation involved k = 10 tasks and 25 observations per task. We varied the number of subjects while keeping other parameters fixed. For each configuration, 30 independent repetitions were performed, and the median and interquartile range of all evaluation metrics (SHD, ARI, and NMI) were reported. Implementation All simulations were implemented in PyTorch using the same optimization settings as described in Section 2.4: Simulation study of the manuscript. 12

S1.10 Real-life evaluation details Cohort selection Patients were included if they met Sepsis-3 criteria: suspected infection and a Sequential Organ Failure Assessment (SOFA) score ≥ 2. Septic shock was defined as sepsis with persistent hypotension requiring vasopressors to maintain a mean arterial pressure (MAP) ≥ 65 mmHg and a serum lactate > 2 mmol·L−1 despite adequate fluid resuscitation. For patients with multiple ICU stays, only the first stay was retained. Data extraction and preprocessing Time t = 0 was set to ICU admission. We extracted the first seven days of the ICU stay and excluded patients with fewer than the 20th percentile of available measurements for creatinine, norepinephrine, or MAP. All time series were resampled into overlapping 24-hour windows for the short-term experiments. Marginal distributions were transformed to standard normal using a quantile transform fitted on the training set. Experimental protocol The dataset was split by subject into 70% training, 15% validation, and 15% test. Each configuration was repeated five times with different initializations for unstructured baseline models. Unless stated otherwise, all hyperparameters and optimizer settings were identical across architectures. Evaluation metrics We report the Root Mean Square Error (RMSE) and empirical 95% predictive coverage on test data. Additional exploratory experiments with time-shifted features and extended covariates were also performed. Extended task set For the large-scale evaluation, we considered 18 tasks comprising vital signs, laboratory biomarkers, and vasoactive drug administration rates: • Hemodynamic variables: mean arterial pressure (MAP), systolic blood pressure (SBP) and heart rate (HR). • Respiratory variables: respiratory rate (RR), oxygen saturation (SpO2 ), partial pressure of oxygen (PO2 ), partial pressure of carbon dioxide (pCO2 ), blood pH, and venous oxygen saturation (sO2 ). • Renal and metabolic biomarkers: serum creatinine (Cr), serum lactate (Lac). • Vasoactive and therapeutic interventions: norepinephrine (NE), epinephrine (Epi), vasopressin (Vaso), phenylephrine (PE), dobutamine (Dobut), and sodium chloride 0.9% (NaCl). • others: body temperature (Temp) All variables were preprocessed and standardized following the same procedure as described above.

S2 Supplementary results S2.1 Missingness

Fig. S2 A patient’s stay observation pattern in Electronic Health Records. Scatter plot of a patient’s stay, with one point per observation. Each line/color represents a different variable (y-axis) as a function of time in days (x-axis).

Fig. S3 Forecasting after uncovering latent pathways.

Toy model illustrating the interpolation between group and individual forecasts after conditioning on observations from t = 0 to t = 5, once latent pathways have been uncovered.

S2.2 Simulation S2.3 PhysioNet Challenge

GraFITi GCS FiO2 Cholesterol TroponinT Lactate TroponinI PaO2 SysABP HCT Bilirubin NISysABP BUN NIMAP Albumin HR HCO3 DiasABP Platelets

MAE ×10−2 [95% CI] StructGP (best)

9.48 [9.30, 9.65] 10.55 [10.25, 10.86] 8.62 [5.75, 13.00] 2.19 [1.59, 2.95] 3.49 [3.24, 3.74] 11.92 [9.10, 15.68] 7.16 [6.94, 7.39] 5.22 [5.17, 5.26] 4.26 [4.16, 4.37] 1.79 [1.54, 2.07] 4.37 [4.33, 4.42] 2.44 [2.35, 2.53] 4.22 [4.17, 4.26] 8.71 [8.03, 9.40] 3.18 [3.15, 3.20] 4.29 [4.18, 4.42] 2.82 [2.79, 2.84] 1.40 [1.35, 1.45]

13.71 [13.44, 14.01] 12.61 [12.27, 12.96] 9.82 [7.18, 13.43] 3.00 [2.20, 3.82] 4.07 [3.72, 4.45] 12.39 [8.27, 16.25] 7.46 [7.23, 7.70] 5.51 [5.46, 5.56] 4.52 [4.41, 4.63] 2.02 [1.75, 2.33] 4.56 [4.51, 4.60] 2.58 [2.47, 2.67] 4.33 [4.29, 4.37] 8.82 [8.14, 9.49] 3.27 [3.25, 3.30] 4.38 [4.24, 4.49] 2.87 [2.85, 2.90] 1.44 [1.39, 1.50]

GraFITi

MSE ×10−3 [95% CI] StructGP (best)

24.67 [23.75, 25.65] 28.88 [27.30, 30.43] 16.41 [4.70, 40.74] 4.07 [1.10, 8.45] 3.79 [3.06, 4.53] 31.36 [16.17, 55.13] 11.88 [10.83, 12.96] 4.72 [4.61, 4.82] 3.14 [2.96, 3.31] 1.38 [0.94, 1.88] 3.39 [3.30, 3.50] 1.35 [1.22, 1.51] 3.03 [2.95, 3.11] 13.46 [11.30, 15.67] 1.91 [1.88, 1.95] 3.29 [3.04, 3.56] 1.47 [1.42, 1.53] 0.43 [0.37, 0.52]

59.73 [57.65, 61.86] 42.17 [40.19, 44.21] 16.62 [6.82, 34.11] 6.72 [2.92, 11.54] 6.24 [4.75, 7.89] 42.85 [18.77, 69.23] 13.48 [12.33, 14.76] 5.40 [5.28, 5.51] 3.55 [3.35, 3.75] 1.84 [1.33, 2.45] 3.74 [3.63, 3.84] 1.54 [1.37, 1.73] 3.24 [3.16, 3.31] 13.78 [11.75, 15.94] 2.07 [2.03, 2.11] 3.50 [3.22, 3.81] 1.63 [1.57, 1.69] 0.43 [0.38, 0.51]

Table S1 Underperforming tasks compared with GraFITi (SOTA) on the PhysioNet challenge Mean absolute error (MAE) and Mean squared error (MSE) with bootstrap 95% confidence intervals. Metrics are reported for up to a 1-day forecasting horizon for 18 tasks, after the first 24 hours of the ICU stay.

PaCO2 MAP Na NIDiasABP RespRate Creatinine K ALP Glucose ALT Temp SaO2 AST Mg Urine WBC Weight pH MechVent

MAE ×10−2 [95% CI] GraFITi StructGP (best)

GraFITi

4.77 [4.65, 4.89] 3.43 [3.39, 3.46] 3.00 [2.91, 3.09] 4.31 [4.26, 4.35] 3.27 [3.22, 3.32] 1.23 [1.16, 1.29] 1.85 [1.79, 1.92] 0.72 [0.58, 0.93] 1.96 [1.88, 2.03] 1.43 [1.17, 1.71] 1.03 [1.01, 1.06] 1.84 [1.74, 1.96] 1.31 [1.09, 1.57] 0.83 [0.80, 0.85] 0.84 [0.83, 0.85] 0.36 [0.36, 0.37] 0.67 [0.67, 0.68] 0.37 [0.36, 0.39] 0.47 [0.46, 0.47]

4.29 [4.02, 4.58] 2.21 [2.13, 2.31] 1.64 [1.52, 1.77] 3.26 [3.18, 3.34] 1.96 [1.87, 2.07] 0.52 [0.39, 0.73] 0.70 [0.61, 0.81] 0.50 [0.08, 1.28] 0.87 [0.73, 1.00] 1.59 [0.90, 2.45] 0.41 [0.27, 0.57] 0.87 [0.61, 1.24] 1.34 [0.71, 2.26] 0.12 [0.10, 0.15] 0.13 [0.12, 0.14] 0.02 [0.02, 0.02] 0.14 [0.12, 0.15] 0.03 [0.03, 0.03] 0.03 [0.03, 0.04]

4.78 [4.65, 4.91] 3.42 [3.39, 3.46] 2.98 [2.89, 3.07] 4.28 [4.24, 4.32] 3.23 [3.19, 3.28] 1.19 [1.12, 1.25] 1.78 [1.72, 1.83] 0.60 [0.46, 0.80] 1.80 [1.73, 1.87] 1.27 [0.98, 1.65] 0.86 [0.84, 0.88] 1.61 [1.50, 1.73] 1.09 [0.83, 1.37] 0.60 [0.58, 0.63] 0.59 [0.58, 0.60] 0.04 [0.04, 0.04] 0.32 [0.31, 0.33] 0.01 [0.01, 0.01] 0.00 [0.00, 0.00]

MSE ×10−3 [95% CI] StructGP (best) 4.50 [4.16, 4.87] 2.38 [2.28, 2.48] 1.64 [1.52, 1.79] 3.31 [3.22, 3.39] 1.98 [1.88, 2.10] 0.52 [0.44, 0.61] 0.60 [0.55, 0.65] 0.49 [0.06, 1.23] 0.81 [0.66, 0.96] 1.85 [0.95, 2.96] 0.35 [0.21, 0.52] 0.99 [0.69, 1.37] 1.52 [0.77, 2.41] 0.08 [0.06, 0.11] 0.13 [0.12, 0.14] 0.00 [0.00, 0.00] 0.11 [0.09, 0.12] 0.00 [0.00, 0.00] 0.00 [0.00, 0.00]

Table S2 Outperforming tasks compared with GraFITi (SOTA) on the PhysioNet challenge Mean absolute error (MAE) and Mean squared error (MSE) with bootstrap 95% confidence intervals. Metrics are reported for up to a 1-day forecasting horizon for 19 tasks, after the first 24 hours of the ICU stay.

High Model Certainty GraFITi StructGP (best) Moderate Model Certainty GraFITi StructGP (best) Low Model Certainty GraFITi StructGP (best)

MAE ×10−2 [95% CI]

MSE ×10−3 [95% CI]

2.53 [2.44, 2.64] 2.46 [2.34, 2.60]

2.03 [1.81, 2.30] 2.20 [1.90, 2.58]

2.81 [2.62, 3.03] 2.96 [2.73, 3.20]

2.97 [2.06, 4.01] 3.79 [2.81, 5.00]

2.80 [2.71, 2.90] 3.05 [2.85, 3.27]

2.27 [2.08, 2.47] 3.59 [2.78, 4.62]

Table S3 Performance Comparison across model confidence levels on the PhysioNet challenge Mean absolute error (MAE) and Mean squared error (MSE) with bootstrap 95% confidence intervals. Metrics are reported for up to a 1-day forecasting horizon for 41 tasks, after the first 24 hours of the ICU stay. Metrics are segmented into terciles based on the StructGP (best) model posterior predictive variance.

References [1] Ormaniec W, Sussex S, Lorch L, Schölkopf B, Krause A. Standardizing structural causal models. arXiv preprint arXiv:240611601. 2024;. [2] Jacobs RA, Jordan MI, Nowlan SJ, Hinton GE. Adaptive mixtures of local experts. Neural computation. 1991;3(1):79–87. [3] Shazeer N, Mirhoseini A, Maziarz K, Davis A, Le Q, Hinton G, et al. Outrageously large neural networks: The sparsely-gated mixture-of-experts layer. arXiv preprint arXiv:170106538. 2017;. [4] Mardia KV, Marshall RJ. Maximum likelihood estimation of models for residual covariance in spatial regression. Biometrika. 1984;71(1):135–146. [5] Zheng X, Aragam B, Ravikumar P, Xing EP. Dags with no tears: Continuous optimization for structure learning. arXiv preprint arXiv:180301422. 2018;. [6] Nemirovsky A. Optimization II. Numerical methods for nonlinear continuous optimization. Technion–Israel Institute of Technology; 1999. [7] Schmidt M, Fung G, Rosales R. Optimization methods for l1-regularization. University of British Columbia, Technical Report TR-2009-19. 2009;. [8] Kingma DP, Ba J. Adam: A method for stochastic optimization. arXiv preprint arXiv:14126980. 2014;. [9] Vowels MJ, Camgoz NC, Bowden R. D’ya like dags? a survey on structure learning and causal discovery. ACM Computing Surveys. 2022;55(4):1–36. [10] Solin A, Särkkä S. Hilbert space methods for reduced-rank Gaussian process regression. Statistics and Computing. 2020;30(2):419–446.

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