零膨胀泊松模型的变量选择及应用

Variable Selection and Application of Zero Inflated Poisson Model

摘要: 零膨胀泊松模型是研究零过多的计数数据的方法。在实际应用中，为避免遗漏，常选择过多变量进行分析，因此需要对模型进行变量选择，本文在零膨胀泊松模型的基础上构造模型进行变量选择，本文通过研究影响德国卫生保健需求的因素，对比岭回归、Lasso、自适应Lasso和加权弹性网方法的效果。

Abstract: The zero inflated Poisson model is a method to study zero and excessive counting data. In practical application, in order to avoid omission, many variables are often selected for statistical analysis, so it is necessary to select variables for the model. This paper compares the effects of ridge regression, Lasso, adaptive Lasso and weighted elastic network by studying the factors affecting the demand for health care in Germany.

1. 引言

泊松分布是研究计数数据中最常用的模型，计数但是面对散度过大的数据时，拟合效果往往不够理想，而这类数据的特点是存在大量的零值，我们把这类数据称为零膨胀数据。对这类数据的研究最早可见于20世纪60年代中关于零膨胀现象的探索，1992年Lamber [1] 提出了零膨胀泊松分布模型，并用于电子制造业的质量控制当中；近些年零膨胀模型得到充分的扩展，Ghosh [2] 等研究了零膨胀模型的贝叶斯方法，Fahrmeir [3] 等提出了一类零膨胀可加模型。当数据变量较多时就需要进行变量选择，本文基于惩罚似然的EM算法 [4] [5] ，用岭回归、Lasso、自适应Lasso和加权弹性网对德国卫生保健数据进行分析。

2. 模型假定和方法分析

2.1. 零膨胀泊松模型

零膨胀泊松分布模型的思想是取值为0的部分和取值为泊松分布的部分按一定比例构成的混合分布，由此我们可以得到零膨胀泊松模型的具体表达：

$y_i~\{\begin{array}{ll}0,\hfill & {\phi }_i\hfill \\ \text{Poisson}\left(\lambda \right),\hfill & 1-{\phi }_i\hfill \end{array}$ (2.1)

其中 ${\phi }_i$ 表示0占的比例，g(y)来自泊松分布，且每一部分的概率都在(0, 1)之间。由(2.1)式可知观测值0由两部分构成，这里把来自额外的零的部分叫做结构零，来自离散分布的零的部分记为分布零，则 $Y=y_i$ 的概率密度为：

$P\left(Y=y_i|{\phi }_i\right)=\{\begin{array}{ll}{\phi }_i+\left(1-{\phi }_i\right){\text{e}}^{-{\lambda }_i},\hfill & y_i=0\hfill \\ \left(1-{\phi }_i\right)\frac{{\lambda }_i^{y_i}}{y_i!}{\text{e}}^{-{\lambda }_i},\hfill & y_i\ge 1\hfill \end{array}$ (2.2)

${\phi }_i$ 为结构零的比例，当 ${\phi }_i=0$ 时模型(2.2)退化为泊松分布。为了考虑零膨胀计数数据中响应变量与协变量之间的关系，Lambert对零膨胀参数 ${\phi }_i$ 和泊松分布参数 ${\lambda }_i$ 引入协变量，二者分别用logistic回归和对数线性回归模型建模，由此得到零膨胀泊松回归模型，连接函数如下所示：

$\{\begin{array}{l}\log \left({\lambda }_i\right)=x_i^{\text{T}}\beta \hfill \\ \text{logit}\left({\phi }_i\right)=\log \left(\frac{{\phi }_i}{1-{\phi }_i}\right)=z_i^{\text{T}}\gamma \hfill \end{array}$ (2.3)

其中 $\beta ={\left({\beta }_1,{\beta }_2,\cdots ,{\beta }_p\right)}^{\text{T}}$ 是协变量 $x_i^{\text{T}}$ 的回归系数， $\gamma ={\left({\gamma }_1,{\gamma }_2,\cdots ,{\gamma }_q\right)}^{\text{T}}$ 是协变量 $z_i^{\text{T}}$ 的回归系数，且 $x_i,z_i$ 分别是p维和q维向量。

2.2. ZIP模型似然函数

引入潜在数据 ${\omega }_i$ ，如果 $y_i$ 来自额外零，记 ${\omega }_i=1$ ，否则 ${\omega }_i=0$ ，故 ${\omega }_i$ 有如下分布：

${\omega }_i=\{\begin{array}{ll}1,\hfill & {\phi }_i\hfill \\ 0,\hfill & 1-{\phi }_i\hfill \end{array}$ (2.4)

令 $Y=\left(Y_0,{\omega }_i\right)$ ，其中 $Y_0=\left(y_i,x_i,z_i\right)$ 为观测数据，则Y为完全数据集，取 $\varphi ={\left({\beta }^{\text{T}},{\gamma }^{\text{T}}\right)}^{\text{T}}$ ，则基于完全数据得到的ZIP模型的似然函数为：

$L\left(\varphi |Y\right)={\phi }_i^{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\omega }_i}}{\left(\left(1-{\phi }_i\right){\text{e}}^{-{\lambda }_i}\right)}^{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{I}_{\left(y_i=0\right)}}-{\displaystyle \underset{k=1}{\overset{n}{\sum }}{\omega }_i}}{\left[\left(1-{\phi }_i\right)\frac{{\lambda }_i^{y_i}{\text{e}}^{-{\lambda }_i}}{y_i!}\right]}^{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{I}_{\left(y_i>0\right)}}}$ (2.5)

则对数似然函数为：

$\begin{array}{c}l\left(\varphi |Y,{\omega }_i\right)=\left({\displaystyle \underset{i=1}{\overset{n}{\sum }}{\omega }_i}\right)\log \left({\phi }_i\right)+\left({\displaystyle \underset{i=1}{\overset{n}{\sum }}{I}_{\left(y_i=0\right)}}-{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\omega }_i}\right)\log \left(\left(1-{\phi }_i\right){\text{e}}^{-{\lambda }_i}\right)+\left({\displaystyle \underset{i=1}{\overset{n}{\sum }}{I}_{\left(y_i>0\right)}}\right)\log \left(\left(1-{\phi }_i\right)\frac{{\lambda }_i^{y_i}{\text{e}}^{-{\lambda }_i}}{y_i!}\right)\\ =\underset{i=1}{\overset{n}{{\displaystyle \sum }}}\left[{\omega }_i{z}^{\text{T}}\gamma -\log \left(1+{\text{e}}^{z^{\text{T}}\gamma }\right)\right]+\underset{i=1}{\overset{n}{{\displaystyle \sum }}}\left(1-{\omega }_i\right)\left(y_i{x}^{\text{T}}\beta -{\text{e}}^{x_i^{\text{T}}\beta }\right)-{\displaystyle \underset{i=1}{\overset{n}{\sum }}\log }\left(y_i!\right)\\ =l_{c,1}\left(\gamma |Y,{\omega }_i\right)+l_{c,2}\left(\beta |Y,{\omega }_i\right)\end{array}$ (2.6)

3. 惩罚似然方法

假设线性回归模型为 $y=X\beta +\epsilon$ ，其中 $y={\left(y_1,\cdots ,y_n\right)}^{\text{T}},X={\left(x_1,\cdots ,x_n\right)}^{\text{T}},\mathcal{A}=\left\{j,{\hat{\beta }}_j\ne 0\right\}$ ，参数 $\hat{\beta }$ 有如下表达：

$\hat{\beta }=\arg \underset{\beta }{\min }\left\{{\Vert y-X\beta \Vert }_2^2+\lambda P\left(\beta \right)\right\}$ (3.1)

其中上式 $P\left(\beta \right)$ 控制了模型的复杂性，在变量选择时 $P\left(\beta \right)$ 可以是岭惩罚、Lasso惩罚、或弹性网惩罚， $\lambda$

为非负的调整参数。其中岭惩罚可以定义为 $P\left(\beta \right)={\displaystyle \underset{j=1}{\overset{p}{\sum }}{\beta }_j^2}$ ；Lasso惩罚 [6] 作为选择变量的正则化方法，可以定义为 $P\left(\beta \right)={\Vert \beta \Vert }_1={\displaystyle \underset{j=1}{\overset{p}{\sum }}\left|{\beta }_j\right|}$ 估计，即：

$\hat{\beta }=\arg {\min }_{\beta }{\Vert Y-X\beta \Vert }_2^2+{\lambda }_n{\displaystyle \underset{j=1}{\overset{p}{\sum }}{\hat{\omega }}_j}\left|{\beta }_j\right|$ (3.2)

其中 ${\hat{\omega }}_j$ 表示惩罚项中参数系数的权重。本节从弹性网 [7] 的思想出发，将惩罚项中的权重用于 $L_1$ 惩罚和 $L_2$ 惩罚中，利用EM算法对加权弹性网的对数似然函数获得最优解，记 $l\left(\varphi ;x_i,z_i\right)$ 作为对数似然函数，且 $\varphi =\left(\beta ,\gamma \right)$ ，则相应的零膨胀计数模型的惩罚似然函数为：

$\hat{\varphi }=\arg \underset{\theta }{\min }\left\{-l\left(\theta ;x_i,z_i\right)+{\lambda }_1{\displaystyle \underset{j=1}{\overset{p}{\sum }}\left({\alpha }_1{\hat{\omega }}_{1j}\left|{\beta }_j\right|+\frac{1}{2}\left(1-{\alpha }_1\right){\hat{\omega }}_{1j}{\beta }_j^2\right)}+{\lambda }_2{\displaystyle \underset{j=1}{\overset{p}{\sum }}\left({\alpha }_2{\hat{\omega }}_{2j}\left|{\gamma }_j\right|+\frac{1}{2}\left(1-{\alpha }_2\right){\hat{\omega }}_{2j}{\gamma }_j^2\right)}\right\}$ (3.3)

其中 ${\hat{\omega }}_{1j}=SE\left({\hat{\beta }}_j\right)/{\hat{\beta }}_{ridge,j},{\hat{\omega }}_{2j}=SE\left({\hat{\gamma }}_j\right)/{\hat{\gamma }}_{ridge,j}$ ，参数 ${\hat{\beta }}_{ridge,j},{\hat{\gamma }}_{ridge}$ 都由岭回归估计得到，调整参数 ${\lambda }_1,{\lambda }_2\ge 0$ ， $SE\left(\hat{\beta }\right),SE\left(\hat{\gamma }\right)$ 是相关系数的标准误。

3.1. EM算法

令 $W_i=I\left(y_i=0\right)$ ，根据零膨胀泊松分布模型的条件期望有：

$P\left(Y_i=y_i,W_i={\omega }_i|x_i,z_i,\beta ,\gamma \right)={\phi }_i^{{\omega }_i}{\left(1-{\phi }_i\right)}^{1-{\omega }_i}{\left(\frac{{\lambda }^{y_i}{\text{e}}^{-{\lambda }_i}}{y_i!}\right)}^{1-{\omega }_i}=\frac{{\left({\text{e}}^{z_i^{\text{T}}\gamma }\right)}^{{\omega }_I}}{1+{\text{e}}^{z_i^{\text{T}}\gamma }}{\left(\frac{{\text{e}}^{y_i{x}_i^{\text{T}}\beta -{\text{e}}^{z_i^{\text{T}}\gamma }}}{y_i!}\right)}^{1-{\omega }_i}$ (3.4)

基于EM算法，首先用无惩罚最大似然估计所得 $\left(\beta ,\gamma \right)$ 作为初始化参数 ${\hat{\varphi }}^0=\left({\beta }^0,{\gamma }^0\right)$ ，首先进行E步，根据完全数据和更新得到参数估计值的条件期望更新 ${\omega }_i^t$ ，假设参数值为 ${\hat{\beta }}^{\left(t\right)},{\hat{\gamma }}^{\left(t\right)}$ ，那么：

${\hat{\omega }}_i^t=E\left({\omega }_i|Y_0,{\hat{\varphi }}^{\left(t\right)}\right)=P\left({\omega }_i|Y_0=y_i,{\hat{\varphi }}^{\left(t\right)}\right)=\{\begin{array}{ll}{\left[1+\exp \left\{-z_i^{\text{T}}{\gamma }^{\left(t\right)}\right\}f\left(y_i;{\hat{\varphi }}^{\left(t\right)}\right)\right]}^{-1},\hfill & y_i=0\hfill \\ 0,\hfill & y_i=1,2,\cdots \hfill \end{array}$ (3.5)

其中 $f\left(y_i;{\hat{\varphi }}^{\left(t\right)}\right)$ 为零膨胀泊松分布，即 $f\left(y_i;{\hat{\varphi }}^{\left(t\right)}\right)={\text{e}}^{-{\lambda }_i}={\text{e}}^{-{\text{e}}^{x_i^{\text{T}}{\beta }^{\left(t\right)}}}$ ；

$\begin{array}{c}Q\left(\varphi |{\hat{\varphi }}^{\left(t\right)}\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}\left[1-E\left({\omega }_i|Y_0,{\hat{\varphi }}^{\left(t\right)}\right)\right]\left(\log \left(f\left(y_i;{\hat{\beta }}^{\left(t\right)},{\hat{\gamma }}^{\left(t\right)}\right)\right)\right)}+E\left({\omega }_i|Y_0,{\hat{\varphi }}^{\left(t\right)}\right)z_i^{\text{T}}{\gamma }^{\left(t\right)}-\log \left(1+{\text{e}}^{z_i^{\text{T}}{\gamma }^{\left(t\right)}}\right)\\ +{\lambda }_1{\displaystyle \underset{j=1}{\overset{p}{\sum }}\left({\alpha }_1{\hat{\omega }}_{1j}\left|{\beta }_j^{\left(t\right)}\right|+\frac{1}{2}\left(1-{\alpha }_1\right){\hat{\omega }}_{1j}{\left({\beta }_j^{\left(t\right)}\right)}^2\right)}+{\lambda }_2{\displaystyle \underset{j=1}{\overset{p}{\sum }}\left({\alpha }_2{\hat{\omega }}_{2j}\left|{\gamma }_j^{\left(t\right)}\right|+\frac{1}{2}\left(1-{\alpha }_2\right){\hat{\omega }}_{2j}{\left({\gamma }_j^{\left(t\right)}\right)}^2\right)}\\ =Q_1\left(\beta |{\hat{\varphi }}^{\left(t\right)}\right)+Q_2\left(\gamma |{\hat{\varphi }}^{\left(t\right)}\right)\end{array}$ (3.6)

M步：对于给定的 ${\hat{\omega }}_i^t$ 将Q函数最小化，即分别最小化 $Q_1\left(\beta |{\hat{\varphi }}^{\left(t\right)}\right)$ 和 $Q_2\left(\gamma |{\hat{\varphi }}^{\left(t\right)}\right)$ ，则：

$\begin{array}{l}{\beta }^{\left(t+1\right)}=\arg \underset{\beta }{\min }{Q}_1\left(\beta |{\hat{\varphi }}^{\left(t\right)}\right),\\ {\gamma }^{\left(t+1\right)}=\arg \underset{\gamma }{\min }{Q}_2\left(\gamma |{\hat{\varphi }}^{\left(t\right)}\right),\end{array}$ (3.7)

其中 ${\hat{\beta }}^{\left(t\right)},{\hat{\gamma }}^{\left(t\right)}$ 分别表示t步迭代时参数 $\beta ,\gamma$ 的估计值；按照上述E步和M步进行迭代最后收敛得到 ${\hat{\beta }}^{\left(n\right)},{\hat{\gamma }}^{\left(n\right)}$ 。

3.2. 调整参数选择

本文根据最小BIC原则确定调整参数 ${\lambda }_1,{\lambda }_2$ ，BIC准则定义如下：

$\text{BIC}=-2l\left(\hat{\theta },{\lambda }_1,{\lambda }_2\right)+\log \left(n\right)k$ (3.8)

其中 $l\left(\hat{\theta },{\lambda }_1,{\lambda }_2\right)$ 为对数似然函数，n为样本个数， $k={\displaystyle \underset{j=1}{\overset{p}{\sum }}I}\left\{{\beta }_j\ne 0\right\}+{\displaystyle \underset{j=1}{\overset{p}{\sum }}I}\left\{{\gamma }_j\ne 0\right\}$ 为零膨胀泊松分布有效参数个数。

4. 实例分析

患者的医疗需求一直是医疗研究中的重要问题之一，本文选择德国卫生需求数据，该用于研究就诊次数与患者状况之间的关系，数据包括响应变量就诊次数数据和13个与患者状况相关数据，共有1812个观测值。本文用岭回归、Lasso、自适应Lasso、和加权弹性网分别对影响就诊次数的相关因素进行变量选择，其中自适应Lasso的权重参数为 $1/\hat{\theta }$ ，加权弹性网的权重参数为 $se\left(\hat{\theta }\right)/\hat{\theta }$ ， $\hat{\theta }$ 是由岭回归估计得到的参数。

Table 1. Coefficient estimation of data model fitting

表1. 数据模型拟合的系数估计

根据表1的结果可以看出加权弹性网相比于其他方法的BIC值更小，说明该方法比岭回归、Lasso和自适应Lasso拟合效果更好；分析加权弹性网的参数估计可知，残疾程度、结婚与否、受教育年限、收入、是否蓝领与就诊次数无关；健康满意度、是否有孩子与就诊次数为零时呈正相关，残疾程度与就诊次数为零时呈负相从泊松部分拟合参数可知，增加就诊次数的影响因素有年龄、是否残疾、是否有公共健康保险和是否有附加保险，其中是否残疾和是否有保险对增加就诊次数的起主要影响，说明身体情况和参保对就诊次数有较强的影响；对就诊次数呈负相关的影响因素为健康满意度、是否结婚、是否自营和是否公务员，其中主要因素为健康满意度、是否自营、是否公务员和是否被雇佣，说明生活和工作越繁忙，就诊次数越少。

5. 结论

本文基于就诊次数数据特征，构建了零膨胀泊松模型，并分别用岭回归、Lasso、自适应Lasso和加权弹性网进行变量选择，得到影响就诊次数的主要变量。但是，目前本文中关注的仅是截面数据，在纵向数据等方面仍有较大的发展空间。