高维线性回归模型稳健变量选择方法综述.pdf

应用概率统计第 40 卷第 1 期2024 年 2 月Chinese Journal of Applied Probability and StatisticsFeb.,2024,Vol.40,No.1,pp.157-181doi:10.3969/j.issn.1001-4268.2024.01.010综述报告高维线性回归模型稳健变量选择方法综述邹航姜云卢(暨南大学经济学院,广州,510632)摘要:随着大数据时代的到来,在经济学、金融学和生物医学等众多研究领域中频繁收集到高维数据.高维数据的特征之一是变量维数 p 随着样本量 n 的增加而变大且通常会超过样本量,同时,异常值也容易出现在高维数据中.因此,如何克服异常值给高维统计推断带来的影响,从而得到更精确的模型,是目前统计学研究的热点问题之一.本文是对高维线性模型下的稳健变量选择方法进行综述.具体地,首先介绍评估稳健性的三个指标:影响函数、崩溃点和最大偏差.其次着重介绍了稳健变量选择方法,包括响应变量含有异常值,响应变量和协变量都含有异常值,高崩溃点且高效的变量选择方法.紧接着介绍相关算法,通过模拟和实例比较不同变量选择方法.最后,简要探讨了高维稳健有效变量选择方法存在的问题及未来的可能发展方向.关键词:高维线性回归模型;稳健性;变量选择;有效性中图分类号:O212.7英文引用格式:ZOU H,JIANG Y L.Overview of robust variable selection methods for high-dimensional linear regression modelJ.Chinese J Appl Probab Statist,2024,40(1):157181.(inChinese)1引言随着科学技术的高速发展,在许多研究领域中,测量和收集到高维数据越来越来普遍,如化学计量学中的光谱数据,生物信息学中的遗传数据和金融投资学中的投资组合数据等.高维数据的本质特征是模型中未知参数的个数 p 远远大于样本量 n.Fan 和 Lv1指出:这类高维数据的存在使得评估模型优劣的三个重要准则统计准确性、模型解释力和计算稳定性均难以满足.当解释变量的个数随着样本容量的增加而增加时,传统的统计方法将不再适用,这就迫切要求统计学家改进传统的统计方法或提出新的统计方法来应对高维数据带来的挑战.国家自然科学基金项目(批准号:12171203)、广东省自然科学基金项目(批准号:2022A1515010045)和中央高校基本科研业务费专项资金项目(批准号:23JNQMX21)资助.通讯作者,E-mail:.本文 2022 年 3 月 30 日收到,2022 年 8 月 7 日收到修改稿.158应用概率统计第 40 卷为了更好地处理高维数据,统计学家们通常会对变量做出稀疏性假设.如何利用稀疏性假设从高维解释变量中选择出与响应变量存在明显关系的真实协变量,从而减少待估参数的个数,即变量选择,是进行高维统计推断的研究热点问题之一.变量选择方法普遍分为两种:最优子集选择法和惩罚正则化方法.最优子集选择法主要包括 Cp准则2、AIC 准则3、BIC 准则4、最优子集回归5等等.这些方法虽然容易理解,但存在一些缺点:例如离散选择过程的不稳定性,计算成本随着变量个数的增加而不可度量,估计的大样本性质难以研究等.基于上述问题,学者们提出了惩罚正则的方法,包括岭回归6、LASSO7、SCAD8、自适应 LASSO9、弹性网10、自适应弹性网11、MCP12、截断 L1函数13、LEP14等等.这些方法大多数可以同时实现参数估计和变量选择,且在一定的条件下具有Oracle 性质,因此被广泛研究.然而,上述方法都是基于平方误差损失函数或者似然函数,容易受到数据中的异常值影响,因此并不稳健.Zhou 等15指出:在高维数据中,异常数据出现的机率越来越大.同时,Hampel 等16指出,一个数据集中一般会包含 1%10%的异常点,在一些如网络数据等更复杂的数据中可能会包含更多.当高维数据中存在异常值时,继续使用传统的变量选择方法将会得到错误的结论.因此,如何应对异常值对估计量造成的不利影响,从而提高模型的预测准确度,是亟待解决的问题.稳健统计方法是一种不需要筛选出数据中的异常值,也能获得合理模型的统计方法.稳健性这一概念最早由 Box17提出.自从 Tukey18发表了他的工作以来,稳健统计越来越受到许多统计学家的重视.Huber19提出了一类位置参数的稳健估计极大似然型估计,并解决了相应的渐近极小极大问题,进一步推动了稳健统计的发展.Hampel20,21在他的博士论文中给出了稳健性的严格定义,并提出了刻画稳健性的两个重要概念:影响函数和崩溃点.目前,稳健统计在线性回归模型中已经有大量研究成果,包括但不限于 M 估计22、分位数回归估计23、R 估计24、LMS 估计25、LTS 估计26、S 估计27、WLS估计27、MM 估计28、估计29、REWLS 估计30、Huber-ESL 估计31等.这些稳健方法可以很好的处理噪声数据,因此陆续有学者开始将稳健方法用于处理高维数据下线性回归模型的变量选择问题.为了便于讨论和理解,我们先对线性模型进行简单阐述.考虑样本量为 n 的随机观测(y1,x1),(y2,x2),(yn,xn),假设它们来自如下的线性回归模型:yi=xi+i,i=1,2,n,其中 xi=(xi1,xi2,xip)是 p 维协变量,=(1,2,p)是未知参数向量,i是随机误差项,E(i)=0,且误差项与协变量是相互独立的.在经典线性模型中,通常假定误差项服从正态分布.Fan 和 Li8提出了惩罚稳健回归估计的一般形式:bn=argminni=1(yi xi)+npj=1pnj(|j|),(1)第 1 期邹航,姜云卢:高维线性回归模型稳健变量选择方法综述159其中,()表示损失函数,pnj(|j|)是关于 j的惩罚函数,协调参数 nj越大,对回归系数的压缩程度就越大,如果某个变量的回归系数被压缩到 0,则该变量未被选入模型中.Fan和 Li32指出,在平方损失函数下,选用的惩罚函数需要使得惩罚最小二乘估计满足稀疏性、无偏性和连续性,从而达到变量选择的目的.Lv 和 Fan33提出了一系列满足这三条性质的惩罚函数.令 A=j:j=0,进一步假定|A|=s p,b=(b1,b2,bp)表示由变量选择过程得到的估计系数.在一定条件下,文献 8,34 研究了变量选择过程的Oracle 性质:(a)变量选择的相合性:PcA=A 1,(b)渐近正态性:n(bA A)N(0,A),其中cA=j:bj=0,A是真实模型下的协方差矩阵.变量选择的相合性表明真实模型被选出的概率随着样本量趋向于无穷时趋向于 1,渐近正态性表明估计出来的非零参数部分具有渐近正态性,且协方差矩阵与真实模型下的协方差矩阵相同.Oracle 性质是变量选择理论研究的主要内容.自从文献 8 提出了惩罚稳健回归估计的一般准则并研究了变量选择的 Oracle 性质后,后续的大多数稳健变量选择方法都沿用了这种形式,它们的区别在于不同损失函数与惩罚函数的选择.比如,Wang 等35提出了 LAD-LASSO 方法,其中损失函数和惩罚函数分别为(t)=|t|,pnj(|j|)=nj|j|.为了使得 LAD-LASSO 在杠杆点存在的情况下更稳健,Arslan36提出了加权 LAD-LASSO 估计.Wang 和 Li37研究了加权 Wilcoxon 类型损失函数的稳健变量选择方法.当损失函数是(t)=t I(t n 的情况下,关于高维线性回归模型的稳健且高效变量选择方法可以查阅文献 39,4951 等.160应用概率统计第 40 卷通过总结上述文献,本文主要从以下几个方面阐述高维线性回归模型稳健变量选择方法的相关工作.文章第 2 部分主要介绍评价稳健性的三个指标.紧接着,文章展开讨论响应变量含有异常值、响应变量和协变量都含有异常值、高崩溃点且高效的稳健变量选择方法,对应第 3 到第 5 部分.第 6 部分介绍惩罚变量选择模型的相关算法,第 7 部分,通过模拟和实证数据,比较不同的变量选择方法.文章的最后部分,我们进行总结和展望.2稳健性评价指标为了从理论上研究模型估计的稳健性,我们首先考虑估计的稳健性度量.目前,已经有不少学者研究了评估稳健性的三个指标:影响函数(influence function,IF)、有限样本崩溃点(the finite sample breakdown point)和最大偏差(maximum bias).Hampel20,21最先给出了稳健性的严格定义,并提出了刻画稳健性的两个重要概念:影响函数和崩溃点.随后,Huber52以及 Li 和 Zhang53研究了一维位置参数稳健估计的稳健性质.Zuo 等54研究了多维位置参数稳健估计的稳健性质,Zuo 和 Cui55研究了稳健散度估计的有限样本崩溃点.在经典的线性回归模型中,也有不少学者从理论上研究了参数估计的稳健性质.更多讨论,可以参考文献 28,30,31,45,56.下面,我们对稳健性的三个评价指标进行介绍.2.1影响函数考虑一个随机变量 X 具有概率分布 P,其中 P P=P:Rd,其分布函数为 F.在很多情况下,参数 常看作分布函数的一个泛函,即 =T(F),例如:均值、方差以及中位数等等.当真实分布稍偏离假设分布 F 时,往往通过度量给定分布 F 关于 T()的效果来测量泛函 T(F)的稳健性.Hampel20探讨了这类问题并引进了影响函数的概念.令 z表示在固定点 z Rd,d 1 的单点概率分布.给定在 Rd上的一个分布 F,对其进行单点污染,污染比例为 (0,1),得到 F 和 z的混合分布,记为 F=(1 )F+z.对于一个统计函数 T:F T(F)Rd,在给定点 z Rd处的影响函数定义为IFz;T,F=lim0+T(1 )F+z)T(F).从上式我们看到,影响函数测量了对分布 F 在点 z 处造成无穷小污染时对统计函数 T(F)的一种相对效果,它度量了 T()的局部稳健性;同时,影响函数的形式同导数相同,表示统计函数 T(F)的变化率.因此,影响函数越大,当分布函数 F 变为 F时,T(F)的变化也越大.所以,我们希望得到一个有界的影响函数.例如,对一元变量均值(F)与中位数 M(F),通过直接计算,我们得到它们的影响函数分别为IFz;(F),F=z (F),IFz;M(F),F=sign(z M(F)2f(M(F),第 1 期邹航,姜云卢:高维线性回归模型稳健变量选择方法综述161其中 sign()是符号函数.从上式我们可以看出,中位数的影响函数是有界的,而均值是无界的,所以中位数是一个稳健的估计.从这里可以看到,影响函数能度量一个估计稳健与否,是一个定性的指标.2.2有限样本崩溃点影响函数从无穷小概念出发,描述单个观测点对估计量的影响,刻画的是统计函数T()的局部稳健性,然而崩溃点衡量的是统计函数 T(F)的全局稳健性,是一个定量的指标.文献 20 最早从渐近角度研究了 T(F)的崩溃点.随后,Donoho 和 Huber57给出了两种有限样本崩溃点概念:有限增加样本崩溃点和有限取代样本崩溃点.接下来,我们首先介绍有限增加样本崩溃点的定义.令 Dn=X1,X2,Xn 是样本量为 n 的一组随机样本,T(Dn)表示基于样本 Dn的一个估计.把 m 个污染样本eDm=Xn+1,Xn+2,Xn+m加到 Dn,得到新的污染样本为 DneDm,其中污染样本的比例为 m/(n+m),则对估计T(Dn)的有限增加样本崩溃点定义为BP(T;Dn)=minmn+m:supeDmT(Dn)T(DneDm)=,(2)其中 表示欧拉范数.对于有限取代样本崩溃点,我们用 m 个污染样本eDm去替代样本 Dn的任意 m 个值,得到新的污染样本为bDn,其中污染样本的比例为 m/n,则对估计 T(Dn)的有限取代样本崩溃点定义为BP(T;Dn/eDm)=minmn:supbDn/DnT(Dn)T(bDn)=.(3)从式(2)和(3)我们可以看到,有限样本崩溃点测量了使得估计变得无穷大之前受污染样本的最小比例.很明显,有限样本的崩溃点越大,所得到的估计克服离群点的能力就越强,也就是说,所得到的估计就越稳健.在一般情况下,上述两个定义都可以使用.但是,Huber52指出,在考虑无结构的问题时,如:位置或尺度参数估计,有限增加样本崩溃点更容易处理.由式(2)可知,有限增加样本崩溃点的取值在 0 与 1 之间.对于一个常数统计量,有限增加样本崩溃点的取值为 1.对于一元随机样本的样本均值,其有限增加样本崩溃点的取值为 1/(n+1),而对样本中位数而言,其值为(n+1)/2,其中 x 表示不大于 x 的最大整数.所以,样本中位数相对于样本均值更稳健.2.3最大偏差在影响函数的定义中,仅仅考虑了单点污染的情况.而在实际情况中,假定分布的偏差可能是由任意分布的污染造成的.接下来,我们将介绍一类更为广泛的全局稳健性的一种162应用概率统计第 40 卷度量工具,即最大偏差.令 G 是污染分布,对一个固定的 0,(1 )F+G 是一个被污染的分布,则统计函数 T(F)的最大偏差定义为B(;T,F)=supGT(1 )F+G)T(F).B(;T,F)测量了最坏情况偏差.对很小的,如果 T(F)有一个适当的最大偏差曲线,则它被认为是稳健的.由于在很多情况下,B(;T,F)很难得到其具体形式,因此,在考虑估计的稳健性时,我们一般用影响函数与有限样本崩溃点度量其稳健性.在后面的讨论中,如果某个估计的的崩溃点达到 1/2,则称它具有高崩溃点,且具有理论上的稳健性质.3响应变量含有异常值的稳健变量选择方法为了在包含异常值的数据中选出对响应变量有贡献的解释变量,Fan 和 Li8提出了稳健变量选择方法的一般框架.经典线性模型中,通常假设误差服从正态分布,这在实际应用中往往很难满足.当误差服从厚尾分布或者样本中包含异常点时,基于最小二乘损失函数或者似然函数的变量选择方法是不稳健的.在这一部分,我们讨论响应变量含有异常值时的稳健变量选择方法.3.1惩罚最小绝对离差估计惩罚最小二乘估计在数据中含有异常值时往往是失效的,为了解决这一问题,Wang等35将惩罚最小二乘中的平方损失函数替换成最小绝对离差(least absolute deviation,LAD)损失函数,提出了 LAD-LASSO 估计,具体形式如下:argminni=1|yi xi|+npj=1j|j|.(4)协调参数 j可以通过 BIC 准则选择.为了求解上述估计,通过增加样本的方法,对原始数据(yi,xi),1 6 i 6 n 简单变换为(yi,xi),1 6 i 6 n+p:(yi,xi)=(yi,xi),1 6 i 6 n;(0,njej),1 6 j 6 p,n+1 6 i 6 n+p,其中 ej是第 j 个元素为 1,其它元素为 0 的 p 维向量.因此式(4)可以转化为argminn+pi=1|yi xi|.这样求解 LAD-LASSO 就转换成了求解传统的 LAD 估计.LAD-LASSO 方法对正负残差施加相同的权重,可以在重尾误差情况下实现稳健变量选择,并且具有 Oracle 性质.另外,第 1 期邹航,姜云卢:高维线性回归模型稳健变量选择方法综述163Gao 和 Huang58研究了在高维稀疏模型下 LAD-LASSO 估计的相合性和变量选择的相合性,Wang59提出的 L1惩罚 LAD 估计(LPLAD)在高维情况下对厚尾误差是有效的.当误差项服从正态分布且数据中不包含异常值时,LAD-LASSO 估计相对于最小二乘损失下的 LASSO 估计效率更低.为了解决这一问题,Lambert-Lacroix 和 Zwald60将惩罚 LAD 损失扩展到惩罚 Huber 函数,当残差较大时,Huber 损失函数等价于 LAD 损失函数;相反,则等价于最小二乘损失函数.与 LAD 损失函数相比,Huber-LASSO 估计保留了稳健性的同时,提高了估计的效率.3.2惩罚分位数估计考虑检验损失函数(u)=u(1(u c,c 0 时,惩罚 M 估计成为惩罚 Huber 估计;当(u)=u+(1)(u)+,0 0,未知参数bj可以用不加惩罚项的加权 LAD估计量替代.借鉴文献 35,在一定条件下,加权 LAD-LASSO 方法也被证明具有n 相合性、稀疏性和 Oracle 性质,具体参看文献 36.166应用概率统计第 40 卷4.2稳健惩罚秩估计为了使得估计量不受响应变量中的异常点和杠杆点的影响,Wang 和 Li37提出了如下的加权 Wilcoxon 型 SCAD(weighted Wilcoxon-type SCAD,WW-SCAD)方法来实现稳健估计和稳健变量选择:argminn1ijbij|ei ej|+npj=1p(|j|),其中 bij是对称的正权重,ei=yixi 是残差,p()是 SCAD 惩罚函数,调整参数 通过BIC 准则获得.当权重 bij是常数时,最小化 n2ijbij|ei ej|等价于最小化 Wilcoxon型散度函数:12ni=1R(yi xi)n+112(yi xi),其中 R(yi xi)表示 yi xi 的秩.为了求解 WW-SCAD 估计,受 LAD-LASSO 思想启发,首先对原始数据做如下变换:(yi,xi)=(bij(yj yi),bij(xj xi),1 6 i 6 n;(0,p(|0j|)j),1 6 j 6 p,n+1 6 i 6 n(n 1)/2+p,i j.这里,p()是惩罚函数 SCAD 的一阶导数,0j是不加惩罚项的加权 Wilcoxon 估计,j表示第 j 个元素为 1,其余元素为 0 的 p 维向量.在对数据简单变换以后,求解 WW-SCAD估计就转化成在新观测样本下求解 LAD 估计:argminn(n1)/2+pi=1|yi xi|.关于 WW-SCAD 方法的理论性质可以具体参看文献 37.4.3稀疏最小截尾二乘估计除了对残差较大的观测值降权来达到稳健目的,对大残差的观测值进行截断也可以作为备选.Alfons 等46提出了稀疏最小截尾平方回归(sparse LTS)估计,并基于崩溃点研究了估计的稳健性质.参数估计通过如下目标函数获得:bsparseLTS=argminhi=1(r2()i:n+hpj=1|j|,(7)其中 h 6 n,(r2()1:n6(r2()2:n6 6(r2()n:n是残差平方的次序统计量.为了求解稀疏最小截尾平方回归估计,式(7)等价的转化为:bH=argminQ(H,)=argminiH(yi xi)2+hpj=1|j|,第 1 期邹航,姜云卢:高维线性回归模型稳健变量选择方法综述167其中集合 H 1,2,n 且|H|=h,最优子集:Hopt=argminH1,2,n:|H|=hQ(H,bH).随后,稀疏截尾最小二乘回归估计由bHopt给出.实际上,稀疏最小截尾平方回归对应于寻找 h 6 n 的观测值的子集,使得其在 LASSO 拟合过程中产生最小的惩罚残差平方和,通常取 h=0.75(n+1),表示取整运算.为了寻找最优子集,文献 46 使用了一个类似于FAST-LTS73的算法.在一些正则条件成立时,文献 46 证明了稀疏截尾最小二乘回归估计量的崩溃点为(n h+1)/n,而 LASSO 和 LAD-LASSO 估计量的崩溃点仅为 1/n.除此之外,稀疏截尾最小二乘回归还具有许多优点,例如,在高维情况下,它可以通过减小方差来提高预测精度,保证了模型的可解释性;避免了传统的稳健统计方法直接推广到高维情况所带来的计算困难问题.4.4其他方法除了前面讨论的几种方法外,还有学者做了相关研究.例如,Gijbels 和 Vrinssen74考虑了三种基于非负绞索(nonnegative garrote)的稳健变量选择方法,即 M、S 和 LTS 非负绞索方法,研究表明,所提方法对响应变量中的异常值和杠杆点都具有稳健性.Loh41研究了高维惩罚稳健 M 估计的相合性和渐近正态性,并提出了一种复合梯度下降算法.Chang 等42基于 Tukey 双权重损失函数提出了 Tukey-LASSO 方法,并利用加速近端梯度算法解决了高维的计算问题.Wang 和 Karunamuni75基于 Lq损失(1 6 q 2)研究了高维线性模型的变量选择问题,并通过实例说明了方法的稳健性.5高崩溃点且高效的稳健变量选择方法正如 Ronchetti44指出,上述大多数方法虽然给出了变量选择的 Oracle 性质,但是很少有学者从理论上研究所提方法的稳健性质:即给出影响函数和有限样本崩溃点的理论结果,它们仅仅使用数值模拟来说明所提方法的稳健性,并且大多数考虑的是重尾误差情况,对杠杆点的情况涉及不多.在这部分,我们关注高崩溃点且高效的变量选择方法,即估计量从理论上被证明具有稳健性质,在数据中没有异常点且误差项服从正态分布时,估计量的渐近方差与惩罚最小二乘估计量的渐近方差之比为 1.为了从理论上研究稳健变量选择方法的稳健性质,Wang 等45提出了一种自适应的指数平方损失函数,考虑了稳健惩罚回归估计方法,并从理论上研究了所提方法的稳健性.文献 45 提出的基于指数平方损失函数的稳健变量选择方法如下:argmaxni=1exp(yi xi)2n npj=1pnj(|j|),(8)168应用概率统计第 40 卷式(8)中的极大化目标函数包含一类特殊的损失函数 :n(u)=1 exp(u2/n),其中n(0,+)是一个协调参数,下标 n 表示依赖于样本量,它用于控制估计的稳健程度.当n很大时,(u)u2/n,指数类型平方损失函数类似于平方损失函数.对于小的 n,残差较大的观测值对 的估计产生的影响很小,这样一来,就使得估计量具备了稳健性质.Wu 和 Ma76指出,文献 45 所提的方法是首次从理论上完全的描述了惩罚过程的稳健性质.至此之后,使得稳健统计理论成为研究的新热点,为稳健统计在高维数据分析提供了理论基础,吸引了一大批学者研究高维数据稳健变量选择的理论分析.Smucler 和 Yohai47提出了惩罚MM回归估计,并从理论上研究了估计的 Oracle 性质和稳健性质.惩罚 MM 回归估计如下:argminRpni=1(ri()sn(r(b1)+pj=1pnj(|j|),其中 是满足某些条件的函数,b1是参数真值的初始相合估计,sn(r(b1)是由 MM 估计方法计算得到的残差尺度估计量.当惩罚函数 pnj(|j|)=n|j|q时,惩罚 MM 估计变为MM 桥回归估计bB.当惩罚函数 pnj(|j|)=n|j|t/|b2,j|时,惩罚 MM 估计变为自适应MM 桥回归估计bA,其中b2是参数真值的初始相合估计,如果b2,j=0,则bA,j=0.当q 0时,文献47证明了MM桥回归估计和自适应MM桥回归估计具有Oracle性质.在惩罚参数 n 0 且固定时,MM 桥回归估计的有限样本崩溃点的下界为 1 1/n,自适应 MM 桥回归估计的有限样本崩溃点将以初始相合估计b2的有限样本崩溃点为下界.Kong 等77针对最小二乘损失下的线性模型提出了一种完全有效的稳健变量选择方法.这种方法在线性模型中增加了一个均值位移参数,来反应离群值的影响,并且对回归系数和均值位移参数都施加了惩罚项,以达到稳健变量选择和离群点检测的目的.如果某个响应变量是异常值,对应的均值位移参数不为 0,如果异常值出现在协变量中,那么对应的均值位移参数将是 0.在响应变量和协变量都含有异常值时,文献 77 从理论上证明了所提出的估计量崩溃点为 1/2,并且具有完全有效性和离群点检测一致性.除了从理论上考虑估计量的稳健性,估计量的有效性也需要考虑.Karunamuni 等48提出了稳健有效的双自适应正则回归方法,该方法对离群点和杠杆点同时施加自适应权重,并以此来作为损失函数的权重.当选取自适应惩罚函数,该估计量具有 Oracle 性质.当误差项服从正态分布,且数据中不包含异常点时,估计量将达到完全有效性.当响应变量和协变量都含有异常值时,估计量的崩溃点可以达到 1/2.关于稳健有效的变量选择方法,还有一些学者做了研究.如 Lozano 等78结合最小距离泛函和 l1正则项,提出了最小距离 LASSO 方法,该方法的稳健性由尺度参数控制.在模型误差分布的适当条件下,文献 78 证明了最小距离 LASSO 估计量能够达到最高的崩溃点和最快的收敛速度.Cohen Freue 等79研究了一种惩罚弹性网 S 估计,并从崩溃点的角度研究了估计量的稳健性.在多元线性回归模型下,Gijbels 等80考虑了稳健的 S 型非负绞索变量选择方法,并基于崩溃点和影响函数研究了该方法的稳健性.Cohen Freue 等81第 1 期邹航,姜云卢:高维线性回归模型稳健变量选择方法综述169基于弹性网提出了一种惩罚稳健估计方法,并在满足一些条件时,证明了所提估计量的稳健性和相合性.Ghosh 和 Majumdar51结合密度功效散度损失函数和惩罚函数研究了超高维稀疏线性模型的稳健且有效的变量选择问题.6算法实现在这一部分,我们主要介绍惩罚稳健变量选择方法的相关算法.惩罚变量选择方法是不同损失函数和惩罚函数的组合,当所优化的目标函数是损失函数和凸惩罚的组合时,如LASSO、自适应 LASSO 等,传统的凸优化算法可以直接应用.然而,当惩罚函数是非凸函数时(如 SCAD、MCP 等),和损失函数组合后的目标函数将不再是凸函数,传统的凸优化算法将不再适用.文献 8 提出了惩罚稳健变量选择方法的一般形式,在式(1)中,nj表示惩罚参数随着样本量和变量的不同而不同,为了表达的简洁,在后面的讨论中,令 pnj(|j|)=p(|j|).当惩罚函数为 SCAD、MCP 等一些非凸函数时,目标函数(1)将不再是关于 的凸函数,极小化将变得十分困难.为此,文献 8 利用二次函数局部近似目标函数,提出了局部二次近似算法(local quadratic approximate,LQA).给定参数的初始值 0,假设 0十分靠近参数真实值.如果 j0非常接近 0,令bj=0,否则,用线性函数局部近似惩罚函数的一阶导数:p(|j|)p(|j0|)+p(|j0|)sign(j0)(j j0)p(|j0|)+p(|j0|)/|j0|j(j j0)p(|j0|)+12p(|j0|)/|j0|(2j 2j0),(9)其中,j j0.通常选不加惩罚项的参数估计值作为初始值.利用式(9),我们可以得到第k 次参数更新的迭代公式:(k+1)=argminn()+npj=1p(|(k)j|)2|(k)j|2j,(10)其中,n()表示损失函数.利用 Newton-Raphson 算法可以求解式(10).当 k收敛时,迭代终止.文献 8 证明了算法的收敛性.LQA 算法虽然使得惩罚回归有了显式解,但它具有后向逐步回归变量选择方法的不足,在迭代过程中,一旦某个变量被筛选出模型,它将不会再被选入模型中.为了解决这一问题,文献 82 提出了在|(k)j|中添加扰动项 的 LQA 算法,但扰动项 的选取十分困难.为了解决 LQA 算法的不足,Zou 和 Li83提出了局部线性近似(local linear approxi-mation,LLA)算法,具体如下:p(|j|)p(|j0|)+p(|j0|)sign(j0)(j j0)170应用概率统计第 40 卷 p(|j0|)+p(|j0|)sign(j0)j sign(j0)j0 p(|j0|)+p(|j0|)(|j|j0|),其中,0为参数初值,j j0.类似于 LQA 算法,参数的迭代更新公式为(k+1)=argminn()+npj=1p(|(k)j|)|j|.当(k)收敛时,迭代终止.利用 LLA 算法,非凸惩罚回归就转变成加权惩罚 L1回归,权重为 p(|j0|).文献83 研究了 LLA 算法的收敛性,并说明了 LLA 算法比 LQA 算法更稳定.实际上 LQA 算法和 LLA 算法都是 MM 算法的特殊情况.MM 算法不是直接极小化原来的非凸函数,而是极小化比原函数更大的替代函数,且这个替代函数近似原来的非凸函数.MM 算法代表极大极小(majorize-minimize)或者极小极大(minorize-maximize)算法,取决于所要优化的目标函数.假设(k+1)是替代函数 ln()+nq()的极小值点,于是有ln(k+1)+np(k+1)6 ln(k+1)+nq(k+1)6 ln(k)+nq(k)=ln(k)+np(k),其中,原函数 p(k+1)和替代函数 q(k+1)的关系为 p()6 q(),p(k)=q(k).文献 82 研究 MM 算法的收敛性质.文献 69 利用 MM 算法求解了惩罚秩回归.为了解决检验损失函数在原点处的不可微问题,文献 61 将带有 LASSO 惩罚的分位数回归模型改写为线性规划形式:min0,ni=1i+(1 )ni=1is.t.pj=1|j|6 s,i6 yi 0pj=1jxij6 i,i,i 0,i=1,2,n.(11)基于线性规划形式(11)得到的拉格朗日原函数,文献 61 提出了一种计算整个解路径的有效算法.文献 39 基于修正的LLA算法利用线性规划问题来求解高维的惩罚分位数回归问题.对于非凸惩罚分位数回归问题,文献 38 将 SCAD 惩罚分解为两个凸函数的差,提出了差分凸算法.实验表明,差分凸算法的计算耗时比 LQA 和 MM 算法更小.对于惩罚最小二乘问题,Fu84和 Daubechies 等85提出了坐标下降算法进行求解.坐标下降法是一种非梯度优化算法,在每次迭代过程中,算法在当前点处沿坐标方向进行一维搜索,以求得函数的局部最小值.Friedman 等86利用坐标下降算法解决了诸如 LASSO和 LAD-LASSO 等凸优化问题,Wu 和 Lange87利用快速贪心坐标下降算法求解带有第 1 期邹航,姜云卢:高维线性回归模型稳健变量选择方法综述171LASSO 惩罚项的惩罚回归模型.此外,对于非凸惩罚最小二乘问题,Breheny 和 Huang88与 Mazumder 等89研究了相关的坐标下降算法并证明了它们的收敛性.在高维情况下,直接优化不光滑的分位数损失函数比较困难.Peng 和 Wang64提出了基于分位数的迭代坐标下降算法(QICD),将惩罚分位数回归模型简化为加权中位数回归问题.除此之外,还有一些学者做了相关研究.如,Tseng 和 Yun90、Wang 等45利用块坐标梯度下降算法来求解惩罚回归问题,Fan 和 Lv91利用路径跟踪坐标优化算法来优化非凸惩罚似然函数.7模拟和实例分析在这一节中,我们通过模拟和实例来比较几种常用的变量选择方法的有限样本性能,包括惩罚最小二乘(LS)、惩罚最小绝对离差(LAD)、惩罚复合分位数回归(CQR)、惩罚指数平方回归(ESL)和惩罚 MM 估计(MM).7.1数值模拟随机样本通过以下方式生成:yi=xTi+i,i=1,2,n,其中,真实参数 =(2,1,2,0,0,0,0,0,0,0),10 个解释变量中只有 3 个是重要变量,协变量 xi=(xi1,xi2,xi10)服从多元正态分布 N(0,),协方差矩阵 =0.5|jk|,(1 6j,k 6 10).在上述模型的基础上,我们考虑如下 4 种情况:(a)数据中没有异常值:误差项服从正态分布 N(0,1);(b)响应变量中含有异常值:误差项服从混合正态分布 0.8N(0,1)+0.2N(5,52);(c)协变量中含有异常值:协变量 xi=(xi1,xi2,xi10)服从混合正态分布 0.8N(0,)+0.2N(5110,52I10),110表示元素全为 1 的 10 维向量,I10表示 10 行 10 列单位矩阵,误差项服从正态分布 N(0,1);(d)协变量和响应变量都含有异常值:协变量和误差项的分布与情况(c)相同,除了在响应变量中添加 10%的杠杆点,这些杠杆点服从分布 N(10,52).分别重复模拟 4 种情况 500 次,每次生成的样本量 n 分别为 200,400,600,1000.为了评估每种变量选择方法,采用如下 4 个评价指标(重复 500 次的平均值):1)被正确识别为 0 系数的真实 0 系数的平均个数(NC);2)被误识别为 0 系数的非 0 真实系数的平均个数(NIC);3)正确选择模型的比例(PC);4)参数估计的均方误差(MSE)和对应的标准差(SD).对于每种方法,我们采用 R 语言中的程序包实现,惩罚项均选用自适应 LASSO 惩罚.LS-ALASSO,LAD-ALASSO,CQR-ALASSO,ESL-ALASSO,MM-ALASSO 分别由172应用概率统计第 40 卷glmnet,rqPen,cqrReg,robustlm,mmlasso 包计算得到.其中,CQR-ALASSO 中的复合分位数水平 =k/10,k=1,2,9,所有方法的参数通过程序包自适应选择得到.模拟结果总结在表 1 至表 4 中.表 1情况一的模拟结果n变量选择方法NCNICPCMSE(SD)LS-ALASSO6.99500.00000.96500.0307(0.0318)LAD-ALASSO6.38500.00000.83500.0484(0.0503)200CQR-ALASSO6.95000.00000.97500.0425(0.0474)ESL-ALASSO6.97500.00000.96500.0465(0.0528)MM-ALASSO6.85290.00000.85290.0366(0.0331)n变量选择方法NCNICPCMSE(SD)LS-ALASSO7.00000.00001.00000.0116(0.0109)LAD-ALASSO6.73000.00000.90000.0233(0.0211)400CQR-ALASSO6.99500.00000.99500.0182(0.0172)ESL-ALASSO7.00000.00001.00000.0245(0.0214)MM-ALASSO7.00000.00001.00000.0154(0.0133)n变量选择方法NCNICPCMSE(SD)LS-ALASSO7.00000.00001.00000.0034(0.0048)LAD-ALASSO6.84500.00000.93500.0062(0.0084)600CQR-ALASSO7.00000.00001.00000.0052(0.0061)ESL-ALASSO7.00000.00001.00000.0068(0.0080)MM-ALASSO7.00000.00001.00000.0049(0.0069)n变量选择方法NCNICPCMSE(SD)LS-ALASSO7.00000.00001.00000.0024(0.0032)LAD-ALASSO6.99500.00000.95000.0041(0.0058)1000CQR-ALASSO7.00000.00001.00000.0040(0.0038)ESL-ALASSO7.00000.00001.00000.0035(0.0044)MM-ALASSO7.00000.00001.00000.0028(0.0041)从表 1 的结果来看,当误差服从正态分布,且数据中不包含异常点时,所有的方法都可以选择出正确的模型和重要变量,随着样本量的增加,所有方法的 MSE 都会越来越小.其中,LS-ALASSO 的 MSE 是最小的,LAD-ALASSO 的 MSE 是最大的,CQR-ALASSO 的MSE 会略小于 MM-ALASSO 和 ESL-ALASSO.当误差项服从混合正态分布,LS-ALASSO虽然正确识别了真实零系数(NC),但是没有识别出真实的非 0 系数(NIC)和真实模型(PC),所得到的 MSE