高维部分线性可加稳健Expectile回归模型

Partially Linear Additive Robust Expectile Regression in High Dimension

摘要: 高维数据一般因具有异方差或非齐次协变量而具有异质性，分位数回归和expectile回归是分析异质高维数据的有力工具，但前者由于损失函数非光滑的特性在计算方面存在较大挑战，而后者会因异常值而不稳健。本文利用一类稳健的非对称损失函数来研究部分线性可加模型的稳健expectile回归，用B样条基函数近似非参数部分，利用加入非凸惩罚的正则化方法来实现变量筛选并进行参数估计。该方法的优势在于：(1) 通过取不同分位水平得到响应变量更完整的条件分布，从而探索数据的异质性分布；(2) 部分线性的模型结构兼顾了线性解释变量和非线性解释变量，一方面增加了模型的灵活性，同时也具有一定的模型可解释性；(3) 稳健expectile回归估计比分位数回归方法计算效率高，比expectile回归稳健。数值模拟和实际数据分析均显示了该方法在模型估计和计算效率上的优势。

Abstract: High-dimensional data are generally heterogeneous due to heteroskedasticity or non-homogeneous covariates. Quantile regression and expectile regression are powerful tools for analyzing heterogeneous high-dimensional data, but the former is a great challenge in calculation due to the non-smooth nature of the loss function, while the latter is unstable due to outliers. In this paper, a class of robust asymmetric loss functions is used to study the robust expectile regression of partial linear additive models, the B-spline basis function is used to approximate the non-parametric part, and the regularization method with non-convex penalty is used to realize variable screening and parameter estimation. The advantages of this method are: (1) A more complete conditional distribution of response variables can be obtained by taking different quantile levels, so as to explore the heterogeneity distribution of data; (2) The partial linear model structure takes into account both linear explanatory variables and nonlinear explanatory variables, which increases the flexibility of the model on the one hand, and has a certain interpretability of the model; (3) The robust expectile regression estimation score digit regression method has higher computational efficiency and is more robust than the expectile regression. Both numerical simulation and actual data analysis show the advantages of the proposed method in model estimation and computational efficiency.

1. 引言

半参数模型拥有线性模型的直观性和非参数模型的灵活性而常被作为处理高维数据的重要工具[1] [2]。本文主要研究半参数模型中的部分线性可加模型(partially linear additive models，简记为PLAMs)，是广义可加非参数回归模型的特殊类型和多元线性回归模型的推广[3]。PLAMs的估计和推断理论已有丰富的研究结果[4]-[7]。

在实际研究中，高维数据常有异质性表现，噪声变量的分布呈现重尾型或非对称尾型。分位数回归(quantile regression)是有效分析异质性的方法，在PLAMs框架下已有很多研究成果产出。Tadao [8]于2014年用核密度法近似非参数部分的部分线性可加分位数回归，使用两步估计的方法并得出估计值的渐近性质；Sherwood和Wang [9]于2016年提出超高维部分线性可加分位数回归，用带惩罚的分位数回归分析超高维数据的异质性并证明了估计量的渐进正态性。

虽然分位数回归受到研究者们的普遍关注，但其损失函数不光滑，在高维情形的计算负担过重。而非对称的最小平方回归(expectile regression) [10]因其具有光滑的损失函数，使得计算效率大幅提升，是分析异质性数据的另一重要工具，同时它也是风险分析领域的常用工具，在金融风险领域中受到关注[11] [12]。Fabian [13]等人于2013年提出部分线性可加expectile回归并得出估计值的渐近性质和Bootstrap置信区间；Zhao [14]等人于2019年提出高维部分线性可加expectile回归，并研究模型的估计与变量选择。

虽然expectile regression比quantile regression计算速度快，却易受异常值的影响造成较大的估计误差，远不如quantile regression稳健。可见expectile regression和quantile regression都是分析数据异质性的典型方法，但都存在不同的缺点。Man等人[15]等人2022年提出一种稳健的expectile回归，用Lipschitz连续局部平方的函数代替expectile regression损失函数的平方项，并调节能随数据自适应的参数$\gamma$ 让偏差与稳健达到平衡，解决了quantile regression计算速度慢和expectile regression不稳健的问题又兼具了二者的优点，于是本文将此种方法推广到更加灵活的PLAMs。

本文主要研究高维部分线性可加稳健expectile回归模型。全文内容安排如下：第2节详细介绍部分线性可加稳健expectile回归模型(retire-PLAMs)，第3节介绍带非凸惩罚的高维部分线性可加expectile回归及算法，第4节是数值模拟，第5节是实证研究，第6节是总结。

2. 高维部分线性可加稳健Expectile回归

2.1. 一类非对称的稳健平方损失

Newey和Powell [10]受quantile regression的启发于1987年提出expectile regression。设$Z\in ℝ$ 是一个具有有限矩的随机变量，同quantile回归的中$\tau \text{-th}$ quantile思想，Z的$\tau \text{-th}$ expectile可被定义为

$e_{\tau }:=\underset{u\in ℝ}{\arg \min }E\left\{{\eta }_{\tau }\left(Z-u\right)-{\eta }_{\tau }\left(Z\right)\right\},\tau \in \left(0,1\right)$ (1)

其中，${\eta }_{\tau }\left(u\right)=\left|\tau -1\left(u<0\right)\right|\cdot \frac{u^2}{2}$ 是非对称的平方损失。$\tau =1/2$ 时，有$e_{1/2}\left(Z\right)=E\left(Z\right)$ ，可见expectile regression是均值回归的不对称类型，expectile的命名也来自“expectation”和“quantile”的组合。expectile regression比quantile regression计算效率快，却不如quantile regression稳健，对重尾分布敏感，尤其是在高维情形。Man等人[15]于2022年提出稳健expectile回归，受文献[16]的启发用Lipschitz连续和局部强凸的函数替换expectile regression损失函数中的平方部分，让$\gamma >0$ 作为平衡偏差与稳健的调节参数。对于$u\in ℝ$ ，令${\ell }_{\gamma }\left(u\right)={\gamma }^2\ell \left(u/\gamma \right)$ ，其中$\ell :ℝ\mapsto \left[0,\infty \right)$ 要满足以下三个条件：

(i) ${\ell }^{\prime }\left(0\right)=0$ 且对于所有的$u\in ℝ$ 有$\left|{\ell }^{\prime }\left(u\right)\right|\le \min \left(a_1,\left|u\right|\right)$ ；(ii) ${\ell }^{″}\left(0\right)=1$ 且对任意的$\left|u\right|\le a_3$ 都有${\ell }^{″}\left(u\right)\ge a_2$ ；(iii) 对于所有的$u\in ℝ$ 有$\left|{\ell }^{\prime }\left(u\right)-u\right|\le u^2$ 成立，其中$a_1,a_2$ 和$a_3$ 都是正常数。

满足上述条件的稳健损失函数有Huber损失：$\ell \left(u\right)=\min \left\{u^2/2,\left|u\right|-1/2\right\}$ ，以及光滑近似Huber损失的函数，例如pseudo-Huber损失：$\ell \left(u\right)=\sqrt{1+u^2}-1$ 和$\ell \left(u\right)=\log \left({\text{e}}^u/2+{\text{e}}^{-u}/2\right)$ ，smoothed Huber损失：$\ell \left(u\right)=\min \left\{u^2/2-{\left|u\right|}^3/6,\left|u\right|/2-1/6\right\}$ 。本文选择以下函数作为稳健expectile回归的损失函数：

$L_{\tau ,\gamma }\left(u\right)=\left|\tau -I\left(u<0\right)\right|\cdot \{\begin{array}{l}{u}^2/2\text{, }\left|u\right|<\gamma ,\\ \gamma u-{\gamma }^2/2,\left|u\right|>\gamma .\end{array}$ (2)

其中，$\tau \in \left(0,1\right)$ ，$\gamma >0$ 是能被数据自适应校准的稳健参数。一阶导和二阶导为：

${L^{\prime }}_{\tau ,\gamma }\left(u\right)=\{\begin{array}{l}-\left(1-\tau \right)\gamma ,u<-\gamma ,\\ \left(\tau -1\right)u,-\gamma \le u<0,\\ \tau u,0\le u\le \gamma ,\\ \tau \gamma ,u>\gamma .\end{array}$ ${L^{″}}_{\tau ,\gamma }\left(u\right)=\{\begin{array}{l}\tau -1,-\gamma \le u<0,\\ \tau ,0\le u<\gamma ,\\ 0,其他.\end{array}$

2.2. 高维部分线性可加稳健Expectile回归

若第i个样本观察值为$\left\{y_i,x_i,z_i\right\}$ ，其中$y_i\in ℝ$ 为响应变量，$x_i$ 和$z_i$ 为解释变量，$x_i={\left(x_{i1},\cdots ,x_{i{p}_n}\right)}^{\prime }$ 是$p_n$ 维的向量，$p_n$ 可随n的变化发散，$z_i={\left(z_{i1},\cdots ,z_{id}\right)}^{\prime }$ 是d维向量，d是固定常数，$i=1,\cdots ,n$ 。对于一个给定的分位水平$\tau \in \left(0,1\right)$ ，本文考虑如下部分线性可加expectile回归模型：

$y_i=x_i^{\text{T}}{\beta }_0+g_0\left(z_i\right)+{\epsilon }_i$ (3)

其中，${\beta }_0$ 是维数为$p_n$ 的线性回归系数向量；$g_0\left(z_i\right)=g_{00}+{\displaystyle {\sum }_{j=1}^d{g}_{0j}\left(z_{ij}\right)}$ 是可加的非参数部分，$g_{00}\in R$ ，并假设$\mathbb{E}\left[g_{0j}\left(z_{ij}\right)\right]=0$ ；${\epsilon }_i$ 为独立的随机噪声变量，满足$e_{\tau }\left({\epsilon }_i\right)=0$ 。模型(1)允许${\beta }_0$ 随不同的$\tau$ 而变化，因此能得出给定$x_i$ 和$z_i$ 后$y_i$ 较全的条件分布情况，进而探究解释变量和反映变量之间更完整的关系[17]。

高维统计推断一般假设系数向量${\beta }_0={\left({\beta }_{01},\cdots ,{\beta }_{0p_n}\right)}^{\prime }$ 是稀疏的，记$A=\left\{j:{\beta }_{0j}\ne 0,1\le j\le p_n\right\}$ 是系数不为零的变量下标集合，个数$q_n=\left|A\right|$ ，A在被估计前不是明确的，不失一般性地记${\beta }_0={\left({{\beta }^{\prime }}_A,{0^{\prime }}_{p_n-q_n}\right)}^{\prime }$ ，其中前$q_n$ 个系数不为零而剩下的$p_n-q_n$ 个系数为零，记$X=\left(X_1,\cdots ,X_{p_n}\right)$ 为$n\times p_n$ 矩阵，与重要线性变量相关的矩阵$X_A$ 是由X前$q_n$ 列组成的子矩阵，设$x_i$ 来自均值为0的分布，$z_{ij}\in \left[0,1\right]$ ，$\forall i,j$ 。值得指出的是，允许$q_n$ 随着n增加，意味有更多的数据被收集时模型可随着信息量的增多而变化以更好地反映事实规律。针对非参数部分$g_0(\cdot )$ ，首先引入如下定义：

定义1 设m为正整数，$v\in \left(0,1\right)$ ，$r\equiv m+v$ 。${\mathscr{H}}_r$ 是$h(\cdot )$ 在$\left[0,1\right]$ 上m次可导且$h^{\left(m\right)}(\cdot )$ 满足阶数为v的Hölder条件的函数集合，即对于任意的$h(\cdot )\in {\mathscr{H}}_r$ ，都存在正整数C使得下式成立：

$\left|h^{\left(m\right)}\left(z^{\prime }\right)-h^{\left(m\right)}\left(z\right)\right|\le C\left|z^{\prime }-z\right|,\forall \text{0}\le z^{\prime },z\le 1$ (4)

假设(1)中非参数部分的可加项$g_{0j}\in {\mathscr{H}}_r$ ，$r\ge 1.5$ ，另外$\pi \left(t\right)={\left(b_1\left(t\right),\cdots ,b_{k_n+l+1}\left(t\right)\right)}^{\prime }$ 是阶数为$l+1$ 且在$\left[0,1\right]$ 上有$k_n$ 个均匀结点的标准化B样条基函数向量。Stone [4]和Schumarker [18]指出若满足上述假设条件，则$g_{0j}(\cdot )$ 可以用B样条基函数的线性组合$\Pi \left(z_i\right)={\left(1,\pi {\left(z_{i1}\right)}^{\prime },\cdots ,\pi {\left(z_{id}\right)}^{\prime }\right)}^{\prime }$ 拟合，且存在一个向量${\xi }_0=\left({\xi }_{00},{\xi }_{01},\cdots ,{\xi }_{0d}\right)\in {ℝ}^{D_n}$ ，其中$D_n=d\left(k_n+l+1\right)+1$ ，使得${\sup }_{z_i}\left|\Pi {\left(z_i\right)}^{\prime }{\xi }_0-g_0\left(z_i\right)\right|=O\left(k_n{}^{-r}\right)$ 成立。

当已知${\beta }_0$ 中重要指标集合A时，给定$\tau \in \left(0,1\right)$ 是，模型(1)中参数的稳健expectile回归估计为：

$\left({\widehat{\beta }}_A,\widehat{\xi }\right)=\underset{\left({\beta }_A,\xi \right)}{\arg \min }\frac{1}{n}{\displaystyle \sum _{i=1}^n{L}_{\tau ,\gamma }\left(y_i-{x^{\prime }}_{A_i}{\beta }_A-\Pi {\left(z_i\right)}^{\prime }\xi \right)}$ (5)

其中${x^{\prime }}_{A_1},\cdots ,{x^{\prime }}_{A_n}$ 代表$X_A$ 的行向量。${\beta }_0$ 的oracle估计值为${\left({\widehat{{\beta }^{\prime }}}_A,{0^{\prime }}_{p_n-q{}_n}\right)}^{\prime }$ ，非参数部分$g_0\left(z_i\right)$ 的估计值为$\widehat{g}\left(z_i\right)={\widehat{g}}_0+{\displaystyle \sum _{j=1}^d{\widehat{g}}_{0j}\left(z_{ij}\right)}$ ，记$\widehat{\xi }={\left({\widehat{\xi }}_0,{{\widehat{\xi }}^{\prime }}_1,\cdots ,{{\widehat{\xi }}^{\prime }}_d\right)}^{\prime }$ ，其中${\widehat{\xi }}_0\in ℝ$ ，${{\widehat{\xi }}^{\prime }}_j\in {ℝ}^{k_n+l+1}$ ，${\widehat{g}}_{0j}\left(z_{ij}\right)=\pi {\left(z_{ij}\right)}^{\prime }{\widehat{\xi }}_j-n^{-1}{\displaystyle \sum _{i=1}^n\pi {\left(z_{ij}\right)}^{\prime }}{\widehat{\xi }}_j$ ，${\widehat{g}}_0={\widehat{\xi }}_0+n^{-1}{\displaystyle \sum _{i=1}^n\pi {\left(z_{ij}\right)}^{\prime }{\widehat{\xi }}_j}$ ，$j=1,\cdots ,d$ 。注意到${\widehat{g}}_{0j}$ 有中心化的步骤，是为了让其满足可识别条件$E\left[g_{0j}\left(z_i\right)\right]=0$ 。

3. 带非凸惩罚的高维部分线性可加Expectile回归及算法

对实际数据进行分析时，通常不清楚维数为$p_n$ 的解释变量里哪几个$q_n$ 变量是重要的，在损失函数后加上惩罚项是筛选重要解释变量的常用手段。惩罚函数Lasso [19]由于其是凸函数便于计算和进行理论分析而倍受欢迎，但不可忽略其会过惩罚较大的系数而导致估计偏差较大的缺点以及设计矩阵需要满足一些较强的条件。相比之下，一些满足条件1的非凸惩罚可以解决Lasso的问题，例如SCAD惩罚[20]。

条件1 惩罚函数$P_{\lambda }\left(\lambda >0\right)$ 对任意$t\ge 0$ 具有表达式$P_{\lambda }\left(t\right)={\lambda }^2{p}_0\left(t/\lambda \right)$ ，其中函数$p_0:{ℝ}_{+}\to {ℝ}_{+}$ 满足：(i) $p_0(\cdot )$ 是在$\left[0,\infty \right)$ 上非降的函数且$p_0\left(0\right)=0$ ；(ii) $p_0(\cdot )$ 在$\left(0,\infty \right)$ 上几乎处处可导且${\lim }_{t\downarrow 0}{p}_0\left(t\right)=1$ ；(iii) 若$t_1\ge t_2>0$ 有${p^{\prime }}_0\left(t_1\right)\le {p^{\prime }}_0\left(t_2\right)$ 。

SCAD惩罚和其一阶导数的表达式为：

$P_{\lambda }\left(\theta \right)=\lambda \theta \cdot {\rm I}\left(\theta \le \lambda \right)+\frac{a\lambda \theta -\left({\theta }^2+{\lambda }^2\right)}{a-1}{\rm I}\left(\lambda \le \theta \le a\lambda \right)+\frac{\left(a+1\right){\lambda }^2}{2}{\rm I}\left(\theta >a\lambda \right)$ (6)

${P^{\prime }}_{\lambda }\left(\theta \right)=\lambda \left\{{\rm I}\left(\theta \le \lambda \right)+\frac{{\left(a\lambda -\theta \right)}_{+}}{\left(a-1\right)\lambda }{\rm I}\left(\theta >\lambda \right)\right\}$ (7)

其中$a>2$ 是固定参数，$\theta >0$ 。于是带非凸惩罚的高维稳健expectile回归的估计值可通过以下最优化问题求

$\left(\widehat{\beta },\widehat{\xi }\right)=\underset{\left(\beta ,\xi \right)}{\arg \min }\frac{1}{n}{\displaystyle \sum _{i=1}^n{L}_{\tau ,\gamma }\left(y_i-{x^{\prime }}_i\beta -\Pi {\left(z_i\right)}^{\prime }\xi \right)}+{\displaystyle \sum _{j=2}^p{P}_{\lambda }\left(\left|{\beta }_j\right|\right)}$ (8)

但直接求解(8)会由非凸性而存在计算不稳定的问题，于是本文采用Zou和Li [21]提出的局部线性近似法，每次迭代的惩罚权重由上一次的迭代结果更新。记${P^{\prime }}_{\lambda }(\cdot )$ 为惩罚函数的一阶导，${\beta }^{\left(0\right)}$ 为迭代的初始值设为0，$S\left(a,b\right)={\left(a-b\right)}_{+}-{\left(-a-b\right)}_{+}$ 是软阀值函数，$\sigma$ 是大于1的正数，${\xi }^{\left(t-1\right)}$ 是第$t-1$ 次迭代的结果，第t次迭代的结果${\beta }^{\left(t\right)}$ 由局部线性近似法转为下面的凸优化问题来求

$\left({\beta }^{\left(t\right)}|{\xi }^{\left(t-1\right)}\right)=\underset{\beta }{\arg \min }\frac{1}{n}{\displaystyle \sum _{i=1}^n{L}_{\tau ,\gamma }\left(y_i-x_i^{\text{T}}\beta -\Pi {\left(z_i\right)}^{\text{T}}{\xi }^{\left(t-1\right)}\right)}\text{+}{\displaystyle \sum _{j=2}^p{P^{\prime }}_{\lambda }\left(\left|{\beta }_j^{\left(t-1\right)}\right|\right)}\left|{\beta }_j\right|$ (9)

其中${P^{\prime }}_{\lambda }\left(\left|{\beta }_j^{\left(t-1\right)}\right|\right)\left|{\beta }_j\right|$ 是非凸惩罚函数$P_{\lambda }\left(\left|{\beta }_j\right|\right)$ 在$\left|{\beta }_j^{\left(t-1\right)}\right|$ 处的局部线性近似。

本文运用两步算法，整体迭代的顺序是先得到不带惩罚的非参数部分的系数$\xi$ ，再求带惩罚的线性系数$\beta$ ，然后求调节参数$\gamma$ 。其中非参数部分的迭代算法为经典的坐标下降法(GDBB) [22]，带惩罚项的线性部分为Fan等人提出的Local Adaptive Majorize-Minimization算法(LAMM) [23]。具体算法描述如表1~3所示。

Table 1. Two-step algorithm

表1. 两步算法

算法1 两步算法

Step 1. Input ${\gamma }^{\left(0\right)}=\sqrt{n/\log \left(np\right)}$ , ${\epsilon }_1$ , ${\epsilon }_2$ , ${\beta }^{\left(0\right)}=0$ , $t=0$ Step 2. 先求非参数部分的系数$\xi$ ，再求带惩罚的线性系数$\beta$ (a) 非参数部分：基于$t-1$ 次的迭代结果${\beta }^{\left(t-1\right)}$ ，通过以下求解最优化问题得到${\xi }^{\left(t\right)}$ ： ${\xi }^{\left(t\right)}=\underset{\xi \in {ℝ}^{D_n}}{\arg \min }\frac{1}{n}{\displaystyle \sum _{i=1}^n{L}_{\tau ,\gamma }\left(y_i-x_i^{\text{T}}{\beta }^{\left(t-1\right)}-\Pi {\left(z_i\right)}^{\text{T}}\xi \right)}$ (b) 线性部分： (b.1) 用第$t-1$ 次的迭代结果${\beta }^{\left(t-1\right)}={\left({\beta }_1^{\left(t-1\right)},\cdots ,{\beta }_p^{\left(t-1\right)}\right)}^{\prime }$ 计算第t次迭代权重${P^{\prime }}_{\lambda }\left(\left|{\beta }_j{}^{\left(t-1\right)}\right|\right)$ (b.2) 把(a)的计算结果${\xi }^{\left(t\right)}$ 带入下面的最优化问题求解${\beta }^{\left(t\right)}$ ： $\left({\beta }^{\left(t\right)}|{\xi }^{\left(t\right)}\right)=\underset{\beta }{\arg \min }\frac{1}{n}{\displaystyle \sum _{i=1}^n{L}_{\tau ,\gamma }\left(y_i-x_i^{\text{T}}\beta -\Pi {\left(z_i\right)}^{\text{T}}{\xi }^{\left(t\right)}\right)}+{\displaystyle \sum _{j=2}^p{P^{\prime }}_{\lambda }\left(\left|{\beta }_j{}^{\left(t-1\right)}\right|\right)}\left|{\beta }_j\right|$ Step 3. 把${\xi }^{\left(t\right)}$ 和${\beta }^{\left(t\right)}$ 带入下面几式中计算$r_i^{\left(t\right)}$ ，${\tilde{r}}_i{}^{\left(t\right)}$ ，$mad\left({\tilde{r}}^{\left(t\right)}\right)$ 和${\gamma }^{\left(t\right)}$ $r_i^{\left(t\right)}=y_i-x_i^{\text{T}}{\beta }^{\left(t-1\right)}-\Pi {\left(z_i\right)}^{\text{T}}{\xi }^{\left(t-1\right)},{\tilde{r}}^{\left(t\right)}=\left(1-\tau \right)r_i^{\left(t\right)}{{\rm I}}_{r_i^{\left(t\right)}\le 0}+\tau r_i^{\left(t\right)}{{\rm I}}_{r_i^{\left(t\right)}>0}$ $mad\left({\tilde{r}}^{\left(t\right)}\right)={\left\{{\Phi }^{-1}\left(0.75\right)\right\}}^{-1}median\left(\left|{\tilde{r}}^{\left(t\right)}-median\left({\tilde{r}}^{\left(t\right)}\right)\right|\right)$ ${\gamma }^{\left(t\right)}=mad\left({\tilde{r}}^{\left(t\right)}\right)\cdot \sqrt{\frac{n}{\log \left(np\right)}}$ Step 4. 如果$\Vert {\beta }^{\left(t\right)}-{\beta }^{\left(t-1\right)}\Vert \le {\epsilon }_1$ , $\left|\Pi {\left(z_i\right)}^{\text{T}}{\xi }^{\left(t\right)}-\Pi {\left(z_i\right)}^{\text{T}}{\xi }^{\left(t-1\right)}\right|\le {\epsilon }_2$ ，那么令$\widehat{\beta }={\beta }^{\left(t\right)}$ ，$\widehat{\xi }={\xi }^{\left(t\right)}$ ，否则令$t=t+1$ ，返回step 2继续进行迭代

Table 2. Step (a) of Algorithm 1

表2. 算法1里的a步骤

算法1.1 Gradient Descent with Barzilai-Borwein Step Size (GDBB)

Step 1. Input ${\xi }_0$ , $0<{\epsilon }_a\le 1$ , $k=0$ Step 2. $d_k=\frac{\partial \left(\frac{1}{n}{\displaystyle \sum _{i=1}^n{L}_{\tau ,\gamma }\left(y_i-x_i^{\text{T}}{\beta }^{\left(t-1\right)}-\Pi {\left(z_i\right)}^{\text{T}}{\xi }_k\right)}\right)}{\partial {\xi }_k}$ . if $\Vert d_k\Vert <{\epsilon }_a$ , then stop, return ${\xi }_k={\xi }^{\left(t\right)}$ ; else Step3 Step 3. 计算${\alpha }_k=\min \left\{\frac{s_{k-1}^{\text{T}}{s}_{k-1}}{s_{k-1}^{\text{T}}{h}_{k-1}},\frac{s_{k-1}^{\text{T}}{h}_{k-1}}{h_{k-1}^{\text{T}}{h}_{k-1}}\right\}$ , where $s_{k-1}={\xi }_k-{\xi }_{k-1}$ , $h_{k-1}=d_k-d_{k-1}$ , Step 4. ${\xi }_{k+1}={\xi }_k-{\alpha }_k{d}_k$ Step 5. $k=k+1$ ,，回到step 2

Table 3. Step (b.2) of Algorithm 1

表3. 算法1里的b.2步骤

算法1.2 Local Adaptive Majorize-Minimization (LAMM)

Step 1. Input ${\beta }^{\left(0\right)}={\beta }^{\left(t-1\right)},0<{\epsilon }_b\le 1,k=1,{\varphi }_0,{\varphi }^{\left(k-1\right)},\sigma$ Step 2. 基于第$k-1$ 次的迭代结果${\beta }^{\left(k-1\right)}={\left({\beta }_1^{\left(k-1\right)},\cdots ,{\beta }_p^{\left(k-1\right)}\right)}^{\prime }$ ，计算第k次迭代相关的权重： $w^{\left(k\right)}={\left({P^{\prime }}_{\lambda }\left(\left|{\beta }_1{}^{\left(k-1\right)}\right|\right),{P^{\prime }}_{\lambda }\left(\left|{\beta }_2{}^{\left(k-1\right)}\right|\right),\cdots ,{P^{\prime }}_{\lambda }\left(\left|{\beta }_p{}^{\left(k-1\right)}\right|\right)\right)}^{\text{T}}={\left(w_1^{\left(k\right)},\cdots ,w_j^{\left(k\right)}\right)}^{\text{T}}$ Step 3. Initialize: ${\varphi }^{\left(k\right)}\leftarrow \max \left\{{\varphi }_0,{\sigma }^{-1}{\varphi }^{\left(k-1\right)}\right\}$ repeat: If $F\left({\beta }^{\left(k\right)},w^{\left(k\right)}\right)\ge {\Psi }_{w^{\left(k\right)},{\varphi }^k}\left({\beta }^{\left(k\right)};{\beta }^{\left(k-1\right)}\right)$ , then ${\varphi }^{\left(k\right)}\leftarrow \sigma {\varphi }^{\left(k\right)}$ until $F\left({\beta }^{\left(k\right)},w^{\left(k\right)}\right)\le {\Psi }_{w^{\left(k\right)},{\varphi }^k}\left({\beta }^{\left(k\right)};{\beta }^{\left(k-1\right)}\right)$ , return ${\beta }_j^{\left(k\right)}\leftarrow T_{w^{\left(k\right)},{\varphi }^{（k）}}\left({\beta }_j^{\left(k-1\right)}\right)$ and ${\varphi }^{(k)}$ Where $F\left({\beta }^{\left(k-1\right)},w^{\left(k\right)}\right)\equiv \Gamma \left({\beta }^{\left(k-1\right)}\right)+{\displaystyle \sum _{j=2}^p{w}_j^{\left(k\right)}}\left|{\beta }_j^{\left(k-1\right)}\right|$ , $\Gamma \left({\beta }^{\left(k-1\right)}\right)=\frac{1}{n}{\displaystyle \sum _{i=1}^n{L}_{\tau ,\gamma }\left(y_i-x_i^{\text{T}}{\beta }^{\left(k-1\right)}-\Pi {\left(z_i\right)}^{\text{T}}{\xi }^{\left(t\right)}\right)}$ , ${\beta }_j{}^{\left(k\right)}\equiv T_{w^{\left(k\right)},{\varphi }^{（k）}}\left({\beta }_j^{\left(k-1\right)}\right)=S\left({\beta }_j{}^{\left(k-1\right)}-\frac{\nabla \Gamma \left({\beta }_j{}^{\left(k-1\right)}\right)}{{\varphi }^{(k)}},\frac{w_j^{(k)}}{{\varphi }^{(k)}}\right)$ ${\Psi }_{w^{\left(k\right)},{\varphi }^k}\left(\beta ;{\beta }^{\left(k-1\right)}\right)\equiv \Gamma \left({\beta }^{\left(k-1\right)}\right)+\langle \nabla \Gamma \left({\beta }^{\left(k-1\right)}\right),\beta -{\beta }^{\left(k-1\right)}\rangle +\frac{{\varphi }^{(k)}}{2}{\Vert \beta -{\beta }^{\left(k-1\right)}\Vert }_2^2+{\displaystyle \sum _{j=2}^p{w}_j^{(k)}}\left|{\beta }_j^{\left(k-1\right)}\right|$ Step 4. If $\Vert {\beta }^{\left(k\right)}-{\beta }^{\left(k-1\right)}\Vert <{\epsilon }_b$ then ${\beta }^{\left(t\right)}={\beta }^{\left(k\right)}$ , otherwise $k=k+1$ ，返回step 2继续进行迭代

4. 数值模拟

本文分别研究了部分线性可加稳健expectile回归(retire-PLAMs)在不带惩罚的低维模型和带惩罚的稀疏高维模型两种情况下数值模拟的表现，在部分线性可加框架下让retire与不同的方法比较：(i) huber；(ii) 非对称最小平方回归(sales)；(iii) quantile回归(qr)。各进行了100次随机模拟，采用以下准则来比较各个方法：

1) MSE：估计${\beta }_0$ 的均方误差，即一百次模拟结果中${\Vert \widehat{\beta }-{\beta }_0\Vert }^2$ 的均值。

2) AE：${\beta }_0$ 的绝对估计误差的均值，即一百次模拟结果中${\displaystyle \sum _{j=1}^p\left|{\widehat{\beta }}_j-{\beta }_j\right|}$ 的均值。

3) AADE：拟合非参数的平均绝对离差的均值，即一百次模拟结果$n^{-1}{\displaystyle \sum _{i=1}^n\left|\widehat{g}\left(z_i\right)-g_0\left(z_i\right)\right|}$ 中的均值。

4) TV：一百次模拟中正确地识别到系数为非零的线性解释变量的次数占比的均值。

5) FV：一百次模拟中错误地识别出系数为非零的线性解释变量的次数占比的均值。

4.1. 不带惩罚的低维模型

首先生成来自多元正态分布$N_{11}\left(0_{11},\Sigma \right)$ 的$\tilde{X}={\left({\tilde{X}}_1,\cdots ,{\tilde{X}}_{11}\right)}^{\prime }$ ，其中$\Sigma ={\left({\sigma }_{jk}\right)}_{11\times 11}$ ，${\sigma }_{jk}={0.5}^{\left|j-k\right|}$ 。让$Z_1=\Phi \left({\tilde{X}}_{10}\right)$ ，$Z_2=\Phi \left({\tilde{X}}_{11}\right)$ ，$\Phi (\cdot )$ 是标准正态分布的分布函数，$X_l={\tilde{X}}_l$ ，$l=1,\cdots ,9$ ，响应变量$Y_i$ 有以下回归模型生成，

$Y_i={X^{\prime }}_i{\beta }_0+\sin \left(2\pi Z_{i1}\right)+Z_{i2}^3+\left(0.5\left|X_{i9}\right|+0.5\right)\left\{{\epsilon }_i-e_{\tau }\left({\epsilon }_i\right)\right\}$ (10)

其中${\beta }_0={\left({\beta }_1,\cdots ,{\beta }_{10}\right)}^{\prime }={\left(2,1.6,1.3,1,0.7,0.4,0.1,-0,1,-0.4,-0.7\right)}^{\prime }$ ，${\beta }_1$ 为截距，${\epsilon }_i$ 可分两种情形，一是来自正态分布$N\left(0,2\right)$ ，二是来自自由度为2.1的t分布，样本量n分为200和400两种情形，$\tau$ 水平分为0.5和0.8两种情况。

表4和表5分别展示了随机噪声来$N\left(0,2\right)$ 和$t\left(2.1\right)$ 的模型(13)在$\tau =\left\{0.5,0.8\right\}$ 时四种不同方法的拟合表现，可以看出随着样本量增加，所有模型的估计效果都有所改进，另外$\tau =0.5$ 时retire与huber估计结果相同，因为$\tau =0.5$ 时二者为同一个模型。表5可见随机噪声来自正态分布时，retire的估计误差是所有模型中表现最好的；表4显示随机误差来自t分布时，retire的估计误差比huber和sales小，比quantile的估计误稍大，原因在于quantile遇到异常值时是最稳健的，但是quantile计算速度比retire慢。总而言之，数值模拟验证retire-PLAMs在拥有稳健性的同时还具有计算效率高的优势。

Table 4. Comparison of the fitting performance of the four methods for model (13) with random noise from t(2.1) at $\tau$ levels of 0.5 and 0.8

表4. 四种方法在随机噪声来自t(2.1)的模型(13)在$\tau$ 水平为0.5和0.8的拟合表现比较

		n = 200	p = 10		n = 400	p = 10
$\tau$	method	AE	MSE	AADE	AE	MSE	AADE
0.5	retire	1.413	0.310	0.303	0.974	0.152	0.205
	huber	1.413	0.310	0.303	0.974	0.152	0.205
	sales	1.784	0.545	0.469	1.437	0.488	0.337
	qr	0.974	0.150	0.270	0.618	0.061	0.192
0.8	retire	1.584	0.402	0.390	1.080	0.192	0.260
	huber	1.431	0.320	0.841	1.002	0.160	0.820
	sales	2.309	1.072	0.623	1.957	1.384	0.485
	qr	1.343	0.290	0.403	0.889	0.130	0.275

Table 5. Comparison of the fitting performance of the four methods for model (13) with random noise from N(0,2) at $\tau$ levels of 0.5 and 0.8

表5. 四种方法在随机噪声来自N(0,2)的模型(13)在$\tau$ 水平为0.5和0.8的拟合表现比较

		n = 200	p = 10		n = 400	p = 10
$\tau$	method	AE	MSE	AADE	AE	MSE	AADE
0.5	retire	1.007	0.163	0.265	0.725	0.082	0.184
	huber	1.007	0.163	0.265	0.725	0.082	0.184
	sales	1.010	0.163	0.272	0.726	0.082	0.189
	qr	1.066	0.181	0.283	0.737	0.087	0.197
0.8	retire	1.088	0.187	0.327	0.822	0.103	0.260
	huber	1.035	0.172	0.612	0.728	0.083	0.608
	sales	1.073	0.182	0.280	0.805	0.099	0.193
	qr	1.348	0.287	0.359	0.870	0.117	0.283

4.2. 带惩罚的稀疏模型

生成来自$N_{p+2}\left(0_{p+2},\Sigma \right)$ 的$\tilde{X}={\left({\tilde{X}}_1,\cdots ,{\tilde{X}}_{p+2}\right)}^{\prime }$ ，其中$\Sigma ={\left({\sigma }_{jk}\right)}_{\left(p+2\right)\times \left(p+2\right)}$ ，${\sigma }_{jk}={0.5}^{\left|j-k\right|}$ 。让$Z_1=\Phi \left({\tilde{X}}_{10}\right)$ ，$Z_2=\Phi \left({\tilde{X}}_{11}\right)$ ，$\Phi (\cdot )$ 是标准正态分布的分布函数；对于$l=1,\cdots ,9$ ，$X_l={\tilde{X}}_l$ ；对于$l=13,\cdots ,p+2$ ，$X_l={\tilde{X}}_{l-2}$ 。响应变量$Y_i$ 有以下回归模型生成：

$Y_i={X^{\prime }}_i{\beta }_0+\sin \left(2\pi Z_{i1}\right)+Z_{i2}^3+\left(0.5\left|X_{i9}\right|+0.5\right)\left\{{\epsilon }_i-e_{\tau }\left({\epsilon }_i\right)\right\}$ (11)

其中${\beta }_0={\left({\beta }_1,\cdots ,{\beta }_{10},0_{p-8}\right)}^{\prime }={\left(2,1.6,1.3,1,0.7,0.4,0.1,-0,1,-0.4,-0.7,{0^{\prime }}_{p-8}\right)}^{\prime }$ ，${\epsilon }_i$ 可分两种情形，一是来自正态分布$N\left(0,2\right)$ ，二是来自自由度为2.1的t分布，样本量$n=400$ ，$p=300,500$ 两种情形，$\tau$ 水平为0.5和0.8两种情况。表6展示了带SCAD惩罚的retire与带Lasso惩罚的另外三种方法的拟合结果表现比较，可见retire-SCAD的线性参数估计误差最小，TP最大，FP最小，是变量筛选最一致的。

Table 6. Comparison of the fitting performance of the four methods for model (14) with random noise from t(2.1) at $\tau$ levels of 0.5 and 0.8

表6. 四种方法在随机噪声来自t(2.1)的模型(14)在$\tau$ 水平为0.5和0.8的拟合表现比较

		n = 400	p = 300			n = 400	p = 500
$\tau$	method	MSE	AADE	TP	FP	MSE	AADE	TP	FP
0.5	retire-SCAD	0.107	0.693	1	0.015	0.132	0.7123	1	0.010
	roer-L1	0.443	0.400	1	0.095	0.510	0.439	1	0.056
	huber-L1	0.443	0.400	1	0.095	0.510	0.439	1	0.056
	sales-L1	0.991	0.519	0.995	0.093	1.092	0.546	0.993	0.577
	qr-L1	0.006	0.018	1	0.017	0.072	0.019	1	0.129
0.8	retire-SCAD	0.157	0.656	1	0.010	0.161	0.607	1	0.007
	roer-L1	0.609	0.463	1	0.099	0.621	0.461	1	0.062
	huber-L1	0.471	0.560	1	0.095	0.497	0.497	1	0.057
	salse-L1	1.002	0.615	0.994	0.091	1.207	0.780	0.996	0.072
	qr-L1	0.263	0.021	1	0.169	0.257	0.019	1	0.135

4.3. 时间对比

笔者增加此部分模拟以展示所提方法retire-L1与qr-L1相比在计算效率上的优势。分别借助R语言的adaHuber和rqPen两个包，稀疏正则参数$\lambda$ 由十则交叉验证选择，retire的稳健参数$\gamma$ 如表1所示由数据自适应调节。100次模拟的数据由模型(15)生成：

$Y_i={X^{\prime }}_i{\beta }_0+\sin \left(2\pi Z_{i1}\right)+Z_{i2}^3+{\epsilon }_i$ (12)

其中$X_i$ ，$Z_{i1}$ 和$Z_{i2}$ 同4.2节的设定，${\epsilon }_i$ 来自自由度为2.1的t分布，$\tau =0.5$ ，样本数量$n=p/2$ ，解释变量维数$p=\left\{100,200,300,400,500\right\}$ 。

模拟结果如表7所示，运行时间的单位为秒，随着n和p的增加估计误差都在减少。虽然在$t_{2.1}$ 的随机噪声下qr-L1比retire-L1较稳健，但计算效率不如retire-L1快，特别是p很大的时候，例$p=500$ 时，qr-L1的计算时间是retire-L1的16倍。

Table 7. Comparison of time and estimation results between retire-L1 and qr-L1 at $\tau =0.5$

表7. Retire-L1与qr-L1在$\tau =0.5$ 的时间和估计结果对比

	p = 100, n = 50	p = 200, n = 100	p = 300, n = 150	p = 400, n = 200	p = 500, n = 250
qr-L1-time/qr-L1-time	1	1	1	1	1
qr-L1-AADE/qr-L1-AADE	1	1	1	1	1
qr-L1-MSE/qr-L1-MSE	1	1	1	1	1
retire-L1-time/qr-L1-time	0.176	0.121	0.092	0.065	0.060
retire-L1-AADE/qr-L1-AADE	1.496	1.351	1.240545	1.261	1.081
retire-L1-MSE/qr-L1-MSE	1.139	1.207	1.340393	1.406	1.432

5. 实证研究

文献[24]指出$\beta$ 胡萝卜素与肺癌、乳腺癌等癌症之间存在直接关系，一些流行病学研究表明其抗氧化特性能清除可能导致癌症的自由基，充足的$\beta$ 胡萝卜素能提高人体免疫力有效地对抗癌症等疾病。我们分析Nierenberg等人[25]收集的血浆$\beta$ 胡萝卜素水平数据集并建立部分线性可加模型，寻找血浆$\beta$ 胡萝卜素水平与个人特征之间的关系，包括年龄、性别、体重指数和其他因素，响应变量Y是以ng/ml为单位的血浆$\beta$ 胡萝卜素，命名为BETAPLASMA，各变量的说明见表8，使用以下部分线性加性模型进行建模：

$\begin{array}{c}Y=u+{\beta }_1\text{SEX}+{\beta }_2\text{SMOK}1+{\beta }_3\text{SMOK}2+{\beta }_4\text{BMI}+{\beta }_5\text{VIT1}+{\beta }_6\text{VIT}2+{\beta }_7\text{CAL}+{\beta }_8\text{FAT}\\ +{\beta }_9\text{ALCOHOL}+{\beta }_{10}\text{BETADIET}+{\eta }_1\left(\text{AGE}\right)+{\eta }_2\left(\text{CHOL}\right)+{\eta }_2\left(\text{FIBER}\right)+\sigma \epsilon \\ =u+{\beta }^{\text{T}}X+{\displaystyle \sum _{j=1}^3{\eta }_1\left(z_j\right)}+\sigma \epsilon \end{array}$ (16)

Table 8. Explanation of the variables for the Plasma $\beta$ Carotene Levels dataset

表8. 血浆$\beta$ 胡萝卜素水平数据集变量的解释说明

变量名	含义
SEX	1 = 男，0 = 女
SMOK1	1 = 以前吸烟，0 = 其他
SMOK2	1 = 现在吸烟，0 = 其他
BMI	身体质量指数
VIT1	1 = 经常补充维生素，0 = 其他
VIT2	1 = 偶尔补充维生素，0 = 其他
CAL	每天消耗的卡路里数
AGE	年龄
FAT	每天消耗的脂肪克数
ALCOHOL	每周饮用的酒精饮料数量
BETADIET	膳食β-胡萝卜素消耗量(微克/天)
CHOL	每日摄入的胆固醇(毫克/天)
FIBER	每日摄入的纤维克数(克/天)
BETAPLASMA	血浆 β-胡萝卜素(纳克/毫升)

在建模前用Min-Max标准化法消除变量量纲不同可能造成的影响，AGE、CHOL和FIBER放入非线性部分，其余解释变量放入线性部分。用带SCAD惩罚的稳健expectile回归筛选影响人体$\beta$ 胡萝卜素含量的重要变量，研究血浆$\beta$ 胡萝卜素在不同expectile水平$\tau =\left\{0.2,0.5,0.9\right\}$ 的条件分布，回归结果见表9。可以看到三个不同$\tau$ 水平的回归系数有差异，说明数据存在异质性，从十个变量里筛选出BMI、胡萝卜消耗、吸烟、维他命摄入和性别对$\beta$ 胡萝卜素含量占主要影响的五个因素。进而推出合理的BMI范围、不吸烟、多摄入维生素能保证人体$\beta$ 胡萝卜素含量充足，能预防癌症等疾病的结论。

Table 9. Regression coefficients of important variables at three different $\tau$ levels

表9. 在三个不同$\tau$ 水平下的重要变量的回归系数

$\tau$	BMI	BETADIET	SMOK1	SMOK2	VIT1	VIT2	SEX
0.2	−143	84	−9	−39	0	−30	−29
0.5	−187	135	0	−42	23	−39	−23
0.8	−222	140	−11	−59	33	−38	−31

6. 总结

本文在PLAMs框架下研究高维稳健expectile回归模型，数值模拟和实证研究显示该法与quantile regression一样能通过取不同的$\tau$ 水平来分析数据异质性，但比quantile regression的计算效率高，而估计效果和quantile regression同样具有稳健性。