一些变量筛选方法——2、《An Introduction to Statistical Learning with R》上的数据降维方法

前面提到，这里介绍的变量筛选的方法全部是基于《An Introduction to Statistical Learning with R》书中的方法，所以这里先开始介绍课本上的变量筛选方法，然后再进行延伸。

---

课本上数据降维方法

标准的回归模型定义为：
$$Y={\beta }_0+{\beta }_1{X}_1+\cdots +{\beta }_p{X}_p+ϵ,$$
其中自变量为$X_1,\cdots ,X_p$，维度为$p$维，因变量为$Y$。

首先介绍子集选择法。

---

子集选择法（Subset Selection）

在子集选择法中，每一个自变量的子集会对应一个模型，而子集选择法就是对这些子集所构成的模型进行比较分析，从中选出最优的模型，以此来进行变量筛选。为了选出最优，我们会用到一些评判指标，包括：Cp，AIC，BIC，以及Adjusted $R^2$等，我们将会用这些评价指标来对测试误差进行估计，下面先介绍这些评判指标（注：这里完全采用书上的符号）：

Cp：
$$C_p=\frac{1}{n}(RSS+2d{\hat{\sigma }}^2),$$

AIC：
$$AIC=\frac{1}{n{\hat{\sigma }}^2}(RSS+2d{\hat{\sigma }}^2),$$

BIC：
$$BIC=\frac{1}{n}(RSS+log(n)d{\hat{\sigma }}^2),$$

Adjusted $R^2$：

$$Adjusted{R}^2=1-\frac{RSS/(n-d-1)}{TSS/(n-1)}.$$

从上式中可以看出，当我们的对于误差为已知方差的正态模型，使用AIC和使用Cp两个准则选择出的模型完全一样。但当方差为未知需要进行估计时，相比于只适用于平方损失的线性模型的Cp准则，AIC可以适用于更多的模型，它对模型的变量个数（复杂度）进行了惩罚。但是在实际使用中，AIC更倾向于选择比真实模型更多的变量，这样就容易导致“过拟合”。

而基于贝叶斯想法得到的BIC，惩罚更重，比AIC更为严苛，尤其是对维度特别高时。虽然它的形式和AIC准则非常类似，但它对变量更多，复杂度更高的模型惩罚更重，所以使用BIC准则进行判别则容易“欠拟合”。

Adjusted $R^2$则是调整过后的离差平方和，也是用来判断模型的优劣。

上面的想法都是从似然函数与模型复杂度角度考虑的。Adjust $R_2$越接近1说明模型拟合得越好，其他三个指标则是越小越好。在实际使用中，就AIC和BIC来进行对比，当样本量非常大时，AIC倾向于选择更多的变量；BIC亲向于选择出正确的那个子模型。而当样本量很小的时候，BIC倾向于选择更少的变量。

交叉验证法：
交叉验证法是另一种对测试误差进行估计的方法，它是从重抽样的角度进行估计。其基本思想是：讲一个训练集随机分为K份，然后选择第K份作为验证集，然后将剩余的K-1份作为训练集训练模型，这样便可以得到K个“测试误差”，然后求其平均值，便可得到估计的测试误差。其公式如下（符号定义完全按照课本）：
$$C{V}_{(K)}=\frac{1}{K}\sum _{i=1}^{K}MS{E}_i,$$
当$K$取样本量$N$时，即为留一交叉验证（Leave-One-Out Cross-Validation，LOOCV）。具体$K$的大小选取根据我们数据的实际情况而定。通常，样本量大时，$K$取5便足够，并且计算快速；当样本量小的时候，$K$通常取10或者$N$。

这种做法可以更加贴近真实的测试误差，所以通常而言，效果会更加出色。但其缺点是计算量较大（LOOCV用于线性的情况，可以化简出显示表达式，计算量不增加）。

下面开始介绍具体的方法：

最优子集法（Best Subset Selection）

最优子集法的思想是将所有的特征组合都进行建模，然后选择最优的模型（最优的判断依据都是前面叙述的几种指标）。

假设数据$p$维，
$\bullet$ 从$p$个变量中任意选择$k$个，共有$C_n^k$种情况，从中选择最优的一个模型$M_k$;
$\bullet$ 再从$M_k,k=1,2,\cdots ,p$中选择最优的模型;
最后选择出的模型，即为最终的模型，其中的变量，即为我们筛选出的变量。

其优点：遍历各种可能的情况，选择出最优的模型；缺点：当维数$p$非常大时，算法复杂度为$O(2^p)$，其呈指数上升。因此这种方法只适用于维数较小的情况。

逐步回归（Stepwise Selection）

逐步回归法分为向前逐步回归与向后逐步回归。其主要思想是：每一次迭代都只能沿着上一次迭代的方向继续进行。向前逐步回归是指初始模型只有常数项，然后逐渐添加变量；向后逐步回归是指初始模型包含了所有变量，然后逐渐删除变量。（注：向前与向后逐步回归筛选出的变量可能不一样，但其思想完全一样，这里只以向前逐步回归法为例。）

假设数据$p$维，
$\bullet$ 对于$p$个变量，$k$从1取到$p$;
$\bullet$ 从$p$个变量中任意选择$k$个，共有$C_n^k$种情况，从中选择最优的一个模型$M_k$;
$\bullet$ 基于上一步最优模型的$k$个变量，再选择加入一个变量，变量从剩下$p-k$个中选取，这样就可以构建$p-k$个模型，从中选取最优的模型;
$\bullet$ 重复以上过程，直到$k=p$迭代完成，然后从$p$个模型中选择最优的模型;
最后选择出的模型，即为最终的模型，其中的变量，即为我们筛选出的变量。

这种方法与最优子集选择法的差别在于，最优子集选择法在第$k$步，可以选择任意$k$个变量进行建模，从而可以保证全局最优；而逐步回归法只能基于之前所选的$(k-1)$个变量进行$k$轮建模。所以逐步回归法不能保证全局最优，因为前面的变量选择选中不重要的变量，但在后面的迭代中也必须加上，从而非最优变量集合。但逐步回归的优势就是计算量大幅度减小，算法复杂度为：$O(\frac{p(p+1)}{2})$，因此针对维数较大的时候更为适用。

---

系数压缩法（Shrinkage）

前面所叙述的方法都是较为传统的方法，从其提出的年代也可以看出，并且子集选择的一个非常大的不足之处是它的不稳定性，即变量选择的结果会由于数据集合的微小变化而发生很大的变化。基于这种情况，后面提出了系数压缩法（Shrinkage），这种方法也就是基于惩罚项变量选择方法。其基本原理是：在极大似然或最小二乘目标函数的基础上，添加一个关于模型复杂度的惩罚函数，构造一个新惩罚目标函数，再对新的惩罚目标函数求最值（最大值或最小值）得到参数估计值。这类方法能大大节省计算时间，降低子集选择法不稳定性带来的风险。

包括前面的最优子集回归，我们都可以叫正则化分析，它都可以统一到下述框架：

$$\underset{{\beta }_0,\beta }{\min }\left\{\frac{1}{N}\sum _{i=1}^N(y_i-{\beta }_0-x_i^T\beta {)}^2\right\}\text{ 限制 }\sum _{j=1}^p{f}_j({\beta }_j)\le t,$$
写成向量形式：
$$\underset{\beta \in {\mathbb{R}}^p}{\min }\left\{\frac{1}{N}{\Vert Y-X\beta \Vert }_2^2\right\}\text{ 限制 }f(\beta )\le t,$$

Lagrangian形式为：
$$\underset{\beta \in {\mathbb{R}}^p}{\min }\left\{\frac{1}{N}{\Vert Y-X\beta \Vert }_2^2+\lambda f(\beta )\right\},$$
其中，$ f(\beta) $为惩罚项，在岭回归和LASSO时分别取不同的范数，后面我们都将针对Lagrangian形式的式子进行讨论，在本小节的最后，会给出直观的理解。

在最优子集回归中，Lagrangian形式的式子变为：
$$\underset{\beta \in {\mathbb{R}}^p}{\min }\left\{\frac{1}{N}{\Vert Y-X\beta \Vert }_2^2+\lambda \Vert \beta {\Vert }_0\right\},$$

其中：
$$\Vert \beta {\Vert }_p={\left(\sum _{i=1}^N|{\beta }_i{|}^p\right)}^{1/p}\text{ 是 }{L}^p\text{ 范数，}{L}^0\text{范数其实就是变量个数。 }$$

前面我们也能发现，这是一个NP-hard，想要找到一个全局最优解非常困难。分支定界算法虽然能避免全局遍历求解，但是速度与稳定性仍然不理想。所以后面有统计学家提出了LASSO的方法。

LASSO

在LASSO中，Lagrangian形式的式子变为：
$$\underset{\beta \in {\mathbb{R}}^p}{\min }\left\{\frac{1}{N}{\Vert Y-X\beta \Vert }_2^2+\lambda \Vert \beta {\Vert }_1\right\},$$

由于$\Vert \beta {\Vert }_1$是凸函数，所以这是一个凸优化问题。并且该式仍可以将一些系数$\beta$收缩到0。现在用的比较多的算法是坐标下降法（Coordinate Descent），它可以在R包glmnet中实现。当然，由于$L^1$范数和$L^0$范数之间的差异性，这两者的结果还是有一些差距的。因此后来的工作是希望找到一种近似$L^0$范数的惩罚项，并且保证该惩罚项有和$L^1$范数一样优秀的性质。这其中，SCAD\cite{article11} 等算法比较出名，但是他们都不是凸函数，并且计算速度慢，一般算法不能保证得到全局最优解。Adaptive Lasso 试图在Lasso的基础上增加一个预先训练出的系数作为权重调整，以使Lasso趋近于$L^0$范数，然而这个权重在高维的情况下并不好确定。因此，Lasso一直保持着很好的实用性。

由于该算法的时间复杂度与特征数$p$是线性关系，因此在处理这种高维特征数据的时候效率很高，并且能够保证最终达到Lasso的最优解。

Lasso由于添加了$L^1$范数的惩罚项，所以其实略微增加了偏差，但却减少了方差，在实际数据中表现非常突出，所以到如今，很多领域仍然在广泛地使用这个方法。

当我们接着将惩罚项变为$L^2$范数惩罚的时候，此时的回归就变为了岭回归。

岭回归（Ridge Regression）

在岭回归中，Lagrangian形式的式子变为：
$$\underset{\beta \in {\mathbb{R}}^p}{\min }\left\{\frac{1}{N}{\Vert Y-X\beta \Vert }_2^2+\lambda \Vert \beta {\Vert }_2\right\},$$

由于$L^2$范数是光滑的凸函数，所以岭回归的显示解为：
$${\hat{\beta }}_j=(1+N\lambda {)}^{-1}{\hat{\beta }}_j^{\text{OLS}},$$
其中
$${\widehat{{\beta }_j}}^{\text{OLS}}=(X^{T}X{)}^{-1}{X}^T{y}_j$$
为标准线性回归解。

其解相当于对标准线性回归的参数估计进行了压缩，这样有效解决了求解过程中矩阵不可逆，或者原本的解非常不稳定的情况，但由于$L^2$范数的特性，这样的解没有稀疏性的特征，达不到变量筛选的目的。

上图解释了为什么岭回归达不到变量筛选的目的，而LASSO可以。因为LASSO的惩罚项是$L^1$范数，在左侧图中表示正方形区域，最小二乘目标函数表示在图上则是椭圆表示。而正方形在坐标轴上有尖点，所以当$\hat{\beta }$发生变化时，很容易碰到尖点。而一旦碰到尖点，则表示有变量的系数被估计为0，从而达到变量筛选的目的。反观右侧的岭回归示意图，由于惩罚项是$L^2$范数，所以在二维平面是是个圆，没有尖点所以相交在坐标轴上的概率为0，因此只能实现变量系数压缩而不能实现筛选。

---

降维法（Dimension Reduction）

课本上提及的降维法是指将变量进行变换后的新变量进行降维，这里主要分为未利用$Y$信息的主成分回归（无监督学习），以及利用了$Y$信息的偏最小二乘回归（有监督学习）。由于这两种方法其实不是属于变量筛选的方法，所以本文只进行简单的叙述。

主成分回归（PCR）

主成分回归的思想非常简单，首先使用主成分分析将数据进行变换（这里只采用一种比较易于理解的形式来进行说明PCA），定义原始数据集： $\underset{n\times p}{X}$ ； 线性变换后的主成分：$\underset{n\times p}{\mathrm{PC}}$ ； 参数矩阵：$\underset{p\times p}{B}$ 。
$$\{\begin{array}{l}\mathbf{p}{\mathbf{c}}_1=b_{11}{x}_1+b_{12}{x}_2+\cdots +b_{1p}{x}_p\\ \mathbf{p}{\mathbf{c}}_2=b_{21}{x}_1+b_{22}{x}_2+\cdots +b_{2p}{x}_p\\ \cdots \\ \mathbf{p}{\mathbf{c}}_{\mathbf{p}}=b_{p1}{x}_1+b_{p2}{x}_2+\cdots +b_{pp}{x}_p\end{array}$$

矩阵形式：
$$\underset{n\times p}{\mathrm{PC}}=\underset{n\times p}{X}\cdot \underset{p\times p}{{\mathbf{B}}^T}$$

首先我们在限制$\stackrel{p}{\sum _{i=1}}{b}_{1i}^2=1$的条件下最大化$\mathbf{p}{\mathbf{c}}_1$的方差。接着，我们再限制$\stackrel{p}{\sum _{i=1}}{b}_{2i}^2=1$的条件下最大化$\mathbf{p}{\mathbf{c}}_2$的方差，并限制它与第一主成分是不相关的。然后第$k$个主成分需要在 $\stackrel{p}{\sum _{i=1}}{b}_{ki}^2=1$的条件在最大化$\mathbf{p}{\mathbf{c}}_{\mathbf{k}}$，并限制它与$\mathbf{p}{\mathbf{c}}_{\mathbf{j}},j=1,\dots ,k-1$不相关。

这个过程会持续到所有主成分被计算出来为止，越靠前的主成分会有越大的方差，并且提供了更多有助于数据分析的信息。

而主成分回归就选取前$k$个主成分作为新的自变量，来进行回归:
$$Y={\beta }_0+{\beta }_1\mathbf{p}{\mathbf{c}}_1+\cdots +{\beta }_k\mathbf{p}{\mathbf{c}}_{\mathbf{k}}.$$

这样做就是对数据先降维再回归，但是数据降维的过程完全没有利用$Y$的信息。相比之下，偏最小二乘回归则利用了$Y$信息降维之后进行回归。

偏最小二乘回归（PLSR）

类似PCR，PLSR也是可以一步一步进行计算来实现。

首先我们同样限制$\stackrel{p}{\sum _{i=1}}{\varphi }_{i1}^2=1$的条件下计算$Z_1={\sum }_{i=1}^p{\varphi }_{i1}{X}_i$，其中${\varphi }_{i1}$就是每个$X_i$与$Y$的Pearson相关系数，然后进行标准化的结果。接着，我们再将每个$X_i$减去$Z_1$得到新的$X_i^{(1)}$，再限制$\stackrel{p}{\sum _{i=1}}{\varphi }_{i2}^2=1$，计算$Z_2={\sum }_{i=1}^p{\varphi }_{i2}{X}_i^{(1)}$，它与$Z_1$是正交的。然后以此类推，直到计算出$Z_m,m\le p$。最后我们再用新得到的变量进行回归，以达到降维的目的：
$$Y={\beta }_0+{\beta }_1{Z}_1+\cdots +{\beta }_m{Z}_m.$$

---

在下一篇博客中，我们将会介绍更多、更先进的降维方法，一些变量筛选方法——3、部分其它变量筛选方法