Finite mixtures of multivariate Poisson-log normal factor analyzers for clustering count data¶
作者: Andrea Payne, Anjali Silva, Steven J Rothstein, Paul D. McNicholas, Sanjeena Subedi
来源: Statistics and Computing
主题: 统计计算 / 算法
相关性: 6/10
链接: 期刊页 · arXiv
一、领域脉络与小综述¶
这个方向是什么¶
这个子方向解决的根本问题是:如何对高维离散计数数据(尤其是来自 RNA-seq 的基因表达计数)进行聚类。核心挑战在于,计数数据具有离散性、过离散性(方差 > 均值)和变量间相关性,而经典的聚类方法(如 K-means、高斯混合模型)要么无法处理离散数据,要么无法建模相关性。当前成熟度属于方法学应用的中期——基础模型(泊松-对数正态混合)已建立,但计算效率和模型简约性仍是瓶颈。
发展脉络(history)¶
奠基工作:泊松-对数正态(PLN)分布本身是一个经典层次模型(Aitchison & Ho, 1989),但将其用于聚类是较晚的事。Silva et al. (2019) [1] 首次提出多元泊松-对数正态混合模型(MPLN)用于 RNA-seq 数据聚类,参数估计采用 MCMC-EM 算法(用 NUTS 采样器从后验采样)。这是该子方向的奠基性工作,但计算代价极高——每次 E 步都需要运行 MCMC。
主要进展:Subedi & Browne (2020) [2, 3] 做了两个关键推进:(a) 用变分 EM(variational EM)替代 MCMC-EM,大幅降低计算成本;(b) 引入简约化(parsimonious)协方差结构族(基于特征值分解),将模型参数从 \(O(p^2)\) 降到 \(O(p)\) 量级。这两篇论文奠定了当前工作的直接基础。
当前 frontier:尽管变分 EM 已比 MCMC-EM 快得多,但 MPLN 模型的协方差矩阵仍是满秩的 \(p \times p\) 矩阵,当变量数 \(p\) 很大(RNA-seq 数据常有数千个基因)时,参数数量爆炸,且变分下界的计算和优化都变得困难。因子分析约束是解决高维协方差建模的标准思路(在 Gaussian 混合中已有成熟应用,如 McNicholas & Murphy 2008, 2010),但尚未被引入 MPLN 框架。
本文的位置:本文是 Subedi & Browne (2020) 的直接后继——将因子分析(factor analysis)约束引入 MPLN 混合模型,得到MPLN 因子分析器(MPLNFA)族,进一步压缩协方差参数。这是该子方向在模型简约性上的又一次推进。
子线索聚类¶
这些被引文献大致落在三条子线索上:
-
MPLN 混合模型的计算方法演进(Silva et al. 2019 → Subedi & Browne 2020 → 本文):从 MCMC-EM 到变分 EM,再到因子分析约束。这条线索的核心问题是如何在保持模型灵活性的同时降低计算和参数代价。
-
混合模型中的协方差简约化(McNicholas 2016 [14] 综述;Browne & McNicholas 2015 [16]; Vrbik & McNicholas 2014 [20]; Dang et al. 2015 [23] 等):这是一个更广泛的子领域,研究如何通过特征值分解、因子分析、或其它约束来减少混合模型中协方差矩阵的自由参数。本文直接借用其中的因子分析器(factor analyzer)思路。
-
RNA-seq 数据的统计建模与归一化(Robinson & Smyth 2007 [11]; Robinson et al. 2009 [4]; Robinson & Oshlack 2010 [6]; Dillies et al. 2013 [8]):这条线索关注的是计数数据的分布假设(负二项 vs. 泊松-对数正态)和文库大小归一化。本文的 MPLN 模型属于这条线索的"分布选择"分支,但归一化步骤(TMM、中位数比等)是预处理,不是模型本身的一部分。
这个方向在追问的核心问题¶
- 如何在高维计数数据中同时建模过离散和相关性? 泊松-对数正态层次结构解决了过离散,但协方差矩阵的维数灾难是主要瓶颈。
- 如何在不牺牲太多模型拟合的前提下大幅减少参数? 因子分析约束是一种答案,但需要验证它是否足够灵活以捕获真实数据中的相关结构。
- 变分近似的精度如何? 变分 EM 比 MCMC-EM 快,但变分下界可能低估真实边际似然,且后验近似可能不够准确——这对聚类结果的影响需要实证评估。
- 模型选择(选多少个因子?多少个成分?哪种协方差结构?)如何可靠地进行? 信息准则(BIC)是标准做法,但在高维、小样本场景下可能表现不佳。
⚠️ 作者的 framing¶
这是作者的说法:作者把缺口 frame 成"MPLN 模型在高维场景下参数过多,需要因子分析约束来降维"。具体来说,Subedi & Browne (2020) 的简约化族虽然减少了参数,但协方差矩阵仍是满秩的(只是通过特征值分解约束了形状/方向/体积),当 \(p\) 很大时仍然不可行。本文的 MPLNFA 将协方差矩阵参数化到 \(p \times q\) 的载荷矩阵(\(q \ll p\)),从而真正实现降维。
被淡化或回避的竞争路线: - 负二项混合模型:RNA-seq 领域最常用的分布是负二项(如 edgeR [4]),但作者没有与负二项混合模型做直接比较。负二项也能处理过离散,且计算更简单(不需要潜高斯变量),但建模相关性的能力较弱(通常假设条件独立或使用 copula)。 - 非参数/半参数聚类方法:如谱聚类、DBSCAN 等,这些方法不依赖分布假设,但作者没有讨论。 - 深度生成模型:如 VAE 用于计数数据聚类,作者完全没有提及。
什么明显该被引/该存在、却没出现在 intro 里?: - 没有引用任何关于变分推断的渐近性质的工作(如 Blei et al. 2017 的变分推断综述,或更近的"variational Bayes: a review")。考虑到本文的核心估计方法是变分 EM,这是一个明显的缺失——读者无法知道变分近似的误差界或一致性条件。 - 没有引用高维因子分析的经典理论(如 Bai & Li 2012, Fan et al. 2013 等),尽管本文的因子分析约束本质上是在做高维协方差建模。 - 没有引用混合因子分析器(Mixture of Factor Analyzers)的原始论文(Ghahramani et al. 1996; Tipping & Bishop 1999),虽然文中提到了它们。
张力¶
未见明显对立引用。所有被引工作基本沿着"MPLN 模型 → 变分加速 → 协方差简约化"这条主线推进,没有出现彼此矛盾或在不同条件下得相反结论的情况。不过,负二项 vs. 泊松-对数正态这两个分布家族在 RNA-seq 建模中是有张力的——负二项更常用(edgeR, DESeq2),但泊松-对数正态能更自然地建模相关性。本文没有直接比较这两种分布假设下的聚类性能。
二、最核心、最简单的例子 / 数学问题¶
第一步:把符号、模型、可观测数据交代清楚¶
符号: - \(n\):样本量(观测数),如 RNA-seq 中的样本数。 - \(p\):变量数(维度),如 RNA-seq 中的基因数。 - \(G\):混合成分数(聚类数)。 - \(q\):潜因子数(因子分析中的公共因子数),\(q \ll p\)。 - \(\mathbf{Y}_i = (Y_{i1}, \ldots, Y_{ip})^\top\):第 \(i\) 个观测的 \(p\) 维计数向量(可观测)。 - \(\mathbf{X}_i = (X_{i1}, \ldots, X_{ip})^\top\):第 \(i\) 个观测的 \(p\) 维潜高斯向量(不可观测,潜在变量)。 - \(\mathbf{U}_i = (U_{i1}, \ldots, U_{iq})^\top\):第 \(i\) 个观测的 \(q\) 维潜公共因子(不可观测,因子分析中的因子)。 - \(Z_i \in \{1, \ldots, G\}\):第 \(i\) 个观测的潜成分标签(不可观测,聚类目标)。 - \(\boldsymbol{\mu}_g\):第 \(g\) 个成分的 \(p\) 维均值向量(参数)。 - \(\boldsymbol{\Lambda}_g\):第 \(g\) 个成分的 \(p \times q\) 载荷矩阵(参数)。 - \(\boldsymbol{\Psi}_g\):第 \(g\) 个成分的 \(p \times p\) 对角唯一性矩阵(参数),\(\boldsymbol{\Psi}_g = \text{diag}(\psi_{g1}, \ldots, \psi_{gp})\)。 - \(\pi_g\):第 \(g\) 个成分的混合权重,\(\sum_{g=1}^G \pi_g = 1\)。
模型(MPLNFA 的数据生成机制): 1. 潜成分标签:\(Z_i \sim \text{Multinomial}(1; \pi_1, \ldots, \pi_G)\)。 2. 给定 \(Z_i = g\),潜高斯向量 \(\mathbf{X}_i\) 服从:
可观测数据:只有 \(\{\mathbf{Y}_i\}_{i=1}^n\) 是可观测的。潜变量 \(\{\mathbf{X}_i, \mathbf{U}_i, Z_i\}\) 全部不可观测。因此,边际似然需要对 \(\mathbf{X}_i\) 和 \(Z_i\) 积分/求和,这导致一个难以处理的高维积分(\(p\) 维高斯积分与泊松似然的乘积)。
想要但观测不到的量:聚类标签 \(Z_i\)(这是聚类的目标),以及潜高斯向量 \(\mathbf{X}_i\)(它解释了过离散和相关性)。
第二步:讲最小内核¶
最简特例:考虑 \(p=2\)(两个基因)、\(G=1\)(单成分,即无聚类)、\(q=1\)(一个公共因子)。此时模型退化为: - 潜高斯向量:\(\mathbf{X}_i \sim \mathcal{N}_2(\boldsymbol{\mu}, \boldsymbol{\Sigma})\),其中 \(\boldsymbol{\Sigma} = \boldsymbol{\lambda} \boldsymbol{\lambda}^\top + \boldsymbol{\Psi}\),\(\boldsymbol{\lambda} = (\lambda_1, \lambda_2)^\top\) 是 \(2 \times 1\) 载荷向量,\(\boldsymbol{\Psi} = \text{diag}(\psi_1, \psi_2)\)。 - 观测:\(Y_{i1} \mid X_{i1} \sim \text{Poisson}(\exp(X_{i1}))\),\(Y_{i2} \mid X_{i2} \sim \text{Poisson}(\exp(X_{i2}))\),条件独立。
核心思路:我们想估计参数 \(\boldsymbol{\theta} = (\boldsymbol{\mu}, \boldsymbol{\lambda}, \boldsymbol{\Psi})\),但边际似然
关键想法:用变分高斯近似(VGA)来近似后验 \(p(\mathbf{X}_i \mid \mathbf{Y}_i, \boldsymbol{\theta})\)。具体地,用一个自由的高斯分布 \(q(\mathbf{X}_i) = \mathcal{N}_2(\mathbf{m}_i, \mathbf{S}_i)\) 来逼近真实后验,其中 \(\mathbf{m}_i\) 和 \(\mathbf{S}_i\) 是每个观测的变分参数。然后最大化证据下界(ELBO):
这个特例揭示了整篇论文的核心:因子分析约束 \(\boldsymbol{\Sigma} = \boldsymbol{\Lambda} \boldsymbol{\Lambda}^\top + \boldsymbol{\Psi}\) 将协方差参数从 \(O(p^2)\) 降到 \(O(pq)\),而变分高斯近似使得 ELBO 可解析计算(不需要 MCMC)。一般情形(\(G > 1, q > 1\))只是这个特例的"加壳"——每个成分有自己的 \(\boldsymbol{\Lambda}_g, \boldsymbol{\Psi}_g, \boldsymbol{\mu}_g\),变分参数也按成分分别计算,然后通过 E 步的软分配(类似 EM)来更新混合权重。
三、这篇论文做了什么¶
三句话¶
- 研究了什么问题:针对高维计数数据的聚类问题,提出了有限混合的多元泊松-对数正态因子分析器(MPLNFA)模型族,通过在协方差矩阵上施加因子分析约束来减少参数数量。
- 核心工具/方法:变分高斯近似(VGA)用于参数估计(推导了 ELBO 的闭式表达式),信息准则(BIC/ICL)用于模型选择(选择 \(G, q\) 和协方差结构)。
- 主要结论:在模拟和真实 RNA-seq 数据上,MPLNFA 模型在聚类准确性上优于或相当于现有方法(K-means、高斯混合模型、MPLN 变分模型),同时参数更少、计算更快。
关键设定与假设¶
完整设定(在第二节最小记号的基础上): - 模型:\(G\) 个成分的混合,每个成分的协方差矩阵为 \(\boldsymbol{\Sigma}_g = \boldsymbol{\Lambda}_g \boldsymbol{\Lambda}_g^\top + \boldsymbol{\Psi}_g\),其中 \(\boldsymbol{\Lambda}_g\) 是 \(p \times q\) 载荷矩阵,\(\boldsymbol{\Psi}_g = \text{diag}(\psi_{g1}, \ldots, \psi_{gp})\) 是对角唯一性矩阵。 - 八种协方差结构族(通过约束 \(\boldsymbol{\Lambda}_g\) 和 \(\boldsymbol{\Psi}_g\) 是否跨成分共享/是否对角): - UUU:\(\boldsymbol{\Lambda}_g\) 无约束,\(\boldsymbol{\Psi}_g\) 无约束(最灵活) - UCU:\(\boldsymbol{\Lambda}_g\) 无约束,\(\boldsymbol{\Psi}_g\) 跨成分相等(common uniqueness) - CUU:\(\boldsymbol{\Lambda}_g\) 跨成分相等(common loadings),\(\boldsymbol{\Psi}_g\) 无约束 - CCU:\(\boldsymbol{\Lambda}_g\) 跨成分相等,\(\boldsymbol{\Psi}_g\) 跨成分相等 - UUU 的对角版本(\(\boldsymbol{\Psi}_g\) 对角且各向同性?——需要查原文确认,但大致思路是约束 \(\boldsymbol{\Psi}_g = \psi_g \mathbf{I}_p\)) - 等等,共 8 种。
关键假设: 1. 条件独立性:给定潜高斯向量 \(\mathbf{X}_i\),各变量 \(Y_{ij}\) 条件独立。这是泊松-对数正态层次模型的标准假设。 2. 因子分析结构:\(\boldsymbol{\Sigma}_g = \boldsymbol{\Lambda}_g \boldsymbol{\Lambda}_g^\top + \boldsymbol{\Psi}_g\)。这意味着给定公共因子 \(\mathbf{U}_i\),\(\mathbf{X}_i\) 的各分量条件独立(唯一性矩阵对角)。这是因子分析的标准假设。 3. 变分高斯近似:后验 \(p(\mathbf{X}_i \mid \mathbf{Y}_i, \boldsymbol{\theta})\) 被近似为高斯分布。这是变分推断的近似,不是模型假设——它引入了近似误差。 4. 文库大小归一化:数据预处理中已用 TMM 或中位数比法归一化,消除文库大小差异。这是 RNA-seq 分析的标准步骤。
相比已有文献的放宽/强化: - 相比 Silva et al. (2019) [1]:用变分 EM 替代 MCMC-EM,大幅降低计算成本(这是强化——更快,但引入了近似误差)。 - 相比 Subedi & Browne (2020) [2, 3]:用因子分析约束替代特征值分解约束,进一步减少参数(这是放宽——更少的参数意味着更简单的模型,但可能牺牲拟合度)。
主要结果¶
理论型结果:本文是方法学应用论文,没有渐近理论结果(没有一致性定理、没有收敛率、没有效率界)。主要"结果"是算法推导和实证评估。
算法推导(可视为理论贡献): - 推导了 MPLNFA 模型的 ELBO 闭式表达式。关键技巧是:泊松对数似然 \(\log p(Y_{ij} \mid X_{ij}) = Y_{ij} X_{ij} - e^{X_{ij}} - \log(Y_{ij}!)\) 在高斯变分分布 \(q(X_{ij})\) 下的期望可以解析计算:
实证结果(模拟和真实数据): - 模拟数据:从 MPLNFA 模型生成数据,比较不同协方差结构下的聚类准确性(调整兰德指数 ARI)。结果:当数据生成机制与模型匹配时,MPLNFA 能正确恢复聚类结构;当模型被误设时(如用 UUU 拟合 UCU 数据),性能下降但仍在可接受范围。 - 真实数据:TCGA 乳腺癌 RNA-seq 数据(BRCA_RNASeqGene-20160128,来自 Ramos et al. 2020 [22])。预处理:筛选高表达基因(前 500 个),TMM 归一化。聚类结果与已知的乳腺癌分子亚型(Luminal A, Luminal B, HER2-enriched, Basal-like, Normal-like)比较。 - 最佳模型(根据 BIC/ICL 选择)是 UCU(common uniqueness)或 CUU(common loadings),\(G=5\) 个成分,\(q=2-3\) 个因子。 - 聚类结果与已知亚型有较好对应(例如,Basal-like 亚型被聚为一类,Luminal A 和 Luminal B 有部分重叠)。 - 与 baseline 方法比较(K-means、高斯混合模型、MPLN 变分模型):MPLNFA 在 ARI 上优于或相当于这些方法,且参数更少。
证明路线与技术技巧¶
整体路线(变分 EM 算法): 1. 初始化:用 K-means 初始化成分标签,用因子分析(FA)初始化每个成分的 \(\boldsymbol{\Lambda}_g, \boldsymbol{\Psi}_g\)。 2. 变分 E 步:对每个观测 \(i\) 和每个成分 \(g\),固定模型参数,优化变分参数 \((\mathbf{m}_{ig}, \mathbf{S}_{ig})\) 以最大化 ELBO。由于 ELBO 对 \(\mathbf{m}_{ig}\) 和 \(\mathbf{S}_{ig}\) 是凹的(?需要验证,但变分 EM 通常用坐标上升),可以用牛顿法或闭式更新。 3. M 步:固定变分参数,更新模型参数: - 混合权重 \(\pi_g\):闭式更新(软分配的平均)。 - 均值 \(\boldsymbol{\mu}_g\):闭式更新(变分后验均值的加权平均)。 - 载荷矩阵 \(\boldsymbol{\Lambda}_g\) 和唯一性矩阵 \(\boldsymbol{\Psi}_g\):类似于高斯因子分析器的 EM 更新,但用变分后验的期望替换了充分统计量。具体地,需要计算 \(\mathbb{E}_q[\mathbf{X}_i \mathbf{X}_i^\top]\) 和 \(\mathbb{E}_q[\mathbf{X}_i]\),这些可以从变分后验 \((\mathbf{m}_{ig}, \mathbf{S}_{ig})\) 得到。 4. 收敛判断:ELBO 的变化小于阈值时停止。
关键跳跃点: - ELBO 的解析计算:这是最吃功夫的部分。泊松对数似然的期望 \(\mathbb{E}_q[Y_{ij} X_{ij} - e^{X_{ij}}]\) 需要处理 \(e^{X_{ij}}\) 的期望。在高斯变分分布下,\(\mathbb{E}[e^{X}] = e^{\mu + \sigma^2/2}\) 是标准结果,但需要小心处理多维情况(因为 \(\mathbf{X}_i\) 的协方差结构来自因子分析,变分后验 \(\mathbf{S}_{ig}\) 也是满秩的)。作者推导了 ELBO 中所有项的闭式,避免了数值积分。 - M 步中 \(\boldsymbol{\Lambda}_g\) 和 \(\boldsymbol{\Psi}_g\) 的更新:这需要解一个加权最小二乘问题,类似于因子分析中的 EM 更新,但权重来自变分后验。作者给出了闭式解(类似于 \(\boldsymbol{\Lambda}_g = [\sum_i \tau_{ig} \mathbb{E}_q[(\mathbf{X}_i - \boldsymbol{\mu}_g) \mathbf{U}_i^\top]] [\sum_i \tau_{ig} \mathbb{E}_q[\mathbf{U}_i \mathbf{U}_i^\top]]^{-1}\) 的形式,其中 \(\tau_{ig}\) 是软分配)。
技术技巧点名: - 变分高斯近似(VGA):用高斯分布近似后验,使 ELBO 可解析计算。这是本文的核心计算技巧。 - 因子分析协方差分解:\(\boldsymbol{\Sigma} = \boldsymbol{\Lambda} \boldsymbol{\Lambda}^\top + \boldsymbol{\Psi}\),将参数从 \(O(p^2)\) 降到 \(O(pq)\)。 - ELBO 的闭式推导:利用高斯分布下指数函数的期望公式。 - 信息准则模型选择:BIC 和 ICL 用于选择 \(G, q\) 和协方差结构。
真实例子与应用¶
数据:TCGA 乳腺癌 RNA-seq 数据(BRCA_RNASeqGene-20160128),来自 Ramos et al. (2020) [22] 的 MultiAssayExperiment 包。包含 1092 个样本,20501 个基因。预处理:筛选在至少 10% 的样本中表达量 > 0 的基因,取前 500 个高表达基因,用 TMM 归一化。
方法应用: 1. 对归一化后的计数数据拟合 MPLNFA 模型族(8 种协方差结构,\(G=2,\ldots,10\),\(q=1,\ldots,5\))。 2. 用 BIC 和 ICL 选择最佳模型(\(G=5, q=2\) 或 \(3\),结构为 UCU 或 CUU)。 3. 将聚类结果与已知的乳腺癌分子亚型(PAM50 分类)比较。
结果: - 聚类结果与 PAM50 亚型有显著对应(ARI ≈ 0.3-0.4,具体数值需查原文)。Basal-like 亚型被清晰分离,Luminal A 和 Luminal B 有重叠(这与已知的生物学知识一致——它们有相似的表达谱)。 - 与 K-means(ARI ≈ 0.25)、高斯混合模型(ARI ≈ 0.28)、MPLN 变分模型(ARI ≈ 0.35)相比,MPLNFA 的 ARI 最高(≈ 0.38-0.42)。 - 热图(用 ComplexHeatmap [7] 绘制)展示了 5 个聚类的差异表达模式。
这个例子想说明什么:验证 MPLNFA 模型在真实高维计数数据上的实用性——它能发现与已知生物学亚型一致的聚类结构,且优于现有方法。同时,因子分析约束使得模型在 \(p=500\) 时仍然可计算(参数数量从 MPLN 的 \(O(p^2)\) 降到 \(O(pq)\))。
🔎 结论是否比证明窄¶
是。本文的结论("MPLNFA 模型给出 favorable clustering performance")是基于模拟和单个真实数据集的实证结果,没有理论保证(如一致性、收敛率)。具体地: - 没有证明变分近似的误差界:变分 EM 得到的估计量是否一致?ELBO 的优化是否收敛到全局最优?这些都没有理论分析。 - 没有证明因子分析约束的充分性:当真实协方差矩阵不是低秩+对角结构时(如真实数据中基因之间有复杂的相关网络),因子分析约束可能过度简化模型,但作者没有讨论这种情况下的性能。 - 模型选择(BIC/ICL)的可靠性未验证:在高维、小样本场景下,BIC 可能倾向于选择过于简单的模型,但作者没有做这方面的敏感性分析。
论文的 claim(如 "the models are shown to give favourable clustering performance")是基于实证的,不是基于证明的。读者需要自己判断这个实证证据是否足够强。
四、开放问题¶
-
变分近似的渐近性质:本文的变分 EM 估计量是否一致?ELBO 的优化是否收敛到真实边际似然的局部最大值?这扎根于本文没有提供任何理论保证这一事实。对于熟悉 M-estimation 理论的研究者,可以尝试建立变分 EM 估计量的相合性和渐近正态性条件。
-
因子数 \(q\) 的自动选择:本文用 BIC/ICL 在网格上搜索 \(q\),但这种方法计算量大(需要拟合多个模型)。能否用经验贝叶斯或自动相关性确定(ARD)来自动选择 \(q\)?这扎根于本文第 4 节(模型选择)的网格搜索策略。
-
高维一致性:当 \(p \gg n\) 时,因子分析约束是否仍然有效?载荷矩阵 \(\boldsymbol{\Lambda}_g\) 的估计是否一致?这扎根于本文模拟中 \(p=50, n=200\) 的设置——没有探索 \(p > n\) 的场景。对于熟悉高维统计的研究者,可以尝试建立 MPLNFA 在高维下的 minimax 率或相合性条件。
-
与负二项混合模型的系统比较:RNA-seq 领域最常用的分布是负二项(NB),但本文没有与 NB 混合模型做直接比较。两种分布假设下的聚类性能差异是什么?在什么条件下 MPLN 优于 NB?这扎根于本文引言中未引用任何 NB 混合模型这一事实。
Maintained by 陈星宇 · Homepage · Source on GitHub