📊 statistics

A scalable Bayesian functional factor model for high-dimensional longitudinal molecular data

Motivated by a longitudinal COVID-19 study, this paper introduces a scalable Bayesian functional factor model that integrates latent factor modeling with functional principal component analysis to jointly uncover coordinated biomarker trajectories and patient heterogeneity in high-dimensional molecular data, utilizing sparsity-inducing priors and an annealed variational algorithm for efficient inference.

Original authors: Salima Jaoua, Daniel Temko, Hélène Ruffieux

Published 2026-03-24

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Salima Jaoua, Daniel Temko, Hélène Ruffieux

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to understand a massive, chaotic orchestra playing a complex symphony. You have 20,000 musicians (variables) playing different instruments (genes, proteins, metabolites), and you are recording them over time (longitudinal data).

The problem is:

Too many musicians: You can't listen to every single instrument individually; it's too much noise.
Messy recordings: Some musicians are recorded every second, others only once a day. The data is "sparse" and irregular.
Hidden patterns: You suspect the musicians aren't playing randomly. They are following a few hidden "conductors" (biological pathways) that tell groups of instruments when to play loud or soft. But you don't know who the conductors are, which musicians they lead, or how the music changes over time.

This paper introduces a new statistical tool called bayesSYNC (a "Bayesian Functional Factor Model") to solve this exact problem. Here is how it works, explained simply:

1. The Core Idea: Finding the "Conductors"

Instead of trying to analyze 20,000 separate songs, the model looks for latent factors. Think of these as the conductors of the orchestra.

Conductor 1 might be the "Inflammation Conductor." They tell the drums (CRP), the brass (cytokines), and the strings (complement proteins) to play a loud, aggressive song.
Conductor 2 might be the "Recovery Conductor." They tell the woodwinds (amino acids) to play a slow, calming melody.

The model's job is to figure out:

Who are the conductors?
Which musicians follow which conductor? (This is the sparsity part: Conductor 1 doesn't talk to every musician, just a specific group).
How does the music change as the concert goes on?

2. The "Functional" Part: Listening to the Melody, Not Just the Notes

Older methods might just look at the average volume of a musician over the whole concert. But in biology, timing is everything. A spike in inflammation at hour 1 is very different from a spike at hour 100.

This model uses Functional Principal Component Analysis (FPCA).

Analogy: Imagine the conductor's baton movement. It's not just a single dot; it's a smooth, flowing curve.
The model breaks this curve down into basic shapes (like a sine wave, a spike, or a slow ramp).
It says: "Conductor 1's movement is mostly made of Shape A (a sharp rise) and Shape B (a slow fall)."
This allows the model to handle messy data where some musicians were recorded at weird times, because it understands the shape of the melody, not just the specific notes caught on tape.

3. The "Bayesian" Part: The Smart Guessing Game

The model is "Bayesian," which means it starts with a set of reasonable guesses and updates them as it hears more music.

The "Spike-and-Slab" Trick: This is a clever way to decide which musicians belong to which conductor.
- Imagine a "Slab" (a flat surface) where a musician could belong to a conductor.
- And a "Spike" (a sharp needle) where a musician definitely does not belong.
- The model tries to push most musicians onto the "Spike" (zero connection) and only keeps the important ones on the "Slab." This keeps the model simple and interpretable, so you don't get a confusing list of 20,000 connections.

4. The "Annealed" Part: Avoiding the Wrong Tune

Finding the best pattern in such a huge mess is like trying to find the highest peak in a foggy mountain range. If you just walk uphill, you might get stuck on a small hill (a local optimum) and think it's the highest point.

The authors use Simulated Annealing.

Analogy: Imagine you are a hiker trying to find the highest peak.
- Standard method: You walk up the nearest hill and stop. You might miss the real mountain.
- Annealing method: You start by shaking the ground violently (high "temperature"). This allows you to jump over small hills and valleys. As you get closer to the top, you slowly calm the shaking down (lowering the temperature) and settle into the true highest peak.
This ensures the model doesn't get stuck on a "good enough" answer but finds the best answer.

5. Why This Matters: The COVID-19 Example

The authors tested this on real data from COVID-19 patients.

The Result: The model found two main "conductors."
1. The Inflammation Conductor: This group of proteins and metabolites went wild in severe patients. The model showed that patients with high scores on this "conductor" had a harder time recovering.
2. The Recovery Conductor: This group helped patients return to normal.
The Benefit: Instead of looking at 57 different blood markers and getting confused, doctors can now look at just two scores per patient. These scores tell them exactly how the patient's immune system is behaving over time, helping to predict who will get better and who might need more help.

Summary

This paper gives scientists a super-powered microscope for time-traveling through biological data.

It handles huge datasets (thousands of variables).
It handles messy data (irregular time points).
It finds hidden patterns (conductors) that explain how groups of biological markers work together.
It tells us who is different (which patients have unique recovery paths).

It's like turning a chaotic, noisy recording of an orchestra into a clear sheet music score, showing exactly who the conductors are and how the music evolves, even if you only have a few scattered recordings of the performance.

1. Problem Statement

The paper addresses the challenge of analyzing high-dimensional longitudinal molecular data (e.g., genes, proteins, metabolites) collected over time in biomedical research. Current statistical tools face a gap in handling these data jointly:

Longitudinal Modeling: Traditional methods often treat variables independently or produce dense, uninterpretable multivariate components when the number of variables ( $p$ ) is large (tens of thousands).
High-Dimensional Analysis: Standard factor models reduce dimensionality but often ignore the functional (time-varying) nature of the data.
Scalability: Existing Gaussian Process (GP) based latent factor models (e.g., MEFISTO) struggle with scalability as the number of subjects ( $N$ ) increases due to the computational cost of inverting large covariance matrices.
Specific Needs: Researchers need a framework that identifies sparse, interpretable biological pathways (latent processes) driving coordinated changes across subsets of biomarkers while capturing subject-specific temporal heterogeneity (how individuals differ in their recovery or disease progression).

2. Methodology: Bayesian Functional Factor Model

The authors propose a unified framework combining Latent Factor Modeling with Functional Principal Component Analysis (FPCA) within a Bayesian hierarchy.

A. Model Structure

Latent Factor Representation:
The observed $p$ -dimensional functional data $y_i(t)$ for subject $i$ is decomposed as:
$y_i(t) = \mu(t) + B h_i(t) + \epsilon_i(t)$
- $\mu(t)$ : Mean function.
- $B$ : A $p \times Q$ loading matrix representing the contribution of $Q$ latent factors to the $p$ variables.
- $h_i(t)$ : Latent factor curves for subject $i$ .
- Sparsity: The loading matrix $B$ is induced to be sparse using a spike-and-slab prior, ensuring that each latent factor influences only a specific subset of biomarkers (interpretable pathways).
Functional Representation of Factors (FPCA):
Each latent factor trajectory $h_i^{(q)}(t)$ is modeled using FPCA:
$h_i^{(q)}(t) = \sum_{\ell=1}^L \zeta_{i\ell}^{(q)} \psi_\ell^{(q)}(t)$
- $\psi_\ell^{(q)}(t)$ : Orthogonal eigenfunctions (temporal modes) estimated via splines.
- $\zeta_{i\ell}^{(q)}$ : Subject-specific scores representing how strongly an individual expresses a specific temporal mode.
- This allows the model to capture shared temporal programs (via eigenfunctions) and individual variation (via scores).
Observation Model:
Data are observed on irregular, sparse grids. The model uses a Gaussian observation model with variable-specific noise variance, allowing for heteroscedasticity.

B. Bayesian Specification

Priors:
- Spike-and-Slab: Used on $B$ to enforce sparsity. The inclusion probability $\omega_q$ follows a Beta prior, allowing the data to determine the effective number of non-zero loadings.
- Gaussian Priors: Applied to spline coefficients and scores.
- Half-Cauchy Priors: Used for variance parameters to ensure robustness.
Identifiability: The model is identifiable up to orthogonal transformations. A post-processing step (eigendecomposition of empirical covariance) is applied to the posterior mean estimates to restore the FPCA interpretation (orthogonal modes ordered by variance).

C. Inference Algorithm: Annealed Variational Inference

Given the high dimensionality and the multimodal nature of the spike-and-slab posterior, standard MCMC is computationally prohibitive. The authors develop a Variational Inference (VI) approach:

Mean-Field Approximation: Assumes independence between parameter blocks (mean functions, eigenfunctions, scores, loadings).
Closed-Form Updates: Due to conditional conjugacy (aided by auxiliary variables for Half-Cauchy priors), all variational updates have closed-form solutions, avoiding numerical optimization at each step.
Simulated Annealing: To escape local optima caused by the multimodal spike-and-slab posterior, the algorithm employs a temperature schedule ( $T \ge 1$ ). The entropy term is inflated at high temperatures to flatten the objective surface, then gradually cooled to $T=1$ . This significantly improves factor selection and convergence.
Scalability: The algorithm scales linearly with $N$ and near-linearly with $Q$ , avoiding the $O(N^3)$ matrix inversions typical of GP methods.

3. Key Contributions

Unified Framework: First scalable Bayesian framework that jointly models high-dimensional variables and temporal dynamics using sparse latent factors and FPCA.
Interpretability: The spike-and-slab prior yields sparse loading matrices, directly identifying which biomarkers belong to which biological pathway.
Scalability: The annealed variational inference algorithm enables analysis of datasets with $p \approx 20,000$ variables and irregular time grids, outperforming Gaussian Process-based alternatives in speed and memory usage.
Robust Inference: The introduction of annealing is shown to be critical for navigating the complex posterior landscape of sparse factor models, leading to better factor selection and higher Evidence Lower Bounds (ELBO) compared to standard VI.
Software: Implementation in the R package bayesSYNC, supporting parallel execution.

4. Results

A. Simulation Studies

Large-Scale Recovery: In simulations with $N=100$ and $p=20,000$ , the model accurately recovered the true sparse loading structure and latent temporal trajectories.
Annealing Benefits: Comparing "Vanilla" VI vs. "Annealed" VI showed that annealing:
- Correctly identified the true number of factors (discarding irrelevant ones).
- Reduced estimation error for mean functions and eigenfunctions.
- Achieved higher ELBO values and faster convergence.
Comparison with MEFISTO:
- Accuracy: bayesSYNC outperformed MEFISTO (both standard and subject-specific versions) in recovering loadings (higher AUC) and reconstructing trajectories (lower Integrated Squared Error).
- Scalability: MEFISTO failed to complete runs for $N=40$ within 36 hours due to memory/computation limits, whereas bayesSYNC scaled efficiently.
- Uncertainty: bayesSYNC provided well-calibrated 95% prediction intervals, whereas MEFISTO-base (shared trajectory) produced overly wide intervals to account for unmodeled subject heterogeneity.

B. Application to Longitudinal COVID-19 Data

Applied to a dataset of 57 biomarkers from 155 patients across 5 severity strata:

Factor 1 (Inflammation/Kynurenine Pathway): Identified a coordinated upregulation of inflammatory markers (CRP, cytokines) and the kynurenine pathway (kynurenine, quinolinic acid), with downregulation of tryptophan and serotonin. This aligns with known immune activation mechanisms in severe COVID-19.
Factor 2 (Amino Acid Metabolism): Dominated by branched-chain amino acids, reflecting metabolic depletion in severe disease.
Clinical Interpretation:
- Severity: The first eigenfunction of Factor 1 acted as a "severity axis." Patients with high positive scores had sustained inflammatory elevation.
- Recovery: The second eigenfunction captured recovery dynamics. Patients with favorable scores returned to baseline faster.
- Heterogeneity: Subject-level scores successfully stratified patients by clinical severity and recovery speed, revealing distinct biological trajectories even within the same clinical strata.

5. Significance

This work provides a critical methodological advance for systems biology and precision medicine. By enabling the analysis of tens of thousands of molecular markers over time, it allows researchers to:

Uncover coordinated biological programs driving disease progression.
Quantify patient heterogeneity in a low-dimensional, interpretable space.
Move beyond static snapshots to understand dynamic disease mechanisms.
The scalability of the approach makes it applicable to modern large-scale cohort studies (e.g., cancer progression, autoimmune diseases, developmental biology) where previous methods were computationally intractable. The open-source implementation (bayesSYNC) ensures reproducibility and accessibility for the broader scientific community.