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

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

---

课本上数据降维方法

标准的回归模型定义为：

Y=β0+β1X1+⋯+βpXp+ϵ, $$Y={\beta }_0+{\beta }_1{X}_1+\cdots +{\beta }_p{X}_p+ϵ,$$

其中自变量为 X1,⋯,Xp $X_1,\cdots ,X_p$ ，维度为 p $p$ 维，因变量为 Y $Y$ 。

首先介绍子集选择法。

---

子集选择法（Subset Selection）

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

Cp：

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

AIC：

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

BIC：

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

Adjusted R2 $R^2$：

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

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

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

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

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

交叉验证法：
交叉验证法是另一种对测试误差进行估计的方法，它是从重抽样的角度进行估计。其基本思想是：讲一个训练集随机分为K份，然后选择第K份作为验证集，然后将剩余的K-1份作为训练集训练模型，这样便可以得到K个“测试误差”，然后求其平均值，便可得到估计的测试误差。其公式如下（符号定义完全按照课本）：

CV(K)=1K∑i=1KMSEi, $$C{V}_{(K)}=\frac{1}{K}\underset{i=1}{\overset{K}{\sum }}MS{E}_i,$$

当 K $K$ 取样本量 N $N$ 时，即为留一交叉验证（Leave-One-Out Cross-Validation，LOOCV）。具体 K $K$ 的大小选取根据我们数据的实际情况而定。通常，样本量大时， K $K$ 取5便足够，并且计算快速；当样本量小的时候， K $K$ 通常取10或者 N $N$ 。

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

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

最优子集法（Best Subset Selection）

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

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

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

逐步回归（Stepwise Selection）

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

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

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

---

系数压缩法（Shrinkage）

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

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

minβ0,β{1N∑i=1N(yi−β0−xTiβ)2} 限制 ∑j=1pfj(βj)≤t, $$\underset{{\beta }_0,\beta }{min}\left\{\frac{1}{N}\underset{i=1}{\overset{N}{\sum }}(y_i-{\beta }_0-x_i^T\beta {)}^2\right\}\text{ 限制 }\underset{j=1}{\overset{p}{\sum }}{f}_j({\beta }_j)\le t,$$

写成向量形式：

minβ∈ℝp{1N‖Y−Xβ‖22} 限制 f(β)≤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形式为：

minβ∈ℝp{1N‖Y−Xβ‖22+λf(β)}, $$\underset{\beta \in {\mathbb{R}}^p}{min}\left\{\frac{1}{N}{\Vert Y-X\beta \Vert }_2^2+\lambda f(\beta )\right\},$$

其中， f(β) $f(\beta )$ 为惩罚项，在岭回归和LASSO时分别取不同的范数，后面我们都将针对Lagrangian形式的式子进行讨论，在本小节的最后，会给出直观的理解。

在最优子集回归中，Lagrangian形式的式子变为：

minβ∈ℝp{1N‖Y−Xβ‖22+λ‖β‖0}, $${\displaystyle \underset{\beta \in {\mathbb{R}}^p}{min}\left\{\frac{1}{N}{\Vert Y-X\beta \Vert }_2^2+\lambda \Vert \beta {\Vert }_0\right\},}$$

其中：

‖β‖p=(∑i=1N|βi|p)1/p 是 Lp 范数，L0范数其实就是变量个数。  $${\displaystyle \Vert \beta {\Vert }_p={\left(\underset{i=1}{\overset{N}{\sum }}|{\beta }_i{|}^p\right)}^{1/p}}\text{ 是 }{L}^p\text{ 范数，}{L}^0\text{范数其实就是变量个数。 }$$

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

LASSO

在LASSO中，Lagrangian形式的式子变为：

minβ∈ℝp{1N‖Y−Xβ‖22+λ‖β‖1}, $${\displaystyle \underset{\beta \in {\mathbb{R}}^p}{min}\left\{\frac{1}{N}{\Vert Y-X\beta \Vert }_2^2+\lambda \Vert \beta {\Vert }_1\right\},}$$

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

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

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

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

岭回归（Ridge Regression）

在岭回归中，Lagrangian形式的式子变为：

minβ∈ℝp{1N‖Y−Xβ‖22+λ‖β‖2}, $${\displaystyle \underset{\beta \in {\mathbb{R}}^p}{min}\left\{\frac{1}{N}{\Vert Y-X\beta \Vert }_2^2+\lambda \Vert \beta {\Vert }_2\right\},}$$

由于 L2 $L^2$范数是光滑的凸函数，所以岭回归的显示解为：

β̂j=(1+Nλ)−1β̂OLSj, $${\hat{\beta }}_j=(1+N\lambda {)}^{-1}{\hat{\beta }}_j^{\text{OLS}},$$

其中

βj^OLS=(XTX)−1XTyj $${\widehat{{\beta }_j}}^{\text{OLS}}=(X^{T}X{)}^{-1}{X}^T{y}_j$$

为标准线性回归解。

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

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

---

降维法（Dimension Reduction）

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

主成分回归（PCR）

主成分回归的思想非常简单，首先使用主成分分析将数据进行变换（这里只采用一种比较易于理解的形式来进行说明PCA），定义原始数据集： Xn×p $\underset{n\times p}{\mathbf{X}}$ ； 线性变换后的主成分： PCn×p $\underset{n\times p}{\mathbf{P}\mathbf{C}}$ ； 参数矩阵： Bp×p $\underset{p\times p}{\mathbf{B}}$ 。

⎧⎩⎨⎪⎪⎪⎪pc1=b11x1+b12x2+⋯+b1pxppc2=b21x1+b22x2+⋯+b2pxp                    ⋯pcp=bp1x1+bp2x2+⋯+bppxp $$\{\begin{array}{l}\mathbf{p}{\mathbf{c}}_{\mathbf{1}}=b_{11}{x}_1+b_{12}{x}_2+\cdots +b_{1p}{x}_p\\ \mathbf{p}{\mathbf{c}}_{\mathbf{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}$$

矩阵形式：

PCn×p=Xn×p⋅BTp×p $$\underset{n\times p}{\mathbf{P}\mathbf{C}}=\underset{n\times p}{\mathbf{X}}\cdot \underset{p\times p}{{\mathbf{B}}^T}$$

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

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

而主成分回归就选取前 k $k$个主成分作为新的自变量，来进行回归:

Y=β0+β1pc1+⋯+βkpck. $$Y={\beta }_0+{\beta }_1\mathbf{p}{\mathbf{c}}_{\mathbf{1}}+\cdots +{\beta }_k\mathbf{p}{\mathbf{c}}_{\mathbf{k}}.$$

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

偏最小二乘回归（PLSR）

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

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

Y=β0+β1Z1+⋯+βmZm. $$Y={\beta }_0+{\beta }_1{Z}_1+\cdots +{\beta }_m{Z}_m.$$

---

在下一篇博客中，我们将会介绍更多、更先进的降维方法，敬请期待～