基于Geman-McClure损失的稳健变量选择

Robust Variable Selection Based on Geman-McClure Loss

摘要: 基于惩罚函数的最小二乘估计或似然估计是变量选择的有效方法。但当数据存在异常值时，罚最小二乘或似然估计的稳健性受到极大挑战。本文提出基于Geman-McClure损失的稳健变量选择方法，该损失函数能够有效抵制数据中异常值的影响。数值模拟和实际数据分析验证了该模型的有效性和稳健性。

Abstract: Least squares or likelihood estimation based on penalty function is an effective method for variable selection. However, the robustness of penalized least squares or likelihood estimation is greatly challenged when there are outliers in the data. In this paper, we propose a robust variable selection method based on the Geman-McClure loss, which is an effective loss function to counteract the influence of outliers in the data. Numerical simulations and analysis of real data validate the validity and robustness of the model.

1. 引言

关于变量选择问题，很多学者做了大量研究。上世纪70年代，人们提出信息准则方法，如1973年Akaike [1] 提出AIC准则，1978年Schwarz [2] 在贝叶斯理论的基础下提出BIC准则。然而，随着变量维数的增加，基于准则选取变量的方法，计算复杂度会急剧增加，效率低下。近年来，变量选择的稀疏正则化方法已逐渐流行起来。1996年，Tibshirani [3] 提出了LASSO (Least Absolute Shrinkage and Selection Operator)，该方法通过L1范数进行惩罚来压缩回归系数的大小，使绝对值较小的回归系数被自动压缩为0。尽管L1范数是凸函数且易于求解，但LASSO为有偏估计，并且变量选择的一致性需要满足一定的不可表示条件 [4] (Irrepresentable Condition)和系数Riesz条件 [5] (Sparse Riesz Condition)。为了解决上述问题，一系列非凸正则化方法被提出。2001年，Fan和Li [6] 提出了SCAD (Smoothly Clipped Absolute Deviation Penalty)罚，是一种近似无偏稀疏估计。2006年，Zou [7] 在LASSO的基础上提出了Adaptive LASSO，该方法是一种LASSO的改进。SCAD和Adaptive LASSO在一定条件下都满足Oracle性质。2010年，T. Zhang [8] 提出Capped L1实现模型的稀疏解。Zhang [9] 提出了MCP (Minimax Concave Penalty)惩罚。很多研究表明，非凸惩罚函数在理论分析以及应用中具有更优的表现 [10] 。

近些年，诸多高维数据变量选择研究工作的前提假设是误差分布为高斯分布或次高斯分布。然而许多实际数据如气候数据，保险理赔数据，电子商务数据等往往服从重尾分布，对于此类数据，上述方法并不适用。模型的稳健性受到极大挑战。统计学家针对具有重尾误差的情形提出了若干稀疏正则化方法，如基于Huber损失的高维M估计 [11] 、基于LAD [12] 或分位数的损失函数的估计、稳健M估计中拟似然估计方法 [13] 、基于梯度下降算法的稳健估计 [14] 、基于指数平方损失 [15] 或t型损失 [16] 的稳健回归、基于Wilcoxon得分函数的秩LASSO估计 [17] 等。

针对误差分布为重尾分布或数据存在异常值的高维模型，本文提出了一种基于Geman-McClure损失的稳健罚估计方法。该方法在X空间或Y空间存在离群值时，依旧能稳健且有效的进行变量选择。

2. 基于Geman-McClure损失的稳健模型

考虑线性回归模型

$y=X^{\text{T}}\beta +\epsilon$ (1)

其中， $y\in R^N$ 为响应变量， $X\in R^{N\times p}$ 为设计矩阵， $\beta \in R^p$ 为回归系数向量， $\epsilon \in R^N$ 为误差向量且服从独立同分布 ${\epsilon }_n\in N\left(0,{\sigma }^2\right)$ ， $n\in \left\{1,2,\cdots ,N\right\}$ 。

高维数据变量选择的稀疏正则化方法的一般框架为：

$\hat{\beta }=\underset{\beta \in R^p}{\arg \min }\left\{{\displaystyle \underset{i=1}{\overset{n}{\sum }}\varphi \left(y_i-x_i^{\text{T}}\beta \right)+n{\displaystyle \underset{j=1}{\overset{p}{\sum }}{p}_{\lambda }\left(\left|{\beta }_j\right|\right)}}\right\}$ (2)

其中， $\varphi \left(\text{·}\right)$ 为损失函数， $p_{\lambda }\left(\left|{\beta }_j\right|\right)$ 为罚函数， $\lambda \ge 0$ 。常见的罚函数如LASSO、SCAD、Adaptive LASSO和MCP等。鉴于Adaptive LASSO的优良统计性质 [7] ，本文选取Adaptive LASSO作为罚函数。Adaptive LASSO形式如下：

$\hat{\beta }=\underset{\beta \in R^p}{\arg \min }\frac{1}{2}{\Vert y-X\beta \Vert }_2^2+\lambda {\displaystyle \underset{p=1}{\overset{P}{\sum }}{\hat{w}}_p\left|{\beta }_p\right|}$ (3)

其中， ${\Vert \cdot \Vert }_2$ 表示L₂范数， $\lambda \ge \text{0}$ ， ${\hat{w}}_p\ge {\text{1}/\left|{\hat{\beta }}_p^{OLS}\right|}^{\gamma }$ 为第p个回归系数的权值， $\gamma >\text{0}$ ， ${\hat{\beta }}_p^{OLS}$ 为普通最小二乘估计的解。可以看出，Adaptive LASSO的实现需通过两步进行：1) 先进行最小二乘估计，将系数估计值的 $\gamma$ 次方的倒数作为第p个变量的权值；2) 对每个变量赋予“量身定做”的权值后，将权值代入(3)式进行求解。

Adaptive LASSO采取的损失函数为平方损失，该损失注定了Adaptive LASSO不适合用于数据存在异常值的情形。本文在Adaptive LASSO以及Geman-McClure损失的基础上，提出一种稳健且有效的变量选择模型，形式如下：

$\hat{\beta }=\underset{\beta \in R^p}{\arg \min }\left\{\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\varphi \left(y_i-x_i^{\text{T}}\beta \right)+\lambda {\displaystyle \underset{p=1}{\overset{P}{\sum }}{\hat{w}}_p\left|{\beta }_p\right|}}\right\}$ (4)

其中，Geman-McClure损失函数为 $\varphi \left(t\right)=t^2/\left(\gamma +t^2\right)$ ，其中 $t=y_i-x_i^{\text{T}}\beta$ 。该损失函数保证正常样本情况下，其自身灵敏度的同时，又降低了对异常样本的敏感程度，提高对变量选择的稳健性。

3. 模拟研究

在本例中，模拟数据由线性回归模型生成，

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

其中变量个数为40，非零系数分别为 ${\beta }_1=3$ ， ${\beta }_2=1.5$ ， ${\beta }_5=2$ ，其余系数均为0，样本大小 $n=200$ ，对于每种情况，重复模拟100次。

情景1： $X_i$ 中存在异常值。每个预测值 $X_i$ 均为40维正态分布的混合样本， $\left(1-c\%\right)N\left(0,1\right)+c\%N\left(3,3^2\right)$ ，其中c为异常样本比例，误差项服从标准正态分布。

情景2： $y_i$ 中存在异常值。每个预测值 $X_i$ 均服从标准正态分布，误差项服从 $\left(1-c\%\right)N\left(0,1\right)+c\%N\left(10,6^2\right)$ ，其中c为异常样本比例。

情景3：误差项为t分布，每个预测值 $X_i$ 均服从标准正态分布，误差项服从 $t\left(k\right)$ ，其中 $k=2,3,5$ 。

为了与本文提出的方法比较，本模拟也利用Adaptive LASSO估计，用10折CV进行参数调节，通过使用R软件包“glmnet”来实现。而基于Geman-McClure损失和自适应LASSO的变量选择方法(以下简记为GM-ALASSO)采用BCGD算法 [18] 对模型进行求解。

为了评价模拟效果，我们计算了估计系数与真实系数之间均方误差的中位数(MSE)；平均模型大小(即非零系数的数量)，MS。更好的模型选择应该产生更准确的预测结果(即较小的MSE值)、更正确的模型大小(即MS)。此外，为了评估变量选择表现，我们还考虑了假阳性率(FPR)和假阴性率(FNR)，定义如下

$\begin{array}{l}FPR=\frac{\#ofselectedunimpo\text{rtant}v\alpha riables}{\#ofselectedvariables},\\ FNR=\frac{\#ofremovedunimpor\text{tant}v\alpha riables}{\#ofremovedvariables}.\end{array}$

当未选择变量时，FPR为0，当所有变量已选择，FNR为0。同时，我们还计算了Hamming距离(HD)，其中HD = FN + FP，其中FN表示非零系数被判为零系数的平均次数，FP表示零系数被判为非零系数的平均次数。在结果上，更好的模型应具备更小的HD。

表格中，在标记为“Under fit”的列中，我们给出了200次重复实验中去除了任何非零系数的比例。同理，“Correct-fit”表示正确选择模型的概率，“Over-fit”列表示选择了一些噪声变量的概率。

从表1可以看出，在X空间数据未被污染的情况下，Adaptive LASSO和GM-ALASSO估计的准确率均在0.9附近，两种方法的估计效果均良好，这证明了GM-ALASSO在X空间无污染情况下估计的有效性。然而，在X空间数据被污染的情况下，随着污染程度的增加，Adaptive LASSO变量选择的能力急剧下降，模型过拟合程度急剧增加，而GM-ALASSO基本稳定在0.7附近，过拟合程度也显著低于Adaptive LASSO，这证明了GM-ALASSO在X空间被污染情况下变量选择方法的稳健性。

Table 1. Scenario 1 estimated results

表1. 情景1估计结果

同理，从表2可看出，在Y空间数据未被污染的情况下，Adaptive LASSO和GM-ALASSO估计的准确率均较为良好，证实GM-ALASSO在Y空间未被污染情况下的有效性。在Y空间被污染的情况下，随着污染程度的加剧，GM-ALASSO的估计效果愈加强于Adaptive LASSO，模型的过拟合程度较Adaptive LASSO也较低，证实GM-ALASSO在Y空间存在污染情况下的有效性。但两者估计效果在Y空间被污染情况下的准确率均会高于在X空间存在污染的情况，说明GM-ALASSO对Y空间存在异常值的抵抗力更为强大。

Table 2. Scenario 2 estimated results

表2. 情景2估计结果

同理，从表3可看出，随着模型重尾程度的加剧，GM-ALASSO的变量选择能力会显著强于Adaptive LASSO，证明了GM-ALASSO在数据服从重尾分布情况下的稳健性。

Table 3. Scenario 3 estimated results

表3. 情景3估计结果

综合而言，GM-ALASSO在X空间或Y空间存在异常值时以及样本服从重尾分布的情况下的变量选择能力均强于Adaptive LASSO，模型整体的稳健性及有效性均较为良好。

4. 实证分析

模拟研究

在本节中，我们将GM-ALASSO变量选择方法应用于波士顿房价数据。该数据是统计数据分析类非常著名的一类数据集，其中包括决定房价的结构因素、环境因素和教育因素。同时，该数据集也属于公开研究的数据，本节采用的数据集来自R软件自带的波士顿数据集，可通过data (“Boston”)命令从R软件自带数据集中调出。包含13个可能影响房价的变量：crim (犯罪率)、zn (高于25000平方英尺房屋比率)、indus (非零售商业区比率)、nox (氮氧化物浓度)、rm (住宅平均房间数)、age (1940年前自住房比率)、dist (与波士顿五大就业中心的加权距离)、rad (高速公路便利指数)、tax (不动产税)、ptratio(学生–老师比例)、black (黑人比例)、lstat (低教育人口比例)、chas (查尔斯河虚拟变量)。响应变量medv (住房价格中位数)。

近年来，有许多方法对波士顿数据进行研究。例如分位数回归算法 [19] 、梯度下降算法 [20] 、SARCH模型、QPLSIM模型和QPLAM模型 [21] 、非参数方法和半参数方法探索数据结构 [22] 。这些方法都可以用来研究波士顿房价数据。波士顿数据集分为两部分：训练集有354个样本，测试集有152个样本。本文分别采用Adaptive LASSO和GM-ALASSO两种方法对波士顿房价数据进行变量选择。为了验证稳健性，将d%的数据进行污染，通过将d%的数据加上0.05×自身值。在本文中，d分别取0，3，5进行实验。将两种算法变量选择后的结果进行线性建模，结果见表4。

由表4可见，当数据不加污染时，虽然GM-ALASSO变量选择方法估计的RMSE以及MAE比Adaptive LASSO方法的要大，但是差别不大，说明在无污染的情况下，GM-ALASSO变量选择方法是有效的。但当数据被污染时，GM-ALASSO变量选择方法的RMSE以及MAE比Adaptive LASSO更小，且R²更大，说明GM-ALASSO在数据由污染时更稳健。从表5可见，在变量选择上，GM-ALASSO选取的变量较为稳定，随着污染的加剧，GM-ALASSO均只在无污染的前提下，新增变量dis，而Adaptive LASSO的波动较为明显，随着污染的加剧，变量选择的改变幅度也在增大，这也说明了GM-ALASSO变量选择在污染情况下的稳健性。

Table 4. Experimental results

表4. 实验结果

Table 5. Estimated regression coefficients for Boston house price data

表5. 波士顿房价数据的估计回归系数

5. 结论

本文在Adaptive LASSO变量选择的框架下提出了一种稳健且有效的变量选择方法，基于Geman-McClure损失提出一种新损失从而达到稳健的效果。模拟结果和实际数据表明，GM-ALASSO方法能够以比较高的概率选择正确的模型且具有较小的模型误差。与传统方法相比，GM-ALASSO方法更适合数据存在异常值的情况。此外，相关的理论性质尚未被讨论，这是未来的研究方向之一。

基金项目

宁波工程学院崇本基金项目(2022014)；宁波市自然科学基金项目(2021J143, 2021J144)。

NOTES

^*通讯作者。