基于SICA罚的变量选择及应用

Variable Selection and Application Based on SICA Penalty

摘要: 高维数据的变量选择一直是统计学领域的热门研究方向。本文研究SICA罚估计在线性模型变量选择中的应用，结合LLA (Local linear approximation)和坐标下降算法给出一种有效的迭代算法，并提出BIC准则选择正则化参数。实际数据的分析表明，与其他变量选择方法相比较，SICA方法在参数估计精度和变量选择方面具有较好的表现。

Abstract: Variable selection of high-dimensional data has always been a hot research direction in the field of statistics. In this paper, we study the application of SICA penalty estimation in variable selection of linear model, give an effective iterative algorithm combined with LLA (local linear approximation) and coordinate descent algorithm, and propose BIC criterion to select regularization parameters. The analysis of actual data shows that SICA method has better performance in parameter estimation accuracy and variable selection compared with other variable selection methods.

1. 引言

在统计建模的过程中，我们总是希望筛选出对响应变量影响较强的变量，剔除没有影响或影响较弱的变量。但面对一些高维数据，往往存在维数灾难(curse of dimension)的现象，即数据的维数要远远多于样本量的大小。这对传统的模型选择方法提出了巨大的挑战。

早在1970年，Hoerl和Kennard提出岭回归 [1]，是对最小二乘估计的改进。岭回归通过放弃最小二乘法的无偏性，以损失部分信息、降低精度为代价获得回归系数更为符合实际、更可靠的回归方法，对病态数据的拟合要强于最小二乘法。但岭回归不具有变量选择的功能，近年来，统计学家门提出过一系列变量选择的方法 [2] [3] [4]。1996年，受NG (Nonnegative Garrot) [5] [6] 方法的启发，Robert Tibshirani首次提出LASSO方法 [7]。该方法是岭回归的一种特殊形式。它通过构造一个惩罚函数来压缩一些回归系数，即强制系数绝对值之和小于某个固定值；同时设定一些回归系数为零，以此达到变量选择的目的。但LASSO估计是有偏估计，2001年，Fan和Li [8] 提出SCAD方法，并理论上证明了SCAD罚估计具有Oracle性质。此后，2010年，Zhang提出MCP估计 [9]，也是一种近似无偏估计。

2009年，Lv and Fan提出了一种新的罚函数估计方法——SICA罚 [10] 用于变量选择和参数估计。本文针对线性模型，研究基于SICA罚的变量选择，结合LLA和坐标下降提出一种迭代算法，并提出BIC准则选择正则化参数。最后通过对实际数据分析，与其他方法进行比较。

2. SICA罚估计原理

考虑线性模型

$y=X\beta +\epsilon$ (1)

其中，y为 $n\times 1$ 的响应变量；X为 $n\times p$ 设计阵； $\beta =\left({\beta }_1,{\beta }_2,\cdots ,{\beta }_p\right)$，为 $p\times 1$ 维未知参数向量； $\epsilon$ 为 $n\times 1$ 独立同分布随机误差向量，均值为0，方差为 ${\delta }^2$。

基于模型(1)的SICA惩罚最小二乘定义如下：

$Q_{\left(n\right)}\left(\beta \right)=\frac{1}{2n}{\Vert y-X\beta \Vert }^2+{\displaystyle \underset{j=1}{\overset{p}{\sum }}{p}_{\lambda ,\tau }}\left(\left|{\beta }_j\right|\right)$ (2)

其中

$p_{\lambda ,\tau }\left(\left|{\beta }_j\right|\right)=\lambda \frac{\left(\tau +1\right)\left|{\beta }_j\right|}{\left|{\beta }_j\right|+\tau }$ (3)

为SICA惩罚函数， $\lambda ,\tau$ 为正则化参数。 $\lambda$ 控制罚函数的惩罚力度， $\tau$ 控制罚函数的凹凸度。当取 $\lambda =1$, $\tau$ 分别取0.01，0.05，0.1，绘制SICA罚函数图像如图1所示。由图可见 $\tau$ 的取值越小，SICA罚越趋近于L0罚。

Figure 1. SICA penalty function

图1. SICA罚函数

对上述(2)式求极小值，得到SICA惩罚最小二乘估计 $\hat{\beta }=\arg {\min }_{\beta }{Q}_n\left(\beta \right)$。记 $\stackrel{⌢}{\beta }=\left\{{\stackrel{⌢}{\beta }}_j;j=1,\cdots ,p\right\}$。极小化 $Q_n\left(\beta \right)$ 的结果使一些系数为0，从而实现了保留对响应变量影响较大的变量，去除对响应变量影响较小的变量。所以我们可以通过对 $Q_n\left(\beta \right)$ 实施极小化同时达到系数估计和变量选择的效果。

3. SICA罚估计的算法及参数选择

3.1. 算法

在本文中使用局部一次逼近(LLA)方法来求解模型。对于固定正则化参数 $\lambda$，LLA方法采用局部线性近似的方法估计 ${\displaystyle \underset{j=1}{\overset{p}{\sum }}{P}_{\alpha }}\left(\left|{\beta }_j\right|\right)$。对于给定的初值 ${\beta }^0={\left({\beta }_1^0,\cdots ,{\beta }_p^0\right)}^{\text{T}}$，LLA将 ${\displaystyle \underset{j=1}{\overset{p}{\sum }}{p}_{\alpha }}\left(\left|{\beta }_j\right|\right)$ 进行一阶泰勒展开：

${\displaystyle \underset{j=1}{\overset{p}{\sum }}\left[p_{\lambda ,\tau }\left(\left|{\beta }_j^0\right|\right)+{p^{\prime }}_{\lambda ,\tau }\left(\left|{\beta }_j^0\right|\left(\left|{\beta }_j\right|-\left|{\beta }_j^0\right|\right)\right)\right]}$ (4)

将上式带入SICA罚函数得到：

$\frac{1}{2}{\Vert y-X\beta \Vert }^2+{\displaystyle \underset{j=1}{\overset{p}{\sum }}\left[p_{\lambda ,\tau }\left(\left|{\beta }_j^{\left(0\right)}\right|\right)+{p^{\prime }}_{\lambda ,\tau }\left(\left|{\beta }_j^{\left(0\right)}\right|\left(\left|{\beta }_j\right|-\left|{\beta }_j^0\right|\right)\right)\right]}$ (5)

其中 $p_{\lambda ,\tau }\left(\left|{\beta }_j\right|\right)=\lambda \frac{\left(\tau +1\right)\left|{\beta }_j\right|}{\left|{\beta }_j\right|+\tau },j=1,2,\cdots ,p$。上式其实是一种加权LASSO模型，然后采用坐标下降算法求最小值即得到各参数的值，模型得以求解。

3.2. 正则化参数选择

使用SICA方法时要确定参数 $\lambda$ 和 $\tau$。对于参数 $\tau$ 可以直接取0.01，相当于SCAD中取a = 3.7 [11]。所以我们只需要固定 $\tau =0.01$，选取适当的参数 $\lambda$ 即可。AIC准则 [12]，GCV准则 [13] 是常用的参数选择方法，但其不具有一致性，易发生过拟合的现象。BIC准则没有这一缺点，所以本文提出基于BIC准则的 的选取：

${\text{BIC}}_{\lambda }=\log \frac{{\Vert y-X{\beta }_{\lambda }\Vert }_2^2}{n-k}+\frac{d{f}_{\lambda }}{n}\log \left(n\right)$ (6)

其中 $d{f}_{\lambda }$ 为广义自由度：

$d{f}_{\lambda }=tr\left\{X{\left(X^{\prime }X+n{\displaystyle \sum \lambda }\right)}^{-1}X\right\}$

其中 ${\displaystyle \sum \lambda }=diag\left({p^{\prime }}_{\lambda }\left(\left|{\stackrel{⌢}{\beta }}_{\lambda 1}\right|\right)/\left|{\stackrel{⌢}{\beta }}_{\lambda 1}\right|,\cdots ,{p^{\prime }}_{\lambda }\left(\left|{\stackrel{⌢}{\beta }}_{\lambda k}\right|\right)/\left|{\stackrel{⌢}{\beta }}_{\lambda k}\right|\right)$，k为模型中非零参数个数，n为样本数量，取 $\stackrel{⌢}{\lambda }=\arg {\min }_{\lambda }\left({\text{BIC}}_{\lambda }\right)$。

4. 实际数据应用

实际数据分析

本文使用Python的Keras库中的数据集Boston_House，该数据集包含美国人口普查局收集的美国马萨诸塞州波士顿住房价格的有关信息，共506个样本，13个变量。通过多种方法对房价与各自变量之间建立回归模型。对数据的具体描述如下：

X₁代表犯罪率；X₂代表住宅用地所占比例；X₃代表非零售地区所占比例；X₄代表是否在河边(0代表不是，1代表是)；X₅代表一氧化氮浓度；X₆代表平均每居民房数；X₇代表建筑年龄；X₈代表与市中心的距离；X₉代表公路可达指数；X₁₀代表物业税率；X₁₁代表城镇师生比例；X₁₂代表黑人比例；X₁₃代表低收入人口所占比例。

首先将数据以7:3的比例划分为训练集与测试集。分别采用经典最小二乘法，岭回归，LASSO回归，SICA方法对测试集进行回归建模，将所得模型运用于测试集做预测。从变量选择的效果与预测准确性来比较不同方法的优劣，准确性以MSE作为判断标准。

MSE称为均方误差(Mean Square Error)，是真实值与预测值的差值的平方然后求平均数，计算公式如下：

$\text{MSE}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(y_i-{\stackrel{⌢}{y}}_i\right)}^2}$ (7)

其中 $y_i$ 表示真实值， ${\stackrel{⌢}{y}}_i$ 表示预测值。

SICA方法首先要确定参数 $\lambda$ 的值，事先给定 $\lambda$ 一个取值范围，查看不同 $\lambda$ 下BIC的值，取使得BIC最小的 $\lambda$ 为参数。BIC随着参数 $\lambda$ 变化而变化的结果如图2所示：

Figure 2. BIC curve

图2. BIC变化曲线

可以看出BIC曲线呈现V型，说明已经找到使得BIC最小的 $\lambda$ 值。模型拟合的系数和MSE如表1，表2所示：

Table 1. MSE under different methods

表1. 不同方法的MSE

Table 2. Coefficient under different methods

表2. 不同方法的变量系数

SICA方法剔除了5个变量，并且所得模型预测结果的MSE较小。这说明相对于最小二乘回归和岭回归，SICA方法准确的剔除了对因变量影响较小的因素，保留了影响较大的因素。LASSO方法虽然剔除了更多的变量，但是MSE却是最大的，说明剔除了应该保留的变量，变量选择的能力上不如SICA方法。逐步回归方法剔除了5个变量，MSE也小于SICA方法，说明对于这份数据，逐步回归的变量选择与模型预测优于SICA方法。综上，SICA方法无论是变量的选择还是模型的预测精度都优于大部分的传统方法。

5. 结论

本文讨论了SICA方法在线性模型的变量选择和参数估计中的应用。通过对实际数据的分析，可以发现相对于大部分的传统方法，本文提出的SCIA方法有更强的变量选择能力，对参数的估计具有更高的精度。

基金项目

全国统计科学研究项目(2019LY06)；浙江省统计研究课题(20TJZZ18)；浙江省自然科学基金资助项目(LY18A010026)；浙江省大学生科技创新活动计划暨新苗人才计划资助项目(2020R475013)。

NOTES

^*通讯作者。