Abstract Effectively stratifying patient risk in chronic diseases like glaucoma is a major clinical challenge. Clinicians need tools to identify patients at high risk of progression from sparse and irregularly-sampled electronic health records (EHRs). We propose a novel deep kernel learning (DKL) architecture that leverages a Gaussian Process (GP) backend. The GP’s kernel is defined by a transformer-based feature extractor applied to clinical-BERT embeddings to model glaucoma patient trajectories from multimodal EHR data. Our method successfully identifies three clinically distinct patient subgroups. Crucially, the model learns to decouple disease progression from current severity, identifying a high-risk group with a worsening trajectory despite having better average visual acuity than a second, stably poor group. This reveals that the model learns to identify progression risk rather than just the current disease state. This ability to stratify patients based on their risk trajectory progression offers a powerful tool for clinical decision support, enabling targeted interventions for high-risk individuals and improving the management of glaucoma care.

1. Introduction Irregular time-series data is a major problem when analyzing longitudinal data sources such as electronic health records (EHRs) (Che et al., 2018; Shukla and Marlin, 2019b; Zhang et al., 2022). Entries with features that vary with visit and that occur at different intervals make it extremely challenging to do forecasting from such data (Karami et al., 2024; Si et al., 2021). This presents even more problems when forecasting models are leveraged to identify different patient clusters; clinicians are often interested in diagnosing when patients might be primed for improvement or downward health trends based on their history of observations, so utilizing models to forecast relevant diagnostic criteria and then cluster patients based on those models’ representations can be extraordinarily useful for clinical practice (Schulam et al., 2015; Feuerriegel et al., 2021). The upshot for medical practitioners on having forecasting methods that can be used to identify patient clinical

Deep Kernel Learning for Glaucoma

trajectories prior to negative trajectories is improved predictive methods that would enable better therapeutic treatments (Placido et al., 2025; Liu et al., 2019). In addition to the sparsity of features and observations over time, there are two further problems that confound time-series forecasts for EHRs. First, the data is highdimensional—especially when including natural language features like clinical notes (Rajkomar et al., 2018; Alsentzer et al., 2019). This can be problematic for traditional methods in irregular time-series forecasts and patient clustering strategies like Gaussian Processes (GP) that suffer high time and memory complexity in dataset size and dimensionality (Williams and Rasmussen, 2006; Schulam et al., 2015). Second, if the goal of the forecast is to identify patient trajectories, then the number of samples of records for most patients will be low. State-of-the-art (SOTA) methods for EHR forecasting that leverage large language models are anything but sample efficient, making them poor candidates for effective trajectory modeling based on limited patient records (Yang et al., 2022; Rasmy et al., 2021). We offer a hybrid solution that employs pre-trained language model embeddings along with a recurrent or transformer-based deep kernel learning (DKL) algorithm to simultaneously address the problems of data irregularity, sparsity, high-dimensionality, and patientrelative low sample size (Wilson et al., 2016; Alsentzer et al., 2019). Notably, our method does not rely upon imputation with temporal process points that are common for addressing data irregularity. Using the SOURCE glaucoma dataset of EHRs, we successfully forecast irregular time-series data for predicting glaucoma patient visual acuity loss. Our method either matches or outperforms most traditional recurrent neural network and transformerbased forecasting methods. We use the GP backend of our DKL algorithm to identify three clinically distinct patient trajectories that reveal a critical insight: patients with moderate vision but high trajectory variance represent high-risk group requiring intensive monitoring, despite those patients’ better vision compared to chronically poor but stable patients. This offers expanded opportunities to identify high-risk patients by decoupling current disease severity from future progression risk, enabling targeted interventions for those most likely to experience vision loss. Our main contributions are as follows: • A novel hybrid architecture that pre-processes EHR natural language features through clinical-BERT embeddings and then trains a DKL algorithm with a transformerencoder feature extractor to predict patient visual acuity loss. • We demonstrate with experiments better hold-out test metrics on the SOURCE dataset using our method compared to standard time-series recurrent neural networks and transformers. • We show that our method identifies three clinically meaningful patient subgroups by decoupling current disease state from progression risk, with the models providing calibrated uncertainty estimates—achieving 53.06% accuracy within 0.1 logMAR. • Our framework provides calibrated uncertainty quantification that enables risk-aware clinical decision support, with high-variance predictions correctly identifying patients with volatile disease courses rather than just reflecting model uncertainty.

Deep Kernel Learning for Glaucoma

Generalizable Insights: We provide probabilistic forecasts that can aid in handling uncertainty in healthcare. By pairing a transformer encoder with a Gaussian Process head, probabilistic deep learning can ingest irregularly sampled, multimodal records, such as unstructured clinical notes, and return calibrated uncertainty that distinguishes irreducible patient variability from model ignorance. Our modeling approach shows that a patient’s current disease severity poorly predicts future trajectory: unsupervised clustering of the learned representations recovers progression subtypes consistent with landmark glaucoma trials, without access to gold-standard structural or functional measurements. We can leverage the uncertainty provided by our method to flag volatile patients for further monitoring. By quantifying uncertainty honestly over messy longitudinal data and decoupling disease state from progression risk, this method can apply to chronic conditions from diabetes to heart failure, where irregular follow-up and heterogeneous trajectories are the rule.

2. Related Work Managing irregular time-series data has been an ongoing project in machine learning, with recent work focusing on leveraging deep learning architectures. While feedforward, convolutional, and recurrent networks commanded the most attention (Lim et al., 2021; Ling et al., 2025), transformers (Vaswani et al., 2017) now predominate. Custom transformer architectures for time-series include modified positional encoders and attention mechanisms (Wen et al., 2023). One prominent line of research involves adjusting the positional encoder with a Hawkes process to integrate interval information (Zhang et al., 2020; Zuo et al., 2020), while sophisticated embeddings can further improve performance (Mei et al., 2022). Recent work converts irregular time-series into line-graphs for vision transformer processing (Li et al., 2023; Dosovitskiy et al., 2021). Machine learning for EHRs has advanced significantly (Si et al., 2021; Shickel et al., 2017). Early Probabilistic Subtyping Models addressed EHR irregularity and sparsity (Schulam et al., 2015), and the subsequent application of deep learning has led to more powerful models. GRU-D (Che et al., 2018) handles irregular data by explicitly modeling informative missingness—missing data salient for predictions. Recurrent networks address high-dimensional EHR data (Rajkomar et al., 2018), despite reproducibility challenges (Si et al., 2021). Attention has increasingly focused on transformers: transformer encoders with point processes can successfully classify irregular EHR data (Karami et al., 2024), and multi-modal EHR data has been integrated to produce better predictions (Zhang et al., 2022). For glaucoma, EHR-based surgery prediction uses regression, decision trees, and neural networks (Wang et al., 2024), including word embedding approaches (Wang et al., 2022). Prior EHR work emphasizes regression or classification. However, clinicians need to group patient clinical history into trajectories for predicting disease progression (Yohannan et al., 2021; Feuerriegel et al., 2021). EHR-based trajectory modeling remains challenging. Gaussian processes address this challenge (Williams and Rasmussen, 2006; Hensman et al., 2013). Disease Trajectory Maps (Schulam and Arora, 2016) use Gaussian processes with low-dimensional clinical markers to model autoimmune progression through latent spaces. More recent work has utilized transformers combined with probabilistic models to measure disease progression and cluster trajectories (Qiu et al., 2025). Their model uses a

Deep Kernel Learning for Glaucoma

transformer-based variational autoencoder with Gaussian mixture priors to jointly model patient clustering and survival times from diagnosis sequences. While this approach successfully identifies disease subgroups with divergent trajectories, it focuses primarily on structured diagnosis codes rather than leveraging the rich semantic information in clinical notes and uses mixture models instead of Gaussian processes for uncertainty quantification.

3. Preliminaries We denote sets as A, random variables and assignments as X and x , random vectors and assignments as X and x, and matrices as A and Σ. Given dataset D = {(xi , yi )}i=1,...,n with features X and targets Y, Bayesian inference computes the posterior predictive distribution for parametric model f : X → Y characterized by parameters w: Z p(y | x, D) = p(y | x; w)p(w | D)dw (1) where p(w | D) ∝ p(D | w)p(w). For Bayesian linear regression, we assume a likelihood N y | x⊺ w, σ 2 and a zero-mean Gaussian prior on the weights, w ∼ N (0, Σp ). Let X and y denote the stacked features and targets from D. For nonlinear mappings, we project inputs ϕ : X → Rd yielding f (x; w) 7→ ϕ(x)⊺ w. With design matrix Φ, the posterior becomes: p(y | x, D) = N y | ϕ⊺x Σp Φ(K + σ 2 I)−1 y, ϕ⊺x Σp ϕx − ϕ⊺x Σp Φ(K + σ 2 I)−1 Φ⊺ Σp ϕx (2) where ϕx = ϕ(x) and K = Φ⊺ Σp Φ. 1/2 Since ϕ and Φ only occur in inner products with Σp , we define kernel k(x, x′ ) = Σp ϕ(x)· 1/2 Σp ϕ(x′ ), yielding a Gaussian process (GP) with mean m(x) 7→ E[f (x)]: f (x) ∼ GP(m(x), k(x, x′ ))

(3)

The common squared exponential kernel is: ∥x − x′ ∥2 kse (x, x ) = exp − 2l ′

(4)

where l is the length parameter. Applying the kernel element-wise over features X and input points X∗ , we get matrix K(·, ·). For noisy observations with zero mean prior, the prior is: y K(X, X) + σ 2 I K(X, X∗ ) ∼ N 0, (5) f∗ K(X∗ , X) K(X∗ , X∗ ) and the posterior predictive is: f∗ | X∗ , D ∼ N f | K(X∗ , X)Λ−1 y, K(X∗ , X∗ ) − K(X∗ , X)Λ−1 K(X, X∗ )

(6)

where Λ = K(X, X)+σ 2 I. Hyperparameter l for our kernel kse is found via empirical Bayes, i.e. maximizing the marginal likelihood log p(y | X, l) = log N (y | 0, Λ) = − 21 y⊺ Λ−1 y − 1 n 2 log|Λ|− 2 log 2π. Computing the posterior requires inverting Λ with O(n3 ) time and O(n2 ) space complexity. This makes utilizing GPs impractical for large datasets with high-dimensional features. 4

Deep Kernel Learning for Glaucoma

Common solutions for this problem include: variational inference with inducing points for large datasets, and neural network dimensionality reduction for high-dimensional features. For scalability, we partition features into training X, inducing point Z, and test X∗ , with joint variational distribution q(fX , fZ , fX∗ ) ≈ p(fX , fX∗ | fZ )q(fZ ). Our optimization target is to minimize the KL divergence between the variational and true posterior by the Evidence Lower Bound (ELBO) (Hensman et al., 2015): L(q) = Eq(fX ) [log p(y | fX )] − DKL (q(fZ )∥p(fZ ))

(7)

To compute the first part of the ELBO, we utilize stochastic gradient descent over minibatches. For high-dimensional feature spaces, we map h : Rd → Rk (k ≪ d) before computing kernels: ∥h(x) − h(x′ )∥2 kdkl se (x, x′ ) = exp − (8) 2l A neural network hθ learns this mapping, yielding deep kernel learning (Wilson et al., 2016). Parameters θ are treated as hyperparameters for the GP and jointly optimized with other hyperparameters via 7. The increased hyperparameters risk overfitting (Ober et al., 2021), which can be mitigated by mini-batching and additional regularization.

4. Disease Trajectories with Deep Kernel Learning EHRs induce challenges for time-series forecasting due to irregular sampling and missing features. GPs address irregular sampling by interpolating through prior and posterior estimates, and they handle missing values naturally by increased posterior variance, providing broader uncertainty estimates. The functions modeled by the multivariate normal distribution in a GP also enable modeling of different types of the underlying generating population. This has made GPs especially useful in the medical domain for tracking population-level disease trajectories (Schulam and Arora, 2016) as the fitted posteriors allow for the clustering of entire trajectory functions. High-dimensional EHR features pose computational challenges for GPs. First, they cannot effectively learn representations from clinical notes. Pre-trained embeddings, such as word2vec (Mikolov et al., 2013), provide usable features, but their high dimensionality create O(n3 ) computational bottlenecks in kernel inversion. Second, time-series architectures like LSTM (Hochreiter and Schmidhuber, 1997), GRU (Cho et al., 2014), and transformers (Vaswani et al., 2017) which could serve as feature extractors for a DKL model, often assume regularly sampled data. Thus, GP-based solutions require both rich features and computationally tractable representations. Combining pre-trained representations with a DKL algorithm provides effective dimensionality reduction. We propose the following pipeline. Clinical-BERT (Huang et al., 2019) first embeds linguistic features from EHRs, such as diagnoses and medications. Features are then aggregated into Rd . Observations are mapped into sequences x1:t ∈ Rt×d for lengths t = 1, . . . , T for a population member with T total observations. An encoder encθ : Rt×d → Rt×k with parameters θ and decoder decθ′ : Rt×k → Rm with parameters θ′ compress features to a latent space with dimension m ≪ d at time t. The resulting feature 5

Deep Kernel Learning for Glaucoma

Deep Kernel Learning Pipeline for Patient Trajectory Clustering

Source EHR Data

Preprocessing Pipeline

Full Dataset De-identification & anonymization

Snellen-to-logMAR conversion

Random Sampling N = 3821 patients Text features

Clinical-BERT Embedding

Patient Sequences

Structured features

Embeddings

Deep Kernel Learning Model

Deep Kernel Learning Transformer Latent Representations

Latent Z

Posterior µ

Posterior σ 2

Agglomerative Clustering

Patient Trajectories

Agglomerative Ward Clustering

(a) Deep Kernel Learning Pipeline

(b) Clustering Pipeline for Latent Space Representations

Figure 1: Deep Kernel Learning Transformer Pipeline for Disease Trajectory Maps. (a) The transformer architecture processes multimodal EHR data through Clinical-BERT embeddings and structured feature extraction. (b) Agglomerative clustering with ward linkage is applied to latent representations to identify distinct patient trajectories, enabling population-level analysis of clinical patterns and outcomes. Model selection was performed across c ∈ {2, 3, 4, 5} clusters, with c = 3 providing optimal separation.

extractor hf e : Rt×d → Rm , is then used in a kernel kf e for a backend GP with a constant mean m(x): f (x1:t ) ∼ GP(m(x1:t ), kf e (x1:t , x′1:t )) ∥hf e (x1:t ) − hf e (x′1:t )∥2 kf e (x1:t , x′1:t ) = exp − 2l hf e (x1:t ) = decθ′ (encθ (x1:t ))

(9) (10) (11)

We require only temporal ordering: t > t′ implies t occurs after t′ . With fixed embeddings, we optimize the feature extractor parameters θ ∪ θ′ via variational inference Wilson et al. (2016); Al-Shedivat et al. (2017). While the GP layer in the DKL model naturally handles irregularly sampled time-series data, temporal information still needs to be processed through the DKL encoder. One such strategy is to use recurrent architectures such LSTM and GRU (Al-Shedivat et al., 2017); another strategy is to utilize transformers (Lyu et al., 2023). For time-series data, we often impute missing entries to improve the regularity of the prediction, but due to the GP layer, we can avoid that strategy and instead adopt explicit temporal features by cyclical encoding (Fisher, 1995): 6

Deep Kernel Learning for Glaucoma

Model

MSE

MAE

+/− 0.1 logMAR

MLE RNN DKL RNN

0.1166 0.1075

0.1795 0.1721

0.2594 0.3172

53.8397 56.9478

MLE LSTM DKL LSTM

0.1156 0.1184

0.1718 0.1712

0.2659 0.2477

57.9248 59.6589

MLE GRU DKL GRU

0.1106 0.1113

0.1708 0.1748

0.2973 0.2934

56.41 53.7741

MLE Transformer DKL Transformer MLE GRU-D MLE GRU ODE MLE RNN CDE

0.1068 0.1061 0.1246 0.1166 0.1246

0.171 0.1703 0.1765 0.1888 0.2017

0.3217 0.3262 0.2084 0.2599 0.2109

53.3605 53.0564 56.8392 43.8709 36.2228

Table 1: Average final hold-out test set performance comparison of DKL algorithms trained via variational inference and baseline models using standard maximum likelihood training with stochastic gradient descent across different neural network architectures. All models were run across 10 random seeds and the final test set scores averaged. DKL models demonstrate competitive performance on glaucoma patient data, with transformer-based architectures achieving the lowest mean squared error, mean absolute error, and highest coefficient of determination. The last column is the percentage of predictions within 0.1 logMAR, which indicates clinical accuracy, where LSTM architectures show superior precision in visual acuity prediction for glaucoma patients.

fsin (t) = sin fcos (t) = cos

2πt P

(12) (13)

where P is the period (31 for days, 12 for months, and 10 for years). This was performed so that we could better compare DKL algorithms with their standard counterparts, since normal architectures struggle with irregular time-series data.

5. Experiments Our methods were tested on electronic health records (EHRs) from the Sight Outcomes Research Collaborative (SOURCE) (Bommakanti et al., 2020; Wang et al., 2024), a multiinstitutional EHR initiative led by the University of Michigan that supports the investigation of ophthalmic outcomes. We restricted our analysis to patients receiving care within the [health system redacted for anonymous review] between January 2016 and December 2023.

Deep Kernel Learning for Glaucoma

Latent Space by Cluster (95% Confidence Ellipses) 1.0

UMAP of Latent Visit Embeddings (2250 sampled patients shown as points)

Cluster

Cluster 0 (n=750) Cluster 1 (n=750) Cluster 2 (n=750)

Cluster 0 (n=821) Cluster 1 (n=5027) Cluster 2 (n=61843)

0.8 15

0.6

UMAP 2

C2 5

0.4

C1 C0

0.2

C1 C2

0.0 0.0

0.2

0.4

0.6

0.8

1.0

(a) Latent Space Representation of Patient Trajectories

UMAP 1

(b) UMAP Projection

Figure 2: Latent space visualizations of patient trajectory clustering. (a) Direct latent space visualization demonstrates nonlinear progression patterns along a curved manifold structure with clear separation between three disease progression archetypes. (b) UMAP dimensionality reduction confirms three unique clinical trajectories rather than a continuous spectrum of disease progression.

We performed data preprocessing and ran experiments on a variety of architectures and used the best architectures for clustering done via functional data analysis. 5.1. Data The dataset comprises 67,691 patients with 402,552 single-day clinical encounters (henceforth, “records”) across 18 fields. Structured features include demographics (age, sex, race, ethnicity), temporal data (encounter date, birth date), and clinical metadata (living status, patient ID). Visit features capture specialty, visit type, and free-text reason. Clinical features encompass diagnoses, procedures, medications, lab results, and surgical history. The outcome variable, ACUITY SCORE, represents best-corrected visual acuity converted from Snellen to logMAR format. Lower logMAR values indicate better visual acuity; higher logMAR values indicate poorer acuity. For each eye, the best-recorded acuity was defined as the minimum logMAR (corresponding to the smallest Snellen denominator), and the worst acuity as the maximum logMAR (largest Snellen denominator). Data was anonymized through the removal of all patient identifiers and personally identifiable information before processing and other identifiers were randomized to protect patient privacy. Preprocessing handled irregular sampling via duplicate removal, temporal alignment on encounter datetime, and standardized missing values (-1 for numeric, empty strings for text). Visual acuity conversion handled standard fractions (e.g., 20/40) and special cases

Deep Kernel Learning for Glaucoma

Visual Acuity Trajectories - Cluster 0

SHAP

1.5 1.0 0.5 0.0

Visual Acuity Trajectories - Cluster 1

Visual Acuity (logMAR) (Higher = Worse Vision)

+-1 SD (between-patient) +-1 GP SD (model uncertainty) Cluster 1 Mean (n=5027)

2.0 1.5 1.0 0.5 0.0

Visual Acuity Trajectories - Cluster 2

Visual Acuity (logMAR) (Higher = Worse Vision)

+-1 SD (between-patient) +-1 GP SD (model uncertainty) Cluster 2 Mean (n=61843)

2.0 1.5 1.0 0.5 0.0 0

100

150

Visit Index (0 = first visit)

High

SURGERIES_743 VISIT_TYPE_366 SURGERIES_669 LAB_RESULTS_69 DX_NAME_499 REASON_FOR_VISIT_559 SURGERIES_13 SURGERIES_729 PROC_NAME_473 SURGERIES_628 REASON_FOR_VISIT_230 SPECIALTY_155 SPECIALTY_142 DX_NAME_759 SURGERIES_295 DX_NAME_474 REASON_FOR_VISIT_538 PROC_NAME_582 LAB_RESULTS_351 SPECIALTY_77

Feature value

Visual Acuity (logMAR) (Higher = Worse Vision)

+-1 SD (between-patient) +-1 GP SD (model uncertainty) Cluster 0 Mean (n=821)

2.0

1.5

200

1.0

0.5

0.0

0.5

1.0

1.5

Low

SHAP value (impact on model output) 1e13

(a) Visual Acuity Trajectories

(b) SHAP Value Analysis

Figure 3: Clinical trajectory analysis and model interpretability. (a) Posterior predictive mean trajectories from DKL transformer demonstrate three distinct patient archetypes with characteristic visual acuity patterns. Mean and standard deviation are conditioned on entire cluster. (b) SHAP analysis reveals surgeryrelated features and specialty codes as primary drivers of model predictions.

(NLP, LP, HM, CF). ACUITY SCORE used the minimum logMAR across 18 measurement types. Clinical text was embedded using Bio ClinicalBERT with mean pooling, yielding 768dimensional representations (batch size: 32, max tokens: 512). Seven text fields were embedded: specialty, visit type, reason, procedures, diagnoses, medications, and lab results. Lab results combined names with values before embedding; medications and procedures were concatenated. This approach preserved clinical semantics while enabling both GP uncertainty quantification and neural feature extraction. 5.2. Methods We compared recurrent neural network (RNN), gated recurrent unit (GRU), long short-term memory (LSTM), and transformer encoder architectures as feature extractors. Baseline models used identical architectures but replaced the GP layer with linear regression trained via maximum likelihood estimation (MLE). Additionally, we trained a GRU-D (Che et al., 2018), GRU ODE (Chen et al., 2018; De Brouwer et al., 2019), and RNN CDE (Kidger et al., 2020) as baselines. DKL models were trained by maximizing the ELBO using stochastic

Deep Kernel Learning for Glaucoma

variational inference. We used a multivariate normal variational distribution with Choleskydecomposed covariances. Our likelihood was a univariate normal distribution. Both model types used cyclical encoding for temporal features and standardized continuous variables. Cyclical encoding provided temporal inductive bias for baselines while enabling GP benefit assessment. Age and embedding features were z-normalized. Hyperparameters were selected via grid search on 100,000 samples. See Appendix A for search details and configurations. We evaluated 1,306 hyperparameter combinations. Models were selected based on the best validation set mean-squared error (MSE). Models were evaluated using a 70/20/10 train/validation/test split. Performance was assessed with MSE, MAE, R2 , and clinical accuracy, defined as the proportion of predictions within ±0.1 logMAR of the observed value. For the DKL models, probabilistic performance was additionally evaluated using 95% credible interval coverage and the continuous ranked probability score (CRPS). All experiments were repeated across 10 random seeds, and we report the mean and standard deviation across seeds. The MLE and DKL RNN, LSTM, GRU, and Transformer models, as well as the MLE GRU-D model, were trained for 200 epochs per seed. Due to longer training times and faster convergence, the GRU-ODE and RNN-CDE models were trained for 25 epochs per seed. The full pipeline, from start-to-finish can be seen in figure 1(a). For clustering, we relied upon functional data analysis and agglomerative ward clustering. The clustering pipeline can be seen in figure 1(b). After training the model, we collected the model’s feature extractor latent representation as well as the posterior predictive mean and variance for each record in the patient’s sequence of records. We clustered all of the patient trajectory representations using agglomerative ward clustering with the L2 norm for distance measurements between trajectories, and we evaluated c = {2, 3, 4, 5} clusters, with c = 3 clusters ultimately selected due to best clustering metrics and clinical interpretation (see appendix C for more details). 5.3. Results DKL models with transformer extractor outperformed baselines on most metrics (see table 1). The DKL LSTM outperformed all other methods on clinical accuracy (59.66%). Compared to MLE, we see that the transformer and RNN feature extractors outperform their MLE counterparts across 10 random seeds, with the RNN improvements being substantial. We found the DKL RNN, GRU, and transformer along with MLE GRU and transformer with encoded temporal features outperformed GRU-D, GRU ODE, and RNN CDE baselines substantially—despite these baselines specialization for irregular time series. See Appendix D for further analyses. We used the best-performing DKL transformer for trajectory analysis. Figure 2(a) shows three distinct patient clusters. The latent space reveals three patient archetypes: worsening, stable, and intermediate trajectories (see Appendix B for 3D visualizations). UMAP visualization (Figure 2(b)) confirms cluster separation in latent space. Figure 3(a) shows cluster-specific posterior means and uncertainties over time. Cluster 1 shows progressive worsening before stabilization and exhibits the highest uncertainty. Higher variance reflects heterogeneous progression patterns within this cluster.

Deep Kernel Learning for Glaucoma

To evaluate clinical markers, we performed SHAP analysis (Lundberg and Lee, 2017) to assess feature relevance for model predictions on the DKL transformer. We found that surgery was most predictive in figure 3(b). To further validate model performance, we experimented with ablations on these features. Table 2 shows our results compared to non-ablated model average performance. We found the largest drop in performance when ablating medicine and surgery names.

6. Discussion The central challenge of applying machine learning to real-world clinical data is not merely managing sparsity or irregularity, but embracing the inherent variability of chronic disease. Patient trajectories are shaped by complex, often unobserved factors, leading to high interand intra-patient variance that is sometimes irreducible, even with complete data. Here, we present a probabilistic deep learning framework that, rather than being confounded by this “messy” data, leverages it to uncover clinically significant insights. We demonstrate that our approach moves beyond static risk assessment by learning to decouple a patient’s current disease state from their future risk of progression. Handling Irreducible Uncertainty in Sequential Medical Data. Our deep kernel learning architecture addresses a fundamental gap in medical ML: most approaches to disease progression either assume regular sampling (violated in real-world EHRs) or provide point estimates without uncertainty quantification (inadequate for clinical decision-making) (Shukla and Marlin, 2019a). By combining a transformer-based feature extractor with a Gaussian Process backend, we achieve three critical capabilities: (1) representation learning from high-dimensional, multimodal data including unstructured text, (2) natural handling of irregular temporal sampling without imputation, and (3) calibrated uncertainty estimates that distinguish between aleatory uncertainty (inherent disease variability) and epistemic uncertainty (model confidence). This distinction is crucial—high variance in a patient’s trajectory may indicate genuinely unpredictable disease behavior rather than model limitations, requiring fundamentally different clinical responses. The architecture’s ability to process sequences of varying lengths with missing features while maintaining uncertainty calibration represents a significant advance over existing methods that typically handle these challenges separately or incompletely. Decoupling Progression Risk from Disease State. Our primary finding validates the model’s ability to learn clinically meaningful representations: clustering in the DKL latent space reveals three distinct patient trajectories that decouple current disease severity from progression risk. Counterintuitively, the group with the worst average visual acuity was not the group with the most severe predicted progression. This demonstrates that our model learns temporal dynamics beyond simple trend extrapolation. This learned separation between current severity and future risk recapitulates a central insight from decades of prospective glaucoma clinical trials. The Canadian Glaucoma Study identified higher baseline IOP, older age, and disc hemorrhages as predictors of open-angle glaucoma progression—but not baseline visual field severity alone (Chauhan et al., 2008). Similarly, the Early Manifest Glaucoma Trial (EMGT) found that baseline IOP, pseudoexfoliation, bilateral disease, and age predicted progression, while 45% of treated patients still progressed despite IOP lower-

Deep Kernel Learning for Glaucoma

Feature Category Removed

MSE

None (Baseline)

0.1081

MED NAME SURGERIES VISIT TYPE PROC NAME DX NAME SPECIALTY REASON FOR VISIT AGE LAB RESULTS

0.1131 0.1126 0.1119 0.1117 0.1100 0.1079 0.1072 0.1071 0.1063

Table 2: Ablation study results showing model performance (MSE) when individual feature categories are removed. Lower MSE indicates better performance. The baseline model (all features) achieves MSE = 0.1081. Removing medication names (MED NAME) or surgeries causes the largest performance degradation, indicating these are the most critical feature categories for prediction accuracy. All experiments on fixed random seed (42).

ing, underscoring that severity at presentation does not determine trajectory (Leske et al., 2003, 2007). The Ocular Hypertension Treatment Study (OHTS) demonstrated that among patients with no baseline glaucomatous damage, a constellation of risk factors — including IOP, central corneal thickness, vertical cup-to-disc ratio, and pattern standard deviation on visual field testing — predicted conversion to glaucoma, establishing that risk is multifactorial and not solely determined by the current disease state (Gordon et al., 2002). A large retrospective analysis combining data from AGIS and a clinical cohort found that after adjusting for IOP, baseline visual field severity did not independently predict the rate of visual field change, further supporting the dissociation between where a patient is and where they are going (De Moraes et al., 2011). Collectively, these trials demonstrate that baseline disease severity is a poor proxy for future progression risk — a principle well understood by glaucoma specialists, but difficult to operationalize in routine clinical practice, particularly at the population level. Our framework provides a computational approach to this clinical problem. In the current research, we identify: (1) a “Stable-Good Vision” group with low vision loss and tight confidence intervals consistent with well-controlled or non-progressive disease, (2) a “High-Progression Risk” group with moderate vision but high trajectory variance, consistent with volatile disease courses requiring intensive monitoring and potential escalation of treatment, and (3) a “Stable-Poor Vision” group with chronically poor but non-worsening vision, potentially representing patients with end-stage but quiescent disease, or patients whose poor acuity reflects non-glaucomatous pathology such as visually significant cataract or macular disease. That an unsupervised model — without access to IOP measurements, visual fields, or OCT data — recovers trajectory subtypes consistent with the progression risk

Deep Kernel Learning for Glaucoma

factors identified by landmark clinical trials provides external validation that our learned representations capture genuine disease heterogeneity, rather than statistical artifacts. Additionally, it suggests that the multimodal EHR features processed through Clinical-BERT embeddings encode clinically meaningful information about disease dynamics that is partially redundant with, and may complement, the structured clinical measurements used in traditional risk calculators. Quantifying and Leveraging Trajectory Uncertainty. The “High-Progression Risk” cluster exemplifies our framework’s novel contribution to uncertainty-aware ML. These patients exhibit high trajectory variance—not due to sparse data or model uncertainty, but reflecting genuinely volatile disease courses. Traditional ML approaches would either smooth over this variability (losing critical information) or produce unstable predictions (limiting clinical utility). Our GP backend quantifies this uncertainty, enabling risk-aware decision support: patients with wide posterior credible intervals trigger more frequent monitoring, while those with narrow bands support extended follow-up. This moves beyond binary classification to provide nuanced, uncertainty-aware recommendations that align with clinical decision-making under uncertainty. The model achieves 53.06% accuracy within 0.1 logMAR (approximately one line on the eye chart), exceeding typical progression prediction tools that achieve AUCs of 0.65-0.75 (Yohannan et al., 2017), while additionally providing calibrated uncertainty estimates. Notably, the clinical significance of a ±0.1 logMAR prediction error varies across the acuity spectrum. In the context of glaucoma trajectory prediction — where the goal is to identify patients trending toward visual impairment (logMAR ≥0.5) or legal blindness (logMAR ≥1.0) — this level of accuracy is sufficient for clinically meaningful risk stratification, even if insufficient for precise acuity prediction. Digital Biomarkers from Multimodal Temporal Learning. The learned patient trajectories represent novel digital biomarkers that integrate the full complexity of longitudinal clinical data (Purushotham et al., 2018). Unlike traditional risk factors (central corneal thickness, cup-to-disc ratio, baseline IOP (Gordon et al., 2002)), our approach leverages Clinical-BERT embeddings to incorporate unstructured clinical notes—capturing subtle observations about disease stability, treatment response, and psychosocial factors typically unavailable to structured models. The transformer architecture learns temporal dependencies across irregular visit sequences, while the DKL framework maps these highdimensional histories to interpretable low-dimensional trajectories. This demonstrates that modern NLP and sequence modeling techniques can extract clinically meaningful signals from the “messy” reality of EHR data, where crucial information is often buried in narrative text rather than structured fields. 6.1. Generalizable Framework for Chronic Disease Trajectories. While demonstrated on glaucoma, our framework addresses challenges universal to chronic disease modeling. The combination of irregular sampling, high-dimensional multimodal data, and irreducible patient-level variability characterizes conditions from diabetes to heart failure to rheumatoid arthritis. Our architecture’s key innovations provide a template for chronic disease analytics. The success in identifying clinically meaningful subtypes without supervision suggests that similar latent structure may exist across chronic conditions,

Deep Kernel Learning for Glaucoma

waiting to be discovered by appropriate ML methods. Future work could explore transfer learning across diseases, potentially sharing learned representations about healthcare utilization patterns, treatment response, and disease progression dynamics. Given glaucoma’s disproportionate impact on African American and Latino populations (prevalence 3-4 times higher than Caucasians (Tielsch et al., 1991)), our work explicitly addresses algorithmic fairness—a critical requirement for clinical ML deployment. Our framework enables systematic fairness audits by: (1) including demographic variables as model inputs, (2) analyzing performance metrics stratified by protected attributes, (3) examining cluster membership distributions across demographic groups, and (4) quantifying whether trajectory uncertainty correlates with social determinants of health. This transparency is essential for responsible AI in healthcare, moving beyond post-hoc bias detection to proactive fairness-aware modeling. The probabilistic nature of our approach allows uncertainty to increase in underrepresented populations, honestly reflecting data limitations rather than producing falsely confident predictions. Preliminary analysis of cluster demographics is ongoing. Given that the SOURCE dataset includes patients from the University of Virginia health system, which serves a diverse catchment area including rural Appalachian communities and urban populations, characterizing whether the identified trajectory subtypes differ across demographic groups is essential for ensuring equitable clinical deployment. Our current limitations highlight important directions for medical ML research. The lack of explicit inter-eye correlation modeling (common in bilateral diseases) motivates development of multi-output GP architectures for correlated outcomes. Reliance on visual acuity rather than gold-standard endpoints (visual fields, OCT measurements) demonstrates the need for methods that can align heterogeneous clinical endpoints with different sampling frequencies and measurement noise characteristics. An important current limitation of this study is the use of visual acuity as the primary outcome measure for glaucoma trajectory modeling. Visual acuity is a late and non-specific indicator of glaucomatous damage— patients can lose substantial peripheral vision and retinal nerve fiber layer thickness while maintaining 20/20 central acuity. VA loss in glaucoma typically occurs only with advanced disease involving central visual field damage, or through comorbid conditions such as retinal disease or visually significant cataract. Consequently, the trajectories identified here likely capture a composite of glaucoma progression, cataract development, and other causes of central vision loss rather than glaucoma-specific progression alone. Gold-standard glaucoma progression measures — including visual field mean deviation (MD) trends, OCT retinal nerve fiber layer (RNFL) thinning rates, and intraocular pressure (IOP) trajectories—were not available for this current analysis, but would substantially improve the clinical specificity of identified patient subtypes. Future iterations of the present research incorporating these modalities will enable earlier detection of progression, as structural (OCT) and functional (VF) changes precede VA loss by years in most glaucoma patients. The inability to distinguish glaucoma-driven visual acuity loss from cataract-driven visual acuity loss represents a significant confounding factor in trajectory interpretation. Cataract is highly prevalent in the age demographic most commonly affected by glaucoma, and cataract surgery (often combined with glaucoma procedures) can dramatically improve VA. Patients in the ’volatile’ cluster may potentially include individuals whose VA fluctuations reflect cataract progression followed by surgical improvement, rather than glaucomatous instability. Additionally, the strong SHAP signal from surgery-related features (Figure 3b) likely 14

Deep Kernel Learning for Glaucoma

captures both glaucoma-specific interventions and combined cataract-glaucoma procedures, which have fundamentally different clinical implications. Future work should ideally stratify patients by cataract status and separate combined cataract-glaucoma procedures from glaucoma-only surgeries to disambiguate these trajectories. Furthermore, the dominance of surgery-related features in the current SHAP analysis warrants careful causal interpretation. In observational EHR data, glaucoma surgery is a treatment response typically triggered by disease progression or failure of medical management, creating an inherent potential confounder: surgery correlates with poor outcomes partly because it is performed on patients already doing poorly. Though this paradigm is being challenged by more proactive ”interventional glaucoma” surgical management paradigms emphasizing early or proactive intervention with (usually) lower-risk minimally invasive glaucoma surgery (MIGS), the complex relationship between decision for surgery and poor outcome is analogous to ’confounding by indication’ in pharmacoepidemiology. The model likely learns surgery as a strong proxy for disease severity and prior progression rather than as a causal driver of outcomes. Disentangling this relationship — and ultimately modeling counterfactual trajectories (e.g., ’what if surgery had been performed earlier?’ or ’what if MIGS were tried rather than a more invasive surgery such as trabeculectomy or glaucoma tube shunt procedure?’) — is a key direction for developing actionable clinical decision support from this framework. The challenge of distinguishing disease-driven from opportunistic interventions (e.g., glaucoma-only vs. combined cataract-glaucoma surgery) represents a broader causal inference problem in observational health data. These limitations are not merely clinical details but fundamental ML challenges in learning from observational data with complex dependencies and confounding (Trottet et al., 2023).

References Maruan Al-Shedivat, Andrew Gordon Wilson, Yunus Saatchi, Zhiting Hu, and Eric P Xing. Learning scalable deep kernels with recurrent structure. Journal of Machine Learning Research, 18(82):1–37, 2017. Emily Alsentzer, John Murphy, William Boag, Wei-Hung Weng, Di Jindi, Tristan Naumann, and Matthew McDermott. Publicly available clinical bert embeddings. In Proceedings of the 2nd Clinical Natural Language Processing Workshop, pages 72–78, 2019. Nikhil K Bommakanti, Yunshu Zhou, Joshua R Ehrlich, Angela R Elam, Denise John, Shivani S Kamat, Jared Kelstrom, Paula Anne Newman-Casey, Manjool M Shah, Jennifer S Weizer, et al. Application of the sight outcomes research collaborative ophthalmology data repository for triaging patients with glaucoma and clinic appointments during pandemics such as covid-19. JAMA ophthalmology, 138(9):974–980, 2020. Balwantray C Chauhan, Frederick S Mikelberg, A Gordon Balaszi, Raymond P LeBlanc, Mark R Lesk, Graham E Trope, Canadian Glaucoma Study Group, et al. Canadian glaucoma study: 2. risk factors for the progression of open-angle glaucoma. Archives of ophthalmology, 126(8):1030–1036, 2008.

Deep Kernel Learning for Glaucoma

Zhengping Che, Sanjay Purushotham, Kyunghyun Cho, David Sontag, and Yan Liu. Recurrent neural networks for multivariate time series with missing values. Scientific reports, 8(1):6085, 2018. Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018. Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio. Learning phrase representations using rnn encoder-decoder for statistical machine translation. arXiv preprint arXiv:1406.1078, 2014. Edward De Brouwer, Jaak Simm, Adam Arany, and Yves Moreau. Gru-ode-bayes: Continuous modeling of sporadically-observed time series. Advances in neural information processing systems, 32, 2019. Carlos Gustavo V De Moraes, Viral J Juthani, Jeffrey M Liebmann, Christopher C Teng, Celso Tello, Remo Susanna Jr, and Robert Ritch. Risk factors for visual field progression in treated glaucoma. Archives of Ophthalmology, 129(5):562–568, 2011. Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby. An image is worth 16x16 words: Transformers for image recognition at scale. In International Conference on Learning Representations, 2021. URL https://openreview.net/forum?id=YicbFdNTTy. Stefan Feuerriegel, Dennis Frauen, Valentyn Melnychuk, Jonas Schweisthal, Konstantin Hess, Alicia Curth, Stefan Bauer, Niki Kilbertus, Isaac Papakonstantinou, and Mihaela van der Schaar. Analyzing patient trajectories with artificial intelligence. Journal of Medical Internet Research, 23(12):e29812, 2021. Nicholas I Fisher. Statistical analysis of circular data. cambridge university press, 1995. Mae O Gordon, Julia A Beiser, James D Brandt, Dale K Heuer, Eve J Higginbotham, Chris A Johnson, John L Keltner, J Philip Miller, Richard K Parrish, M Roy Wilson, et al. The ocular hypertension treatment study: baseline factors that predict the onset of primary open-angle glaucoma. Archives of ophthalmology, 120(6):714–720, 2002. James Hensman, Nicolò Fusi, and Neil D. Lawrence. Gaussian processes for big data. In Uncertainty in Artificial Intelligence, volume 29, pages 282–290. AUAI Press, 2013. URL http://auai.org/uai2013/prints/papers/244.pdf. James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational gaussian process classification. In Artificial intelligence and statistics, pages 351–360. PMLR, 2015. Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural computation, 9(8):1735–1780, 1997. 16

Deep Kernel Learning for Glaucoma

Kexin Huang, Jaan Altosaar, and Rajesh Ranganath. Clinicalbert: Modeling clinical notes and predicting hospital readmission. arXiv preprint arXiv:1904.05342, 2019. Hojjat Karami, Tyler Guo, David Svoboda, Ahmadreza Argha, Aditya Mahadevan, Richard Dobson, and Behraad Firouzi. Tee4ehr: Transformer event encoder for better representation learning in electronic health records. Artificial Intelligence in Medicine, 153:102909, 2024. Patrick Kidger, James Morrill, James Foster, and Terry Lyons. Neural controlled differential equations for irregular time series. Advances in neural information processing systems, 33:6696–6707, 2020. M Cristina Leske, Anders Heijl, Mohamed Hussein, Bo Bengtsson, Leslie Hyman, Eugene Komaroff, and Early Manifest Glaucoma Trial Group. Factors for glaucoma progression and the effect of treatment: the early manifest glaucoma trial. Archives of Ophthalmology, 121(1):48–56, 2003. M Cristina Leske, Anders Heijl, Leslie Hyman, Boel Bengtsson, LiMing Dong, Zhongming Yang, and EMGT Group. Predictors of long-term progression in the early manifest glaucoma trial. Ophthalmology, 114(11):1965–1972, 2007. Zekun Li, Shiyang Li, and Xifeng Yan. Time series as images: Vision transformer for irregularly sampled time series. Advances in Neural Information Processing Systems, 36: 49187–49204, 2023. Bryan Lim, Sercan Ö Arık, Nicolas Loeff, and Tomas Pfister. Time-series forecasting with deep learning: a survey. Philosophical Transactions of the Royal Society A, 379(2194): 20200209, 2021. Chen Ling, Xujiang Dong, Yangli Zhao, Yucheng Chen, Huan Liu, and Yanzhi Wei. A comprehensive survey of deep learning for time series forecasting: architectural diversity and open challenges. Artificial Intelligence Review, 57(1):1–47, 2025. Zhenyu Liu, Chen Chen, Luyu Li, Jiayu Zhou, Xiaoqian Li, Le Song, and Yuan Qi. Representation learning for clinical time series prediction tasks in electronic health records. BMC Medical Informatics and Decision Making, 19(1):1–14, 2019. Scott M Lundberg and Su-In Lee. A unified approach to interpreting model predictions. Advances in neural information processing systems, 30, 2017. Wenlong Lyu, Shoubo Hu, Jie Chuai, and Zhitang Chen. Efficient bayesian optimization with deep kernel learning and transformer pre-trained on multiple heterogeneous datasets. arXiv preprint arXiv:2308.04660, 2023. Hongyuan Mei, Chenghao Yang, and Jason Eisner. Transformer embeddings of irregularly spaced events and their participants. In International Conference on Learning Representations, 2022. URL https://openreview.net/forum?id=Rty5g9imm7H. Tomas Mikolov, Kai Chen, Greg Corrado, and Jeffrey Dean. Efficient estimation of word representations in vector space. arXiv preprint arXiv:1301.3781, 2013. 17

Deep Kernel Learning for Glaucoma

Sebastian W Ober, Carl E Rasmussen, and Mark van der Wilk. The promises and pitfalls of deep kernel learning. In Uncertainty in Artificial Intelligence, pages 1206–1216. PMLR, 2021. Davide Placido, Boya Yuan, Jonas X Hjaltelin, Chunlei Zheng, Amalie D Haue, Piotr J Chmura, Søren Brunak, et al. Vadesc-ehr: a transformer-based variational autoencoder for clustering longitudinal survival data from electronic health records. Nature Communications, 16(1):1–14, 2025. Sanjay Purushotham, Chuizheng Meng, Zhengping Che, and Yan Liu. Benchmarking deep learning models on large healthcare datasets. Journal of biomedical informatics, 83:112– 134, 2018. Jiajun Qiu, Yao Hu, Li Li, Abdullah Mesut Erzurumluoglu, Ingrid Braenne, Charles Whitehurst, Jochen Schmitz, Jatin Arora, Boris Alexander Bartholdy, Shrey Gandhi, et al. Deep representation learning for clustering longitudinal survival data from electronic health records. Nature Communications, 16(1):2534, 2025. Alvin Rajkomar, Eyal Oren, Kai Chen, Andrew M Dai, Nissan Hajaj, Michaela Hardt, Peter J Liu, Xiaobing Liu, Jake Marcus, Mimi Sun, et al. Scalable and accurate deep learning with electronic health records. NPJ digital medicine, 1(1):18, 2018. Laila Rasmy, Yang Xiang, Ziqian Xie, Cui Tao, and Degui Zhi. Med-bert: pretrained contextualized embeddings on large-scale structured electronic health records for disease prediction. NPJ digital medicine, 4(1):86, 2021. Peter Schulam and Raman Arora. Disease trajectory maps. Advances in neural information processing systems, 29, 2016. Peter Schulam, Fredrick Wigley, and Suchi Saria. Clustering longitudinal clinical marker trajectories from electronic health data: Applications to phenotyping and endotype discovery. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 29:1, pages 2956–2964, 2015. Benjamin Shickel, Patrick James Tighe, Azra Bihorac, and Parisa Rashidi. Deep ehr: a survey of recent advances in deep learning techniques for electronic health record (ehr) analysis. IEEE journal of biomedical and health informatics, 22(5):1589–1604, 2017. Satya Narayan Shukla and Benjamin M Marlin. Interpolation-prediction networks for irregularly sampled time series. arXiv preprint arXiv:1909.07782, 2019a. Satya Narayan Shukla and Benjamin M Marlin. Modeling irregularly sampled clinical time series. arXiv preprint arXiv:1906.04716, 2019b. Yuqi Si, Jingcheng Du, Zhao Li, Xiaoqian Jiang, Timothy Miller, Fei Wang, W Jim Zheng, and Kirk Roberts. Deep representation learning of patient data from electronic health records (ehr): a systematic review. Journal of biomedical informatics, 115:103671, 2021.

Deep Kernel Learning for Glaucoma

James M Tielsch, Alfred Sommer, Joanne Katz, Richard M Royall, Harry A Quigley, and Jonathan Javitt. Racial variations in the prevalence of primary open-angle glaucoma: the baltimore eye survey. Jama, 266(3):369–374, 1991. Cécile Trottet, Manuel Schürch, Ahmed Allam, Imon Barua, Liubov Petelytska, Oliver Distler, Anna-Maria Hoffmann-Vold, Michael Krauthammer, et al. Modeling complex disease trajectories using deep generative models with semi-supervised latent processes. arXiv preprint arXiv:2311.08149, 2023. Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. Advances in neural information processing systems, 30, 2017. Sally Yu-Yun Wang, Brian Tseng, and Tina Hernandez-Boussard. Deep learning approaches for predicting glaucoma progression using electronic health records and natural language processing. Ophthalmology Science, 2(2):100127, 2022. Sally Yu-Yun Wang, Kuldev Singh, Sophia C Lin, Tina Hernandez-Boussard, Robert Ritch, Leopold Rosen, and Joshua D Stein. Prediction models for glaucoma in a multicenter electronic health records consortium: The sight outcomes research collaborative. Ophthalmology Glaucoma, 7(2):162–172, 2024. Qingsong Wen, Tian Zhou, Chaoli Zhang, Weiqi Chen, Ziqing Ma, Junchi Yan, and Liang Sun. Transformers in time series: A survey. In International Joint Conference on Artificial Intelligence(IJCAI), 2023. Christopher KI Williams and Carl Edward Rasmussen. Gaussian processes for machine learning. MIT press Cambridge, MA, 2006. Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P Xing. Deep kernel learning. In Artificial intelligence and statistics, pages 370–378. PMLR, 2016. Xi Yang, Aokun Chen, Nima PourNejatian, Hoo Chang Shin, Kaleb E Smith, Christopher Parisien, Colin Compas, Cheryl Martin, Anthony B Costa, Mona G Flores, et al. A large language model for electronic health records. NPJ digital medicine, 5(1):194, 2022. Jithin Yohannan, Jiangxia Wang, Jamie Brown, Balwantray C Chauhan, Michael V Boland, David S Friedman, and Pradeep Y Ramulu. Evidence-based criteria for assessment of visual field reliability. Ophthalmology, 124(11):1612–1620, 2017. Joel Yohannan, Michael V Boland, Valentin Bolduc, Jacky M Sivak, Teodora Sajovic Chen, David S Friedman, Pradeep Y Ramulu, Carlos Gustavo Vde De Moraes, et al. Assessing glaucoma progression using machine learning trained on longitudinal visual field and clinical data. Ophthalmology, 128(7):1016–1026, 2021. Qiang Zhang, Aldo Lipani, Omer Kirnap, and Emine Yilmaz. Self-attentive hawkes process. In International conference on machine learning, pages 11183–11193. PMLR, 2020.

Deep Kernel Learning for Glaucoma

Xinlu Zhang, Shiyang Xu, Luoyao Zhang, Jingran Qiao, and Jing Chen. Improving medical predictions by irregular multimodal electronic health records modeling. arXiv preprint arXiv:2210.12156, 2022. Simiao Zuo, Haoming Jiang, Zichong Li, Tuo Zhao, and Hongyuan Zha. Transformer hawkes process. In International conference on machine learning, pages 11692–11702. PMLR, 2020.

Deep Kernel Learning for Glaucoma

Appendix A. Our feature extractor architectures consisted of recurrent neural networks (RNNs), gated recurrent units (GRUs), long short-term memory (LSTM), and transformer encoders. Our grid search initially started with transformer encoders for prototyping and we expanded to the other architectures. After we found that the GRU and transformer encoder feature extractors had the best performance, per validation mean squared error, we expanded the parameter configuration searches for those architectures to considerably more than the RNN and LSTM models. The total number of parameter configurations search was 1306 and the different configurations are given in the tables. Optimal parameters for each architecture are in bold: Table 3: Optimal Hyperparameters (bold) for All Feature Extractor Architectures Hyperparameter

RNN

GRU

LSTM

Transformer

CDE

ODE

GRUD

Num. Inducing Points Batch Size Model/Hidden Dim. Num. Heads Feedforward Dim. Num. Layers Decoder Dim. Latent Dim. Learning Rate

64, 128 32 16, 64, 256, 512 – – 1, 2, 3, 4, 6 64, 128, 256 2, 3, 4 5e-5, 1e-4, 2e-4

64, 128 32 16, 64, 256, 512 – – 1, 2, 3, 4, 6 64, 128, 256 2, 3, 4 1e-4, 2e-4, 3e-4, 5e-4

64, 128 32 16, 64, 256, 512 – – 1, 2, 3, 4 64, 128, 256 2, 3, 4 5e-5, 1e-4, 2e-4

64, 128 32, 64 128, 256, 512, 1024 4, 8, 16, 32, 64, 128 512, 1024, 2048 3, 4, 6, 8, 12 16, 32, 64, 128, 256 2, 3, 4, 5 1e-4, 2e-4, 5e-4

– 32 8 – – – 64 2 1e-4, 2e-4, 3e-4

– 32 256 – – – 32 3 1e-4, 2e-4

– 32 512 – – – 128 4 1e-4, 2e-4

Appendix B. We plot the raw and clustered trajectories with respect to the latent representations of a DKL-transformer.

Appendix C. Here, we provide additional cluster statistics for our experiments. We compared agglomerative ward, agglomerative average linkage, functional k-means, and GMM across k ∈ {2, 3, 4, 5} (Table 4). While agglomerative average linkage achieved marginally higher Silhouette scores at k=3, it produced severely imbalanced clusters (n=67,033, 508, 150), rendering the solution clinically uninterpretable. Agglomerative ward with k=3 achieved a Silhouette score of 0.743 and Davies-Bouldin score of 0.576, with balanced, clinically meaningful cluster sizes, and was selected as the final configuration.

Deep Kernel Learning for Glaucoma

Patient Trajectories (Z1 / Z2 / Mean)

Patient Trajectories (Z1 / Z2 / Variance)

0.9

2.0

0.8

1.5 Mean

0.6 0.5

1.0

0.4 0.3 0.2 0.1

0.5 0.0 1.0

1.0

0.8 0.0

0.2

0.8

0.6 0.4 Z2

0.4 Z1

0.6

0.0

0.2

0.2 0.8

1.0

Variance

0.7

0.6 0.4 Z2

0.4 Z1

0.6

0.0

(a) Means

0.2 0.8

1.0

0.0

(b) Variances

Figure 4: Posterior predictive mean (4(a)) and variance (4(b)) trajectories with respect to DKL transformer latent dimensions, renormalized on 3,821 patients. The Z-axis represents mean in logMAR units where higher indicates worse glaucoma. 3D Trajectories by Cluster (Mean)

3D Trajectories by Cluster (Variance)

Cluster 0 (n=821) Cluster 1 (n=1500 sampled) Cluster 2 (n=1500 sampled)

2.0

0.8

1.5 Mean

0.6

0.9

0.5

1.0

0.4 0.3 0.2 0.1

0.5 0.0 1.0

1.0

0.8 0.0

0.2

0.8

0.6 0.4 Z2

0.4 Z1

0.6

0.0

0.2 0.8

1.0

Variance

0.7

0.2

0.6 0.4 Z2

0.4 Z1

0.0

(a) Clustered Means

0.6

0.2 0.8

1.0

0.0

(b) Clustered Variances

Figure 5: Posterior predictive mean (5(a)) and variance (5(b)) trajectories with respect to DKL transformer latent dimensions, renormalized on 3,821 patients. The Z-axis represents mean in logMAR units where higher indicates worse glaucoma. Each trajectory represents a patient and three clusters represent three distinct patient trajectories, with purple indicating worst, blue indicating moderate, and yellow indicating stable.

Table 4: Clustering Method Comparison Across k ∈ {2, 3, 4, 5} K

Silhouette ↑

Davies-Bouldin ↓

Cluster Sizes

Agglomerative Ward Agglomerative Avg Functional K-Means GMM

2 2 2 2

0.8983 0.9047 0.8682 0.6077

0.3925 0.3582 0.5184 1.0027

[66870, 821] [67033, 658] [65882, 1809] [55437, 12254]

Agglomerative Ward Agglomerative Avg Functional K-Means GMM

3 3 3 3

0.7428 0.8863 0.7349 0.4880

0.5756 0.4565 0.5795 1.0280

[61843, 5027, 821] [67033, 508, 150] [60513, 6474, 704] [52738, 11438, 3515]

Agglomerative Ward Functional K-Means GMM

4 4 4

0.5366 0.6380 0.3682

0.6918 0.6172 1.3219

[48155, 13688, 5027, 821] [54369, 10994, 1882, 446] [46946, 11142, 7556, 2047]

Agglomerative Ward Functional K-Means GMM

5 5 5

0.5369 0.5353 0.2978

0.6456 0.6531 1.3383

[48155, 13688, 5027, 665, 156] [43890, 18026, 4349, 1103, 323] [37764, 16749, 8761, 3415, 1002]

Note: Silhouette scores range from −1 to 1; higher is better. Davies-Bouldin scores measure intra-cluster similarity relative to inter-cluster distance; lower is better. Agglomerative average linkage at k = 3 achieves a higher Silhouette score but produces highly imbalanced clusters (n=67,033, 508, 150), limiting clinical interpretability. Agglomerative ward with k = 3 (bold) was selected as the final configuration based on its strong cluster quality metrics and clinically meaningful, balanced cluster sizes. Agglomerative average (k = 2, 3) results excluded from k ≥ 4 as cluster sizes collapse to near-singleton groups.

Deep Kernel Learning for Glaucoma

Method

Deep Kernel Learning for Glaucoma

To assess robustness, we evaluated agglomerative ward stability across 100 random seed initializations, yielding high consistency (NMI: 0.950 ± 0.120, ARI: 0.957 ± 0.103; Table 5), confirming the presence of genuine cluster structure rather than algorithmic artifacts. Table 5: Clustering Stability Metrics for Transformer + DKL (k = 3) Metric

Mean

Std

NMI ARI

0.950 0.957

0.120 0.103

Note: NMI = Normalized Mutual Information; ARI = Adjusted Rand Index. Stability metrics computed across 100 clustering runs with different random initializations for k = 3 clusters.

Table 6 describes the final agglomerative ward (k = 3) cluster metrics including Silhouette and Calinski-Harabasz scores. Table 6: Clustering Quality Metrics and Overall Dataset Statistics Metric

Value

Total Patients Total Data Points Number of Clusters

67,691 402,552 3

Cluster Distribution: Cluster 0 Cluster 1 Cluster 2 Quality Metrics: Silhouette Score Calinski-Harabasz Score

821 patients 5027 patients 61843 patients 0.5373 23910.874

Appendix D. To assess model stability, we trained the ODE and CDE architectures for 25 epochs and all other architectures for 200 epochs using 10 random seeds (42, 123, 456, 789, 1024, 2048, 3141, 5926, 8765, and 4321), with each seed producing a different train-validation-test split. Table 8 reports the mean and standard deviation of test performance across seeds. Table 7 summarizes probabilistic performance for the DKL models using the continuous ranked probability score (CRPS), reported as the mean and standard deviation across seeds.

Deep Kernel Learning for Glaucoma

Table 7: CRPS performance across 10 random seeds for DKL models. Lower values indicate better probabilistic calibration and sharpness. Architecture DKL GRU DKL LSTM DKL RNN DKL Transformer

CRPS Mean

CRPS SD

0.0839 0.0936 0.0785 0.0122

0.0516 0.0093 0.0091 0.0291

Table 8: Experimental Variance and Performance. All models were tested across 10 random seeds. All models were trained to validation loss convergence. We observe that DKL models tended to outperform alternatives. Mean MSE

MSE SD

Mean MAE

MAE SD

Mean (%) ± 0.1

SD (%) ± 0.1

MLE GRU DKL GRU MLE LSTM DKL LSTM MLE RNN DKL RNN MLE Transformer DKL Transformer GRU-D GRU ODE RNN CDE

0.1106 0.1113 0.1156 0.1184 0.1166 0.1075 0.1068 0.1061 0.1246 0.1166 0.1246

0.0036 0.0168 0.0042 0.0027 0.0062 0.0025 0.0015 0.0026 0.0047 0.0029 0.0033

0.1708 0.1748 0.1718 0.1712 0.1795 0.1721 0.171 0.1703 0.1765 0.1888 0.2017

0.0023 0.0197 0.0037 0.0018 0.0071 0.0021 0.0015 0.0026 0.0033 0.0028 0.0055

56.41 53.7741 57.9248 59.6589 53.8397 56.9478 53.3605 53.0564 56.8392 43.8709 36.2228

1.6192 7.9965 1.3425 0.4508 8.3422 1.9967 1.2419 1.3695 0.8086 2.1627 2.2172

Deep Kernel Learning for Glaucoma

Architecture

Deep Kernel Learning for Glaucoma

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