Zero-inflation in the multivariate poisson lognormal family¶

作者: Bastien Batardière, Julien Chiquet, François Gindraud, Mahendra Mariadassou
来源: Statistics and Computing
主题: 统计计算 / 算法
相关性: 6/10
链接: 期刊页 · arXiv

一、领域脉络与小综述¶

这个方向是什么¶

本方向关注的是高维计数数据的统计建模，特别是当数据中存在大量“零”值（零膨胀）时，如何构建一个既能捕捉变量间复杂依赖关系、又能处理零膨胀、且具有可解释性和计算可扩展性的概率模型。当前成熟度较高，已有如泊松-对数正态（PLN）模型等有效框架，但零膨胀问题仍是活跃的研究前沿。

发展脉络（history）¶

奠基工作：Akaike (1974) 提出的AIC准则为模型选择提供了基础工具，至今仍被用于比较不同复杂度的模型。O’Hara & Kotze (2010) 则明确指出对计数数据进行对数变换的弊端，为直接建模计数数据（而非变换后数据）的必要性提供了关键论据。
主要进展（PLN模型及其变分推断）：Chiquet, Robin, Mariadassou 等人 的一系列工作奠定了PLN模型作为高维计数数据分析通用框架的地位。Chiquet et al. (2018) 将PLN模型与稀疏网络推断结合，Chiquet et al. (2017) 提出了概率性泊松PCA的变分推断，而 Chiquet et al. (2020) 则系统总结了PLN模型在群落生态学中的多功能性（样本比较、聚类、网络推断等）。这些工作的核心是使用变分推断（VI） 来近似处理PLN模型中因潜变量积分导致的难解似然。Hui et al. (2017) 和 Blei et al. (2016) 的综述为VI在广义线性潜变量模型中的应用提供了理论基础和方法论。
当前Frontier（零膨胀与单细胞数据）：在单细胞RNA测序（scRNA-seq）领域，零膨胀（或称“dropout”）是核心挑战。Risso et al. (2017) 提出了ZINB-WaVE模型，使用零膨胀负二项分布。Lopez et al. (2018) 的scVI和 Wang & Gu (2018) 的VASC则利用深度生成模型（VAE）来处理零膨胀。Choudhary & Satija (2021) 的大规模比较研究指出，对于稀疏数据泊松模型合适，但对于有足够测序深度的基因，负二项模型（即存在过离散）是必要的。Lee et al. (2017) 则从贝叶斯变量选择角度处理微生物组数据中的零膨胀问题。
本文的位置：本文作者（Chiquet, Mariadassou等）是PLN模型的主要推动者。本文是他们在PLN框架内工作的自然延伸：将零膨胀成分显式地、以多元方式整合进PLN模型，提出ZIPLN模型。它区别于ZINB-WaVE（使用负二项）和scVI/VAE（使用深度网络），保留了PLN模型的潜变量高斯结构和可解释性，同时通过变分推断保持计算可扩展性。

子线索聚类¶

PLN模型及其变分推断：核心是Chiquet, Mariadassou, Robin等人的工作（2017, 2018, 2020）。他们建立了PLN模型作为高维计数数据分析的通用框架，并发展了基于变分推断的高效估计算法。这是本文的直接基础。
零膨胀计数模型：包括ZINB-WaVE (Risso et al., 2017)、贝叶斯零膨胀模型 (Lee et al., 2017)、双变量ZINB模型 (Cho et al., 2020) 等。这些工作证明了零膨胀建模的必要性，但大多基于负二项分布或处理特定场景。本文则是在泊松-对数正态框架下处理零膨胀。
深度生成模型（VAE）：scVI (Lopez et al., 2018)、VASC (Wang & Gu, 2018)、siVAE (Choi et al., 2023) 等。这些方法利用深度神经网络进行非线性降维和零膨胀建模，计算能力强，但可解释性较弱。本文的ZIPLN模型与之形成对比，强调可解释性。
变分推断的理论与计算：Blei et al. (2016) 的综述和 Westling & McCormick (2015) 关于变分估计渐近性质的工作，为VI提供了理论支撑。Liu & Zhong (2024) 提出的拉普拉斯-泰勒近似方法，被本文引用作为改进变分E步的潜在方向。

这个方向在追问的核心问题¶

如何在高维、零膨胀的计数数据中，同时估计变量间的依赖结构（协方差矩阵）和零膨胀机制？
如何在保持模型可解释性（如潜变量结构）的同时，实现可扩展至数千变量的高效计算？
对于零膨胀的建模，是采用“额外伯努利过程”（如本文）还是“负二项分布”（如ZINB-WaVE）更合适？ 前者将零分为“结构零”和“采样零”，后者则通过过离散来吸收额外零。
变分推断的近似误差对参数估计和后续推断（如假设检验）有何影响？ 这是Westling & McCormick (2015) 试图回答的问题，也是本文未深入探讨的。

⚠️ 作者的 framing¶

作者的缺口：作者将缺口frame为“PLN模型不处理零膨胀，而真实数据（如微生物组）中零膨胀普遍存在”。因此，ZIPLN模型是PLN模型的“显然的下一步”扩展。
被淡化的竞争路线：作者淡化了深度生成模型（如scVI）的路线，强调ZIPLN的可解释性（“preserves explainability”）。他们未直接与scVI等模型进行性能比较，而是与PLN模型本身对比，以展示零膨胀成分的贡献。
值得研究者去查的问题：为什么本文没有引用或讨论ZINB-WaVE (Risso et al., 2017)？ ZINB-WaVE同样是一个处理零膨胀的多元潜变量模型，且同样使用变分推断。它与ZIPLN的核心区别在于：ZINB-WaVE使用负二项分布（通过一个额外的离散参数处理过离散），而ZIPLN使用泊松-对数正态（通过潜变量方差处理过离散）。比较这两种零膨胀建模策略（伯努利-泊松 vs. 负二项）的优劣，是一个值得探索的张力点。 此外，Westling & McCormick (2015) 关于变分估计渐近性质的工作，本文仅引用但未深入讨论其对ZIPLN参数估计标准误和假设检验的影响，这也是一个潜在缺口。

张力¶

未见明显对立引用。不同工作（如PLN vs. ZINB-WaVE vs. scVI）更多是建模哲学和计算效率的权衡，而非结论性矛盾。

二、最核心、最简单的例子 / 数学问题¶

第一步：把符号、模型、可观测数据交代清楚¶

符号：
- \( n \)：样本量（如微生物组研究中的样本数）。
- \( p \)：特征数（如微生物物种或OTU的数量）。
- \( Y_{ij} \)：可观测的计数数据，表示第 \( i \) 个样本中第 \( j \) 个特征的计数。\( Y \) 是一个 \( n \times p \) 的矩阵。
- \( Z_{ij} \)：不可观测的伯努利潜变量，表示第 \( i \) 个样本中第 \( j \) 个特征是否处于“结构零”状态。\( Z_{ij} = 1 \) 表示结构零（计数必然为0），\( Z_{ij} = 0 \) 表示非结构零（计数可能为0或正数）。
- \( \pi_{ij} \)：零膨胀概率，即 \( P(Z_{ij} = 1) \)。它是模型参数，可以依赖于协变量。
- \( X_i \)：可观测的协变量向量（\( d \) 维），用于解释均值或零膨胀概率。
- \( W_i \)：不可观测的 \( q \) 维潜变量（通常 \( q \ll p \)），服从多元高斯分布。它捕捉了特征间的依赖结构。
- \( B \)：\( p \times d \) 的系数矩阵，将协变量 \( X_i \) 映射到对数均值。
- \( C \)：\( p \times q \) 的载荷矩阵，将潜变量 \( W_i \) 映射到对数均值。
- \( \Sigma \)：\( q \times q \) 的协方差矩阵，描述潜变量 \( W_i \) 的依赖结构。通常假设 \( \Sigma = I_q \) 或 \( C C^\top \) 来捕捉特征间依赖。
- \( \Theta \)：模型参数集合，包括 \( B, C, \Sigma \) 以及零膨胀参数（如 \( \pi_{ij} \) 中的系数）。
模型：
1. 零膨胀机制：\( Z_{ij} \sim \text{Bernoulli}(\pi_{ij}) \)。如果 \( Z_{ij} = 1 \)，则 \( Y_{ij} = 0 \) 是确定的。
2. 计数生成机制：如果 \( Z_{ij} = 0 \)，则 \( Y_{ij} \) 服从泊松分布，其对数均值由一个线性潜变量模型决定：
  \[Y_{ij} | (Z_{ij}=0, W_i) \sim \text{Poisson}(\exp(\eta_{ij}))\]
  其中 \( \eta_{ij} = X_i^\top B_j + W_i^\top C_j \)，\( B_j \) 是 \( B \) 的第 \( j \) 列，\( C_j \) 是 \( C \) 的第 \( j \) 行。
3. 潜变量分布：\( W_i \sim \mathcal{N}_q(0, \Sigma) \)。
可观测数据：
- 可观测：计数矩阵 \( Y \)（\( n \times p \)）和协变量矩阵 \( X \)（\( n \times d \)）。
- 不可观测 / 潜在：零膨胀指示变量 \( Z \)（\( n \times p \)）和潜变量 \( W \)（\( n \times q \)）。我们只能通过模型假设和观测到的 \( Y \) 来推断它们。

第二步：讲最小内核¶

最简特例：假设我们只有一个样本（\( n=1 \)）和一个特征（\( p=1 \)），且没有协变量（\( d=0 \)）。那么模型退化为一个单变量零膨胀泊松-对数正态模型。

记号简化：去掉下标 \( i, j \)。我们有：
- 可观测：\( Y \)（一个计数）。
- 潜变量：\( Z \in \{0,1\} \)（伯努利），\( W \sim \mathcal{N}(0, \sigma^2) \)（一维高斯）。
- 参数：零膨胀概率 \( \pi \)，对数均值截距 \( \mu \)，潜变量方差 \( \sigma^2 \)。
模型：
\[P(Y = y) = \begin{cases} \pi + (1-\pi) \cdot P_{\text{PLN}}(Y=0), & \text{if } y = 0 \\ (1-\pi) \cdot P_{\text{PLN}}(Y=y), & \text{if } y > 0 \end{cases}\]
其中 \( P_{\text{PLN}}(Y=y) = \int_{-\infty}^{\infty} \frac{e^{-e^{\mu + w}} (e^{\mu + w})^y}{y!} \cdot \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{w^2}{2\sigma^2}} dw \)。
核心思路：这个模型要解决的核心问题是：观测到的零值，有多少是“结构零”（由 \( Z=1 \) 导致），有多少是“采样零”（由泊松过程恰好抽到0导致）？
- 如果没有零膨胀（\( \pi=0 \)），模型退化为标准的单变量PLN。此时，零值只能来自泊松抽样。
- 如果引入零膨胀，模型允许一部分零值由独立的伯努利过程产生。这为数据中“过多”的零值提供了一个额外的解释。
为什么这个特例是核心：整篇论文的多元ZIPLN模型，本质上就是将这个单变量特例推广到 \( p \) 个特征和 \( n \) 个样本，并引入协变量和潜变量依赖结构。所有技术难点（如变分推断、ELBO推导）都源于这个推广。理解了这个最简特例，就理解了ZIPLN模型的核心思想：用一个额外的伯努利潜变量来“吸收”多余的零，从而更准确地估计泊松过程的参数（\( \mu, \sigma^2 \)）和变量间的依赖关系。

三、这篇论文做了什么¶

三句话¶

研究了什么问题：针对高维计数数据中常见的零膨胀现象，在多元泊松-对数正态（PLN）模型基础上，引入一个额外的伯努利潜变量来显式建模零膨胀机制，提出了ZIPLN模型。
核心工具/方法：采用变分推断（VI） 进行参数估计，比较了两种变分近似策略：独立高斯与伯努利变分分布（mean-field）和以伯努利为条件的高斯变分分布（structured）。推导了相应的变分下界（ELBO）和梯度，并开发了R包ZIPLN。
主要结论：模拟实验表明，即使在零膨胀比例高达90%时，ZIPLN仍能有效恢复参数。在牛微生物组数据集（零占比90.6%）上的应用显示，考虑零膨胀显著提升了对数似然并降低了潜空间离散度，从而改善了组别判别。

关键设定与假设¶

模型设定：在第二节最小记号的基础上，完整设定如下：
- \( Y_{ij} | Z_{ij}, W_i \sim \text{Poisson}(\exp(\eta_{ij})) \)，其中 \( \eta_{ij} = X_i^\top B_j + W_i^\top C_j \)。
- \( Z_{ij} \sim \text{Bernoulli}(\pi_{ij}) \)，其中 \( \pi_{ij} \) 可以是常数（固定）、依赖于特征 \( j \)（特征特异）、依赖于样本 \( i \)（位点特异）或依赖于协变量 \( X_i \)。
- \( W_i \sim \mathcal{N}_q(0, \Sigma) \)，通常假设 \( \Sigma = I_q \) 以简化，特征间依赖由载荷矩阵 \( C \) 和潜变量 \( W \) 共同捕捉。
假设：
- 条件独立性：给定潜变量 \( W_i \) 和零膨胀指示 \( Z_{ij} \)，不同特征的计数 \( Y_{ij} \) 是条件独立的。这是PLN模型的核心假设，也是其可解释性的来源。
- 潜变量高斯性：\( W_i \) 服从多元高斯分布。这是一个建模假设，便于计算和解释。
- 变分分布族假设：变分推断假设一个近似的后验分布 \( q(Z, W) \) 来逼近真实后验 \( p(Z, W | Y) \)。本文考虑了两种形式：
  1. Mean-field (MF)：\( q(Z, W) = \prod_{i,j} q(Z_{ij}) \prod_i q(W_i) \)，其中 \( q(Z_{ij}) \) 是伯努利，\( q(W_i) \) 是高斯。
  2. Structured (STR)：\( q(Z, W) = \prod_{i,j} q(Z_{ij}) \prod_i q(W_i | Z_i) \)，其中 \( q(W_i | Z_i) \) 是条件高斯，其均值和方差依赖于 \( Z_i \)。这个近似更精确，但计算更复杂。
相比已有文献：相比标准PLN模型（Chiquet et al., 2020），本文增加了零膨胀成分。相比ZINB-WaVE（Risso et al., 2017），本文使用泊松-对数正态而非负二项分布，保留了PLN的潜变量结构。相比深度VAE模型（Lopez et al., 2018），本文的变分分布更简单（高斯/伯努利），计算更快，可解释性更强。

主要结果¶

理论结果：本文主要是方法论文，没有提供渐近理论（如估计量的相合性或渐近正态性）。核心“结果”是推导了两种变分近似下的ELBO和梯度更新公式。
- MF-ELBO：一个封闭形式的ELBO，包含与 \( Z \) 和 \( W \) 相关的项。梯度可以解析计算。
- STR-ELBO：由于 \( q(W_i | Z_i) \) 依赖于 \( Z_i \)，ELBO中包含一个对 \( Z \) 的期望项，无法解析计算。作者使用一阶泰勒展开来近似 \( \mathbb{E}_q[\exp(-W_i^\top C_j)] \) 这一项，从而得到一个可计算的近似ELBO。
模拟实验：
- 设定：生成 \( n=100 \) 个样本，\( p=50 \) 或 \( 100 \) 个特征，潜变量维度 \( q=2 \)。零膨胀比例从 \( 0\% \) 到 \( 90\% \) 变化。
- 核心量化结论：
  - 参数恢复：ZIPLN模型（MF和STR）能够有效恢复载荷矩阵 \( C \) 和回归系数 \( B \)。随着零膨胀比例增加，估计误差（如Frobenius范数）略有增加，但即使在90%零膨胀下，ZIPLN仍显著优于忽略零膨胀的PLN模型。
  - 模型选择：AIC和BIC能够正确选择零膨胀模型（ZIPLN）而非PLN模型，尤其是在零膨胀比例较高时。
  - 变分近似比较：STR近似在参数恢复上通常略优于MF近似，但计算时间更长。两者差距不大，MF近似在计算效率上更具优势。
真实数据例子：
- 数据：来自 Mariadassou et al. (2023) 的牛微生物组数据集。包含45头奶牛，4个时间点，5个身体部位（口腔、鼻腔、阴道、牛奶），共 \( n=900 \) 个样本，\( p=200 \) 个最丰富的ASV。数据中零的比例高达 90.6%。
- 方法应用：将ZIPLN和PLN模型应用于该数据，潜变量维度 \( q=5 \)。模型包含协变量（时间、部位、个体）。
- 结果：
  - 对数似然：ZIPLN模型的对数似然显著高于PLN模型（例如，MF-ZIPLN的ELBO为 -1.42e5，而PLN为 -1.56e5），表明ZIPLN更好地拟合了数据。
  - 潜空间离散度：ZIPLN模型估计的潜变量 \( W_i \) 在潜空间中的离散度（即样本间的变异）显著小于PLN模型。作者解释为：PLN模型为了解释大量零值，不得不将样本在潜空间中“拉开”，而ZIPLN通过零膨胀成分吸收了这些零，使得潜变量能更真实地反映样本间的生物学差异。
  - 组别判别：基于潜变量 \( W_i \) 的样本聚类（如按身体部位）在ZIPLN模型中比在PLN模型中更清晰、分离度更好。这直接验证了上述“降低离散度、改善判别”的论点。
- 这个例子想说明什么：该例子旨在证明，在真实的高零膨胀数据中，忽略零膨胀（使用PLN）会导致对潜变量结构的扭曲估计，而ZIPLN模型能够更准确地揭示数据背后的生物学结构（如不同身体部位的微生物群落差异）。

证明路线与技术技巧¶

本文是方法论文，没有传统意义上的“定理证明”。其“证明路线”是变分推断算法的推导。

整体路线：
1. 写出完全数据对数似然：\( \log p(Y, Z, W | \Theta) \)。
2. 写出证据下界（ELBO）：\( \mathcal{L}(q, \Theta) = \mathbb{E}_q[\log p(Y, Z, W | \Theta)] - \mathbb{E}_q[\log q(Z, W)] \)。目标是最大化ELBO来同时估计变分参数和模型参数。
3. 选择变分分布族：选择MF或STR近似。
4. 推导ELBO的封闭形式：对于MF近似，ELBO可以解析计算。对于STR近似，需要处理一个棘手的期望项 \( \mathbb{E}_q[\exp(-W_i^\top C_j)] \)。
5. 处理STR-ELBO中的棘手项：这是关键跳跃点。
  - 难点：\( \mathbb{E}_q[\exp(-W_i^\top C_j)] \) 涉及对条件高斯分布 \( q(W_i | Z_i) \) 的期望，而 \( Z_i \) 本身是随机的，导致无法得到封闭形式。
  - 作者的技巧：使用一阶泰勒展开来近似 \( \exp(-x) \approx 1 - x \)。这使得期望变为 \( 1 - \mathbb{E}_q[W_i^\top C_j] \)，而 \( \mathbb{E}_q[W_i^\top C_j] \) 可以解析计算。作者也提到，可以使用二阶泰勒展开（如Liu & Zhong, 2024所建议）来获得更精确的近似。
6. 坐标上升法优化：交替更新变分参数（\( q(Z), q(W) \) 的参数）和模型参数（\( B, C, \Sigma, \pi \) 的参数）。这构成了一个变分EM算法。
技术技巧点名：
- 变分推断（VI）：核心工具，用于近似后验分布。
- 泰勒展开：用于近似STR-ELBO中难解的期望项，是处理结构化变分近似的关键。
- 坐标上升法：用于优化ELBO。
- 自动微分：利用TMB或PyTorch等工具自动计算梯度，简化了算法实现。

🔎 结论是否比证明窄¶

是。本文的“结论”主要基于模拟实验和单一真实数据案例。作者没有提供任何关于ZIPLN估计量的渐近性质（如相合性、渐近正态性、效率）的理论证明。文中提到“variational inference lacks theoretical guarantees on the estimates”（引用Stoehr & Robin, 2024），并承认这是VI方法的普遍缺点。因此，ZIPLN模型的有效性目前仅由有限的经验证据支持，其统计推断（如置信区间、假设检验）的理论基础尚不牢固。作者在结论部分也明确将“theoretical guarantees for the variational estimator”列为未来工作。

四、开放问题¶

变分估计的渐近理论：ZIPLN模型的变分估计量是否具有相合性和渐近正态性？其渐近方差是多少？如何构建有效的置信区间和假设检验？这扎根于本文引言中引用的 Westling & McCormick (2015) 和 Stoehr & Robin (2024) 的工作，以及作者在结论中明确指出的未来方向。
零膨胀机制与过离散的区分：ZIPLN模型使用伯努利-泊松机制处理零膨胀，而ZINB-WaVE使用负二项分布。在何种数据生成机制下，一种模型会优于另一种？是否存在模型选择准则来区分这两种不同的零膨胀/过离散来源？这扎根于 Choudhary & Satija (2021) 关于不同误差模型比较的工作，以及本文未与ZINB-WaVE进行直接比较的现状。
结构化变分近似的改进：本文使用的STR近似依赖于一阶泰勒展开。使用二阶泰勒展开（如 Liu & Zhong (2024) 所建议）是否能带来显著的性能提升？或者是否存在其他更精确的近似方法（如重要性采样、重参数化技巧）？这扎根于本文第4.2节关于STR近似的讨论。
高维协方差估计：本文假设潜变量协方差 \( \Sigma = I_q \)，特征间依赖由载荷矩阵 \( C \) 捕捉。如果直接对 \( p \times p \) 的协方差矩阵进行稀疏估计（如Chiquet et al., 2018），ZIPLN模型的计算和理论会如何变化？这扎根于PLN模型在网络推断中的应用，以及本文未探索高维协方差结构的现状。

Maintained by 陈星宇 · Homepage · Source on GitHub