1. 引言
在数字化发展的今天,金融、气象、生物医学等诸多领域的数据中各指标之间的关系变得愈加错综复杂,简单地运用线性模型已经不足以分析各指标的联系。而半参数模型不仅具有强大的拟合能力,且不需要对数据做出过多的假设,十分适用于处理分布未知的复杂数据。部分线性模型 [1] 则是常见的一类半参数模型,因其同时具有线性模型和非参数模型的优点,吸引了广大统计学者的关注,成为了各领域处理数据的重要工具。
而随着大数据时代的到来,高维数据频繁出现在各个领域。高维数据的变量个数过多,往往在建模分析过程中容易导致“维数灾祸”,传统的变量分析方法已无力应对这种问题。对此,可同时进行变量选择和系数估计的正则化方法被学者们提出,其一般形式为“损失函数 + 惩罚函数”。当采用岭回归 [2] 作惩罚函数时,其通过有偏估计很好地解决了变量间多重共线性的问题,但这种方法只能限制模型回归系数尽可能地小,不具备稀疏性。LASSO [3] 惩罚虽然能实现稀疏变量的目的,但对于具有群组效应的特殊数据结构,它只能勉强选出该变量组的一个变量,容易忽略其他重要变量,对于高维数据中存在多重共线性时,LASSO的估计效果会变得很差。SCAD [4] 惩罚是对称且非凸性的,它虽然可以产生稀疏解,但其迭代算法却极其缓慢,难以适应高维数据情形。而Zou和Hastie [5] 提出了弹性网络法,因其结合了LASSO和岭回归惩罚项的优势,不仅可以解决多重共线性的问题,而且还能对变量进行有效筛选。
高维数据下的部分线性模型因其需要考虑线性协变量的选择和系数估计、光滑参数的选择、非参数函数的估计等方面的问题,因此如何在高维数据下对部分线性模型进行稳健的变量选择成为了一大挑战。目前,关于高维部分线性模型的变量筛选问题研究还不多。Xie和Huang [6] 运用SCAD惩罚来实现线性部分的稀疏性,并用多项式样条来估计非参数分量。杨宜平等 [7] 结合样条方法和Dantzig或Lasso变量选择方法,同时进行变量选择和未知参数估计,并证明了估计误差的非渐近界。Chen等 [8] 采用自适应弹性网对部分线性模型进行变量选择,并证明了其Oracle性质。Wang等 [9] 利用组Bridge和最小二乘估计讨论了高维部分线性模型的估计问题。赖秋楠等 [10] 提出了profile贪婪向前回归变量筛选方法,证明了在一定正则条件下所提PGFR方法具有筛选相合性。但这些文献大多都是基于最小二乘估计作为损失函数,而在实践中最小二乘估计对噪音较敏感,稳定性较差,特别是在面对高维数据中含有极端值、离群点或是呈现尖峰、厚尾分布时,其预测效果往往不佳。为了实现更稳健的估计,Koenker和Basset [11] 提出了分位数回归,因其不需要对误差分布作特定假设仍可获得全局分布的特征,且保持较高的稳健性,进而被不同领域广泛研究运用。
本文就将利用将弹性网络法和分位数回归相结合的正则化方法,对高维部分线性模型进行稳健变量选择。此方法不需要对模型中误差项的分布情况做任何假设,在面对高维数据中出现异常值或极端值时仍保持较好的稳健性,能够丰富具体地反应响应变量的分布特征。除此之外,这种组合的正则化方法还可以灵活处理高维数据中多重共线性与群组效应问题,能对变量间的强相关性信息进行变量筛选,产生稀疏解。理论上将证明该模型估计量的相合性和稀疏性,并通过数值模拟,与其他变量选择方法进行比较,表明该方法更有效更稳健。
2. 稳健变量选择过程
考虑如下部分线性模型:
, (1)
其中,
为n维响应变量;
为
设计矩阵,且
;
为p维未知参数向量;
为n维随机误差向量,与
独立且期望为0,方差为
;
为未知光滑函数,为避免维数灾难,假定Z为一维随机变量,且在闭区间
上取值。
对模型(1)两边关于Z求条件期望,可得
, (2)
由(2)式可得
的表达式,再将其代入(1)式,经整理得
, (3)
若
和
已知,则模型(3)即为标准的线性模型。
记
,
,令
和
分别为
和
的估计量,本文采用核估计,即
;
,
其中
为核函数,h为窗宽。
记
,
;
,
,
则式(3)可改写为
. (4)
在面对高维数据时(
),为了确保模型的可识别性,提高模型拟合的精确度,假定真实的系数向量
为稀疏的,即仅小部分系数非零。在不失一般性的情况下,我们定义
,用s表示真实回归系数的非零元素的数量,假定真实模型表示为
,
其补集
表示噪音协变量,即
的前s个分量是非零的,设计矩阵可写为
,
其中,
,
;
表示非零系数的信号协变量矩阵,
表示系数为零的噪音协变量矩阵。本文将设计矩阵
的每一列标准化使得
,
。
结合(4)式和线性模型的弹性网惩罚分位数估计 [12] ,可得求解
的模型如下:
,
其中,
,
为分位数回归的损失函数,
为
的
范数,
为
的
范数,
为弹性网惩罚项,其中
为正则化系数。
3. 统计性质
本节将建立高维部分线性模型参数
的弹性网估计量
的性质。参照文献Fan [13] ,定义预先知道协变量位置的Oracle正则化估计
,并且满足
,并做出如下假设条件:
(A1)
独立于
,且满足
,
。
(A2) 最佳窗宽
。
(A3)
是一个正定矩阵,
,
。
(A4) 核函数
的支撑集为
,存在常数
,核函数满足
,且
,
,
。
(A5)
,
和
的二阶导数都是一阶Lipschitz连续的。
(A6) 存在常数
和
,对任意满足
的u,
一致有界且不为0和
,有
,
其中,
和
分别为误差
的密度函数和分布函数。
(A7) 定义
,则
的特征值介于
与
之间,且
.
(A8) 令
,
为大于零的常数,有
,
其中,对于矩阵
与向量
,
,
且对某些常数
,有
。
条件(A1)~(A5)是研究半参数模型常用的假设条件,其中条件(A2)为核估计最优窗宽,条件(A4)是对核函数的一般假定。条件(A6)对噪音分布进行假设,任意
在0点周围是Lipschitz连续的,如常用的柯西分布就满足这一条件。条件(A7)对信号协变量矩阵
和设计矩阵
的规模进行限制,当设计矩阵
产生某种分布时,
的界以渐近1的概率满足。条件(A8)对设计矩阵
与子矩阵
的列向量的相关性进行控制。
下面给出定理说明可用Oracle信息来确定信号协变量的位置,正则化估计
以趋于1的概率可估计出真实系数向量
。
定理1 若
,
,且满足条件(A1)~(A7),则存在常数
,使得
,
其中,
,
为大于零的常数。
定理2 假设条件(A1)~(A8)成立,若
,
,
,
且
,c为大于零的常数。那么,存在目标函数
的全局最小值
,至少以
的概率满足:
1)
,
2)
。
4. 定理证明
在证明上述定理前,本节先给出几个引理。
引理1 (Chen [14] )假设条件(A1)~(A5)成立,令
,则有
.
引理2 (Fan [13] )假设条件(A7)成立,则对任意
,都有
.
引理3 (Bühlmann [15] )霍夫丁不等式:设独立随机变量
来自某个空间
,
为
上的实值函数,若存在正实数
,使得对任意的i满足
,
,则对所有的
,有
.
引理4 考虑以
为中心的球形领域
内:
,
其中序列
,假设
,
,
,
则当条件(A6)~(A8)成立时,有
,
其中,
,
为大于零的常数。
定理1的证明:
令
,则目标函数可写作
.
对于给定的一个常数
,定义集合
,
然后,定义函数
.
首先,结合引理1可知,对任意
,当n充分大时,有
, (5)
其中,
,
是
在0的某个领域内的一个下界。
接着证明不等式(5),令
,则对于任意
,根据Hölder不等式有
.
此时,再针对
的符号进行讨论,当
时,依次根据
,Fubini定理,泰勒展开及条件(A6)可得
(6)
其中,
对于
一致成立。
当
时,可得式(6)右侧同样的结果。进一步,由条件(A7)得
(7)
因此,结合式(6)和式(7)以及
的定义,对于任意
,可得式(5)成立。但所设定的估计器
可能不在
集合中,因此,令
,其中
,
,
此时
,再根据目标函数为凸函数的性质以及
的定义,可得
. (8)
结合式(8)与三角不等式,有
(9)
定义事件
,
根据引理2,可得
.
又由Cauchy-Schwarz不等式,有
,
,
其中,
。因此,在事件
上,根据不等式(9)有
.
取
,由假设
,
和条件(A7)可知
,再结合式(5),则在事件
上有
,
进一步,得
.
根据
,可得
,故在事件
上,有
,
可知事件
发生的概率大于等于
,定理1得证。
定理2的证明:
此处先给出引理4的证明,其是定理2证明的基石。对于固定的
,
。定义
,
其中,
是设计矩阵
的第i行。
首先,需对引理4做如下分解:
根据上式可知,若能证明以下各式至少以
的概率成立,那么引理4可得证。
, (10)
, (11)
, (12)
首先证明
成立,将式(10)改写为
. (13)
根据条件(A6),有
,
其中,
是
的累计分布函数,
,
对任意
,有
.
结合式(13)及Cauchy-Schwarz不等式,可推出
. (14)
依次考虑不等式(14)右边的两项,由条件(A8)可得第一项的上界被限制为
, (15)
再根据条件(A6)~(A7)及
,可得第二项的上界被限制为
.
又因为
,由假设条件
有
, (16)
再将上述不等式(16)与式(15)代入式(14),可得
成立。接下来证明
,根据引理3,如果
,则有
,
因此,
至少以概率
成立。
现将利用Bühlmann [15] 中的推论14.4证明
成立。首先,对每个固定的j,定义函数空间
,则对任意的
,有
。又因函数
有界,可得
.
此外,对任意的
和
,由条件(A6)和均值定理可得
(17)
其中,
介于
与
之间。又因
一致有界,结合条件(A7)和式(17),有
(18)
其中,C是大于零的常数。
由Bühlmann [15] 中引理14.27可知
空间的球
可被
个半径为
的球所覆盖,又因
只取3种结果,结合式(18)可得
的覆盖数为
。因此,
对任意的
,有
.
又由Bühlmann [15] 中的推论14.4可知,对任意的
,有
,
令
,其中
为充分大的常数,则可得
.
因此,若
,那么
至少以
的概率成立,从而引理4得证。
下面接着证明定理2,由于
是
的最小值点且满足KTT条件。若要证明
是目标函数
的全局最小值,只需验证如下条件
, (19)
其中,对任意n维向量
,有
,则根据KKT条件及
为凸函数的性质,可知
是
的全局最小值点。定义事件
,
,
其中,
同定理1中所定义,并且
,
那么,结合定理1和引理4可知
.
又因为
属于事件
,故在事件
上不等式(19)成立,定理2得证。
5. 数值模拟
本节将通过数值模拟来考察所提出方法的有限样本性质。从如下部分线性模型产生数据
,
其中,真实的回归参数被固定为
,参数数量
,样本容量
,协变量
,其中
;非参数函数
,其中T服从
上的均匀分布。
为了验证所提方法关于厚尾分布的稳健性,考虑如下误差分布:标准正态分布
;混合正态分布
;自由度为4的t分布
;Laplace分布;Cauchy分布。本节将在这五种误差分布下,将所提的弹性网惩罚分位数回归(Q-EN)与弹性网惩罚最小二乘回归(EN)、Lasso惩罚分位数回归(Q-lasso)和Ridge惩罚分位数回归(Q-ridge)进行比较,其中
。通过计算以下4个值来评估所提方法:
1) L1 loss:即
的L1范数。
2) L2 loss:即
的L2范数。
3) FP:选入模型的噪音协变量的数量。
4) FN:未选入模型的信号协变量的数量。
程序重复运行100次,模拟结果如表1所示,其中各项数值为其均值。

Table 1. Simulation results of relevant covariates
表1. 相关协变量模拟结果
比较Q-EN和EN:在这五种误差分布情形下,Q-EN的L1损失和L2损失都较EN更低。当误差分布为MN1时,Q-EN比EN的FN更小。当误差分布为t4和Cauchy时,Q-EN较EN模型具有更小的FP。当误差分布为Laplace和
时,EN的FP和FN都优于Q-EN。总体而言,在面对具有厚尾分布和异常值,如t4等情形时,分位数回归较最小二乘估计具有更稳健的特性。
比较Q-EN和Q-lasso、Q-ridge:在前四种误差分布情形下,Q-EN的L1损失和L2损失都明显更低。当误差分布为Cauchy时,虽然Q-ridge具有比Q-EN更低的L2损失,但它选择了模型中的所有变量,并没有进行有效的变量选择。当误差分布为
、MN1、Laplace和t4时,Q-lasso具有更小的FP,表明Lasso惩罚项倾向于将真实有效的信号协变量选择到模型中。当误差分布为
、MN1和Laplace时,相比Q-lasso,Q-EN具有更小的FN,说明Q-EN剔除大量不相关的预测因子,表现更好。
通过在不同的误差分布情形下,与其他变量方法进行比较,表明Q-EN的总体性能优于其他四种方法。Q-EN可以有效的进行变量筛选产生稀疏解,且在面对厚尾分布或异常值时仍表现优异。
6. 总结
针对高维数据下的部分线性模型变量选择问题,提出了将弹性网络法和分位数回归相结合的正则化方法。该方法可以有效地处理高维数据中的多重共线性与群组效应问题,还可以在面对数据含有极端值、离群点或是呈现尖峰、厚尾分布时,仍保持较好稳健性。理论结果表明,在一定条件下模型估计量具有相合性和稀疏性,可将不相关变量压缩至零,从而进行有效的变量选择。在数值模拟方面,通过比较不同惩罚函数和不同损失函数在五种误差分布下的估计效果,证实了所提估计方法优于其他变量选择方法,具有更强的稳健性和有效性。在实践应用中,所提模型方法可为生物医学、基因学等领域进行高维数据下的实证分析提供理论支撑。利用此模型方法对高维基因学数据进行基因检测分析,探究与疾病相关程度高的协变量基因,为药物研究提供方向,可推动相关药物研发领域的发展。
基金项目
重庆市自然科学基金(CSTB2022NSCQ-MSX1370)。
NOTES
*通讯作者。