在第11讲中，我们介绍了响应变量是二分类变量而采取的逻辑回归模型，以及多分类情况下的广义逻辑回归模型。然而，除了考虑二分类的响应变量，我们还有一类重要的变量——计数变量。这类变量的特点是，其取值是离散的，但是又具有次序大小之分，同时我们研究的目标不是概率而是响应变量的具体计数，尤其是一段时间内事件发生的平均次数。要处理计数变量，建立回归模型，首先我们需要解析泊松分布。

一、泊松回归的基本模型

1.1 计数随机变量与泊松分布

计数随机变量，是指随机变量的取值是离散的可数整数，即X=0,1,2,⋯X=0,1,2,\cdots。要描述计数随机变量的分布情况，就需要使用泊松分布。泊松分布是一个含有参数λ\lambda的分布：

P(Y=y)=λyy!e−λ≡f(y;λ)P(Y=y)=\frac{\lambda^y}{y!}e^{-\lambda}\equiv f(y;\lambda)

这里的λ\lambda就表示事件发生的平均次数。

泊松分布也属于指数分布族，因为如果令θ=log⁡λ\theta=\log\lambda，那么还是可以写成如下形式：

exp[yθ−eθ−log⁡y!]exp\left[y\theta-e^\theta-\log y!\right]

不难验证：

EY=∂b(θ)∂θ=eθ=λVar Y=ϕ∂b2(θ)∂θ2=eθ=λ\begin{aligned} EY&=\frac{\partial b(\theta)}{\partial\theta}=e^\theta=\lambda\\ Var\ Y&=\phi\frac{\partial b^2(\theta)}{\partial\theta^2}=e^\theta=\lambda \end{aligned}

这与泊松分布的性质相同：期望和方差相等且都等于参数λ\lambda。

1.2 泊松回归及其拟合

确定好响应变量的分布之后，接下来应该选择联系函数。再次回顾，广义线性模型是要找这样的联系函数：

g(EY)=Xβg(EY)=X\beta

由于响应变量服从泊松分布，所以EY=λEY=\lambda，也就是说，我们其实要找λ\lambda和线性预测子之间的联系函数。首先，肯定不能直接λ=Xβ\lambda=X\beta，这是因为λ\lambda代表的是事件发生的平均次数，一定是正值，而线性预测子有正有负。因此，我们需要找一个这样的联系函数：

  • 在定义域[0,+∞)[0,+\infty)上值域为R\mathbb R
  • 线性单增函数
  • 最好能和θ\theta产生关联（逻辑回归就是这样，性质很好）

那不难想到，我们可以取g(λ)=log⁡λg(\lambda)=\log\lambda作为联系函数，也就是对响应变量取自然对数，此时有：

λ=g−1(Xβ)=eXβ\lambda=g^{-1}(X\beta)=e^{X\beta}

如果把前后联系起来，那么θ\theta可以改写为：

θ=log⁡λ=g−1(Xβ)=Xβ\theta=\log\lambda=g^{-1}(X\beta)=X\beta

也就是说，θ\theta就等价于线性预测子。

但是，这样做还不够完整。泊松回归是根据事件发生的平均次数来建模的，这是它最大的特色，但也是它最大的问题——没有考虑样本总量/观测时间。极端情况下，假设两个样本的自变量的取值完全相同，但是两组样本的样本量/观测时间不同，结果发现响应变量（事件发生的次数）不同，这是不合理的。所以，泊松回归需要考虑样本量/观测时间长短，我们需要引入nin_i：

λini=exiTβ\frac{\lambda_i}{n_i}=e^{x_i^T\beta}

或者写成：

λi=niexiTβ\lambda_i=n_ie^{x_i^T\beta}

这里的nin_i可以是样本量，那么λi/ni\lambda_i/n_i代表人均发生次数；nin_i还可以是观测时间，例如”天“，那么λi/ni\lambda_i/n_i代表每天发生次数。

进一步地，θ\theta参数就可以改写为：

θi=log⁡λi=log⁡ni+xiTβ\theta_i=\log\lambda_i=\log n_i+x_i^T\beta

这里的log⁡ni\log n_i被称为偏移量(offset)。偏移量相当于线性预测子中加入一个已知常数项。

接下来就是拟合泊松回归了，我们可以计算得分统计量：

Uj=∑i=1n[yi−μivar(yi)xij∂μi∂(xiTβ)]=∑i=1n[yi−λiλixijniexiTβ]=∑i=1nxij(yi−λi)\begin{aligned} U_j&=\sum_{i=1}^n\left[\frac{y_i-\mu_i}{var(y_i)}x_{ij}\frac{\partial\mu_i}{\partial(x_i^T\beta)}\right]\\ &=\sum_{i=1}^n\left[\frac{y_i-\lambda_i}{\lambda_i}x_{ij}n_ie^{x_i^T\beta}\right]\\ &=\sum_{i=1}^nx_{ij}(y_i-\lambda_i) \end{aligned}

令Uj=0U_j=0，可以发现每一个yiy_i就是λi\lambda_i的极大似然估计，也就是说λ^i=yi\hat\lambda_i=y_i。

同理，我们可以计算信息阵：

Ijk=∑i=1n[xijxikvar(yi)(∂μi∂(xiTβ))2]=∑i=1n[xijxikλi(niexiTβ)2]=∑i=1nxijxikλi\begin{aligned} \mathfrak I_{jk}&=\sum_{i=1}^n\left[\frac{x_{ij}x_{ik}}{var(y_i)}\left(\frac{\partial\mu_i}{\partial(x_i^T\beta)}\right)^2\right]\\ &=\sum_{i=1}^n\left[\frac{x_{ij}x_{ik}}{\lambda_i}(n_ie^{x_i^T\beta})^2\right]\\ &=\sum_{i=1}^nx_{ij}x_{ik}\lambda_i \end{aligned}

之后就可以使用Fisher Scoring算法了。或者，我们直接使用IRWLS算法，先求出工作变量：

zi=log⁡λi+yi−λiλi,λi=niexiTβ(m)z_i=\log \lambda_i+\frac{y_i-\lambda_i}{\lambda_i},\lambda_i=n_ie^{x_i^T\beta^{(m)}}

然后求权重矩阵：

W=diag(λi)W=diag(\lambda_i)

之后迭代更新回归系数：

β(m+1)=(XTWX)−XTWZ\beta^{(m+1)}=(X^TWX)^-X^TWZ

其中Z=(z1,⋯ ,zn)TZ=(z_1,\cdots,z_n)^T。

1.3 泊松回归模型检验与诊断的基本结论

前面我们在学习逻辑回归时已经学习过如何对模型进行检验和诊断，这一系列方法可以迁移到泊松回归模型当中。这里我们不再展示推导过程，直接给出迁移的结论。

首先给出泊松回归的Deviance统计量，用于检验模型结构：

D=2∑i=1N[yilog⁡yiλ^i−(yi−λ^i)]D=2\sum_{i=1}^N\left[y_i\log\frac{y_i}{\hat\lambda_i}-(y_i-\hat\lambda_i)\right]

接着计算Deviance残差：

di=sign(yi−λ^i)Dd_i=sign(y_i-\hat\lambda_i)\sqrt D

其中sign()sign()是符号函数。标准化的Deviance残差为：

rD(i)=di1−hiir_D^{(i)}=\frac{d_i}{\sqrt{1-h_{ii}}}

其中hii=X(XTX)−1XTh_{ii}=X(X^TX)^{-1}X^T是帽子矩阵的对角元。

泊松回归也有Pearson残差，即：

ri=yi−λ^iλ^ir_i=\frac{y_i-\hat\lambda_i}{\sqrt{\hat\lambda_i}}

标准化的Pearson残差为：

rP(i)=ri1−hiir_P^{(i)}=\frac{r_i}{\sqrt{1-h_{ii}}}

上述两种残差都可以用来绘制残差图。

拟合优度检验需要构造Pearson卡方统计量，也就是：

∑i=1nri2=∑i=1n(yi−λ^i)2λ^i∼χ2\sum_{i=1}^nr_i^2=\sum_{i=1}^n\frac{(y_i-\hat\lambda_i)^2}{\hat\lambda_i}\sim\chi^2

泊松回归无法使用HL检验来评价拟合优度，但可以使用伪决定系数，例如McFadden伪决定系数，计算方法和逻辑回归是一致的，这里不再赘述。

二、泊松回归的过度分散与解决方案

2.1 拟泊松分布

我们已经知道，要考察模型是否具有过度分散的问题，只需要计算

ϕ=Pearson χ2n−p\phi=\frac{Pearson\ \chi^2}{n-p}

如果大于1就是出现过度分散，小于1就是欠分散。

在逻辑回归中，如果出现过度分散，我们就会放宽方差，即Var(yi)=ϕπi(1−πi)Var(y_i)=\phi \pi_i(1-\pi_i)，相当于使用拟二项分布(quasi-binomial)。那么同理，泊松回归也可以放宽方差，令Var(yi)=ϕλiVar(y_i)=\phi\lambda_i，也就是让方差稍大于期望，这样就相当于使用拟泊松分布(quasi-poisson)。

使用quasi型的分布确实简单易懂，但是这并不是一个真正的分布模型，AIC值是不可用的。对于泊松回归来说，我们还有更好的办法。

2.2 负二项回归模型

解决过度分散问题的基本原则不变，即放松对方差的限制，允许Var(yi)>E(yi)Var(y_i)>E(y_i)。quasi-poisson做到了，但是很遗憾不是一个分布模型。我们希望能找到一个分布模型，并且满足方差大于期望，这个模型就是：负二项回归模型。

要谈论负二项回归，首先要讨论负二项分布。负二项分布(Negative Binomial, NB)描述了在一项重复试验中，直到第xx次试验时事件才成功发生rr次的概率分布。设事件成功发生一次的概率为π\pi，那么负二项分布可以表示为前x−1x-1次试验成功r−1r-1次且第xx次试验成功，即：

Pr⁡(X=x)=f(x,r,π)=Cx−1r−1πr−1(1−π)x−r×π=Cx−1r−1πr(1−π)x−r\Pr(X=x)=f(x,r,\pi)=C_{x-1}^{r-1}\pi^{r-1}(1-\pi)^{x-r}\times\pi=C_{x-1}^{r-1}\pi^r(1-\pi)^{x-r}

负二项分布的期望为：

EX=∑x=r∞xCx−1r−1πr(1−π)x−r=rπ∑x=r∞Cxrπr+1(1−π)x−r=rπ\begin{aligned} EX&=\sum_{x=r}^\infty xC_{x-1}^{r-1}\pi^r(1-\pi)^{x-r}\\ &=\frac{r}{\pi}\sum_{x=r}^\infty C_x^r\pi^{r+1}(1-\pi)^{x-r}\\ &=\frac{r}{\pi} \end{aligned}

再根据：

EX2=∑x=r∞x2Cx−1r−1πr(1−π)x−r=rπ∑x=r∞[(x+1)Cxrπr+1(1−π)x−r−Cxrπr+1(1−π)x−r]=rπ[r+1π−1]=r(r−π+1)π2\begin{aligned} EX^2&=\sum_{x=r}^\infty x^2C_{x-1}^{r-1}\pi^r(1-\pi)^{x-r}\\ &=\frac{r}{\pi}\sum_{x=r}^\infty \left[(x+1)C_x^r\pi^{r+1}(1-\pi)^{x-r}-C_x^r\pi^{r+1}(1-\pi)^{x-r}\right]\\ &=\frac{r}{\pi}\left[\frac{r+1}{\pi}-1\right]\\ &=\frac{r(r-\pi+1)}{\pi^2} \end{aligned}

得到负二项分布的方差为：

Var(X)=EX2−(EX)2=r(1−π)π2Var(X)=EX^2-(EX)^2=\frac{r(1-\pi)}{\pi^2}

不过，在广义线性回归中，我们不这样直接表达负二项分布，而是定义第rr次试验成功前的“失败”次数，即：

Y=X−rY=X-r

此时负二项分布可以表达为：

Pr⁡(Y=y)=f(y,r,π)=Cy+r−1r−1πr(1−π)y\Pr(Y=y)=f(y,r,\pi)=C_{y+r-1}^{r-1}\pi^r(1-\pi)^y

这样做相当于把负二项分布进行了平移，YY仍然服从负二项分布。那么接下来的问题是，为什么要这样去表达负二项分布？有两个原因。

第一，均值和方差的表达很自然。可以计算得到：

EY=E(X−r)=rπ−r=r(1−π)π=μVar(Y)=Var(X−r)=Var(X)=r(1−π)π2=μ+αμ2\begin{aligned} &EY=E(X-r)=\frac{r}{\pi}-r=\frac{r(1-\pi)}{\pi}=\mu\\ &Var(Y)=Var(X-r)=Var(X)=\frac{r(1-\pi)}{\pi^2}=\mu+\alpha\mu^2 \end{aligned}

其中α=1/r\alpha=1/r，π=1/(1+αμ)\pi=1/(1+\alpha\mu)。

这样处理之后可以发现负二项分布天然具有Var(Y)>EYVar(Y)>EY的属性，多出来的αμ2\alpha\mu^2项可以解决过度分散的问题。

第二，负二项分布和泊松分布之间有紧密联系。事实上，负二项分布可以看作是一个复合分布，即负二项分布=泊松分布+Gamma分布。泊松分布是一个关于参数λ\lambda的分布，是一个条件分布，它假定参数λ\lambda是一个未知但确定的值：

Y∣λ∼Poisson(λ)Y|\lambda\sim Poisson(\lambda)

然而，λ\lambda本身是有可能存在波动的，可以认为它服从Gamma分布：

λ∼Gamma(r,π1−π)\lambda\sim Gamma(r,\frac{\pi}{1-\pi})

所以通过贝叶斯公式，我们就可以求得关于YY的边缘分布，就是负二项分布。

如果是单纯的泊松分布，那么EY=Var(Y)=μEY=Var(Y)=\mu；而负二项分布中方差多出来的αμ2\alpha\mu^2就是由Gamma分布（也就是λ\lambda的随机性）产生的。当r→∞r\rightarrow\infty时，波动αμ2\alpha\mu^2近乎消失，就退化成泊松分布了。总结来说，负二项分布=泊松分布+λ\lambda波动，负二项分布中的μ\mu等价于泊松分布中的EλE\lambda。

这一讲我们先不展开Gamma分布，只需要了解负二项分布是由泊松分布和Gamma分布复合即可。

有了负二项分布，就可以引入负二项回归了。负二项回归使用的联系函数和泊松回归一样，都是取：

g(μ)=log⁡μg(\mu)=\log\mu

因此，我们不难写出各种统计量了，这里我们直接给出结果：

Uj=∑i=1nyi−μi1+αμixijIjk=∑i=1nμi1+αμixijxikzi=log⁡μi+yi−μiμiwi=μi1+αμi\begin{aligned} &U_j=\sum_{i=1}^n\frac{y_i-\mu_i}{1+\alpha\mu_i}x_{ij}\\ &\mathfrak I_{jk}=\sum_{i=1}^n\frac{\mu_i}{1+\alpha\mu_i}x_{ij}x_{ik}\\ &z_i=\log\mu_i+\frac{y_i-\mu_i}{\mu_i}\\ &w_i=\frac{\mu_i}{1+\alpha\mu_i} \end{aligned}

最后，部分软件会估计α\alpha的值，也就是波动系数。估计方法同样采用极大似然估计，这里不再赘述。

三、零截断与零膨胀

在泊松回归的最后，我们来看一种特殊情况——0的出现。在泊松分布当中是有一定概率出现0的，即：

Pr⁡(Y=0)=f(0,λ)=e−λ\Pr(Y=0)=f(0,\lambda)=e^{-\lambda}

但是真实数据可能不尽如人意。有时候0在数据集中不出现，比如我们要统计病人住院天数，那么没住院的就不在数据集中；有时候0在数据集中大量出现，比如我们要统计保险索赔的次数，这个次数往往都是0，否则保险公司一定是亏的。对于这些情况，使用纯泊松回归就无法进行描述，因为0出现的频率和模型给出的概率完全不同。

因此，我们需要对泊松回归模型进行一定的改进。根据0出现的情况，我们把模型分为两类：零截断模型和零膨胀模型。

3.1 零截断模型

零截断(Zero Truncate, ZT)是指0在数据集中不出现（≠很少出现）。既然0不出现，那么概率分布的所有取值的概率之和就不为1了，这就不是概率分布了。解决办法也很简单，就是把0从概率分布函数中去掉，也就是取条件分布：

f(y∣y>0)=f(y)1−f(0)f(y|y>0)=\frac{f(y)}{1-f(0)}

对于泊松分布来说，由于其概率分布为：

f(y,λ)=exp{ylog⁡λ−λ−log⁡y!}f(y,\lambda)=exp\{y\log\lambda-\lambda-\log y!\}

同时不难求出：

f(0,λ)=e−λf(0,\lambda)=e^{-\lambda}

因此可以求得：

f(y,λ∣y>0)=exp{ylog⁡λ−λ−log⁡y!}1−e−λf(y,\lambda|y>0)=\frac{exp\{y\log\lambda-\lambda-\log y!\}}{1-e^{-\lambda}}

这就是零截断泊松(Zero Truncate Poisson, ZTP)。对应的对数似然函数为：

l(y,λ)=∑i=1n[ylog⁡λ−λ−log⁡y!−log⁡(1−e−λ)]l(y,\lambda)=\sum_{i=1}^n\left[y\log\lambda-\lambda-\log y!-\log(1-e^{-\lambda})\right]

可以看到，零截断泊松比纯泊松多了−log⁡(1−e−λ)-\log(1-e^{-\lambda})项。

ZTP回归仍然取g(λ)=log⁡λg(\lambda)=\log\lambda，但是已经不属于广义线性回归的范畴，只能使用极大似然法或者梯度法完成。

同理，我们在2.2小节学习过负二项分布，该分布也可以转换为零截断模型。负二项分布的分布函数为：

f(y)=exp{ylog⁡(1−π)−rlog⁡1π+log⁡Cy+r−1r−1}=exp{ylog⁡αμ1+αμ−1αlog⁡(1+αμ)+log⁡Cy+r−1r−1}\begin{aligned} f(y)&=exp\left\{y\log(1-\pi)-r\log\frac{1}{\pi}+\log C_{y+r-1}^{r-1}\right\}\\ &=exp\left\{y\log\frac{\alpha\mu}{1+\alpha\mu}-\frac{1}{\alpha}\log(1+\alpha\mu)+\log C_{y+r-1}^{r-1}\right\} \end{aligned}

并且取到0的概率为：

f(0)=(1+αμ)−1/αf(0)=(1+\alpha\mu)^{-1/\alpha}

因此可以求得零截断负二项(Zero Truncate Negative Binomial, ZTNB)为：

f(y∣y>0)=exp{ylog⁡αμ1+αμ−1αlog⁡(1+αμ)+log⁡Cy+r−1r−1}1−(1+αμ)−1/αf(y|y>0)=\frac{exp\left\{y\log\frac{\alpha\mu}{1+\alpha\mu}-\frac{1}{\alpha}\log(1+\alpha\mu)+\log C_{y+r-1}^{r-1}\right\}}{1-(1+\alpha\mu)^{-1/\alpha}}

ZTNB回归同样取g(λ)=log⁡λg(\lambda)=\log\lambda，但同样也不属于GLM的范畴，需要使用极大似然法或梯度法求解。

3.2 零膨胀模型

零膨胀(Zero Inflate, ZI)是指0在数据集中频繁出现。与零截断相反，这个时候原始的概率分布无法充分地描述0的出现，因此我们需要加入别的分布来描述。零膨胀模型的基本思想是，把数据集中的0分为两类：

  1. 第一类，是因为计数产生的0，即计数模型中的0
  2. 第二类，是因为这个事件不发生而产生的“结构化”的0，即二分类模型中的0

也就是说，零膨胀模型认为数据集中的0有两个来源：一个是事件不发生产生的0，另一个是事件发生但计数产生的0。因此，零膨胀模型可以这样写：

f(y)={Binary(y=0,π)+[1−Binary(y=0,π)]Count(y=0,λ)[1−Binary(y=0,π)]Count(y>0,λ)f(y)= \begin{cases} Binary(y=0,\pi)+[1-Binary(y=0,\pi)]Count(y=0,\lambda)\\ [1-Binary(y=0,\pi)]Count(y>0,\lambda) \end{cases}

对于观测值为0的数据，要么通过二分类模型Binary(y=0)Binary(y=0)产生，剩余部分1−Binary(y=0)1-Binary(y=0)通过计数模型Count(y=0)Count(y=0)产生。二分类模型可以是我们之前学过的逻辑回归模型，也可以是probit模型、log-log模型等；计数模型可以是泊松回归模型，当然还可以是负二项回归模型。

如果计数模型是泊松回归，那么就变成了零膨胀泊松(Zero Inflate Poisson, ZIP)。我们设二分类模型的联系函数为g(x)g(x)，反函数为g−1(x)=h(x)g^{-1}(x)=h(x)，定义数据集中0对应的样本集合为D0D_0，那么ZIP模型可以写作：

ZIP:f(y)=f(yi=0)+f(yi>0)f(yi=0)=h(xiTβ)+[1−h(xiTβ)]e−λ,i∈D0f(yi>0)=[1−h(xiTβ)]exp{yilog⁡λi−λi−log⁡yi!},i∉D0\begin{aligned} &ZIP:f(y)=f(y_i=0)+f(y_i>0)\\ &f(y_i=0)=h(x_i^T\beta)+[1-h(x_i^T\beta)]e^{-\lambda},i\in D_0\\ &f(y_i>0)=[1-h(x_i^T\beta)]exp\left\{y_i\log\lambda_i-\lambda_i-\log y_i!\right\},i\notin D_0 \end{aligned}

如果计数模型是负二项回归，那么就变成了零膨胀负二项(Zero Inflate Negative Binomial, ZINB)。同上述假定下，ZINB模型可以写作：

ZINB:f(y)=f(yi=0)+f(yi>0)f(yi=0)=h(xiTβ)+[1−h(xiTβ)](1+αμ)−1/α,i∈D0f(yi>0)=[1−h(xiTβ)]exp{yilog⁡αμ1+αμ−1αlog⁡(1+αμ)+log⁡Cy+r−1r−1},i∉D0\begin{aligned} &ZINB:f(y)=f(y_i=0)+f(y_i>0)\\ &f(y_i=0)=h(x_i^T\beta)+[1-h(x_i^T\beta)](1+\alpha\mu)^{-1/\alpha},i\in D_0\\ &f(y_i>0)=[1-h(x_i^T\beta)]exp\left\{y_i\log\frac{\alpha\mu}{1+\alpha\mu}-\frac{1}{\alpha}\log(1+\alpha\mu)+\log C_{y+r-1}^{r-1}\right\},i\notin D_0 \end{aligned}