负二项回归 | Codeboy's Blog

两个窗口的平均计数都为3，实际波动却可能很不同。一种原因是窗口背后的事件速率本身在变：给定速率时仍是Poisson，跨窗口平均后却不再满足方差等于均值。负二项回归把这种异质性写成完整概率分布，因此除了平均数，也能预测零次概率和尾部。

形式陈述

先固定参数的含义

给定固定协变量 $x_i\in\mathbb R^d$ 、已知暴露量 $e_i>0$ ，设响应 $Y_1,\ldots,Y_n$ 独立，取值为0、1、2、……。本页使用NB2参数化：

\mu_i=e_i\exp(x_i^\mathsf T\beta),\qquad k>0, \qquad p_k(y;\mu_i)=\frac{\Gamma(k+y)}{\Gamma(k)y!} \left(\frac{k}{k+\mu_i}\right)^k \left(\frac{\mu_i}{k+\mu_i}\right)^y.\tag{1}

β是回归系数，k是全体观察共用的形状参数。均值和方差分别为

E(Y_i\mid x_i,e_i)=\mu_i,\qquad \operatorname{Var}(Y_i\mid x_i,e_i)=\mu_i+\frac{\mu_i^2}{k}.\tag{2}

“共用k”只是同一个参数，不表示共享同一次随机抽出的潜在速率。k可以事先固定，也可以与β一起估计。越小的k允许越大的相对异质性；另一常见记号α=1/k给出方差μ+αμ²。NB1的方差与均值成线性比例，是另一种参数约束，不能混用本页得分。

当k固定时，式(1)属于广义线性模型。其自然参数为 $\vartheta=\log\{\mu/(k+\mu)\}<0$ ，配分项为 $-k\log(1-e^\vartheta)$ ；本页采用的log均值链接不是这个自然链接。若再估k，就要额外处理形状参数，不能把一次固定族的GLM拟合称为全部工作。

Gamma–Poisson构造与归一化

令各 $U_i$ 独立，密度为

g_k(u)=\frac{k^k}{\Gamma(k)}u^{k-1}e^{-ku},\qquad u>0.

这是shape为k、rate为k的Gamma分布；scale则是1/k。利用Gamma函数的正实积分与递推，代换t=ku得到积分为1，并得EU_i=1、Var(U_i)=1/k。再令给定全部U_i时，各计数独立且

Y_i\mid U_i\sim\operatorname{Pois}(\mu_iU_i).

对Poisson质量积分，

\int_0^\infty\frac{e^{-\mu u}(\mu u)^y}{y!} \frac{k^ku^{k-1}e^{-ku}}{\Gamma(k)}du =\frac{\mu^yk^k}{y!\Gamma(k)} \frac{\Gamma(k+y)}{(k+\mu)^{k+y}},

恰为式(1)。非负混合的总质量等于被混合分布总质量的平均，所以归一化为1。再用全期望与全方差，

EY=E(\mu U)=\mu,\qquad \operatorname{Var}Y=E(\mu U)+\operatorname{Var}(\mu U)=\mu+\mu^2/k.

这个构造是式(1)的一种生成解释。观察到负二项边缘分布，并不能单凭它证明真实系统恰好存在一个Gamma隐变量。

固定k的得分与信息

令η_i=log e_i+x_iᵀβ。忽略与β无关的项，对数似然为

\ell(\beta;k)=\sum_i\{y_i\eta_i-(y_i+k)\log(k+e^{\eta_i})\}+C(k,y).

逐项求导得到得分和Fisher信息：

U_\beta=\sum_i x_i\frac{k(y_i-\mu_i)}{k+\mu_i},\qquad I_{\beta\beta}=\sum_i\frac{k\mu_i}{k+\mu_i}x_ix_i^\mathsf T.\tag{3}

观测负Hessian中的权重为 $k\mu_i(k+y_i)/(k+\mu_i)^2$ ；将y_i换为其期望μ_i才得到式(3)。这些权重都严格为正，所以固定k、满列秩设计下ℓ对β严格凹。但严格凹不保证极大点在有限处存在：只有截距且所有y_i=0时，令β→−∞，似然一直增加并趋于1。

直觉

Poisson模型把每个窗口的强度固定住；Gamma混合让相同已知条件下的窗口还有未观察到的速率差异。总波动由两部分相加：给定速率仍有Poisson计数噪声，速率本身也会波动。第二部分与μ²成比例，说明暴露量增大时，异质性可以比计数噪声增长得更快。

k大时U集中在1附近，模型接近Poisson。固定μ和每个整数y，利用 $\Gamma(k+y)/\Gamma(k)=\prod_{j=0}^{y-1}(k+j)$ ，式(1)在k→∞时趋于 $e^{-\mu}\mu^y/y!$ 。这是极限关系；任何有限k都仍有正的额外方差。

回归中的e_i控制已知暴露量。若两个对象协变量相同，暴露量加倍会使μ加倍；在同一暴露量下，某协变量增加1而其他条件不变，均值比为对应系数的指数。这个解释针对计数均值，不是零次概率的固定比值。

例子与边界

均值相同，零概率不同

取μ=3、k=2，得到方差15/2、零次概率4/25。递推比为

\frac{p_k(y+1;\mu)}{p_k(y;\mu)} =\frac{k+y}{y+1}\frac{\mu}{k+\mu},

所以P(Y=1)=24/125、P(Y=2)=108/625。均值同为3的Poisson分布零次概率约0.04979，而这里为0.16。多零可以由连续速率异质性出现，不必另设一类“永不发生事件”的对象。

NB2仍不是任意计数分布。它在μ>0、有限k时必定方差大于均值，无法匹配欠离散总体。仅凭某份有限样本方差略低于均值，也不能断言总体必定欠离散；样本矩与总体模型约束要分开。

共享一个U会改变联合模型

设两个窗口的边缘均值均为3、k=2。如果各自抽独立U_i，二者独立，总方差为15。如果共享同一个U，给定U后仍条件独立，却有

\operatorname{Cov}(Y_1,Y_2)=\operatorname{Var}(3U)=9/2,

因此总方差为24。共享情形的联合零概率为 $E(e^{-6U})=(2/8)^2=1/16$ ，独立情形却是 $(4/25)^2=16/625$ 。每一个单窗的NB分布都没有变。

更一般地，共享U时，两个计数的联合质量为

P(Y_1=a,Y_2=b)=\frac{\Gamma(k+a+b)}{\Gamma(k)a!b!} \frac{k^k\mu_1^a\mu_2^b}{(k+\mu_1+\mu_2)^{k+a+b}}.

这是对同一个潜变量积分的结果，不等于两个边缘质量的乘积。重复窗口的分析若忽略这一点，连似然对应的实验都会改变；只需要边际均值而允许簇内相关时，可回到GEE的独立簇合同。

下面的图把两个边缘均值改为3、6，共享倍数仍保持k=2。每个面板使用相同的边缘分布，颜色却记录不同的联合质量；热图外的尾部没有被重新归一化进来。

均值3和6、形状2的两个NB边缘，在独立与共享Gamma机制下形成不同联合质量；红框是联合零事件 — 相同边缘与不同联合

一次固定k更新可以手算

只有截距、暴露量都为1，计数为0、1、5。取k=2，起点μ=1、β=0。式(3)给U=2、I=2，因此一次Fisher scoring取β新=1、μ新=e。真正得分根满足μ=样本均值=2，即β=log2；一次更新不应被当成最终答案。

若只按两个矩设定k，常用 $k=\bar y^2/(v-\bar y)$ ，其中v是约定好分母的样本中心二阶矩。它只在v>ȳ时给正值，而且通常不是联合似然极大点。报告时必须说明这是矩匹配，而不能把它标成MLE。

推论与应用

一轮拟合做哪些运算

固定k、给定当前β，先计算所有μ_i，再形成式(3)，解 $I\Delta=U$ ，尝试β+tΔ。回溯t=1、1/2、……直到实际对数似然不下降；若设计退化，先处理不可识别列，不直接求逆。也可用观测Hessian作Newton步。

在d≥1、稠密设计且指数/对数按数值原语计算时，形成信息矩阵需O(nd²)，分解并解线性方程需O(d³)，额外工作存储O(d²)；逐行读数据时不必另存所有响应。每次回溯检验另需O(nd)。这些是每轮代价，不是达到某精度的迭代次数保证。

若k未知，可以在β更新后优化ξ=log k以保证k>0。整数响应允许把形状得分写成无需另定义特殊导数的有限和：

\frac{\partial\ell}{\partial k} =\sum_i\left\{ \sum_{j=0}^{y_i-1}\frac1{k+j} +\log\frac{k}{k+\mu_i} +\frac{\mu_i-y_i}{k+\mu_i} \right\}.\tag{4}

y_i=0时内层和为空。它来自Gamma递推后对数求导；直接计算内层和总共需O(Σy_i)项。较大计数也可使用经过数值验证的log-Gamma及其导数实现。交替步骤只要求所选目标增加，不能据此声称全局最优、有限k极大点或已达收敛；Poisson极限k→∞也可能是最优边界。

完整概率与均值合同各自能交付什么

NB2给出每个y的概率，因而可以讨论零频数、预测尾部和完整似然。准Poisson的均值加比例方差合同只规定部分矩，不会自动给出式(1)。反过来，若实际目的只需要正确的均值估计和稳健协方差，旧有sandwich推断不要求把负二项质量当成真分布。

计数观测机制终点进一步把独立与共享异质性放到不同暴露量的两窗中，并比较零截断、门槛与零膨胀的观测合同。

参考资料

A. Colin Cameron、Pravin K. Trivedi，Essentials of Count Data Regression，作者稿，1999，§3.1（pp.6–7）Gamma混合、NB1/NB2，§4.1（pp.9–10）均值与准似然。
Achim Zeileis、Christian Kleiber、Simon Jackman，Regression Models for Count Data in R，Journal of Statistical Software 27(8),2008,pp.1–25，§2.1，固定与估计形状参数的区别。

形式陈述 ​

先固定参数的含义 ​

Gamma–Poisson构造与归一化 ​

固定k的得分与信息 ​

均值相同，零概率不同 ​

共享一个U会改变联合模型 ​

一次固定k更新可以手算 ​

一轮拟合做哪些运算 ​

完整概率与均值合同各自能交付什么 ​

反向链接