Spectral structure learning for clinical time series

Venue: NeurIPS_cc_2024_Workshop_BDU
Published: 2024-08-24 13:47:34

Abstract: We develop and evaluate a structure learning algorithm for clinical time series.
Clinical time series are multivariate time series observed in multiple patients and irregularly sampled, challenging existing structure learning algorithms.
We assume that our times series are realizations of StructGP, a $k$-dimensional multi-output or multi-task stationary Gaussian process (GP), with independent patients sharing the same covariance function.
StructGP encodes ordered conditional relations between time series, represented in a directed acyclic graph.
We implement an adapted NOTEARS algorithm, which based on a differentiable definition of acyclicity, recovers the graph by solving a series of continuous optimization problems.
Simulation results show that up to mean degree 3 and 20 tasks, we reach a median recall of 0.93% [IQR, 0.86, 0.97] while keeping a median precision of 0.71% [0.57-0.84], for recovering directed edges.
We further show that the regularization path is key to identifying the graph.
With StructGP we proposed a model of time series dependencies, that flexibly adapt to different time series regularity, while enabling us to learn these dependencies from observations.

Main paper: 
 introduction: 
 Structure learning is the task of learning the dependency structure of either time-independent variables or, in this study, time series [1, 2] . The structure of dependency is usually represented as a directed acyclic graph (DAG), in which nodes represent variables, and links between nodes represent relations between variables. These relations encode conditional or marginal independence relations [3, 4] , or under certain strong additional assumptions (e.g., no hidden cofounders), have a causal interpretation and the structure can be used for causal modeling [1, Chapter 6] .Structure learning algorithms are typically classified into constraint-based methods and score-based methods [2] . Score-based methods assume a parametric model in which the parameters support the graph, and then search for the graph that best fits the data [5] . The main limitation of this approach is that the acyclicity of the graph requires an iterative search over graph structures that satisfy the constraint. And that the number of acyclic graphs grows superexponentially with the number of nodes in the graph [6] . Formulating the acyclicity of a directed acyclic graph (DAG) as a differentiable function of its adjacency matrix allowed Zheng et al. [7] to frame the problem as a continuous optimization problem. Following this line of work, DYNOTEARS [8] uses the acyclicity constraint to identify the structure of linear structural vector autoregressive models (SVAR) [9] .In this study, we aim to develop a structure learning algorithm for clinical time series collected from Electronic Health Records (EHRs). EHRs time series are collections of irregularly sampled multivariate time series, collected in multiple patients. To learn graphical models of dependence between those, we work in continuous time and develop StructGP, a structured multi-task Gaussian process (GP) [10] .Submitted to Workshop on Bayesian Decision-making and Uncertainty, 38th Conference on Neural Information Processing Systems (BDU at NeurIPS 2024). Do not distribute. 
structured gaussian process: 
 We consider the k-dimensional multi-output Gaussian process Yptq defined as the filtration by Hptq of white noise wptq:Yptq " pH ˚wqptq,where Hptq is a sparse k ˆk lower triangular matrix-valued impulse response function, and wptq is a k-dimensional white noise vector. Taking the Fourier transform of Yptq, we find Zpωq a non-stationary complex white noise multi-output process [11, p. 418] :EQUATION)where Hpωq is the Fourier transform of Hptq, and Wpωq complex independent white noise processes.In addition, the covariance of Zpωq or spectral density is:Cov `Zpωq, Zpω 1 q ˘" " Hpωq HT pωq if ω " ω 1 0 if ω ‰ ω 1 .Thus, it identify Hpωq as the Cholesky factor of the covariance matrix of Zpωq. Given Equation 5 in Appendix A6.3, we find that the sparsity pattern of H parameterizes ordered conditional relations between time series:H vu ptq " 0, @t P R ô Hvu pωq " 0, @ω P p´π, πq ô Z u K K Z v | Z t1,2,...,u´1u ô f uv|C pωq " f uv pωq ´fuC pωqf ´1 CC pωqf Cv pωq " 0, @ω P p´π, πq ô Y u K K Y v | Y t1,2,...,u´1u, where f uv|C is the partial cross spectrum of Y u and Y v given C " Y t1,2,...,u´1u .We therefore parameterize Hptq as follows:Hptq " pI ´Sq ˝Lptq,where L vu ptq " exp `´t 2 Lvu ˘, S is a sparse lower triangular matrix up to permutation, L is a positive square matrix, I is the identity matrix, and ˝the Hadamard product. The support of S can then be interpreted as the adjacency matrix of a directed acyclic graph (DAG) G that encodes ordered conditional independence relation:S vu ‰ 0 (3) ôY u G Ý Ñ Y v ôY u K KY v | Y t1,2,...,u´1u .As our time series are irregularly sampled from multiple patients we switch here to a set-based indexing, thanks to the marginalization properties of GP. Our data is then a collection of n scalar observations y from r individuals and k tasks, typxq : x P pN, N, R `qu. Each observation is indexed by the input vector x, a triplet composed of the patient index i, task index j, and time t, such that x " pi, j, tq, with i P t1, ..., ru, j P t1, ..., ku and t P R `. We observe multiple independent patients, and the intra-patient covariance of the process for patient i can be written:CovrYpi, u, tq, Ypi, v, t 1 qs " ph u ˚hT v qpt ´t1 q. With classical GP, the set of free parameters θ " tS, Lu is learned by maximizing log ppy|X, θq, the marginal likelihood of the training observations y given inputs X [12] . With structured GP, we follow the NOTEARS algorithm to learn the graph, order of tasks, and sparsity pattern [7] and impose an acyclicity constraint on S. NOTEARS leverages the trace of the matrix exponential of the adjacency matrix as a differentiable acyclicity constraint. Our objective is therefore to solve the constrained optimization problem below:θ ˚" argmin θ ´log ppy, X, θq `λ}S} 1 (4) s.t. TrpexppS ˝Sqq ´k " 0,where λ is the penalty strength.The above is solved by dual ascent following the augmented Lagrangian method [13] , such that the constrained problem is equivalent to solving a series of unconstrained problems, the primal and the dual. The primal is the penalized objective function augmented with a quadratic penalty term, and is solved with a proximal gradient method (see Appendix A6.5), the dual is solved by gradient ascent (see Appendix A6.4). Finally, we find by grid search λ ˚, the sparsity penalty that minimizes an equivalent of the Akaike information criterion (AIC): AIC " 2}S} 0 ´2 log Lpy, X, θq.Grid search is conducted from λ max to λ min on a log-scale, with warm-start. Because solving the augmented Lagrangian problem to small error is computationally expensive, and leads to numerical instability when ρ becomes large, we choose to only loosely solve it, typically with a large tolerance for the acyclicity constraint (ϵ " 0.1, see Appendix A6.4). Thus, to ensure dagness, we apply a hard-threshold operation, i.e. we mask elements in S lower than the minimal threshold that ensures dagness.We empirically assess the accuracy of the algorithm to identify a graph from observations through a series of simulations. For each simulation, we sample a graph, the covariance parameters θ, and observations from the sampled prior GP. The definition of the GP follows that of section 2.1, with the additional constraint that lengthscales parameters are tied for each task and exponentiated ( L vu " exppℓ v q for all v P t1, 2, . . . , ku, u P t1, 2, . . . , ku). We then fit the model using the overall algorithm from section 2.2 (including the grid search), and compare the predicted graph Ĝ with the true simulated graph G. The comparison is made with the structural hamming distance (SHD).SHD is the "edit distance" of graphs, it counts the number of edge modifications (insertion, deletion, inversion) necessary to transform a predicted graph into the true simulated graph. We also compare the root mean square error (RMSE) between the predicted Ŝ and true (simulated) S parameters. Each simulation is repeated 100 times and we report the average metric along with its bootstrapped 95% confidence intervals. For comparison purposes, we also report the same metrics for a random graph from the same distribution as simulated.In all simulations, the support of S is sampled from a random (Erdős-Rényi) graph in which the sparsity level is controlled by the mean degree of the graph md. S parameters are uniformly sampled in r´2, ´0.5s Y r0.5, 2s. ℓ u parameters are uniformly sampled in r´0. 3 Results 
simulation study: 
 We report in Figure 1 the results of 'EXP1' which shows that the algorithm from section 2.2 identifies the true graph almost perfectly in this 'ideal' regime of parameters (large weights, sparse graphs, low noise). Indeed, with 100 patients, most of the errors that contribute to the SHD are extra links with small weights. We then show in supplementary Figure S2 , that the number of grid search steps is heavily impacting the capacity to recover the true graph (red dots), but not because it allows us to precisely infer the optimal λ ˚(green dots). Indeed, fitting the model from a random initialization with the optimal λ ˚produces very poor performances, and improving the number of grid search steps does not improve the solutions (green dots). In supplementary Figure S3 , we see that up to mean degree 3 and 20 tasks, most errors are spurious links. Precision (P) and recall (R) median and interquartile for varying mean degree and number of task parameters, from 100 replications of 'EXP3' simulations with md " 3. False positives include extra links and reversed links. Metrics are computed from models learned (Predicted) and compared with metrics for random graphs of the same distribution as the simulated data (Random). 
conclusion: 
 We developed a structure learning algorithm that learns ordered conditional independence relations between irregularly sampled time series. These relations are parametrized by StructGP, a GP model built upon the convolution between a sparse lower-triangular matrix-valued impulse response function and independent white noises. It corresponds to assuming a linear additive Gaussian SCM for the Fourier representation of the time series, whose structure is invariant over all frequencies. Based on a differentiable definition of acyclicity, this algorithm recovers the true graph by solving a series of continuous optimization problems with high sensitivity and good precision on simulated data. Further work will be needed to bridge the gap between simulated data and real-world data, regarding the sensitivity to scaling and unmeasured cofounders. These sections adapt proofs from [4, 14] with the same notations as the remainder of the article for simplicity. 
sparsity pattern of the precision matrix encodes conditional independence: 
 Let Y " pY 1 , Y 2 , . . . , Y k q J be a multivariate normal random vector with mean vector µ and covariance matrix K, i.e., Y " N pµ, Kq.Consider the two sets of variables A " tu, vu and B " t1, 2, . . . , kuztu, vu. The covariance matrix can be partitioned as:K " ˆKA K AB K T AB K B ˙.After conditioning on B, the conditional covariance matrix is given by:K A|B " K A ´KAB K ´1 B K T AB .Now, define the precision matrix Ω " K ´1 with the same partition:Ω " ˆΩA Ω AB Ω T AB Ω B ˙.Using the inversion of a partitioned matrix, we recognize that:Ω A " pK A ´KAB K ´1 B K T AB q ´1.Or otherwise stated that the precision matrix is invariant to conditioning:Ω A|B " K ´1 A|B " Ω A .Thus, let:Ω A " ˆωuu ω uv ω uv ω vv ˙.Then, we have:K A|B " 1 detpΩ A q ˆωvv ´ωuv ´ωuv ω uu ˙.This shows that:Y u K K Y v | Y ztu,vu ô ω uv " 0.Let Y " pY 1 , Y 2 , . . . , Y k q J be a multivariate normal random vector with mean vector µ and covariance matrix K, i.e., Y " N pµ, Kq.Consider the two sets of variables A " t1, 2, . . . , u ´1u and B " tu, . . . , ku. The covariance matrix K can be partitioned and factorized as follows:K " ˆKA K AB K T AB K B " ˆI 0 K J AB K ´1 A I ˙ˆK A 0 0 K B ´KJ AB K ´1 A K AB ˙ˆI K ´1 A K AB 0 I ˙.Hence, the Cholesky factor L of K is given by:L " ˆLA 0 K J AB L ´1 A L S ˙, where L A " cholpK A q, and L S " cholpK B ´KJ AB K ´1 A K AB q.The element of the v-th row and u-th column, L vu , can be found in L S in the first column of at the row pv ´u `1q: L vu " L S,pv´u`1q1 .The conditional covariance given A (i.e., all the variables preceding u) is:K B|A " L S L T S .This shows that:EQUATIONô K B|A,v1 " 0ô k ÿ u"1L S,pv´u`1qu L S,1u " 0 ô L S,pv´u`1q1 " 0 ô L vu " 0.Given:• Objective function: f pθq " ´log Lpy, X, θq `Pλ pSq • Constraint: gpθq " trpexppS ˝Sqq ´k " 0 • Convergence criteria: constraint tolerance ϵ and ρ max • Primal solver: Solver Define:• Lagrange multiplier α pkq at step k• Augmented Lagrangian: Lpθ, α pkq , ρq " f pθq `αpkq gpθq `ρ 2 gpθq 2 Do:1. Choose initial guess θ p0q , Lagrange multipliers α p0q " 0 and ρ " 1 2. For k " 0, 1, 2, . . . (a) While ρ ă ρ max update θ pk`1q .• Minimize the augmented Lagrangian function:θ pk`1q " Solver `Lpθ, α pkq , ρq • Break if gp.q sufficiently decreases:gpθ pk`1q q ă 0.25 gpθ pkq q• Else augment ρ: ρ " 10 ρ (b) Update Lagrange multipliers:α pk`1q " α pkq `ρ ¨gpθ pk`1q q (c) Check convergence criteria. If satisfied, stop.gpθ pk`1q q ă ϵ or ρ ě ρ max 6.5 Proximal gradient method PGM [15] is a generalization of gradient descent for objective functions that can be split between a differentiable and non-differentiable part. It requires defining a proximal operator for the nondifferentiable part:prox λ pSq uv " $ & % s uv ´λ, if s uv ą λ 0, if s uv ď |λ| s uv `λ, if s uv ă ´λ.Parameters are found by taking a classic gradient step followed by a proximal step at each iteration k:(i) θ k`1 2 " θ k ´αk ∇p´log Lpy, X, θ k qq (ii) θ k`1 " prox α k λ pθ k`1 2 q with α k a learning rate determined by line search at each iteration, following Beck and Teboulle [16] . Reports of 'EXP2', an increasing number of grid search steps for 50 patients, 10 tasks, 10 observations per task, and random graph of mean degree 3. The predicted graph (red dots) is compared with a random graph from the same graph distribution (blue dots), and also with a predicted graph directly fitted from a random initialization with the optimal λ˚found from grid search (green dots). Error bars represent bootstrap 95% confidence intervals. 
 
 Appendix:  



Reviews: review: title: title: Spectral structure learning for clinical time series
main_review: main_review: This is a well-written paper that proposes a novel method for learning the dependency structure of clinical time series data collected from Electronic Health Records (EHRs).

Quality:

The paper is well-structured and easy to follow.
The authors clearly define the problem and their approach.
The methodology is described in detail, including relevant formulas and references.
The evaluation methodology is sound, with simulations to assess the accuracy of the algorithm.

Clarity:

The writing is clear and concise.
The authors use appropriate technical language.
Figures and tables are informative and well-explained.

Originality:

The proposed method, StructGP, is a novel approach to learning the dependency structure of clinical time series data.
The use of a differentiable definition of acyclicity for structure learning is innovative.

Significance:

This work has the potential to improve the understanding of relationships between different clinical time series, which could lead to better diagnosis and treatment of diseases.
The method is applicable to a wide range of clinical time series data, making it potentially impactful.
Pros:

The proposed method can handle irregularly sampled data, which is common in EHRs.
The method is able to learn both the direction and strength of dependencies between time series.
The simulations show promising results in terms of accuracy.
Cons:

The method is computationally expensive due to the grid search and proximal gradient steps.
The performance of the method may be sensitive to the specific scaling of the data and the presence of unmeasured confounders.
The paper does not address how the method would scale to very large datasets.

Overall:

This is a well-written and innovative paper that proposes a promising method for learning the dependency structure of clinical time series data. Further research is needed to address the limitations mentioned above, but this work has the potential to be a valuable tool in clinical research.
title: title: This paper introduces StructGP, a structure learning algorithm for irregularly sampled clinical time series based on a multi-task Gaussian process (GP) model.
main_review: main_review: -The StructGP model is specifically designed for clinical time series, which are often irregularly sampled. Many existing algorithms struggle with non-uniform data, but StructGP effectively models these challenges using continuous-time Gaussian processes.
-The incorporation of a Directed Acyclic Graph structure allows StructGP to model ordered conditional dependencies, offering potential for causal inference in time series data. This is valuable in clinical applications, where understanding cause-effect relationships can influence decision-making.
Limitations:
-While recall is strong, precision is lower at 0.71%, indicating a higher number of false positive dependencies. This can lead to overfitting in real-world clinical data where spurious connections may be identified, potentially introducing misleading conclusions.
-StructGP assumes that the observed time series capture the full dependency structure, but clinical data often contain unmeasured confounders or latent variables. The model might not be robust in the presence of such confounders, leading to biased or incorrect dependency structures.
-While the NOTEARS-based optimization is efficient for small to medium graphs (up to 20 tasks), scaling it up to larger graphs or datasets could lead to significant computational costs. The grid search for optimal sparsity parameters also adds to the computational burden.