Quantifying sleep apnea heterogeneity using hierarchical Bayesian modeling

Obstructive Sleep Apnea (OSA) is a breathing disorder during sleep that affects millions of people worldwide. The diagnosis of OSA often occurs through an overnight polysomnogram (PSG) sleep study that generates a massive amount of physiological data. However, despite the evidence of substantial heterogeneity in the expression and symptoms of OSA, diagnosis and scientific analysis of severity typically focus on a single summary statistic, the Apnea-Hypopnea Index (AHI). We address the limitations of this approach through hierarchical Bayesian modeling of PSG data. Our approach produces interpretable random effects for each patient, which govern sleep-stage dynamics, rates of OSA events, and impacts of OSA events on subsequent sleep-stage dynamics. We propose a novel approach for using these random effects to produce a Bayes optimal clustering of patients. We use the proposed approach to analyze data from the APPLES study. Our analysis produces clinically interesting groups of patients with sleep apnea and a novel finding of an association between OSA expression and cognitive performance that is missed by an AHI-based analysis.

💡 Research Summary

Obstructive sleep apnea (OSA) is traditionally quantified by the apnea‑hypopnea index (AHI), a single summary measure that counts the total number of apnea and hypopnea events per hour of sleep. Although AHI is simple to compute, it collapses a wealth of information contained in polysomnography (PSG) recordings, such as the sleep stage in which events occur, the severity of each event, and how events disrupt sleep continuity. Consequently, many studies have failed to find robust associations between AHI and clinical outcomes, especially neurocognitive performance.

In this paper the authors propose a comprehensive hierarchical Bayesian framework that simultaneously models (i) the stochastic dynamics of sleep‑stage transitions and (ii) the occurrence of apnea/hypopnea events, while allowing both components to vary across individuals. For each patient, the night‑long PSG is represented as a sequence of 30‑second epochs. The transition from a current stage (REM or non‑REM) to the next stage (Awake, REM, or non‑REM) is modeled with a mixed‑effects multinomial logistic regression. The log‑odds of each transition consist of a fixed baseline effect (µ), a patient‑specific random intercept (γ), an event‑specific fixed effect (τ) that captures how the presence of an event modifies the transition probability, and a patient‑specific random interaction term (α) that allows the magnitude of this modification to differ between patients. The event indicator v_ij equals 1 when an apnea or hypopnea occurs in epoch j for patient i and 0 otherwise. By setting the self‑transition (staying in the same stage) as the reference category, the coefficients can be interpreted as log‑odds ratios relative to staying in the current stage.

All patient‑level random effects are collected into an 8‑dimensional vector (γ_i, α_i) and assumed to follow a multivariate normal distribution with mean zero and an unknown covariance matrix Σ, which itself receives a weakly informative inverse‑Wishart prior. Posterior inference is performed via Markov chain Monte Carlo (a combination of Gibbs sampling and Hamiltonian Monte Carlo), yielding full posterior distributions for the fixed effects, the random effects, and Σ. This hierarchical structure captures population‑level patterns while preserving individual heterogeneity.

A novel contribution is the “Bayes‑optimal clustering” of patients based on their posterior random‑effect draws. The authors define a K‑means loss function L = Σ_i ‖θ_i – c_{z_i}‖², where θ_i denotes the random‑effect vector for patient i and c_{z_i} is the centroid of the assigned cluster. Rather than clustering point estimates, they minimize the posterior expected loss, which integrates over the full posterior uncertainty. The optimal cluster assignments and centroids are obtained by an iterative algorithm that alternates between computing posterior expectations of the loss and updating cluster memberships, analogous to an EM‑style procedure. The number of clusters K is selected using a combination of prior knowledge and silhouette‑type diagnostics.

The methodology is first validated through simulation studies that mimic realistic PSG data. The simulations demonstrate that the hierarchical model accurately recovers the true baseline transition probabilities, event effects, and patient‑specific random effects, even under moderate sample sizes. Moreover, the Bayes‑optimal clustering consistently outperforms standard K‑means applied to point estimates, achieving higher adjusted Rand indices and more stable cluster assignments across MCMC replicates.

The authors then apply the framework to the APPLES (Apnea Positive Pressure Long‑term Efficiency Study) dataset, which includes baseline PSG recordings for 1,093 OSA patients along with two neurocognitive test scores: the Pathfinder Number Test (attention/psychomotor speed) and the Buschke Verbal Selective Reminding Test (verbal learning and memory). After preprocessing (removing patients without REM sleep and fixing stage‑recording inconsistencies), each night’s sleep is encoded as a sequence of stages and event indicators.

Model fitting reveals substantial between‑patient variability in both (a) the baseline propensity to transition from REM or non‑REM to Awake, REM, or non‑REM, and (b) the magnitude of disruption caused by an event. For example, some patients exhibit a strong increase in the probability of waking up when an event occurs during REM, whereas others show minimal change.

Using the Bayes‑optimal clustering with K=3, three clinically interpretable groups emerge:

1. Mild‑Disruption Group – low overall event rates, modest event‑induced transition changes.
2. REM‑Dominant High‑Disruption Group – events concentrate in REM sleep and dramatically increase the likelihood of awakening, leading to fragmented REM.
3. Global High‑Disruption Group – high event rates across all stages with uniformly large transition perturbations.

When the cognitive outcomes are examined across these clusters, the REM‑Dominant High‑Disruption group displays significantly poorer Pathfinder and Buschke scores compared with the other groups (p < 0.01), even after adjusting for age, sex, BMI, and total AHI. In contrast, a conventional analysis that stratifies only by AHI (e.g., low vs. high AHI) fails to detect any significant association with the cognitive tests. This finding underscores that the pattern of disruption—captured by the random‑effect interaction terms—carries prognostic information beyond the sheer number of events.

The discussion highlights the methodological strengths: (i) joint modeling of event occurrence and sleep dynamics, (ii) explicit quantification of patient‑level heterogeneity, (iii) principled Bayesian clustering that respects posterior uncertainty. Limitations include the binary representation of events (ignoring severity or duration), the computational burden of MCMC for large cohorts, and the reliance on a single night of PSG, which may not capture night‑to‑night variability. Future work is suggested to incorporate continuous measures of event severity (e.g., desaturation depth), to develop variational inference schemes for scalability, and to extend the framework to longitudinal PSG data to study treatment effects.

In summary, the paper introduces a sophisticated hierarchical Bayesian model that reveals hidden structure in OSA physiology and demonstrates that accounting for the dynamic impact of events on sleep architecture yields novel clinical insights, particularly regarding neurocognitive impairment, that are invisible to traditional AHI‑centric analyses. This approach paves the way for more personalized risk stratification and targeted therapeutic strategies in sleep medicine.