Making Multi-Axis Models Robust to Multiplicative Noise: How, and Why?

1. Problématique et Contexte

Le problème :
Les modèles multi-axes (ou structurés par somme de Kronecker) visent à apprendre simultanément des réseaux de dépendance (graphes) pour les lignes et les colonnes d'une matrice de données (par exemple, les cellules et les gènes en séquençage ARN à cellule unique ou scRNA-seq). Ces modèles décomposent la matrice de précision $\Omega$ en une somme de Kronecker $\Psi_{cols} \oplus \Psi_{rows}$, réduisant ainsi le nombre de paramètres de $O(d_{rows}^2 d_{cols}^2)$ à $O(d_{rows}^2 + d_{cols}^2)$.

Cependant, les méthodes existantes (comme GmGM) supposent généralement que les données suivent une distribution gaussienne ou une copule gaussienne. Or, dans des domaines comme le scRNA-seq, les données sont corrompues par un bruit multiplicatif naturel. Ce bruit provient de biais techniques (biais GC, incomplétude de la capture des gènes) où la valeur observée est le produit d'une valeur latente et de facteurs de bruit strictement positifs ($r_{rows} \times r_{cols}$).

Les limites des approches actuelles :

1. Incapacité à modéliser le bruit : Aucune méthode n'existait pour ajuster des modèles multi-axes sous un régime de bruit multiplicatif.
2. Hypothèses restrictives : Les copules gaussiennes supposent une dépendance asymptotiquement nulle pour les événements extrêmes (pas de dépendance de queue), ce qui est souvent faux dans les données biologiques réelles.
3. Manque de flexibilité : Ignorer ce bruit réduit la capacité expressive du modèle et fausse l'estimation des réseaux de dépendance conditionnelle.

2. Méthodologie : MED-MAGMA

Les auteurs proposent un nouvel algorithme, MED-MAGMA (Matrix Elliptical Distribution's Multi-Axis Graphical Modelling Algorithm), basé sur l'estimation du maximum de vraisemblance (MLE) via l'algorithme Expectation-Maximization (EM).

A. Modèle Probabiliste

Le modèle suppose que les données observées $X$ sont le résultat d'un bruit multiplicatif appliqué à des données latentes $Z$ suivant une distribution normale multivariée structurée :
$$X_{ij} = r_{rows, i} \cdot r_{cols, j} \cdot Z_{ij}$$
$$\text{vec}(Z) \sim \mathcal{N}(0, (\Psi_{cols} \oplus \Psi_{rows})^{-1})$$
où $r$ suit une distribution strictement positive arbitraire. Ce modèle généralise les distributions elliptiques au cas matriciel.

B. Étape de "Dénouage" (Denoising) - Transformation $f$

Comme la distribution exacte du bruit $r$ est inconnue, la vraisemblance directe est intractable. Les auteurs définissent une fonction de transformation $f(X)$ qui élimine l'information corrompue par le bruit tout en préservant la structure de dépendance.

• Fonction $f$ : Elle normalise les données en utilisant des moyennes géométriques sur les lignes, les colonnes et l'ensemble de la matrice.
  $$f(X)_{ij} = \frac{x_{ij}}{(\prod x_{i'j'})^{\frac{1}{d_{rows}d_{cols}}} (\prod x_{ij'})^{\frac{1}{d_{cols}}} (\prod x_{i'j})^{\frac{1}{d_{rows}}}}$$
• Propriété clé : $f(X) = f(Y)$ si et seulement si $X$ et $Y$ diffèrent uniquement par un facteur multiplicatif de type $r_{rows} \otimes r_{cols}$. Cela projette les données dans un "espace quotient" où le bruit est éliminé.
• Gestion des zéros : Une version modifiée de $f$ est proposée pour gérer les zéros (fréquents en scRNA-seq) en ignorant les termes nuls dans les produits.

C. Algorithme EM et Approximation de Laplace

L'algorithme itère entre deux étapes pour estimer $\Psi_{rows}$ et $\Psi_{cols}$ :

1. E-Step (Estimation de la statistique suffisante) :

   • L'objectif est de calculer l'espérance de la statistique suffisante ($Z Z^T$ et $Z^T Z$) conditionnée aux données observées transformées $y = f(x)$.
   • L'intégrale sur l'espace des fibres (l'ensemble des $Z$ tels que $f(Z)=y$) est intractable. Les auteurs utilisent l'approximation de Laplace pour approximer cette intégrale autour du point $z^*$ qui minimise l'énergie quadratique $z^T \Omega z$ sous la contrainte $f(z)=y$.
   • La recherche de $z^*$ est formulée comme un problème d'optimisation quadratique sous contrainte, résolu par un algorithme "flip-flop" (alternance entre les variables de ligne et de colonne).

2. M-Step (Maximisation) :

   • Une fois les statistiques suffisantes approchées ($S_{rows}$ et $S_{cols}$), l'étape de maximisation revient exactement au problème résolu par l'algorithme GmGM (GmGM: a fast multi-axis Gaussian graphical model).
   • Cela permet de réutiliser l'efficacité de GmGM pour mettre à jour les graphes $\Psi$.

3. Contributions Clés

1. Premier algorithme robuste au bruit multiplicatif : Introduction de MED-MAGMA, capable d'estimer des graphes conditionnels dans des matrices corrompues par un bruit multiplicatif, un scénario omniprésent en biologie.
2. Généralisation des distributions elliptiques : Le modèle ne suppose pas une distribution gaussienne stricte ni une copule gaussienne, permettant de capturer des dépendances de queue (tail dependence) et des modes de dépendance plus complexes.
3. Traitement théorique de l'identifiabilité : Résolution des problèmes d'identifiabilité inhérents aux modèles elliptiques et aux produits de Kronecker via des contraintes de moyenne géométrique.
4. Implémentation logicielle : Mise à disposition du package Python MED-MAGMA.

4. Résultats Expérimentaux

Les auteurs ont évalué MED-MAGMA de trois manières :

• Validation sur données synthétiques :

  • Sur des matrices $100 \times 150$ et $200 \times 200$ avec des niveaux de bruit multiplicatif variables.
  • Résultat : MED-MAGMA maintient des performances stables (AUPR, assortativité, AMI) quelle que soit la force du bruit, tandis que GmGM et les estimateurs de base (MLE simple) se dégradent sévèrement.

• Validation systématique sur données réelles (scRNA-seq) :

  • Test sur 30 datasets publics du "Single Cell Expression Atlas" (moins de 1000 cellules).
  • Résultat : MED-MAGMA surpasse systématiquement GmGM (avec et sans "nonparanormal skeptic") et les méthodes uniaxes en termes d'assortativité (cohérence des clusters de cellules) et d'AMI (accord avec les types cellulaires connus).

• Étude de cas (E-MTAB-7249) :

  • Analyse d'un dataset de carcinome ovarien.
  • Les modules de gènes appris correspondent biologiquement aux sous-types cellulaires (stromal vs néoplasique).
  • Tests de robustesse avec l'outil Robin : Les réseaux appris par MED-MAGMA montrent une grande stabilité face à la perturbation des arêtes (facteurs de Bayes élevés), confirmant la robustesse des structures découvertes.

5. Signification et Limites

Signification :
Ce travail démontre que la prise en compte explicite du bruit multiplicatif, souvent ignoré, est cruciale pour l'inférence de réseaux en biologie. En passant d'une hypothèse gaussienne à une distribution elliptique généralisée, MED-MAGMA offre une flexibilité accrue pour modéliser les dépendances extrêmes et les biais techniques inhérents aux technologies de séquençage, produisant des réseaux biologiques plus fiables.

Limites :

• Évolutivité : La méthode gère actuellement des matrices de quelques milliers de lignes/colonnes mais pourrait rencontrer des difficultés avec des dizaines de milliers (coût cubique de l'algorithme).
• Hypothèses de dépendance : Bien que plus flexible que le gaussien, le modèle suppose toujours des dépendances linéaires et radialement symétriques.
• Approximation : L'utilisation de l'approximation de Laplace pour l'intégrale sur les fibres est une approximation, bien que jugée efficace en pratique.

En conclusion, MED-MAGMA représente une avancée méthodologique majeure pour l'analyse de données matricielles complexes et bruitées, en particulier dans le domaine de la génomique à cellule unique.