高维纵向数据的惩罚expectile 估计①

樊梅红， 李婷婷

西南大学 数学与统计学院，重庆 400715

在统计学中, 有很多方法可以寻找数据间的潜在关系, 刻画数据结构． 文献[1]提出的广义估计方程(generalized estimating equaiton, GEE)是很常见的一种分析纵向数据的统计方法, 在研究数据内部关系及预测等方面有重要的作用． 此方法的优点是即使相关结构被误判, 所得估计仍然是相合的． 在大数据时代, 高维纵向数据能比时序数据和横截面数据提供更多的信息． 然而数据的高维性使模型变得复杂, 降低了模型的估计精度． 带惩罚项的正则化估计是解决高维数据的常用方法． 文献[2]提出的SCAD惩罚和文献[3]提出的MCP惩罚是常见的非凸惩罚方法, 具有Oracle性质． 文献[4-5]将GEE与惩罚函数相结合, 提出了惩罚广义估计方程(penalized generalized estimating equation, PGEE), 模拟研究表明该方法在筛选出重要变量的同时得到模型回归系数的无偏估计．

在实际应用中, 数据往往会呈现异质性． 文献[6]首次提出分位数回归(quantile regression, QR)方法, 可以捕捉整个条件分布的特征． 文献[7]基于独立的数据结构提出纵向数据的线性分位数回归模型． 这不可避免地会损失估计效率． 进一步, 文献[8]考虑纵向数据重复观测样本间的相关性, 建立分位数GEE回归模型, 提高了估计效率． 文献[9]对纵向数据的分位数回归模型添加惩罚项, 提出了惩罚分位数回归模型．

QR方法对应的损失函数具有不可微性, 这给数值计算带来了很大的难度, 尤其对于高维复杂数据来说, 该问题变得更加突出． 受分位数回归的启发, 文献[10]将分位数回归中的非对称绝对值损失函数替换为非对称最小平方损失函数, 提出了期望分位数(expectile)估计量． Expectile方法不仅继承了QR方法可以处理异质性的优点, 且具有连续可微的损失函数, 相较QR方法在计算上也有很大的优势． 在独立同分布的截面数据中, 文献[11-12]将expectile回归与惩罚函数相结合, 提出带有惩罚项的expectile回归模型, 建立了Oracle性质, 同时实现了变量选择和异方差识别． 文献[13]将expectile应用到纵向数据, 提出了广义expectile估计方程(generalized expectile estimating equation, GEEE)． 模拟结果显示, GEEE估计量可以识别出异方差, 在保留分位数优点的同时, 降低了计算难度． 近年来, 作为QR方法的替代, expectile方法受到部分学者的关注, 但在纵向数据变量选择方面的研究还不多见． 本文将截面数据的惩罚expectile回归模型扩展到纵向数据, 提出PGEEE(penalized generalized expectile estimating equation)估计量． 模拟结果和实证分析显示, PGEEE估计量不仅可以实现高维数据的变量选择, 并且同时为重要变量的回归系数进行估计． 更重要的是, PGEEE方法可以得到一系列τ水平下的变量选择和模型估计结果, 详细地刻画了数据的异质结构, 能够比GEE提供更多的信息．

1 纵向数据的expectile回归和PGEEE模型

1.1 Expectile回归模型和PGEEE

定义随机变量Y的τ-expectile值为

(1)

其中τ∈(0, 1),ρτ(θ)=|τ-I(θ≤0)|·θ2是非对称平方损失函数, I是示性函数． 由τ-expectile的定义易知, 当τ=0.5时,ρτ(·)等价于经典的最小二乘损失函数, 则模型(1)对应经典的均值回归模型,μτ(Y)为随机变量Y的数学期望．

假设有纵向样本数据(yij,Xij),i=1,…,n,j=1,…,mi, 满足如下的expectile线性回归模型

(2)

yi=Xiβn+εi

(3)

对βn的估计可以通过求解如下目标函数的最小值来获得, 即

(4)

考虑重复观测时个体内的相关性, 文献[13]在纵向数据协变量数固定的情况下提出了GEEE模型, 即通过求解如下估计方程

(5)

进一步地, 本文在协变量维数pn发散的情况下, 提出纵向数据的惩罚非对称最小二乘PGEEE估计, 即通过求解如下估计方程

Q(βn)=S(βn)-nP′λn(|βn|)Sign(βn)=0

(6)

获得系数βn的PGEEE估计． 其中,P′λn(|βn|)=(p′λn(|βn1|), …,p′λn(|βnpn|))T,pλn(t)是一个含有调节参数λn的非负惩罚函数,p′λn(t)为pλn(t)的导数． Sign(βn)=(sign(βn1), …, sign(βnpn))T, sign(t)=I(t>0)-I(t<0)为符号函数．P′λn(|βn|)Sign(βn)定义为对应元素相乘得到的向量． 本文考虑MCP和SCAD两种惩罚方法． MCP惩罚函数的数学表达式为

(7)

为简化模型, 参考文献[14], 取γ=3． SCAD惩罚函数的数学表达式为

(8)

根据文献[2]建议取γ=3.7． 此时模型(6)中需要选择的参数只有λn, 本文使用BIC准则来选取, 表达式见算法过程．

1.2 求解算法

(10)

Step4: 重复Step2-Step3直至收敛, 并计算λn对应的BIC值, 其表达式为

(11)

其中,df表示λn对应模型所选择的变量个数．

2 Oracle性质

注定理1表明所提出的方法可以选出正确的模型, 同时实现对重要变量回归系数的参数估计, 称为Oracle性质[2]．

定理1的证明:

(12)

成立即可． 根据表达式, 有

(13)

(βn-βn0)TS(βn)=(βn-βn0)TS(βn0)+(βn-βn0)T[S(βn)-S(βn0)]=I1+I2

(14)

其中

(15)

考虑I11, 有

(16)

(17)

且

(18)

|I1|=Op(pn)‖u‖

(19)

将I2分为两部分计算, 有

(20)

记

(21)

其中由(A3)知

又因为

(24)

(25)

其中

又

(28)

I2=-Op(pn)‖u‖2

(29)

由(19),(29)式可得, (14)式的值由(29)式控制, 小于0． 易知(13)式中的第二项以nαn2‖u‖ +nbnαn2‖u‖2为界, 因此可以找到一个足够大的D, 使得(13)式的值完全由(29)式决定． (12)式得证．

(30)

(31)

由(A7)可知, (31)式的符号完全由βj的符号决定． (30)式得证．

即

(32)

(33)

定理证毕．

3 模拟研究

为了研究所提方法的有限样本性质, 本文比较了不同的惩罚方法及相关结构下所提出方法的效果． 数据来源于以下模型

(34)

情形1pn=10,k=9,mi=4,n=50, 100, 200,βn=(-3, 5, 0, 0, 4, 0, 0, 2, 0, 0)T．Ri是参数为0.9的等相关结构矩阵．

情形2k=2,mi服从参数为(3, 6)的均匀分布,Ri是参数为0.9的AR(1)结构矩阵． 其余设置和情形1一样．

情形3pn=30,n=100, 200．βn=(-3, 5, 0, 0, 4, 0, 0, 2, 0, 0, …, 0)T． 其余设置和情形1一样．

表1 情形1模拟结果

表2 情形2模拟结果

表3 情形3模拟结果

(i) SCAD和MCP两种惩罚方法并无明显的优劣之分． FN均为0, 表示所有重要的变量都被识别, FP接近0, 表明噪音变量被选择的可能性很小;

(ii) 在情形1和情形3中,τ=0.9时, Prob等于1, 而τ=0.5时, Prob的值接近0． 这表明所提出的估计量PGEEE可以在不同的τ水平下, 有效识别出正确的模型, 刻画数据中的异方差结构;

(iii) 在不同的τ水平下, 即使选择的变量相同, 参数估计值也可能不同(情形2)． 在此情形下, 估计量的MSE和MAE随着样本量增大而减小, 表示该方法可以在识别出异方差的同时实现回归参数的一致估计;

(iv) 对比情形1和情形3, 协变量维数pn从10增加至30, 结果显示模型中噪音变量数量增加时, PGEEE估计表现依然较好, 且估计量MSE减小, 表明该方法可以用于分析高维数据, 排除无关变量, 识别出重要变量．

(v) 考虑相关结构时估计量的表现总体上优于独立(IND)的情形． 即使相关结构被误判后, 参数估计效果依然很好, 尤其使用UN结构时．

4 实证分析

数据来自1976年至1982年间对美国经济收入动态的面板研究, 包含了连续7年595名民众的工资水平, 属于平衡数据, 更多详细信息参考文献[15]． 该研究中, 协变量包括工作经历E, 工作时间W, 工作职业O(蓝领取1, 否则0), 工作行业I(制造业取1, 否则0), 居住地S(居住在南部取1, 否则0), 种族B(黑人取1, 否则0), 是否住在都市统计区A(如果是取1, 否则0), 是否结婚M(结婚取1, 否则0), 性别F(女性取1, 否则0), 劳动保障U(签合同取1, 否则0) 及受教育程度D, 响应变量为对数变换后的工资水平．

表4给出了τ=0.01,0.5,0.95下参数的PGEEE估计, 其中τ=0． 5对应经典的均值回归估计． 分析结果可知, 不同的惩罚方法和不同的相关结构选出的变量基本一致． 可以看到, 在3个水平下均被选择的变量有O,B,F,D; 均未被选择的变量有W． 截距项,B,F的系数估计随着τ不同而变化, 图1a,b为不同种族及性别对应的工资随时间变化的箱线图． 男性的工资明显高于女性, 白人的工资明显高于黑人． 在τ=0.01时,E被认为是噪音变量, 而在τ=0.5和0.95时被认为是重要变量． 在τ=0.95时, 除了独立结构下MCP估计外, 工作行业I, 居住地S, 是否结婚M, 劳动保障U均被剔除在模型外; 而在τ=0.01 和0.5时则被认为是重要变量． 图1c,d,e,f为这些变量对应的工资分布箱线图． 以变量S为例, 可以看到, 在低分位点时, 居住在北部的工资要明显高于南部, 但是在高分位点时, 两者的区别并不明显, 这与PGEEE的估计结果相吻合． 由此可见, 该方法比采用普通最小二乘估计(τ=0.5)挖掘出了更多的信息．

图1 工资箱线图

表4 工资数据参数估计结果

5 结语

本文基于expectile提出了高维纵向数据的PGEEE估计量, 在实现模型变量选择的同时, 对模型的回归系数进行估计． 在正则条件下本文建立了PGEEE估计量的Oracle性质． 数值模拟结果显示, MCP与SCAD惩罚及不同的协方差结构在变量选择方面并无明显差异． 相较于独立结构, 考虑相关结构时回归系数的估计效率更高． 多数情况下, 不确定结构(UN)的PGEEE估计量具有较好的估计精度． 最后建立工资数据的PGEEE模型, 可以看到在不同的τ水平下, 影响工资的因素有所区别, 同一个因素影响程度也可能不同． 这表明PGEEE可以有效识别数据中的异质结构, 比经典的惩罚估计方程估计(PGEE)挖掘出更丰富的信息, 更合理地分析了工资的影响因素．