676 微生物组零膨胀模型¶

# 零膨胀模型的数学直觉

# 普通负二项分布（NB）：
# P(Y=0) 只来自NB分布本身
# 问题：真实数据的零值比NB预测的多得多

# 零膨胀负二项分布（ZINB）：
# P(Y=0) = π + (1-π) × P_NB(Y=0)
# π = 结构零的概率（"真的没有"）
# (1-π) × P_NB(Y=0) = 抽样零的概率（"有但没测到"）
# P(Y=k) = (1-π) × P_NB(Y=k), k > 0

# Hurdle模型：
# 第一部分（二项）：P(Y=0) vs P(Y>0) → logistic回归
# 第二部分（截断NB）：P(Y=k | Y>0) → 只对非零值建模
# 区别：hurdle模型不区分两种零，把所有零归为一类

# 零膨胀模型分析
library(pscl)         # 零膨胀模型
library(MASS)         # 负二项分布
library(glmmTMB)      # 广义线性混合模型（支持ZI）
library(DHARMa)       # 模型诊断

# ============ 准备数据 ============
# 假设有OTU表和元数据
otu <- read.csv("otu_table.csv", row.names = 1)  # OTU计数表
meta <- read.csv("metadata.csv", row.names = 1)  # 元数据

# 选择一个物种分析
species <- "Bacteroides_fragilis"        # 目标物种
y <- otu[, species]                      # 该物种的计数
zero_pct <- mean(y == 0) * 100           # 零值比例
cat(sprintf("%s: 零值比例=%.1f%%\n", species, zero_pct))

# ============ 1. 模型比较：选择最佳分布 ============
# 1a. 泊松模型（最简单，通常不适合）
m_pois <- glm(y ~ Group + Age, data = meta, family = poisson)

# 1b. 负二项模型（处理过度离散）
m_nb <- glm.nb(y ~ Group + Age, data = meta)

# 1c. 零膨胀泊松模型（ZIP）
m_zip <- zeroinfl(y ~ Group + Age | 1,  # 计数部分 | 零膨胀部分
                   data = meta,
                   dist = "poisson")     # 泊松分布

# 1d. 零膨胀负二项模型（ZINB）——最常用
m_zinb <- zeroinfl(y ~ Group + Age | Group,  # 零膨胀部分也可包含协变量
                    data = meta,
                    dist = "negbin")     # 负二项分布

# 1e. Hurdle模型
m_hurdle <- hurdle(y ~ Group + Age | Group,
                    data = meta,
                    dist = "negbin")     # 负二项分布

# 模型比较（AIC越小越好）
aic_table <- data.frame(
  Model = c("Poisson", "NB", "ZIP", "ZINB", "Hurdle"),
  AIC = c(AIC(m_pois), AIC(m_nb), AIC(m_zip), AIC(m_zinb), AIC(m_hurdle))
)
aic_table <- aic_table[order(aic_table$AIC), ]  # 按AIC排序
print(aic_table)
cat("最佳模型:", aic_table$Model[1], "\n")

# ============ 2. ZINB模型详细分析 ============
summary(m_zinb)
# 输出包含两部分：
# Count model coefficients (negbin): 计数部分的系数
#   → 在物种"存在"时，分组对丰度的影响
# Zero-inflation model coefficients (binomial): 零膨胀部分
#   → 分组对"结构零概率"的影响

# 提取系数
count_coefs <- coef(m_zinb, model = "count")   # 计数部分系数
zero_coefs <- coef(m_zinb, model = "zero")     # 零膨胀部分系数
cat("\n计数部分系数:\n")
print(round(count_coefs, 3))
cat("\n零膨胀部分系数:\n")
print(round(zero_coefs, 3))

# ============ 3. 使用glmmTMB（支持随机效应） ============
# 对于纵向/重复测量数据
library(glmmTMB)
m_glmmtmb <- glmmTMB(
  y ~ Group + Age + (1 | Subject_ID),   # 固定效应 + 随机效应
  ziformula = ~ Group,                   # 零膨胀部分
  family = nbinom2,                      # 负二项分布
  data = meta
)
summary(m_glmmtmb)

# ============ 4. 模型诊断 ============
# 使用DHARMa包进行残差诊断
sim_res <- simulateResiduals(m_zinb)     # 模拟残差
plot(sim_res)                            # 诊断图
testZeroInflation(sim_res)               # 检验零膨胀是否显著
# 如果p<0.05 → 确实存在零膨胀 → 需要ZI模型
# 如果p>0.05 → 普通NB模型就够了

# ============ 5. 批量分析多个物种 ============
analyze_all_species <- function(otu, meta, formula_count, formula_zero) {
  """对所有物种批量拟合ZINB模型"""
  results <- list()                      # 存储结果
  for (sp in colnames(otu)) {
    meta$y <- otu[, sp]                  # 当前物种计数
    zero_pct <- mean(meta$y == 0)        # 零值比例

    tryCatch({
      if (zero_pct > 0.3 & zero_pct < 0.95) {  # 零值30-95%用ZINB
        m <- zeroinfl(
          as.formula(paste("y ~", formula_count, "|", formula_zero)),
          data = meta, dist = "negbin"
        )
      } else if (zero_pct <= 0.3) {      # 零值<30%用普通NB
        m <- glm.nb(
          as.formula(paste("y ~", formula_count)),
          data = meta
        )
      } else {                            # 零值>95%跳过
        next
      }
      # 提取p值
      s <- summary(m)
      coefs <- s$coefficients$count       # 计数部分
      results[[sp]] <- data.frame(
        species = sp,
        zero_pct = zero_pct,
        coef = coefs[2, 1],              # Group系数
        pvalue = coefs[2, 4],            # Group p值
        model = ifelse(zero_pct > 0.3, "ZINB", "NB")
      )
    }, error = function(e) {
      cat("  跳过", sp, ":", e$message, "\n")
    })
  }
  result_df <- do.call(rbind, results)   # 合并结果
  result_df$padj <- p.adjust(result_df$pvalue, method = "BH")  # FDR校正
  return(result_df[order(result_df$padj), ])  # 按校正p值排序
}

# ZINQ-L: 零膨胀分位数回归（纵向版本，2025年）
# 适用于：纵向微生物组数据中零膨胀物种的差异分析

# 安装
# devtools::install_github("xxx/ZINQ-L")  # 从GitHub安装

# ZINQ-L优势：
# 1. 不需要假设特定分布（分位数回归是非参数的）
# 2. 处理零膨胀
# 3. 支持纵向/重复测量数据
# 4. 对异常值更稳健

# 使用示例（伪代码，实际包名以发表为准）
# library(ZINQL)
# result <- zinql(
#   y = abundance_vector,       # 物种丰度
#   X = covariates_matrix,      # 协变量矩阵
#   Z = zero_covariates,        # 零膨胀部分协变量
#   id = subject_id,            # 个体ID（纵向数据）
#   tau = 0.5                   # 分位数（0.5=中位数）
# )

# Python零膨胀模型
import numpy as np   # 数值计算
import pandas as pd  # 数据处理
import statsmodels.api as sm  # 统计模型
from statsmodels.discrete.count_model import (
    ZeroInflatedPoisson,       # ZIP模型
    ZeroInflatedNegativeBinomialP  # ZINB模型
)

# 1. 诊断零膨胀
def diagnose_zero_inflation(counts, plot=True):
    """诊断数据是否需要零膨胀模型"""
    zero_pct = (counts == 0).mean()      # 零值比例
    mean_val = counts[counts > 0].mean() # 非零均值
    var_val = counts[counts > 0].var()   # 非零方差
    dispersion = var_val / mean_val      # 离散指数

    print(f"零值比例: {zero_pct:.1%}")
    print(f"非零均值: {mean_val:.2f}")
    print(f"非零方差: {var_val:.2f}")
    print(f"离散指数(方差/均值): {dispersion:.2f}")

    if dispersion > 2:
        print("→ 过度离散: 需要负二项分布")
    if zero_pct > 0.3:
        print("→ 零膨胀: 需要ZI模型")

    # 与泊松分布的零值比较
    from scipy.stats import poisson
    expected_zeros = poisson.pmf(0, counts.mean())  # 泊松预期零值
    print(f"泊松预期零值: {expected_zeros:.1%}")
    print(f"实际零值: {zero_pct:.1%}")
    print(f"超额零值: {zero_pct - expected_zeros:.1%}")

    if plot:
        import matplotlib.pyplot as plt
        fig, axes = plt.subplots(1, 2, figsize=(12, 4))
        # 直方图
        axes[0].hist(counts, bins=50, edgecolor='black', alpha=0.7)
        axes[0].set_title(f"分布 (零值={zero_pct:.0%})")
        axes[0].set_xlabel("计数值")
        # 非零值分布
        axes[1].hist(counts[counts > 0], bins=30, edgecolor='black',
                     alpha=0.7, color='orange')
        axes[1].set_title("非零值分布")
        axes[1].set_xlabel("计数值")
        plt.tight_layout()
        plt.savefig("zero_inflation_diagnosis.png", dpi=150)

# 2. 拟合ZINB模型
def fit_zinb(y, X, Z=None):
    """拟合零膨胀负二项模型"""
    X = sm.add_constant(X)               # 添加截距
    if Z is None:
        Z = X.copy()                     # 零膨胀部分用相同变量
    else:
        Z = sm.add_constant(Z)

    model = ZeroInflatedNegativeBinomialP(
        y, X,                            # 因变量和自变量
        exog_infl=Z,                     # 零膨胀部分的变量
        inflation='logit'                # 零膨胀用logistic
    )
    result = model.fit(disp=False)       # 拟合模型
    print(result.summary())             # 打印结果

    # 提取关键信息
    print("\n== 计数部分（物种存在时的丰度差异）==")
    for name, coef, pval in zip(X.columns, result.params[:len(X.columns)],
                                 result.pvalues[:len(X.columns)]):
        sig = "***" if pval < 0.001 else "**" if pval < 0.01 else \
              "*" if pval < 0.05 else ""
        print(f"  {name}: coef={coef:.4f}, p={pval:.4f} {sig}")

    print("\n== 零膨胀部分（结构零的概率差异）==")
    inflate_params = result.params[len(X.columns):]
    inflate_pvals = result.pvalues[len(X.columns):]
    for name, coef, pval in zip(Z.columns, inflate_params, inflate_pvals):
        sig = "***" if pval < 0.001 else "**" if pval < 0.01 else \
              "*" if pval < 0.05 else ""
        print(f"  {name}: coef={coef:.4f}, p={pval:.4f} {sig}")

    return result

# 3. 使用示例
# 模拟数据
np.random.seed(42)
n = 200                                  # 样本数
group = np.random.binomial(1, 0.5, n)   # 随机分组
age = np.random.normal(50, 10, n)        # 年龄

# 模拟零膨胀负二项数据
from scipy.stats import nbinom
pi = 1 / (1 + np.exp(-(0.5 - 1.0 * group)))  # 结构零概率
mu = np.exp(2.0 + 0.5 * group - 0.01 * age)  # 期望丰度
r = 2                                    # 离散参数
y = np.where(np.random.binomial(1, pi, n) == 1,
             0,  # 结构零
             nbinom.rvs(r, r / (r + mu)))  # NB分布

# 拟合
X = pd.DataFrame({"Group": group, "Age": age})
result = fit_zinb(y, X)

# 零膨胀模型选择
零值<30%: 普通NB或泊松
零值30-90%: ZINB或hurdle
零值>90%: 考虑转为presence/absence二分类

# ZINB vs Hurdle
ZINB: 区分结构零和抽样零，生物学解释更强
Hurdle: 不区分零的来源，数学上更简单
通常选择: AIC/BIC比较决定

# R包选择
pscl: zeroinfl/hurdle（最经典）
glmmTMB: 支持随机效应（纵向数据）
brms: 贝叶斯ZINB（最灵活）
DHARMa: 模型诊断（必用）

# 诊断步骤
1. 看零值比例（>30%？）
2. 看离散指数（方差/均值>2？）
3. 比较泊松预期零值vs实际零值
4. 拟合多个模型，AIC比较
5. DHARMa残差诊断

# 结果解读
计数部分系数: 物种"存在"时的丰度变化方向和大小
零膨胀部分系数: 物种"不存在"概率的变化方向和大小
→ 两部分可能方向不同！（疾病组丰度高但出现率低）