基于两部模型的组合惩罚似然估计方法研究及其应用

Research and Application of Combined Penalty Likelihood Estimation Method Based on Two-Part Model

摘要: 在统计学中，多借助零膨胀模型研究零膨胀数据潜在的模型结构及变量选择问题。然而，在多数情况下，响应变量的非零部分为定量数据，简单的零膨胀模型无法刻画这类数据的模型结构，对应的参数估计方法也不再适用。鉴于此，学者提出处理零膨胀半连续数据的两部模型。本文将组合惩罚似然估计方法引入两部模型，研究其变量选择问题。提出一种新的处理高维统计分析问题的惩罚似然估计方法：NCPM (New Combined Punishment Method)，并将该方法应用于太原市降水量数据，分析其影响因素。模拟及实例分析结果均表明本文的方法行之有效，较传统的惩罚似然估计方法具有更高的预测精度。

Abstract: In statistics, the potential model structure and variable selection problems of zero expansion data are often studied by means of zero expansion model. However, in most cases, the non-zero part of the response variable is quantitative data. A simple zero expansion model cannot describe the model structure of such data, and the corresponding parameter estimation method is no longer applicable. In view of this, scholars proposed a two-part model to deal with zero expansion semi-continuous data. In this paper, the combined penalty likelihood estimation method is introduced into the two-part model to study the problem of variable selection. A new penalty likelihood estimation method, NCPM (New Combined Punishment Method), is proposed to deal with high-dimensional statistical analysis problems. The method is applied to Taiyuan precipitation data and its influencing factors are analyzed. The results of simulation and case analysis show that the proposed method is effective and has higher prediction accuracy than the traditional penalty likelihood estimation method.

1. 引言

太原市由于地理环境的影响，形成了北温带大陆性气候，四季分明、冬无严寒、夏无酷暑。但其昼夜温差较大，年际气候变化较大，季风环流交替明显，气象灾害频发。尤其近几年，随着城市化进程的不断加快，公共绿地面积减少，道路等不透水面积不断增大，使得渗透力下降，排水压力加重。2005年8月16日，太原市最大降水量达到28.2 mm/h，导致道路积水严重，交通瘫痪数小时；2016年7月连日暴雨袭击给太原市带来了巨大损失：城市道路淹没，造成交通瘫痪；街边店铺、停泊车辆等不同程度的灌入雨水，给人民带来了不可逆转的财产损失；甚至对市民生命安全产生了威胁。研究降水量的影响因素，不仅可以为农作物的播种、培育提供便利，而且可以给气象部门提前预警，以便人们在暴雨来临之前能够做出有效的预防措施，避免造成人员损伤、财务损失等。

降水量数据是典型的零膨胀半连续数据。这类数据集中“零”值占很大比例，数据中的非零部分服从某一连续分布，所以使用传统的数据模型不能很好地解释这类数据。针对这种特殊的数据类型，专家学者提出一种行之有效的方法——两部模型 [1]。两部模型的第一部分用来判断响应变量是否为零，第二部分则用于描述非零响应变量的分布。该模型可以更合理、准确地研究此类特殊数据的内部规律，并在数据预测、检验等方面有着非常重要的作用。针对两部模型，学者们提出多种估计方法，包括：极大似然估计(MLE)、拟似然估计 [2] (McCulloch and Searle, 2001)，惩罚似然估计 [3] (Yau and Lee, 2001)和贝叶斯估计 [4] (Ghosh et al., 2006)等。其中，惩罚似然估计方法使用最为广泛，学者提出多种惩罚似然估计方法，如：Bridge [5] (Frank和Friedman, 1993)，MCP [6] (Minimax Concave Penalty)等，已广泛应用到多种领域，并取得良好效果。但当 $p\gg n$ 或者解释变量之间存在较强的相关性时，上述方法的性能会有一定的局限性。为此，Zou [7] 等提出了弹性网(Elastic net, Enet)，该方法将Lasso和 $L_2$ 混合形成一种新的惩罚，改善了变量间相关性对预测性能造成的影响；Wang [8] 等将SCAD (Smoothly Clipped Absolute Deviation)与 $L_2$ 组合，构成组合惩罚(Combined Penalization, CP)。但针对后来提出的MCP函数的组合惩罚问题，目前国内并没有具体阐述。本文将该方法引入两部模型对其进行具体说明，并将其应用于太原市降水量数据中，分析其影响因素。

2. 模型及研究方法

2.1. 模型建立

假设Y表示事件，根据两部模型的基本思想，Y的分布可表示为如下形式：

$f\left(y\right)=\pi \ast {\mathbb{I}}_{\left(y=0\right)}+\left[\left(1-\pi \right)\ast f\left(y\right)\right]{\mathbb{I}}_{\left(y>0\right)},y\ge 0,0\le \pi \le 1$ (1)

式中 ${\mathbb{I}}_{\left(\mathcal{A}\right)}$ 为示性函数， $\pi$ 为零事件发生的概率， $f\left(y\right)$ 为非零事件服从的分布。描述零事件发生可能性的函数包括：logit、probit、log-log等 [9]。本文选用logit作为模型第一部分的连接函数进行研究，记为：

$\pi =Pr\left(y=0|X\right)=\frac{1}{1+\exp \left({\alpha }_0+{\alpha }^{\text{T}}X\right)}$ (2)

其中， $\alpha ={\left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_p\right)}^{\text{T}}$ 表示未知的回归系数， ${\alpha }_0$ 表示截距项， $X={\left\{x_1,x_2,\cdots ,x_p\right\}}^{\text{T}}$ 表示p维的解释变量。

根据参考文献 [10] 令 $f\left(y\right)$ 为Gamma分布，建立对应的概率密度函数为：

$\begin{array}{l}f\left(y\right)=\frac{{\left(\frac{\mu }{{\sigma }^2}\right)}^{\left(\frac{{\mu }^2}{{\sigma }^2}\right)}}{\Gamma \left(\frac{{\mu }^2}{{\sigma }^2}\right)}{y}^{\frac{{\mu }^2}{{\sigma }^2}-1}{\text{e}}^{-\frac{\mu }{{\sigma }^2}y},\\ \sigma =\exp \left({\beta }_0+{\beta }^{\text{T}}X\right),y>0,\mu >0,{\sigma }^2>0\end{array}$ (3)

其数学期望和方差分别为： $\mu$， ${\sigma }^2$， $\beta ={\left({\beta }_1,{\beta }_2,\cdots ,{\beta }_p\right)}^{\text{T}}$ 表示未知的回归系数， ${\beta }_0$ 表示截距项。

对每个 $i,i=1,2,\cdots ,n$ 有：

${\pi }_i=\frac{1}{1+\exp \left({\alpha }_0+{\alpha }^{\text{T}}{x}_i\right)},{\sigma }_i=\exp \left({\beta }_0+{\beta }^{\text{T}}{x}_i\right)$ (4)

可得对应的对数似然函数为：

$\begin{array}{l}L\left(\theta \right)=\underset{i=1}{\overset{n}{{\displaystyle \sum }}}\{I_{\left(y_i=0\right)}\log \left({\pi }_i\right)+I_{\left(y_i>0\right)}[\log \left(1-{\pi }_i\right)-\log \left(gamma\left(\frac{{\mu }^2}{{\sigma }_i^2}\right)\right)\\ +\frac{{\mu }^2}{{\sigma }_i^2}\ast \log \left(\frac{\mu }{{\sigma }_i^2}\right)+\left(\frac{{\mu }^2}{{\sigma }_i^2}-1\right)\ast \log \left(y_i\right)-y_i\ast \frac{\mu }{{\sigma }_i^2}]\},\theta =\left({\alpha }_0,\alpha ,{\beta }_0,\beta \right)\end{array}$ (5)

2.2. 研究方法

对2.1节所述模型中的未知参数采用基于组合惩罚函数的极大似然方法进行估计，即对目标函数：

$Q\left(\theta \right)=-L\left(\theta \right)+J_{\lambda ,v}\left(\theta \right)$ (6)

$J_{\lambda ,v}\left(\theta \right)=p_{\lambda ,\gamma }\left(\theta \right)+\frac{v}{2}{\theta }^2$ (7)

求最小值。(6，7)式中的 $\lambda ,v$ 称为调整参数， $\gamma$ 为正则化参数。惩罚函数中只考虑系数 $\alpha$ 和 $\beta$。 $p_{\lambda ,\gamma }\left(\theta \right)$ 选取MCP函数，其形式如(8)。对正则化参数 $\gamma$，Breheny and Huang (2011)的模拟中建议 $\gamma =3$，所以本文取 $\gamma =3$，并且尝试令 $\gamma$ 取了几个不同的值，得到的结果基本一致 [11]。一般来说，调整参数的选择方法有很多，包括：AIC，BIC，GCV和CV等。本文利用10次5折交叉验证法确定调整参数 $\lambda$ 的值。另外，还可以选用其他惩罚函数代替MCP。例如Hard、Lasso (Least Absolute Shrinkage and Selection Operator)、Ridge和SCAD等。

MCP函数形式如下：

$\begin{array}{l}{p}_{\lambda ,\gamma }\left(\alpha ,\beta \right)=\underset{j=1}{\overset{p}{{\displaystyle \sum }}}\left(\rho \left(\left|{\alpha }_j\right|;\lambda ,\gamma \right)+\rho \left(\left|{\beta }_j\right|;\lambda ,\gamma \right)\right),\lambda >0,\gamma >0\\ \rho \left(\left|\theta \right|;\lambda ,\gamma \right)=\{\begin{array}{l}\lambda \left|\theta \right|-\frac{{\theta }^2}{2\gamma },\left|\theta \right|\le \lambda \gamma \\ \frac{1}{2}{\lambda }^2\gamma ,\left|\theta \right|>\lambda \gamma \end{array},\lambda \ge 0,\gamma >1,\theta =\left(\alpha ,\beta \right)\end{array}$ (8)

图1展示了几种惩罚函数的阈值函数图。结果表明，MCP、CP和提出的NCPM均具有稀疏性和连续型。组合惩罚函数对于 $\forall v>0$ 都不会产生近似无偏估计量，而当 $v\to 0$ 时会产生渐近无偏估计量。

(a) Lasso (b) MCP (c) CP(SCAD + L₂) (d) NCPM (MCP + L₂)

Figure 1. Threshold Function Diagram of Lasso, MCP, CP and NCPM

图1. Lasso、MCP、CP和NCPM的阈值函数图

2.3. 求解算法

3. 数值模拟及结果分析

3.1. 评价标准

本文选取一些常用的评价指标来衡量模型的泛化能力，具体指标如下：

1) 预测均方误差： $\text{PMSE}=\frac{1}{m}{\displaystyle {\sum }_{i=1}^m{\left(f\left(x_i\right)-y_i\right)}^2}$ (9)

2) 标准误差： $\text{SE}=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^m{\left(f\left(x_i\right)-y_i\right)}^2}}{m}}$ (10)

3) Freq：重要变量被优先选择的平均次数；

4) $N_{nz}$ ：最终模型中非零系数个数；

5) 平均误差： $\text{ME}={\Vert \theta -\hat{\theta }\Vert }_2^2$ (11)

6) AUC：指ROC曲线下的面积

7) 错误发现率： $\text{FDR}=\frac{\text{FP}}{\text{FP}+\text{TP}}$ (12)

8) 假反例率： $\text{FNR}=\frac{\text{FN}}{\text{FN}+\text{TP}}$ (13)

上式m表示样本量； $y_i$， $f\left(x_i\right)$ 分别表示第i个样本的预测值和真实值； $\theta$ 表示参数真实值， $\hat{\theta }$ 表示参数的估计值；TP为真正例；FP为假正例；FN为假反例；TP为真正例。PMSE、SE、ME均是衡量样本真实值与预测值之间偏差的综合指标，FDR、FNR是性能度量指标，这些指标的值越小，表明模型描述数据的精确度越高。ROC曲线是以假正例率为横坐标，真正例率为纵坐标作图得到。若一种方法对应的ROC曲线被另一条曲线完全“包住”，则后者的性能优于前者 [12]。

3.2. 模拟及结果

根据前述模型随机生成模拟数据集。下面给出的四个不同场景，仅通过改变解释因子的维数以及设计矩阵的相关结构来实现，每个场景重复模拟100 ( $N=100$ )次。在场景1和2中，分别生成训练集和测试集，且每生成一个训练集，在相同设置下独立生成相应的测试集，用于评估所得模型的预测性能。每个训练集由100个独立的观察值组成，而每个测试有500个独立的样本；在场景3和4中，整体生成模拟数据，不区分训练集和测试集。

模拟1 在该例中，令 $p=50$， ${\alpha }_j=-0.5,-1,3,1,0.2$， $j=1,11,31,41,50$。 ${\alpha }_j=0$， $j\ne 1,11,31,41,50$ ； ${\beta }_j=1.49,0.4,-0.01,0.39,-0.06$， $j=1,\cdots ,5$。 ${\beta }_j=0$， $j\ne 1,\cdots ,5$。另外，假设： $X=\left(X_1,X_2,\cdots ,X_p\right)~N_p\left(0,\Sigma \right)$， ${\Sigma }_{i,j}={\rho }^{\left|i-j\right|}$， $\rho$ 分别取0.1、0.4、0.7。

模拟1旨在解释因子相关程度不同的情况下，将各惩罚方法的预测精度和变量选择性能进行比较。作为参考，还考虑了Oracle估计结果，即所有重要变量都事先已知的岭回归。

Table 1. Simulation results of example 1

表1. 例1模拟结果

从表1中报告的预测精度结果可以看出，在 $\rho$ 的取值不同的情况下，Ridge是四种惩罚方法中效率最低的。显然，当真实模型中有多余变量时，简约的回归模型比非简约的回归模型具有更强的预测能力。此外，Enet的预测精度略低于MCP和NCPM，这可能是Lasso惩罚的偏差性引起的。可以看出，这些方法中MCP和NCPM是最优的。至于变量选择的性能：岭回归使用所有变量预测，Enet倾向于包含太多冗余变量，而MCP和NCPM通过对未知参数的值进行压缩，实现变量选择。因此，MCP和NCPM是这四种方法中最好的变量选择方法。

模拟2 在该例中，令 $p=50$， ${\alpha }_j=-0.5,-1,3,1,0.2,-0.14$， $j=1,\cdots ,6$。 ${\alpha }_j=0$， $j\ne 1,\cdots ,6$ ； ${\beta }_j=1.49,0.5,0.4,-0.01,0.39,-0.06$， $j=1,\cdots ,6$。 ${\beta }_j=0$， $j\ne 1,\cdots ,6$。令 $X_j,j=1,\cdots ,50$ 为一组 的随机变量，且 $X_j~N\left(0,1\right)$。本例设置两种情况：

1) 当 $j=7,\cdots ,12$ 时，令 $X_j=X_{j-6}+{\eta }_j$，其中， ${\eta }_j\left(i.i.d\right)$ 且 ${\eta }_j~N\left(0,0.01\right)$ ；

2) 当 $j=7,\cdots ,28$ 时，令 $X_j=X_{j+22}+{\eta }_j$，其中， ${\eta }_j\left(i.i.d\right)$ 且 ${\eta }_j~N\left(0,0.01\right)$ ；

(1)与(2)的区别在于：情况(1)中，重要变量间存在相关关系，情况(2)中，非重要变量间具有相关关系。

Table 2. Simulation results of example 2

表2. 例2模拟结果

模拟2旨在比较当解释变量间存在强相关性时，各惩罚方法性能的优越性。表2展示了实例2的模拟结果。MCP在场景1的变量选择中是最好的，但是在场景2中，比其他变量选择要差得多，得到的结果完全是误导性的。在场景2中，使用组合惩罚(Enet和NCPM两种方法)的性能比使用单一惩罚要好得多；NCPM对这两种情况都是最好的。从这个例子中，我们可以看出，在MCP中增加 $L_2$ 惩罚可以显著降低解释因子之间的高共线性带来的风险。

接下来的两个例子旨在比较两种组合惩罚方法(CP和NCPM)在 $\rho$ 较大情况下的效果。

模拟3在该例中，令 $p=200$，剩余条件与模拟1相同。

模拟4在该例中，令 $X_j\left(i.i.d\right)$， $j=1,\cdots ,p$ 且 $X_j~N\left(0,1\right)$。分别，令 $p=100,300,500$。 ${\alpha }_j=-0.5,-1,3,1,0.2$， $j=1,11,31,41,50$。 ${\alpha }_j=0$， $j\ne 1,11,31,41,50$ ； ${\beta }_j=1.49,0.4,-0.01,0.39,-0.06$， $j=1,\cdots ,5$。 ${\beta }_j=0$， $j\ne 1,\cdots ,5$

表3总结了组合惩罚方法的模拟结果。在表格中，我们分别列出了系数 $\alpha$ 和 $\beta$ 在不同衡量指标下的具体情况。模拟结果表明：两种方法的变量选择均有很好的效果；增加解释变量间的相关性对两种方法的变量选择性能影响不大，这表明它们均有很好的处理解释变量间共线性的能力；NCPM比CP给出了更精确的变量选择结果。例如：例3 $\rho =0.4$ 的情境下(表3中第3、4行)，系数 $\alpha$ 在使用NCPM方法时，得到的FDR、FNR、ACU的值分别为：0.03、0.39、0.93，使用CP方法得到的值分别为0.05、0.43、0.92。并且发现，NCPM方法在ME、PMSE上均有所改善。另外，对系数 $\beta$ 也有类似的发现。

Table 3. Simulation 3 and Simulation 4 Results

表3. 模拟3、模拟4结果

4. 实例分析

4.1. 降水量定义

降水量是指从天空落到地面的液态或固态(经融化后)水，未经流失，在水平面上积聚的深度(https://baike.baidu.com/item)。1 mm的降水量是指：在666.7 m²上，降水总深度达到1 mm。

4.2. 数据来源及预处理

降水量相关数据源于气象数据库https://www.aqistudy.cn/historydata/index.php，https://www.wunderground.com/history。本文主要收集了山西省太原市2017年1月~2018年12月期间风速、能见度、气压、日照时长、相对湿度、AQI、PM2.5以及PM10的日平均值。具体说明见表4。

Table 4. Description of relevant variables

表4. 相关变量说明

其中，2018年中有65天降水量大于0；2017年中有58天降水量大于0。

对数据收集整理后，共得到730组观测值。其中607 (约占83.2%)组内降水量的取值为0。另外，为了消除量纲的影响，需要对变量进行标准化处理。本文使用z-score标准化法对原始数据预处理，且经过z-score标准化处理的数据服从 $N\left(0,1\right)$。

z-score标准化公式为：

${\tilde{x}}_i=\frac{x_i-\bar{x}}{\sigma },\bar{x}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i},\sigma =\sqrt{\frac{1}{n-1}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(x_i-\bar{x}\right)}^2}}$ (14)

4.3. 模型建立

根据参考文献 [13] 得到：Gamma分布对降水量非零部分的数据拟合效果最佳。因此，本节使用前述的Logit-Gamma模型进行研究。

4.4. 结果分析及模型比较

运用NCPM方法对实例数据估计，并与Enet方法对比，结果见表5。

Table 5. Example analysis results (Estimated value (Standard error))

表5. 实例分析结果(估计值(标准误差))

由表5新提出方法得到的结果可知:气温与降水量成正相关，这一现象主要受全球变暖的影响。气候变暖，气温升高，水循环加快，大气中的水蒸气增多，降水量也随之增大。风速与降水量成正比。这是因为风速越大，单位时间内进入空气中的水分子越多，蒸发量就越大，导致降水量越大。日照时长与降水量成负相关；空气相对湿度与降水量成正相关，这些结论与实际相符。PM2.5、PM10的浓度以及AQI均与降水量呈负相关。这表明：降水能有效去除大气中的颗粒物，降低空气中PM2.5、PM10的浓度，从而起到净化空气的作用，导致AQI的值降低。

另外发现，新提出的方法得到的模型更简约，模型可解释性更高。如：表5中PM10的Logit部分，用两种方法获得的估计值分别：−2.089、0.077；PM10的Gamma部分，用两种方法得到的估计值分别为：−0.539、1.943，显然利用Enet方法得到的估计与实际意义不符。

5. 结束语

本文以变量选择切入，从理论和数值模拟两方面系统地研究组合惩罚函数的极大似然估计方法在Logit-Gamma两部模型中的表现。具体总结如下：

利用L₂在高度相关解释变量间的良好表现，将基于SCAD + L₂惩罚函数的极大似然估计方法扩展到MCP函数，提出一种新的处理高维统计分析问题的惩罚似然估计方法，NCPM极大似然估计方法。该方法改善了变量间相关性对模型稳定性、精确度的影响。模拟研究表明，当 $p\gg n$ 或解释变量间的相关性较强时，该方法高效便捷且易于实现。

关于组合惩罚似然估计两部模型的变量选择，本文采用LLA-CGD算法。该算法解决了目标函数非线性问题，同时实现计算可行性。数值模拟结果显示，该算法选择效果良好，为两部模型的变量选择提供了新思路。

将提出的NCPM方法应用于Logit-Gamma两部模型，分析太原市降水量的影响因素。结果显示，是否降水主要受露点温度、风速、日照时长、空气相对湿度、PM2.5及PM10浓度等的影响；当降水产生时，降水量多少更易受日平均气温、风速、日照时长、空气相对湿度、PM2.5及PM10浓度、AQI等的影响。最后与Enet方法对比，进一步证实了提出方法具有估计的稳定性、模型可解释性等优势。