引言

随着科技的发展，超高维数据越来越多地出现在诸如基因表达、信号处理、金融分析等领域中。在这类数据中，协变量的维数远大于样本量，且随着样本量的增加呈指数级增长，然而只有少量的协变量对响应变量有影响，呈现稀疏性的特点。现有的基于惩罚的变量选择方法都面临着计算成本、统计精度及算法稳定性等挑战^[1-3]，不能很好地处理超高维数据的降维问题。为了解决这一问题，近年来许多学者提出了各种超高维数据的特征筛选方法。Fan等^[4]针对线性模型提出一种基于Pearson相关系数的确定独立性筛选(SIS)方法，随后他们将SIS方法进一步推广到广义线性模型和非参数可加模型^[5-6]。在无模型假定的条件下，Zhu等^[7]基于响应变量的条件分布和协变量的边际相关性提出了确定独立性排序筛选(SIRS)方法；Li等^[8]基于距离相关系数提出了距离相关性确定独立性筛选(DC-SIS)方法；Li等^[9]基于Kendall τ相关系数提出了稳健秩变量筛选(RRCS)方法。

为了解决超高维数据中异方差的问题，结合分位数回归的稳健性和全面性，He等^[10]通过样条函数逼近边际分位数回归的方式，提出了分位数自适应确定独立筛选(QaSIS)方法。Wu等^[11]提出了条件分位数特征筛选(Q-SIS)方法和条件分布函数特征筛选(DF-SIS)方法，该方法计算简单快捷，无须进行非参数估计，且对协变量没有有限矩条件的限制。Chang等^[12]创造性地将经验似然方法用于超高维数据，通过对零点处的边际似然比进行排序，提出了线性模型下基于边际经验似然的特征筛选(EL-SIS)方法，该方法只涉及单变量优化问题，便于计算，对分布的假设较宽松。随后，EL-SIS方法又被进一步推广到半参数和非参数模型中^[13-14]。

然而上述基于边际经验似然的筛选方法都是在一定的模型框架下，为了避免特定的模型假设，同时为了有效解决超高维数据中异方差的问题，本文将EL-SIS与条件分位数筛选方法相结合，提出了基于边际经验似然的分位数特征筛选(EL-QSIS) 和分布函数特征筛选(EL-DFSIS)方法。沿袭经验似然方法的特征，所提方法具有计算简单快捷、无须参数估计、对分布的假设较宽松、不依赖于模型假定等优点。通过理论证明、数值模拟和实例分析进一步验证了所提方法满足确定筛选性质且具有良好的有限样本性质。

1 变量筛选方法 1.1 基于边际经验似然的分位数特征筛选

令Y和X=(X₁, X₂, …, X_pn)^T分别表示响应变量和p_n维协变量，其中维数p_n随着n的增加呈指数级增长。不失一般性，假定E(X_k)=0，E(X_k²)=1，k=1, 2, …, p_n。假设Y与X之间满足稀疏性原则, 即只有少部分协变量对响应变量有影响。对某给定的τ∈(0, 1)，记条件τ分位数Q_τ(Y|X)=inf {y: P(Y≤y|X)≥τ}，为了筛选出其中的重要变量，定义活跃指标集为M_τ={k: Q_τ(Y|X)依赖于X_k，k=1, 2, …, p_n}。记s_n=|M_τ|，其中|M_τ|表示M_τ中元素的个数，根据稀疏性原则，s_n＜n。

由文献[11]，若Q_τ(Y|X_k)=Q_τ(Y)，k=1, 2, …, p_n，则对任意t∈R，d_k(t)=0，其中d_k(t)=E{[τ－I(Y＜Q_τ(Y))]I(X_k＜t)}。因此可以用d_k(t)衡量X_k和Y之间的边际相关性。记g_k=E{d_k²(X_k)}, 则当Q_τ(Y|X_k)=Q_τ(Y)时，g_k=0，反之g_k＞0。由此可见, g_k越大，则越说明X_k是影响Y的条件τ分位数的重要变量。因此，可以通过度量边际效用g_k是否等于0来进行特征筛选。

设有独立同分布的样本{(X_i, Y_i)}_i=1ⁿ, 且样本协变量已经标准化, 即

$ \frac{1}{n} \sum\limits_{i=1}^n X_{i k}=0, \frac{1}{n} \sum\limits_{i=1}^n X_{i k}^2=1, k=1, 2, \cdots, p_n $

令g_jk=d_k²(X_jk), j=1, 2, …, n, 定义如下边际经验似然。

$ \begin{array}{l} \;\;\;\;E L_k=\sup \left\{\prod\limits_{j=1}^n \omega_j: \omega_j \geqslant 0, \sum\limits_{j=1}^n \omega_j=1, \sum\limits_{j=1}^n \omega_j g_{j k}=\right. \\ 0\} \end{array} $

(1)

利用拉格朗日乘子法求解式(1)，得到边际经验似然比为

$ l_k=-2 \ln \left\{E L_k\right\}-2 n \ln n=2 \sum\limits_{j=1}^n \ln \left(1+\lambda g_{j k}\right) $

其中，拉格朗日乘子λ满足$ \sum\limits_{j=1}^n \frac{g_{j k}}{1+\lambda g_{j k}}=0$。根据Chang等^[12]的理论分析，当k∉M_τ时，l_k值不会太大，而当k∈M_τ时，l_k将以很大的概率发散，故可将l_k作为度量指标进行排序从而筛选出重要变量。

为了给出l_k的估计，定义g_jk(j=1, 2, …, n)的经验估计为

$ \begin{aligned} & \;\;\;\;\hat{g}_{j k}=\left\{\frac { 1 } { n } \sum _ { i = 1 } ^ { n } [ \tau - I ( Y _ { i } < \hat { Q } _ { \tau } ( Y ) ) ] I \left(X_{i k}<\right.\right. \\ \left.X_{j k}\right) & \}^2 \end{aligned} $

式中，$ \hat{Q}_\tau(Y)$是基于Y₁, Y₂, …, Y_n的样本τ分位数。相应地，可估计l_k为

$ \hat{l}_k=2 \sum\limits_{j=1}^n \ln \left(1+\hat{\lambda} \hat{g}_{j k}\right) $

式中，$ , \hat{\lambda} $满足$ \sum\limits_{j=1}^n \frac{\hat{g}_{j k}^\tau}{1+\hat{\lambda} \hat{g}_{j k}^\tau}=0$。从而M_τ可估计为

$ \hat{M}_\tau=\left\{k: \hat{l}_k \geqslant \gamma_n, k=1, 2, \cdots, p_n\right\} $

式中，γ_n是预先设定的阈值。在实际应用中，常将$ \hat{l}_k$, k=1, 2, …, p_n按降序排列，此时可得M_τ的估计为

$ \hat{M}_\tau=\left\{k: \hat{l}_k \text { 为前 } d_n \text { 个最大的, } k=1, 2, \cdots, p_n\right\} $

其中指定的模型大小d_n可仿照文献[4]取为[n/ln n]，这里[a]表示取不大于a的整数。

1.2 理论性质

为了证明EL-QSIS方法满足确定筛选性质，假设下列正则化条件成立。

1) 存在常数c＞0, 使得

$ \min \limits_{k \in M_\tau} g_k \geqslant c n^{-\kappa} \text {, 对某 } \kappa \in[0, 1 / 2) \text { 。} $

2) 在Q_τ(Y)附近, F(y)二阶可微，对于f(y), 存在正数c₁、c₂, 使得0＜c₁＜f(y)＜c₂＜∞一致成立，且f′(y)一致有界, 其中F(y)和f(y)分别为Y的分布函数和密度函数。

3) 在给定X_Mτ下，X_Mτ^c与I(Y＜Q_τ(Y))条件独立，且X_Mτ与X_Mτ^c相互独立，其中X_Mτ={X_j: j∈M_τ}, X_Mτ^c={X_j: j∉M_τ}。

文献[11]有相同的条件1)~3)。条件1)要求重要变量对应的g_k中的最小值不能太小，这也意味着重要的协变量的信号不能太弱，这个条件被广泛应用于超高维数据的特征筛选中。条件2)是分位数回归的常见条件。利用条件3)可以把重要变量和非重要变量区分开，从而保证筛选排序的一致性。注意在文献[12]的条件A.2中要求协变量与响应变量的尾部满足指数衰减速率，而本文对协变量没有任何限制条件，因此，本文所提出的筛选方法对重尾分布更加稳健。

引理1  在条件1)、2)下，有

$ \max \limits_j\left|\hat{g}_{j k}-g_{j k}\right|=O_p\left(n^{-\kappa}\right) $

由$ \left|\hat{Q}_\tau(Y)-Q_\tau(Y)\right|=O\left(n^{-1 / 2}(\ln n)^{1 / 2}\right) a . s .$和文献[11]中定理1的证明步骤，容易得到引理1的结论。

定理1  在条件1)~3)下，存在常数C₁＞0，对任意α∈(0, 1/2－κ)，有

$ \max \limits_{k \in M_\tau} P\left\{\hat{l}_k<c^2 n^{2 \alpha}\right\} \leqslant \exp \left\{-C_1 n^{1-2 \kappa}\right\} $

证明：对∀k∈M_τ，由文献[15]得

$ \hat{l}_k=2 \max \limits_{\hat{\lambda} \in {\mathit{\Lambda}}_n} \sum\limits_{j=1}^n \ln \left\{1+\hat{\lambda} \hat{g}_{j k}\right\} $

其中$ {\mathit{\Lambda}}_n=\left\{\hat{\lambda} \text { : 对任 } j=1, 2, \cdots, n, 1+\hat{\lambda} \hat{g}_{j k} \geqslant n^{-1}\right\}$。注意0≤g_jk≤1，且仿照文献[12]中定理1、命题2的证明步骤, 对某ε＞0，取$ \hat{\lambda}=\left(n^{\varepsilon} \max _j \hat{g}_{j k}\right)^{-1}$，则对某t＞0, 得

$ \begin{array}{l} \;\;\;\;P\left\{\hat{l}_k<2 t\right\} \leqslant P\left\{\sum\limits_{j=1}^n \frac{\hat{g}_{j k}}{n^{\varepsilon} \max _j \hat{g}_{j k}}<t+n^{1-2 \varepsilon}\right\}= \\ P\left\{\sum\limits_{j=1}^n \hat{g}_{j k}<\left(t n^{\varepsilon}+n^{1-\varepsilon}\right) \max _j \hat{g}_{j k}\right\} \leqslant P\left\{\frac { 1 } { \sqrt { n } \sigma _ { k } } \sum \limits_ { j = 1 } ^ { n } \left(g_{j k}-\right.\right. \\ \left.g_k\right)<\frac{1}{\sigma_k}\left[\left(t n^{\varepsilon-1 / 2}+n^{1 / 2-\varepsilon}\right) \max _j \hat{g}_{j k}-\sqrt{n} g_k+\sqrt{n} \max _j \mid \hat{g}_{j k}-\right. \\ \left.\left.g_{j k} \mid\right]\right\} \leqslant P\left\{\frac{1}{\sqrt{n} \sigma_k} \sum\limits_{j=1}^n\left(g_{j k}-g_k\right)<\frac{1}{\sigma_k}\left[\left(t n^{\varepsilon-1 / 2}+\right.\right.\right. \\ \left.n^{1 / 2-\varepsilon}\right) \max _j g_{j k}-\sqrt{n} g_k+\left(t n^{\varepsilon-1 / 2}+n^{1 / 2}+n^{1 / 2-\varepsilon}\right) \cdot \\ \left.\left.\max _j\left|\hat{g}_{j k}-g_{j k}\right|\right]\right\} \leqslant P\left\{\frac{1}{\sqrt{n} \sigma_k} \sum\limits_{j=1}^n\left(g_{j k}-g_k\right)<\frac{1}{\sigma_k} \cdot\right. \\ {\left[\left(t n^{\varepsilon-1 / 2}+n^{1 / 2-\varepsilon}\right)-\sqrt{n} g_k+\left(t n^{\varepsilon-1 / 2}+n^{1 / 2}+n^{1 / 2-\varepsilon}\right) \cdot\right.} \\ \left.\left.\max _j\left|\hat{g}_{j k}-g_{j k}\right|\right]\right\} \end{array} $

式中，g_k=E(g_jk)，σ_k²=Var(g_jk)。对于L→∞，取ε使得n^ε=L/g_k，且令2t=ng_k²/L²，则有

$ \begin{aligned} & \frac{t n^{\varepsilon}}{n g_k}=\frac{1}{2 L} \\ & \frac{n^{1-\varepsilon}}{n g_k}=\frac{1}{L} \\ & \frac{\left(t n^{\varepsilon-1 / 2}+n^{1 / 2-\varepsilon}\right)-\sqrt{n} g_k}{\sigma_k}=\frac{\left(\frac{3}{2 L}-1\right) \sqrt{n} g_k}{\sigma_k} \end{aligned} $

由引理1，有

$ \begin{aligned} & \;\;\;\;\frac{\left(t n^{\varepsilon-1 / 2}+n^{1 / 2}+n^{1 / 2-\varepsilon}\right) \max _j\left|\hat{g}_{j k}-g_{j k}\right|}{\sigma_k}=O_p\left(n^{1 / 2} \cdot\right. \\ & \left.\max \nolimits_j\left|\hat{g}_{j k}-g_{j k}\right|\right)=O_p\left(n^{1 / 2-\kappa}\right) \end{aligned} $

这一项可以被忽略或被O_p(n^1/2g_k)代替，这是因为在条件1)下，对∀k∈M_τ，n^1/2g_k≥cn^1/2－κ, 从而

$ \begin{array}{l} \;\;\;\;\;\;P\left\{\hat{l}_k<\frac{c^2 n^{1-2 \kappa}}{L^2}\right\} \leqslant P\left\{\hat{l}_k<\frac{n g_k^2}{L^2}\right\} \leqslant \\ P\left\{\frac{1}{\sqrt{n} \sigma_k} \sum\limits_{j=1}^n\left(g_{j k}-g_k\right)<\frac{(3 /(2 L)-1) \sqrt{n} g_k}{\sigma_k}\right\} \end{array} $

进一步由文献[12]的引理1、命题2可知, 存在常数C₁＞0，使得

$ P\left\{\hat{l}_k<\frac{c^2 n^{1-2 \kappa}}{L^2}\right\} \leqslant \exp \left\{-C_1 n^{1-2 \kappa}\right\} $

最后, 对某α∈(0, 1/2－κ)，取L=n^{1/2－κ－α}，则定理1成立。

定理2(确定筛选性质)  在条件1)~3)下，存在常数C₁＞0，对任α∈(0, 1/2－κ)和γ_n=c²n^2α，有

$ P\left\{M_\tau \subset \hat{M}_\tau\right\} \geqslant 1-s_n \exp \left\{-C_1 n^{1-2 \kappa}\right\} $

证明：由定理1及$ P\left\{M_\tau \not \subset \hat{M}_\tau\right\}$=P{存在k∈M_τ, 使得$ \left.\hat{l}_k<c^2 n^{2 \alpha}\right\} \leqslant s_n \max \limits_{k \in M_\tau} P\left\{\hat{l}_k<c^2 n^{2 \alpha}\right\}$，定理2显然成立。

从定理2可知，协变量维数p_n随样本量n的增加呈指数级增长，且满足

$ \ln p_n=O\left(n^{1-2 \kappa}\right) $

则当n→∞时, 有$ P\left\{M_\tau \subset \hat{M}_\tau\right\} \rightarrow 1$，说明估计的重要变量指标集$ \hat{M}_\tau$以概率1包含真实的重要变量指标集M_τ，即所提出的EL-QSIS方法满足确定筛选性质。

1.3 基于边际经验似然的分布函数特征筛选

若关注的活跃指标集为

$ M=\left\{k: F(y \mid \boldsymbol{X}) \text { 依赖于 } X_k, k=1, 2, \cdots, p_n\right\} $

其中F(y|X)=P(Y≤y|X), 则可将所提出的EL-QSIS方法进行推广，得到一种基于边际经验似然的分布函数特征筛选(EL-DFSIS)方法，且该方法不依赖于模型假定。这里令

$ h_k(y, t)=E\left\{[F(y)-I(Y \leqslant y)] I\left(X_k<t\right)\right\} $

则可通过度量E{h_k²(Y, X_k)}是否等于0来进行特征筛选。类似地，令$ {\widetilde g_{jk}} = \left\{ {\frac{1}{n}\sum\limits_{i = 1}^n {\left[ {\hat F\left( {{Y_j}} \right) - I\left( {{Y_i} \le } \right.} \right.} } \right.{\left. {\left. {\left. {{Y_j}} \right)} \right]I\left( {{X_{ik}} < {X_{jk}}} \right)} \right\}^2}$，j=1, 2, …, n, $ \hat F(y) = \frac{1}{n} \cdot \sum\limits_{l = 1}^n I \left( {{Y_l} \le y} \right)$，边际经验似然比$ {\tilde l_k} = 2\sum\limits_{j = 1}^n {\ln } \{ 1 + \left. {\widetilde \lambda {{\widetilde g}_{jk}}} \right\}$, 其中$ {\widetilde \lambda }$满足$ \sum\limits_{j=1}^n \frac{\tilde{g}_{j k}}{1+\tilde{\lambda} \widetilde{g}_{j k}}=0$。基于$ {\tilde l_k}$, M可估计为

$ \hat{M}=\left\{k: \tilde{l}_k \geqslant \tilde{\gamma}_n, k=1, 2, \cdots, p_n\right\} $

其中$ \tilde{\gamma}_n$是预先设定的阈值。

类似地，在一定的正则条件下，仿照文献[12]和上述定理1、2的证明步骤，可以证明EL-DFSIS方法也满足确定筛选性质。

2 数值模拟与实例分析

本节通过数值模拟和实例分析来验证所提出的EL-QSIS、EL-DFSIS筛选方法的有限样本性质，并且分别将它们与QaSIS^[10]、Q-SIS^[11]和SIRS^[7]、DF-SIS^[11]、EL-SIS^[12]等方法进行比较。

在数值模拟中考虑样本量n为150或300，协变量维数p_n=3 000, 筛选出的变量个数d_n=n/ln n, 对每种情形重复300次试验。评价指标包括：p₀—真实的模型大小；P_ALL—在给定模型尺寸d_n下，300次重复试验中所有重要预测变量被选中的比例；Median—300次重复试验中包含所有重要预测变量的最小模型尺寸的中位数；IQR—300次重复试验中包含所有重要预测变量的最小模型尺寸的四分位差。

例1  考虑异方差线性模型

$ Y=X_1+0.8 X_2+0.6 X_3+0.4 X_4+0.2 X_5+\sigma(\boldsymbol{X}) \varepsilon $

式中，X=(X₁, X₂, …, X_pn)^T~N_pn(0, Σ)，Σ=(0.8^|i－j|)(i, j=1, 2, …, p_n)，σ(X)=X₂₀+X₂₁+X₂₂，且误差ε~N(0, 1)或t(4)。考虑分位点τ=0.5或τ=0.75，此时真实的重要预测变量的个数分别为5和8。模拟结果见表 1。

例2  考虑异方差非线性模型

$ Y=X_1^2 \sin X_2+X_3^3+\left(\cos X_4\right)^3+X_5+\sigma(\boldsymbol{X}) \varepsilon $

其他设置条件与异方差线性模型相同，模拟结果见表 2。特别地，在给定模型尺寸d_n下, 表 3给出了协变量X₂₀、X₂₁、X₂₂在300次重复试验中被选中的比例(除去τ=0.5的情形)P₂₀、P₂₁、P₂₂。

从表 1和表 2的结果来看，EL-DFSIS与DF-SIS方法表现相当，并且都要优于SIRS方法，表现为具有更小的最小模型尺寸和更高的重要变量覆盖率。与DF-SIS方法相比，EL-DFSIS方法具有更小的最小模型尺寸IQR；当模型为线性模型时，EL-DFSIS方法与EL-SIS方法表现相当，但当模型为非线性模型时，EL-DFSIS方法能以更高的比例筛选出所有重要预测变量，并且具有更小的最小模型尺寸。随着样本量n的增加，所有方法筛选出的包含所有重要预测变量的比例均显著增加。

当τ=0.5时，EL-QSIS、QaSIS和Q-SIS方法均有良好的筛选表现，最小模型尺寸基本都为5，且重要预测变量被选出的比例都接近于1。但当τ=0.75时，由于存在异方差的影响，QaSIS方法在所有结果中均表现最差，EL-QSIS和Q-SIS方法筛选表现良好，且EL-QSIS方法比Q-SIS方法有更高的重要变量覆盖率和更小的最小模型尺寸IQR，说明EL-QSIS方法对真实模型中重要预测变量的筛选更加准确有效。

表 3展示了例2中对误差产生影响的重要变量被选中的比例。可以发现，与其他方法相比，EL-QSIS和EL-DFSIS方法都能以较高的比例筛选出这3个重要变量。

例3(实例分析)  将EL-DFSIS方法应用到实际数据分析中。考虑文献[16]中的白血病数据集，该数据集包含47例急性淋巴细胞白血病(ALL)患者和25例急性髓系白血病(AML)患者的7 129个基因表达水平数据。本文感兴趣的是识别可能影响白血病种类ALL和AML的基因。Golub等^[16]、Barut等^[17]发现Zyxin、hSNF2b和TCRD(对应的基因为X95735_at、D26156_s_at、U29175_at和M21624_at)为识别两类白血病的重要变量。在给定模型尺寸d=25下，表 4展示了EL-DFSIS与SIRS、DF-SIS、EL-SIS方法对这4个重要基因的筛选结果。

从表中可见这4种方法均能筛选出X95735_at、D26156_s_at和U29175_at这3个重要基因，而EL-DFSIS和EL-SIS还筛选出了重要基因M21624_at，说明基于经验似然的EL-DFSIS筛选方法在实际应用中也具有良好的变量筛选能力。

综上所述，与其他几种方法相比，本文提出的EL-QSIS、EL-DFSIS筛选方法具有良好的筛选降维效果，表现在可以更精确地筛选出重要变量，且能更稳定地用较小的模型尺寸识别出所有的重要变量。

3 结束语

针对超高维异方差数据，本文提出了基于边际经验似然的分位数特征筛选(EL-QSIS)和分布函数特征筛选(EL-DFSIS)方法，两种方法均不依赖于模型假定，无须进行复杂的参数估计和迭代计算，且对分布的假设较宽松。分析了所提方法的确定筛选性质，数值模拟与实例分析的结果表明与QaSIS、Q-SIS和SIRS、DF-SIS、EL-SIS等方法相比，EL-QSIS和EL-DFSIS方法可以更精确地筛选出重要变量，且能更稳定地用较小的模型尺寸识别出所有的重要变量，具有良好的筛选降维效果。未来可考虑将该方法进一步推广到纵向数据、缺失数据的情况。