零膨胀计数模型（zero inflated count model） - ricket-sjtu/bi028 GitHub Wiki

对于常规的计数资料（count data），我们一般会采用泊松（Poisson）或负二项（negative binomial）分布进行拟合。但在某些特殊情境下，会出现“零”出现的概率被低估，造成泊松或者负二项分布都无法正确预测的情况。

早在1960年代，就有学者注意到了这种所谓的“零膨胀（zero-inflated）”现象（Johnson & Kotz, 1969）。然而，直到1986年，才有学者提出了解决针对经济学领域的“零膨胀”现象的模型——Hurdle模型（Murllahy, 1986）。1992年，Lambert对于零和非零计数值，提出了混合分布模型——零膨胀泊松分布模型（zero-inflated poisson, ZIP），同时纳入协变量，应用于电子制造业中的质量控制（Lambert, 1992）。在此基础上，Greene将ZIP模型的思想应用到负二项分布，提出了零膨胀负二项模型（zero-inflated negative binomial, ZINB），并用BHHH方法估计模型参数的标准误差（Greene, 1994）。近年来，这些模型被广泛应用在高通量组学数据的分析。

零膨胀模型的数学原理

零膨胀模型的原理是将事件发生的计数看作两个不同的过程：

  1. 事件发生次数只能取零的贝努利（Bernoulli）过程；
  2. 事件发生次数可能为零或者正整数的泊松（Poisson）或者负二项(Negative Binomial)过程。

也就是说，将模型泛化为一个混合模型——事件发生与否的logit模型以及事件发生次数的Poisson或Negative-Binomial模型。零膨胀计数模型实际上就是一种针对零较多且符合等离散（epi-dispersed）的泊松分布模型或过离散（over-dispersed）的负二项分布模型的数据进行的复合计数模型。

零膨胀计数模型的复合概率分布

其中$p_i$表示贝努利过程中事件发生的概率，也就是“过多零（extra zero）”发生的概率；而$g(y_i)$表示观察来源于第二个过程，服从泊松或者负二项分布。观察值中的零一部分来自那些始终不可能发生事件的个体，概率为$p_i$；而另一部分来自于泊松分布或者负二项分布理论的计数值0，概率为$1-p_i$。因此，$Y=y_i$的概率密度为

若$p_i$的取值受个体自身协变量的影响，也就是 $$ p_i = F(w_i^T \gamma) $$ 这里的$F(\cdot)$称为零膨胀连接函数（zero-inflated link function），可以选择用probit或者logit函数：

  • logit连接函数：
  • probit连接函数：

这里$w_i \in \mathbb{R}^q$为零膨胀自变量向量，$\gamma \in \mathbb{R}^q$为零膨胀参数。

零膨胀泊松模型（zero-inflated Poisson model, ZIP）

零膨胀泊松模型中，

  • ZIP模型中的期望值为和方差为
  • $x_i$与$w_i$可以一致，也可以不同。当$w_i$仅包含常数项时，ZIP模型比Poisson模型需要多估计一个参数；如果$w_i = x_i$，则ZIP模型相比Poisson模型需要多估计一倍的参数。

零膨胀负二项（zero-inflated negative binomial, ZINB）模型

对于负二项分布模型，

  • ZINB模型中$y_i$的期望值$E(y_i) = \mu(1-p)$，方差为$\operatorname{Var}(y_i) = E(y_i)(1+\mu(p+\alpha))$
  • 当$\alpha=0$时，ZINB模型等同于ZIP模型。

如何选择模型

  1. 观察计数资料是否存在零膨胀（zero-inflation）现象；
  2. 观察计数资料是否存在过度分散（overdispersion）的情况。观察计数资料的均值与方差是否相等，alpha检验是否显著。若两者基本相等，且alpha检验不显著，则为等分散，宜采用零膨胀泊松模型；若方差明显大于均值，且alpha检验显著，则属于过度分散的情况，宜采用零膨胀负二项模型；
  3. 以Vuong检验决定模型的选择，并用图形比较观测数据的分布与截距回归、泊松模型、负二项模型、零膨胀负二项模型拟合曲线的差异。

R语言分析

  • R语言的glm(){stats}可进行泊松回归分析（family=poisson）
  • R语言的glm.nb(){MASS}可进行负二项回归分析；
  • R语言中的zeroinfl(){pscl}可进行零膨胀泊松回归分析（dist="poisson"）或零膨胀负二项分析（dist="negbin"）；
  • 比较嵌套模型（nested models），可用似然比检验（likelihood-ratio test）lrtest(){lmtest}进行；
  • 用vuong(){pscl}可进行相同数据集非嵌套模型的比较分析；

嵌套模型的比较分析

可用似然比检验进行（Likelihood ratio test, LRT）:

  • 比较两个嵌套模型，需要分别计算两个模型的对数似然值（logLik(){stats}）；
  • 两模型对数似然值的差符合自由度为两模型参数个数之差的卡方分布（$\Chi^2$ distribution）；
  • 这样就算可以计算p值；
  • lrtest(){lmtest}或者``

也可用Wald检验进行（Wald test）

  • Wald检验只能用于嵌套模型，默认将模型与平凡模型（trivial model，仅包含截距的模型）进行比较；
  • Wald检验利用的是卡方检验或者F检验：一般来说，大样本采用卡方检验；而对小样本数据一般采用F检验；
  • 一般来说，卡方检验的统计量$\chi^2$在数值上是F检验统计量$F$的$k$倍，其中$k$为两模型参数个数的差异；
  • 卡方检验的的自由度是两模型的参数个数差异，$k$；
  • F检验的自由度是$(k, df1)$，其中$df1$是简单模型的参数个数；
  • R语言的waldtest(..., test=c("F", "Chisq")){lmtest}函数；

非嵌套模型比较

如果两个非嵌套模型$M_1, M_2$进行比较，

  • 模型1的最大似然估计的概率$p_i = \hat{Pr}(y_i | M_1)$（predprob(){pscl}）
  • 模型2的最大似然估计的概率$q_i = \hat{Pr}(y_i | M_2)$（predprob(){pscl}）
  • 计算$m_i = \log(p_i / q_i)$
  • Vuong统计值为$z = \sqrt{N} \bar{m} / s_m$，其中$\bar{m} = \sum m_i / N$，$s_m$为$m_i$的standard deviation。
  • 当两个模型没有显著差异时，$v$逼近标准正态分布；当统计值显著大于0时，模型1优于模型2；如统计值显著小于0，则模型2优于模型1；
  • Vuong检验可用AIC或BIC进行校正，同时考虑模型的复杂度。

这是一个检验的实例：

data("bioChemists")
## compare Poisson GLM and ZIP
glm1 <- glm(art ~ ., data = bioChemists, family = poisson)
zip <- zeroinfl(art ~ . | ., data = bioChemists, EM = TRUE)
vuong(glm1, zip)

这是返回结果：

Vuong Non-Nested Hypothesis Test-Statistic: 
(test-statistic is asymptotically distributed N(0,1) under the
 null that the models are indistinguishible)
-------------------------------------------------------------
              Vuong z-statistic             H_A    p-value
Raw                   -4.180398 model2 > model1  1.455e-05
AIC-corrected         -3.638468 model2 > model1 0.00013713
BIC-corrected         -2.332709 model2 > model1 0.00983172

这里的p-value使用pnorm(z, 0, 1)计算得到。

The zero-inflated Poisson model predicting *** from , *** and *** was statistically significant (chi-square=, df=, p-val=*). The Vuong test suggests that our zero-inflated model is a significant improvement over standard Poisson one.

⚠️ **GitHub.com Fallback** ⚠️