高维数据变量选择中MCP正则化参数选择研究

MCP Regularization Parameter Selection in High Dimensional Data Variable Selection

摘要: 大数据时代，高维数据的变量选择是现代统计的研究热点问题之一。MCP正则化方法是常用的变量选取方法，但MCP正则化方法的优劣取决于能否选取出最优的正则化参数。本文在BIC准则的基础上，提出适用于MCP正则化参数选择的MBIC准则。通过数据模拟及实际应用表明，MCP方法在MBIC准则下能够以更高的概率选择正确的模型，MBIC准则明显优于其它参数选择方法。

Abstract: In the era of big data, variable selection of high-dimensional data is one of the hot topics in modern statistics. The MCP regularization method is a commonly used variable selection method, but the merits of the MCP regularization method depend on whether the optimal regularization parameter can be selected. Based on the BIC criterion of regularization parameter selection, an MBIC criterion is proposed for MCP regularization parameter selection. Through data simulation and practical application, the MCP method with MBIC criterion can select the correct model with higher probability, which is obviously superior to other regularization parameter selection methods.

1. 引言

随着大数据时代的来临，庞大的数据资源吸引了越来越多领域的关注。各行各业都希望通过数据处理与挖掘发现数据隐含的信息，为相关决策提供现实依据。特别在分析建模中，为了全面而准确地反应信息的特征及其内在规律，常常引入多个指标，进而形成高维数据。然而并不是高维数据中的所有信息都是有效的，过多的变量反而会导致模型复杂度提升，以至于模型拟合效果和预测精度的降低。因此，如何从海量的高维数据中提取有用特征是一个亟待解决的问题。变量选择就是这样一种从大量信息中提取相关变量从而建立稀疏稳健模型的技术。

传统的变量选择方法如最佳子集选择或逐步向前向后回归，需要结合AIC [1]，BIC [2] 等准则。但在高维数据下，容易出现难以克服的NP-Hard问题。为了克服传统方法的缺陷，统计学家们提出了众多基于惩罚函数的变量选择方法 [3]，如：Lasso估计(Least Absolute Shrinkage and Selection Operator) [4]，SCAD估计(Smoothly Clipped Absolute Deviation) [5]，MCP估计(Minimax Concave Penalty) [6] 等。MCP由Zhang等提出用于高维数据的变量选择，MCP估计满足变量选择的Oracle性质，即一致地选择出正确的模型，且参数的估计满足渐进正态性。在实际应用中，MCP方法优劣取决于能否选择合适的正则化参数λ，正则化参数越小，模型复杂程度越高；正则化参数越大，模型精确程度越低。因此，如何选择合适的正则化参数λ是一个至关重要的问题。常见的正则化参数选择方法有交叉验证(CV)，广义交叉验证(GCV) [4] 以及AIC，BIC [7] 等信息准则。Wang et al. [8] 考虑到GCV方法的过拟合性，提出了参数选择的BIC准则，并从理论上证明了模型选择的一致性。SCAD等常借助BIC准则选择正则化参数 [9]，但该准则对于MCP估计未必能选出最优的模型。

鉴于上述原因，本文通过对BIC准则进行改进，提出一种更适合于MCP正则化参数选择的修正BIC准则(MBIC)。通过数据模拟，比较MBIC准则与BIC准则在MCP方法中的效果。最后，讨论不同方法在实际数据中的应用，分析了1986和1987年赛季美国职业棒球大联盟的棒球运动员收入数据，探究与美国棒球运动员收入相关的影响因素。

2. 罚估计方法

考虑线性模型

$y_i=x_i^{\text{T}}\beta +{\epsilon }_i,i=0,1,\cdots ,n$

其中 $y_i$ 是第 个响应变量， $x_i$ 是 $p\times 1$ 阶的协变量， ${\epsilon }_i$ 是均值为0，方差为 ${\sigma }^2$ 的i.i.d的随机误差项。为了同时进行变量选择和参数估计，常采用很多基于罚函数的稀疏正则化方法，其一般框架为

$L\left(\beta ;\lambda \right)={\Vert y-X{\beta }_{\lambda }\Vert }_2^2+n{\displaystyle {\sum }_{j=1}^{p}p\left(\left|{\beta }_j\right|;\lambda \right)}$ (1)

其中 $p\left(\left|{\beta }_j\right|;\lambda \right)$ 表示惩罚函数。

2.1. MCP估计

2010年，CUN-Hui Zhang提出MCP [6]，MCP是一种非凸罚函数，在 $\left[0,\infty \right)$ 的定义为

$p\left({\beta }_j;\lambda \right)=\lambda \underset{0}{\overset{{\beta }_j}{{\displaystyle \int }}}{\left(1-\frac{x}{\gamma \lambda }\right)}_{+}\text{d}x$

其一阶导数为

$p_{\lambda ,\gamma }\left({\beta }_j\right)=\{\begin{array}{l}\lambda -\frac{{\beta }_j}{\gamma },\text{if}0<{\beta }_j<\gamma \lambda \\ 0,\text{if}{\beta }_j>\gamma \lambda \end{array}$ (2)

其中 $\lambda \ge 0$ 和 $\gamma >1$ 为正则化参数。

结合MCP罚函数的一阶导数的形式，可以看出MCP从0到 $\gamma \lambda$ 惩罚力度呈线性下降趋势，当 ${\beta }_j>\gamma \lambda$ 是惩罚力度为0，即不惩罚。MCP罚函数满足近似连续性，稀疏性和无偏性。

2.2. 其他罚函数

Lasso方法对参数的L1范数进行惩罚，Lasso的惩罚项为 $p_{\lambda }\left(\left|{\beta }_j\right|\right)=\lambda \left|{\beta }_j\right|$，估计形式为 ${\min }_{\beta }{\Vert y-X{\beta }_{\lambda }\Vert }^2,\text{s}\text{.t}\text{.}{\Vert \beta \Vert }_1\le t$ 其中 ${\Vert \beta \Vert }_1={\displaystyle {\sum }_{i=1}^p\left|{\beta }_i\right|}$，上式也等价于

$\hat{\beta }=\arg \underset{\beta }{\min }\left\{\frac{1}{n}{\Vert y-X{\beta }_{\lambda }\Vert }_2^2+\lambda {\Vert \beta \Vert }_1\right\}$

但是，Lasso对较大系数的估计是有偏估计，并且Lasso估计也不满足变量选择的Oracle性质。2001年Fan and Li提出SCAD方法 [5]，同时证明了其满足变量选择的Oracle性质。与Lasso相比，SCAD是无偏估计，因而受到广泛关注。SCAD罚函数的惩罚项为

$p_{\lambda }\left(\left|{\beta }_j\right|\right)=\{\begin{array}{l}\lambda \left|{\beta }_j\right|,\left|{\beta }_j\right|\le \lambda \\ -\frac{{\left|{\beta }_j\right|}^2-2\alpha \lambda \left|{\beta }_j\right|+{\lambda }^2}{2\left(\alpha -1\right)},\lambda \le \left|{\beta }_j\right|\le \alpha \gamma \\ \frac{{\left(\alpha +1\right)}^2\lambda }{2},\left|{\beta }_j\right|>\alpha \gamma \end{array}$ (3)

其中 $\lambda \ge 0$ 和 $\alpha >1$ 为正则化参数，在实际应用中常取 $\alpha =3.7$。

3. 正则化参数选择方法

在实际应用中，正则化模型(1)的优劣与正则化参数 $\lambda$ 取值密切相关，不同的参数 $\lambda$ 会导致不同的惩罚力度，进而影响最终的模型。因此，参数λ的选择至关重要。常见的选择参数 $\lambda$ 的方法有CV，GCV和各种信息准则，如AIC及BIC等。

针对LASSO估计，Zou H. et al. [8] 给出了估计的自由度，并提出了适用于Lasso估计的BIC准则，定义如下

$\text{BIC}\left(\lambda \right)=\frac{{\Vert y-X{\beta }_{\lambda }\Vert }_2^2}{n{\sigma }^2}+\frac{\log \left(n\right)}{n}{\widehat{df}}_{\lambda }$ (4)

其中 ${\widehat{df}}_{\lambda }={\displaystyle {\sum }_{i=1}^{n}cov\left(x_i{\beta }_{\lambda }-y_i\right)/{\sigma }^2}$， ${\sigma }^2=Var\left(\epsilon \right)$。

此外，Wang et al. [7] 证明了GCV方法易出现模型选择的过拟合现象，针对SCAD估计提出了BIC准则。BIC准则是在有限的模型集合中的模型选择准则，BIC准则认为具有最小BIC值的模型是模型集合中最优良的模型。BIC具体定义如下

$\text{BIC}\left(\lambda \right)=\log \frac{{\Vert y-X{\beta }_{\lambda }\Vert }_2^2}{n}+\frac{{\widehat{df}}_{\lambda }}{n}\log \left(n\right)$ (5)

其中 ${\widehat{df}}_{\lambda }$ 为广义自由度，

${\widehat{df}}_{\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|{\hat{\beta }}_{\lambda 1}\right|\right)/\left|{\hat{\beta }}_{\lambda 1}\right|,\cdots ,{p^{\prime }}_{\lambda }\left(\left|{\hat{\beta }}_{\lambda d}\right|\right)/\left|{\hat{\beta }}_{\lambda d}\right|\right)$。

但在实际应用中，MCP估计在BIC准则下选择了较为复杂的模型，故本文提出MBIC准则，定义如下

$\text{MBIC}\left(\lambda ,\alpha \right)=\log \left(\frac{{\Vert y-X{\beta }_{\lambda ,\alpha }\Vert }_2^2}{n-p_0}\right)+\frac{p_0}{n}\log \left(n\right)$ (6)

其中， $p_0$ 表示非零变量个数。

4. 模型研究和实际数据分析

4.1. 模拟研究

小节通过模拟实验比较LASSO，SCAD，MCP变量选择方法的性能。

考虑线性模型

$y=X^{\text{T}}\beta +\sigma \epsilon$

进行随机模拟，从而产生数据x和y。在模拟实验中， $n=200$， $\epsilon ~N\left(0,1\right)$， $\sigma =2$，变量个数p分别取8，12，20，且 $\beta ={\left(3,1.5,0,0,2,0,0,0,\cdots \right)}_{1\times p}^{\text{T}}$， $x_i,x_j$ 之间的相关系数为 $cor\left(j_1,j_2\right)={0.5}^{\left|j_1-j_2\right|}$。

算法上，Lasso估计，SCAD估计和MCP估计均采用坐标下降算法 [10]。MCP估计分别利用BIC准则(5)，MBIC准则(6)选择正则化参数，SCAD估计采用BIC准则(5)，而LASSO估计采用BIC准则(4)选择正则化参数。所有模拟实验重复进行100次，模拟结果如表1所示。

为比较Lasso、SCAD、MCP估计精确性，给出模型误差公式

$ME\left(\hat{\mu }\right)={\left(\hat{\beta }-\beta \right)}^{\text{T}}E\left(X{X}^{\text{T}}\right)\left(\hat{\beta }-\beta \right)$

其中，“MME”表示100次重复实验中模型误差ME的中位数；“SD”表示100次重复实验中模型误差ME的标准差；“C”表示100次重复实验中非零系数被正确估计为非零个数的均值；“IC”表示100次重复实验中零系数被错误估计为非零个数的均值；“Underfit”表示欠拟合，即在100次模拟实验中将非零系数错误估计为零的比例；“Correctfit”表示正确拟合，即在100次模拟实验中将非零系数正确估计为非零的比例；“Overfit”表示过拟合，即100次模拟实验中选择了3个重要变量并且包含了非零系数的比例。

从表1可以看出，在BIC正则化参数选择方法下，MCP估计和SCAD估计方法在变量选择和模型误差方面优于LASSO方法。在模型误差中，所有变量选择的方法均能减小模型误差，而MCP方法在MBIC准则下具有最小的模型误差，而且能够以更高概率选择真实模型。综上，在MBIC正则化参数选择方法下，MCP估计在变量选择能力和模型误差方面均最优。

Table 1. Simulation results

表1. 模拟结果

4.2. 实际数据分析

本文利用来自R语言的ISLR包中的数据集Hitters，该数据集包含20个变量，322次样本，描述关于1986和1987赛季的棒球大联盟中的棒球运动员收入的相关信息。数据集各个变量描述如下：

X1：1986年击球的次数；X2：1986年的点击次数；X3：1986年的本垒打数量；X4：1986年的运行次数；X5：1986年击败的次数；X6：1986年的散步次数；X7：联赛的年份；X：职业生涯中击球的次数；X9：职业生涯中的点击次数；X10：职业生涯中的本垒打数量；X11：职业生涯中的跑步次数；X12：职业生涯中击球的次数；X13：职业生涯中的散步次数；X14：表示1986年底的球员联赛A级和N级的因素；X15：表示1986年底的分裂E和W等级的因素；X16：1986年的罢工数量；X17：1986年的助攻数量；X18：1986年的错误数量；X19：表示1987年初的球员联赛A级和N级的因素；Y：1987年开业日的年薪数(千美元)。

假设棒球运动员在各个赛季的表现与棒球运动员收入呈线性关系，即有如下线性模型

$y_i=\underset{j=1}{\overset{20}{{\displaystyle \sum }}}{x}_{ij}{\beta }_j+{\epsilon }_i,i=0,\cdots ,n$

其中 $y_i$ 表示第i个运动员的收入， $x_{ij}$ 是他的第j个变量， ${\epsilon }_i$ 是均值为0，方差为 ${\sigma }^2$ 的i.i.d的随机误差项。

利用最小二乘估计(OLS)、Lasso、SCAD和MCP估计分析该数据。变量选择结果如表2所示。从表2可以看出，无罚的最小二乘估计(OLS)选择了所有的变量，Lasso选择了15个变量，SCAD选择了6个变量，基于BIC准则的MCP估计选择了8个协变量，而基于MBIC准则的MCP估计选择了7个协变量，选择了相对稀疏的模型。从参数估计的结果看，基于MBIC的MCP估计的结果更接近于最小二乘估计值。

Table 2. Parameter estimation under different methods

表2. 不同方法下的参数估计

5. 结论

本文讨论了MCP方法在变量选择和参数估计的应用，提出了更适合MCP估计的MBIC准则。数据模拟以及实际数据分析中都表明在MBIC准则下MCP估计的结果更优于BIC准则估计结果。

基金项目

浙江省自然科学基金资助项目(LY18A010026)；全国统计科学研究项目(2019LY06)；宁波市自然科学基金资助项目(2017A610143)；国家级大学生创新创业训练计划项目(201911058025)；王伟明助创基金资助项目(2018020)；浙江省大学生科技创新活动计划暨新苗人才计划资助项目(2018R428027)。

NOTES

^*通讯作者。