一些变量筛选方法——3、部分其它变量筛选方法

由于《An Introduction to Statistical Learning with R》书中的方法书中的方法都是一些比较基础的方法，在做模拟实验以及真实超高维数据时，会出现很多局限性。因此本文后半部分介绍了课本上未提及到的一些方法。这些方法会在后面的模拟实验以及真实数据中进行应用，并且比较书上传统的方法与下述三种方法的真实变量筛选效果。

首先介绍将L0$L^0$范数与L1$L^1$范数相结合的SCAD方法。

---

SCAD(Smoothly Clipped Absolute Deviation)

与岭回归相比，SCAD降低了模型的预测方差，与此同时与Lasso相比，SCAD又缩小了参数估计的偏差，同时它还有很多前面算法所不具备的优秀性质，因而受到了广泛的关注。

SCAD将前面博客提到的gλ(β)=λf(β)$g_{\lambda }(\beta )=\lambda f(\beta )$，变为如下形式：

gλ(β)=⎧⎩⎨⎪⎪⎪⎪λ|βj|,−(|βj|2−2aλ|βj|+λ2)/(2a−2),(a+1)λ2/2,0≤|βj|<λ,λ≤|βj|<aλ,|βj|≥aλ.$$\begin{array}{r}{g}_{\lambda }(\beta )=\{\begin{array}{ll}\lambda |{\beta }_j|,& 0\le |{\beta }_j|<\lambda ,\\ -(|{\beta }_j{|}^2-2a\lambda |{\beta }_j|+{\lambda }^2)/(2a-2),& \lambda \le |{\beta }_j|<a\lambda ,\\ (a+1){\lambda }^2/2,& |{\beta }_j|\ge a\lambda .\end{array}\end{array}$$

其中，λ≥0,a>2$\lambda \ge 0,a>2$，Fan和Li\cite{article11} 建议a取3.7。特别地，若设计矩阵X$X$正交时，SCAD法参数估计显式表达式如下：

β^SCAD=⎧⎩⎨⎪⎪⎪⎪sign(βj^)|βj^|−λ),((a−1)βj^−sign(βj^)aλ)/(a−2),βj^,0≤|βj|<2λ,2λ≤|βj^|<aλ,|βj^|≥aλ.$$\begin{array}{r}{\hat{\beta }}^{SCAD}=\{\begin{array}{ll}sign(\widehat{{\beta }_j})|\widehat{{\beta }_j}|-\lambda ),& 0\le |{\beta }_j|<2\lambda ,\\ ((a-1)\widehat{{\beta }_j}-sign(\widehat{{\beta }_j})a\lambda )/(a-2),& 2\lambda \le |\widehat{{\beta }_j}|<a\lambda ,\\ \widehat{{\beta }_j},& |\widehat{{\beta }_j}|\ge a\lambda .\end{array}\end{array}$$

上图说明了L0$L^0$惩罚，L1$L^1$惩罚，与SCAD三者惩罚之间的差别。可以看出，L0$L^0$方法只会进行变量筛选，不会进行压缩，L1$L^1$（LASSO）既会进行变量筛选，也会对系数继续一定的调整。而SCAD可以从图中很明显的其结合了两种方法，对系数较大的变量不进行惩罚，对系数较少的进行压缩或者删去，因此这种方法既可以筛选变量，也有着Oracle的性质，使其预测效果和真实模型别无二致。

SCAD虽然有相应的迭代算法，但是由于其复杂度高，所以计算速度相对较慢。另外老师上课讲过的将L1$L^1$与L2$L^2$范数相结合的Elastic Net方法\cite{article16}，也是基于前面的一种衍生方法，本文不再进行阐述。

---

SIS（Sure Independence Screening）

当今大数据时代，维数远大于样本量的情况已经非常多见。尽管前面所提出的方法，而且也能一定程度上解决髙维数据问题。但当遇到超高维数据，即维数P无穷大时，上述的算法也会出现问题。针对这类超高维问题，Fan和Lv\cite{article12} 提出了SIS的方法。

针对线性回归模型(2)，按照SIS的思想，首先Y$Y$为中心化向量，计算Y$Y$与每一个自变量xi$x_i$的相关系数，记为

ω=XTY,$$\omega =X^{T}Y,$$

其中 ω=(ω1,⋯,ωp)T$\omega =({\omega }_1,\cdots ,{\omega }_p{)}^T$,若ωi${\omega }_i$越大，说明xi$x_i$与Y$Y$相关性越强。所以，可以根据|ωi|$|{\omega }_i|$的大小来进行变量选择。对任意的γ∈(0,1)$\gamma \in (0,1)$，对|ωi|$|{\omega }_i|$进行从大到小排序，然后取其一个子集

Mγ={1≤i≤p:|ωi|是前[γn]个最大的},$$M_{\gamma }=\{1\le i\le p:|{\omega }_i|\text{是前}[\gamma n]\text{个最大的}\},$$

其中，n$n$是样本数，[γn]$[\gamma n]$是γn$\gamma n$的整数部分，进而保证了[γn]<n$[\gamma n]<n$，与之对应的自变量则入选模型。如果觉得选择[γn]$[\gamma n]$不便于确定，可以选择n−1$n-1$或n/logn$n/\log n$。

而关于相关系数，可以选用自己认为合适的。本文后面的模拟选用传统的Pearson相关系数，以及近几年比较火的可用于检验独立的无参数假设的距离相关性（Distance Covariance），下面其计算公式：

距离相关性（Distance Covariance）

aj,kbj,k=∥Xj−Xk∥,j,k=1,2,…,n,=∥Yj−Yk∥,j,k=1,2,…,n,$$\begin{array}{rl}{a}_{j,k}& =\Vert X_j-X_k\Vert ,j,k=1,2,\dots ,n,\\ b_{j,k}& =\Vert Y_j-Y_k\Vert ,j,k=1,2,\dots ,n,\end{array}$$

其中：||⋅||$||\cdot ||$表示Euclidean范数（欧几里得距离），有：

Aj,k:=aj,k−a¯¯¯j⋅−a¯¯¯⋅k+a¯¯¯⋅⋅,Bj,k:=bj,k−b¯¯j⋅−b¯¯⋅k+b¯¯⋅⋅,$$A_{j,k}:=a_{j,k}-{\overline{a}}_{j\cdot }-{\overline{a}}_{\cdot k}+{\overline{a}}_{\cdot \cdot },B_{j,k}:=b_{j,k}-{\overline{b}}_{j\cdot }-{\overline{b}}_{\cdot k}+{\overline{b}}_{\cdot \cdot },$$

其中：a¯¯¯j⋅${\overline{a}}_{j\cdot }$ 表示由aj,k$a_{j,k}$组成的矩阵，第j$j$行均值，a¯¯¯⋅k${\overline{a}}_{\cdot k}$ 表示第k$k$列均值，以及a¯¯¯⋅⋅${\overline{a}}_{\cdot \cdot }$是X$X$样本中所有数取平均。b$b$的符号标记同a$a$一样，则样本的距离相关性定义为：

dCov2n(X,Y):=1n2∑j=1n∑k=1nAj,kBj,k.(2)$$\begin{array}{}\text{(2)}& {\text{dCov}}_n^2(X,Y):=\frac{1}{n^2}\sum _{j=1}^n\sum _{k=1}^n{A}_{j,k}{B}_{j,k}.\end{array}$$

---

利用随机森林进行变量筛选

其实使用随机森林进行变量筛选是一个比较小众的方法，但其实代表了一类方法。模型本身是用于预测的模型，但在预测过程中，可以对变量重要性进行排序，然后通过这种排序来进行变量筛选。这类方法其实还适用于最近比较火的xgboost，lightgbm等一些非常流行的基于树的机器学习算法，在实际应用中，效果都非常突出。

本文只以较为基础的随机森林中的变量筛选为例：

变量重要性评判用Gini指数为标准，针对一棵树中的每个节点k$k$，我们都可以计算一个Gini指数：

Gk=2p^k(1−p^k),$$G_k=2{\hat{p}}_k(1-{\hat{p}}_k),$$

其中p^k${\hat{p}}_k$表示样本在节点k$k$属于任意一类的概率估计值。

一个节点的重要性由节点分裂前后Gini指数的变化量来确定：

I△k=Gk−Gk1−Gk2,$$I_{\mathrm{△}k}=G_k-G_{k1}-G_{k2},$$

Gk1$G_{k1}$和Gk2$G_{k2}$分别表示Gk$G_k$产生的子节点。针对森林中的每棵树，都用上述的标准来递归产生，最终随机抽取样本和变量，产生森林，假设森林共产生T$T$棵树。

森林中，如果变量Xi$X_i$在第t$t$棵树中出现M$M$次，则变量Xi$X_i$在第t$t$棵树的重要性为：

Iit=∑j=1MI△j.$$I_{it}=\sum _{j=1}^M{I}_{\mathrm{△}j}.$$

则Xi$X_i$在整个森林中的变量重要性为：

I(i)=1n∑t=1TIit.$$I_{(i)}=\frac{1}{n}\sum _{t=1}^T{I}_{it}.$$

最终我们根据变量重要性来选择变量，选择的个数可以用SIS中的方法，选取n−1$n-1$或n/logn$n/\log n$个。

至此，变量筛选的一些方法已进行了简要的概述，包括课本中的以及一些延伸的方法。下面将用模拟实验以及真实数据，来对这些方法进行比较分析。

---

原始对偶**集算法(PDAS)

原始对偶**集算法（Primal Dual Active Set，PDAS）是一个非常新的方法，但做的事情是最优子集选择的事情。其主要思想是引入**集，对所有的β$\beta$进行批量迭代更新。这个方法的优势在于，可以处理超高维数据（上万维），而最优子集选择一旦超过了50维，基本就完全没办法进行运算。后面我们也将采用PDAS来进行模拟。

其算法如下：

1. 给定某固定的T$T$，初始的β0${\beta }^0$，d0=−g(β0)h(β0)$d^0=-{\displaystyle \frac{g({\beta }_0)}{h({\beta }_0)}}$，根据β0${\beta }^0$和d0$d^0$得出A0${\mathcal{A}}^0$、I0${\mathcal{I}}^0$。令k=0$k=0$
2. For k=0,1,2,…,Kmax$k=0,1,2,\dots ,K_{max}$, do
   (2.a) 更新(βk+1,dk+1)$({\beta }^{k+1},d^{k+1})$:
   
   ⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪βk+1Ik=0dk+1Ak=0βk+1Ak=argminl(βAk|Y,XAk)dk+1Ik=−g(βkIk)h(βkIk)$$\{\begin{array}{l}{\beta }_{I^k}^{k+1}=0\\ d_{A^k}^{k+1}=0\\ {\beta }_{A^k}^{k+1}=argminl({\beta }_{A^k}|\mathit{\text{Y}},{\mathit{\text{X}}}_{A^k})\\ d_{I^k}^{k+1}=-{\displaystyle \frac{g({\beta }_{I^k}^k)}{h({\beta }_{I^k}^k)}}\end{array}$$
   
   (2.b) 通过以下方式计算新的**集Ak+1${\mathcal{A}}^{k+1}$和非**集Ik+1${\mathcal{I}}^{k+1}$:
   
   Ak+1={j:−h(βk+1j)−−−−−−−−√|βk+1j+dk+1j|⩾−h(βk+1j)−−−−−−−−√∥βk+1j+dk+1j∥T,∞},$$A^{k+1}=\{j:\sqrt{-h({\beta }_j^{k+1})}|{\beta }_j^{k+1}+d_j^{k+1}|⩾\sqrt{-h({\beta }_j^{k+1})}\Vert {\beta }_j^{k+1}+d_j^{k+1}{\Vert }_{T,\mathrm{\infty }}\},$$
   
   
   Ik+1=(Ak+1)c$$I^{k+1}=(A^{k+1}{)}^c$$
   
   (2.c) 如果Ak+1=Ak${\mathcal{A}}^{k+1}={\mathcal{A}}^k$，则停止迭代；否则令k=k+1$k=k+1$，继续(2.a)和(2.b)步。
   (2.d) 输出β=βk+1$\beta ={\beta }^{k+1}$。

---

后面我们将对前面提及的一些算法进行simulation，以及真实案例操作，敬请期待。