ISL 真是好书！

资料

An Introduction to Statistical Learning，下文简称 ISL

概述

标准的线性模型如下：

线性模型在接近实际问题中依然富有竞争力，通过使用别的拟合过程替换最小二乘法，可以改善线性模型的效果（prediction accuracy 和 model interpretability ）。

• Prediction Accuracy : 假设因变量与自变量之间的真实关系就是几乎线性的，那么使用最小二乘法得到的估计值的 bias 会很低。 如果 n >> p ，也就说观测值数目远大于自变量数目，那么最小二乘估计值的 variance 也会很低。因此，如果满足这两个条件，最小二乘方法表现会很好。然而，如果 n 不是远大于 p， 那么最小二乘拟合结果的 variance 很高，容易造成 overfitting，然后在预测新的样本时结果很差。如果 p > n ，那么根本无法得到一个唯一的最小二乘参数估计值，最小二乘法根本无法使用。如果此时我们人为地降低因变量的个数，那么我们的代价就是在降低了 variance 的同时 提高了 bias。
• Model Interpretability : 很多情况下，多重回归的一些自变量其实压根和因变量没关系。将这些无关变量添加到模型中会使得模型增加不必要的复杂度。通过移除这些变量 （即将它们的系数估计值设为0），我们可以得到一个更容易解释的模型。但是，使用最小二乘法估计参数是不可能得到正好是0的参数的（可以接近0，但不是0）。

ISL 这里提到了三类替代最小二乘的方法

• Subset Selection : 这种方法包括确认我们认为与因变量有关的一组自变量（从总共 p 个自变量中挑子集），然后我们用减少了变量个数的这组自变量预测因变量。
• Shrinkage : 这个方法同时使用所有的 p 个变量来拟合模型。但是，所有的最小二乘预测值均向0收缩（称为 regularization）。这种方法可以减少 variance ，有一些系数可能会正好估计为0，从而降低了自变量个数，可以实现变量选择。
• Dimension Reduction : 这个方法包含了将 p 个自变量投射到 (projecting) 一个M维的子空间中，其中 M < p 。

Subset Selection

Best Subset Selection

为了找到 best subset， 我们需要使用 p 的自变量的所有组合来进行一个个单独的最小二乘拟合。因此我们需要进行p次包含一个变量的模型拟合，进行 p(p-1)/2 次包含两个变量的模型拟合，以此类推（总共 2^p 次）。然后我们查看所有模型拟合的结果，找出我们认为是最好的模型。

这个过程的算法表述如下：这里我们挑出所有固定因子数目时的最佳组合，然后最后在比较不同因子数目中的最佳模型。

但是，这里需要注意一点, RSS是会随着变量数目的增加而递减的，无论增加的自变量有没有用。

**因此，如果你使用 RSS 或 R² 作为评价指标，最终最好的模型一定是包含了所有自变量的模型。**如下图所示

这里的问题在于，这里计算得到的 RSS 和 R² 实际衡量地是 training error ，但是我们要的模型应该是需要 test error 最低（training error 倾向于比 test error 更低，trainging error 很低无法保证 test error 很低）。因此在上图中的第三步，我们需要使用交叉验证的预测误差, C_p , BIC , adjusted R² 来选择后 M₀ , M₁ , ……， M_p 。

我们这里提到的思路虽然是用在最小二乘回归中，其他模型也是可以这么做的，比如逻辑回归。不过逻辑回归用的不是 RSS， 而是 deviance 。deviance 等于 -2 times maximized log-likelihood 。值越小，说明拟合程度越好。

虽然 best subset seleciton 方法很简单，逻辑也很清晰，但是它的问题就在于计算量太大。如果 p 大于40，这个算法就基本无法实现了。

Stepwise Selection

出于计算量的原因， 当 p 很大时，best subset selection 无法实现；除此之外，当 p 很大时这个方法还有其他问题。当 p 很大时，很大概率容易出现过拟合现象（模型仅适用于训练集，对于其他数据不适用，参见维度诅咒）

处于上面的原因，逐步回归是一个很好的备选方法。

Forward Stepwise Selection

向前逐步回归指模型从没有自变量出发，然后每次新增一个自变量，至到所有的自变量都加入到模型中。在实际中，每一步拟合到模型中的变量是 gives the greatest additional improvement 的变量。这个算法具体过程如下：

这个算法比上面的 best subset 方法计算量小很多。虽然一般来说向前逐步回归效果很好，但是这种方法并不能保证一定能找到所有组合的最优组合。

例如，加入一个p=3的数据集，最好的一个变量的模型包含 X₁，最好的两个变量的模型包含 X₂ 和 X₃ 。那么向前逐步回归在第二步找最优的二变量模型程序会失败，因为 M₁ 包含 X₁ ，因此 M₂ 必须包括 X₁ 。

这里举了个例子，比对 best subset 方法 和 forward stepwise 方法， 在四变量的时候，best subset 将 rating 替换为了 cards，而 forward stepwise 方法必须要包含 rating 。

Backward Stepwise Selection

与向前迭代回归相反，向后迭代回归从包含所有的变量开始，然后每次剔除一个最没有用的变量。具体算法见下：

同样类似于向前迭代，向后迭代同样不保证能够挑出最好的模型。但是向后迭代要求 n > p 。如果 p 很大，只能使用向前迭代。

Hybrid Approaches

我们可以同时杂合使用向前迭代和向后迭代。我们先按照类似于向前迭代回归的方法来添加因子，但是当我们添加了一个新的变量后，这个方法会剔除不能使得这个模型拟合更好的变量（就是走一步，再退一步）。这种方法最接近于 best subset selection ，并且同时利用了向前迭代回归和向后迭代回归的优势。

其实意思就是说，基本用的方法是向前迭代，但是向前迭代是一种贪婪算法，之前的变量会一直保留。但是在不断新增变量的过程中，可能发生一种现象，那就是最开始加入的变量变得不重要了（变量之间的相关性），这个时候剔除不显著的初始变量更为合理。因此这种混合迭代的方法，每新增一个变量，就会看看旧的变量有没有无用变量，有的话就剔除。这种做法更加合理。

Choosing the Optimal Model

上面的这些方法都会构建很多模型，我们需要一种方法来决定哪个模型是最优的。就像我们提到的，R² 和 RSS 会随着自变量数目的增加而降低，因此不能作为决定自变量数目不同的模型之间比较的指标。这里我们需要能估计 test error 的指标作为评判标准。这里有两种方法：

1. 我们可以i通过对 training error 做出调整，间接估计 test error
2. 我们可以通过对 validation set 或使用交叉验证方法，直接估计 test error

下面会详细介绍这两种方法

C_p, AIC, BIC, and Adjusted R²

这些参数都是通过调整 training error 来估计 test error ，主要是避免了自变量数目会降低training error 的影响。这里主要考虑四种方法：C_p, AIC, BIC, and Adjusted R² 。下图显示了 best subset selection 方法的四个参数的值的变化。

对于包含了 d 个因子的最小二乘模型， 使用 C_p 估计 test MSE 的公式如下（下面是MSE的估计公式）：

这里 σ_hat² 是对因变量预测的残差方差的估计值。

我们可以 C_p 添加了一个惩罚项，如下式，来调整 training error 容易低估 test error 的趋势。而且这个惩罚项会随着自变量数目的增加而增加（含有d），从而趋于弥补 training RSS 随着变量增加而降低的趋势。

由于 σ_hat² 是 σ² 的无偏估计值，因此 C_p 是 test MSE 的一个无偏估计值（?）。

AIC 应用在使用最大似然法来拟合参数的模型中，但是在线性模型中，最大似然法和最小二乘是一回事。在这种情况下，AIC 计算公式如下

我们可以看到 AIC 和 C_p 是一个比例关系，AIC 乘以 σ_hat² 就得到了 C_p ，因此这两个是等价的，所以在上图中只显示了 C_p 的结果。

BIC 来源于贝叶斯的观点，但是最终的公式和 AIC/ C_p 差不多，对于最小二乘模型，公式如下，这里相比较 AIC ，就是把 2 改为了 log(n) （这里其实应该写成 ln(n) , 底为 e, 只要 n > 7, ln(n) 就大于2 ），因此 BIC 对自变量数目很多的模型的惩罚效应更重。因此上面的图中 BIC 的最佳变量数目为 4，而 AIC/ C_p是 6 ，其实在图中我们看到变量数目为4后面的曲线就很平了，变量数为4还是6没啥区别。

校正R2 公式如下，也是会对变量数目过多有惩罚效应

这四种方法均有严谨的理论证明，但是ISL没讲，它们的证明均有一定的前提或假设（如 n 非常大）。虽然校正R² 很流行，但是从理论上校正R² 的效果不如另外三种方法。

Validation and Cross-Validation

除了通过校正 training error 间接估计 test error ，我们也可以通过验证集或交叉验证直接估计 test error ，好处就是需要满足的假设更少，而且可以应用在非常广泛的模型选择任务中（不一定是线性模型）。

在过去，由于执行交叉验证需要很大的计算量无法实现，从而导致 AIC ，BIC 等方法更加流行。但是如今随着计算机算力的增加，交叉验证的计算量已经不再是问题了，因此交叉验证成为了模型选择更好的方法。

下图展示了 BIC， validation set errors , CV （k=10）三种方法的结果，我们可以看后面两种方法的最佳变量数目均是6，但是我们看到 变量数是4 还是 6 差不多。validation set errors 是通过随机挑选3/4的个体作为训练集，剩下的个体就作为验证集。

其实我们可以看到，变量数为 3到11之间的差距很多，如果计算 valiation set errors 时换一种区分参考集和验证集的方式，或者 CV 方法换一个 k 值，那么可能找到的最佳的变量数目会变化。其实，实际我们选模型的方式是 one-standard-error rule ：我们需要同时计算 test MSE 估计值的标准误，然后挑 test MSE 最低点 的 一倍标准误内的最简单的模型。 后面的逻辑就是：如果有很多模型效果差不多好，那么我们应该挑最简单的模型。 在上面的例子中，我们最终挑选的模型就是 3个变量的模型。

Shrinkage Methods

上面提到的那些方法是用来挑选变量的方法。作为另外一种备选方法，我们可以同时将 p 个预测因子均包含在模型中，然后 constrain 或者说是 regularize 系数估计值，或者说是，shrinks the coefficient estimates towards zero。

在这当中，最出名的两个方法就是 ridge regression 和 lasso 。

Ridge Regression

一般的线性模型估计参数的标准是，使得RSS最小，公式如下：

岭回归（Ridge regression）类似于最小二乘，但是它拟合参数用的需要最小化的指标不同，公式如下，其中 λ ≥ 0 ， 称之为 tuning parameter 。每设置一个不同的 λ ，就能得到一套不同的模型参数。因此，选择一个合适的 λ 是至关重要的。

其实就是新增了一部分，称为 shrinkage penalty ， 作用是将所有的参数尽量往0压缩。

注意到， shrinkage penalty 仅仅针对 β₁, …… ， β_p ， 但是不会针对截距 β₀ ，理由显而易见。

岭回归的效果如下图，左边是估计系数随着 λ 的变化的趋势，可以看到随着 λ 的增加，所有因变量的系数均趋向于 0。

右图的横坐标改了，分子/分母的计算方式如下：

因此当 λ = 0，岭回就是普通的线性回归，右图横坐标为1；当 λ 趋于 ∞ 时， 右图横坐标趋于0。

标准的线性回归的系数估计值是 scale equivariant ：如果将自变量 X_j 乘以一个常数 c ，那么这个自变量的系数的估计值就会乘以常数 1/c 。用另外一句话说，无论自变量的值如何缩放，下面这个式子的值永远不变。

但是，与之相反，岭回归的系数估计值会随着自变量系数乘以一个常数发生实质性变化。例如工资这个变量，单位是美元，如果单位换成千美元，那么标准的线性回归的系数估计就会缩小1000倍。这个时候由于岭回归的代价函数中含有系数估计值的平方项，岭回归系数就不会是正好缩小1000倍。换句话说，下式不仅仅取决于 λ 的大小，而且也取决于自变量的缩放比例，甚至还会受到其他自变量缩放比例的影响。

因此，在进行岭回归前，需要对自变量进行标准化。公式如下，分母为估计的标准差。这样所有的自变量有相同 的 scale（不知道啥叫 scale） 。经过下面这个公式，所有自变量均拥有值为1的标准差。

所以上图6.4 的横坐标是 standardized ridge regression coefficient estimates, 也就是说自变量已经经过标准化了。

为什么岭回归比最小二乘更好

这还是来自于 bias-variance trade-off 。随着 λ 增加，岭回归的 flexibility 会降低，导致 variance 降低，但是 bias 增加。如下图所示，当 λ = 0 时， variance 很高，但是没有bias；最佳的 λ 值接近于30 。

一般来说，当因变量与 p个自变量之间的真实关系接近线性关系，那么最小二乘估计值拥有很低的 bias ，但是 variance 很高。这意味着只要 training data 发生一个很小的变化，参数估计值就会发生很大的变化。在实际情况中，如果 p 非常大，接近于 n，那么最小二乘估计的参数会是 extremely variable 。 如果当 p > n， 最小二乘甚至无法给出一个唯一解。这个时候岭回归依然可以发挥作用，仅仅会增加一点 bias，但是会极大地降低 variance 。因此，岭回归适用于 p 很大的情况（最小二乘的 variance 很高）。

The Lasso

岭回归的一个明显的缺点是，它最终的模型中会包含所有p个自变量。虽然它会将所有的系数向0压缩，但是它不会真的把某些系数设置为0，仅仅是接近0。这对于预测没有影响，但是对于模型解释很有影响。

lasso 方法是岭回归的一个替代方法，它克服了上面的这个劣势。lasso 的 代价函数如下：

我们可以看到， lasso 和 岭回归公式相似，仅仅是把系数估计值平方改为绝对值。同样，类似于岭回归，lasso 也同样将系数估计值 shrinks towards zero 。但是，lasso 可以使得某些系数估计值正好是0。因此，类似于 best subset selection 方法， lasso 方法同时可以实现变量选择的目的。我们一般说，the lasso yields sparse models ，也就是说模型仅包含部分向量。同样地，选择一个合适的 λ 值对 lasso 方法也很重要。

下图为 lasso 方法的效果，我们可以只要选择不同的λ值，就可以得到不同的因子数目的模型。

岭回归和 lasso 方法的另一种理解

我们认为 lasso 和 岭回归的系数估计过程，可以视为接近下面这两个问题

也就是说，对于任意一个 λ 值，都会有有一个 s 值作为系数估计值绝对值/平方和的最大值，在满足小于 s 的前提下，RSS 最小的系数就是估计的系数。

best subset selection 可以理解为下式，也就是说条件是系数数目要低于 s 。

当 p 很大时，best subset selection 方法不可用，这时候必须用岭回归或者 lasso 方法。

为什么lasso可以实现变量选择？

下图为p=2时一个解释说明， 红色的圈为RSS等高线，中间的黑点为最小二乘估计值；蓝色图形为 lasso/ 岭回归限制条件的区域，lasso 是一个菱形，岭回归是一个圆形。如果 s 很大，那么lasso/岭回归的限制条件区域就可能包含最小二乘估计值，那么此时解就等于最小二乘估计值。但是我们图中的限制条件区域并不大，左图为 lasso 结果，估计值为RSS等高线与菱形限制区域外边缘刚刚好接触的那个交点，岭回归同理。由于岭回归的限制区域是一个圆形，没有突出点，因此这个交点很难正好在某个轴上；而 lasso 的四个突出点均在轴上，因此交点很有可能就是轴上的突出点，如图中所示。

在更高的维度上，可能就会有很多的系数估计值都是0。

比较 lasso 和 岭回归

我们很容易看到 lasso 方法相比于岭回归有很大的优势，产生的模型更容易解释，但是哪种方法的预测准确性更好呢？下图中右图为比较岭回归和lasso的结果，横坐标均为 training data 的 R² ，这是一种很好的比较不同的正则方法的方式。我们可以岭回归和 lasso 的 bias 差不多，但是岭回归的 variance 更低。这个图是由于所有的45个自变量都和因变量有关系，所以岭回归效果更好。

如果我们用的模拟数据，因变量只和45个变量中的2个变量有关系，这时我们可以从下图看到，lasso 预测效果更好。

这说明**这两种方法之间在预测准确性上没有绝对的优劣之分。**如果因变量实际与大量的预测因子均有关系，那么岭回归效果更好。但是，与因变量有关的因子数目在实际项目中是一个未知数，但是我们可以通过交叉验证来决定哪种方法更好。

一个说明岭回归和lasso的简单例子

为了更好地理解岭回归和 lasso 的特性，我们考虑一个特殊的例子，其中 n = p， 然后 X 是一个单位矩阵。为了进一步简化问题，假设我们正在执行没有截距的回归。根据上面这些假设，最小二乘问题简化为寻找使得下式最小的一组系数

在这个例子中，最小二乘解为

在这种情况下，岭回归的解为下式最小的一组系数估计值

lasso 的公式如下

岭回归的解为

lasso 的解为

下图展示了这种情况。我们可以看到 岭回归和 lasso 的 shrink 策略不一样。在岭回归中，所有的系数的压缩比例都是相同的；与之相反，lasso 将正负 λ/2 之间的系数都直接压缩为0，其他范围的系数压缩比例相同。

岭回归和 lasso 的贝叶斯解释

我们现在可以用一个 贝叶斯方法的视角来重新审视岭回归和 lasso 。从贝叶斯方法的角度出发，首先我们需要知道回归系数的先验分布，我们称为 p(β) 。 数据的似然性可以写作 f(Y|X, β) 。 先验分布乘以似然值 ，我们就得到了后验分布，形式如下（?）

一般的线性模型如下：

我们假设残差之间是独立的，而且来自于一个相同的正态分布。还有，我们假设下式成立，其中g() 为一个密度函数，岭回归和 lasso 服从两个特殊的密度分布。

我们可以发现岭回归和 lasso 服从两个不同的 g() 函数。

• 如果 g 是一个均值为0，标准差为 λ 函数的正态分布，那么 β 的 posterior mode （基于给定数据的最可能的 β值）就是岭回归的解。换句话说，岭回归的解就是这个时候的后验分布的均值。
• 如果 g 是 一个 双指数分布 (double-exponential) (也叫拉普拉斯分布) ， 均值为0，scale parameter 是 λ 的函数，那么这时的 posterior mode 就是 lasso 的解。

这两种分布如下图所示。

选择 λ 的值

岭回归和 lasso 均需要选择一个合适的值。交叉验证提供一个简单的解决办法，通过计算不同的 λ 的 CV errors ，来选择一个 CV errors 最低的 λ 值，如下图所示，左图为不同的λ值计算得到的 CV errors，右图为不同的 λ 值计算得到的岭回归标准化系数。我们通过左图可以看到，最佳的 λ 值的 CV error 降低幅度并不明显，左图的左半边基本都是平的，这种情况一般我们会直接用最小二乘法，而不用岭回归/lasso。

下图为 lasso 的十倍交叉验证与系数估计结果。左图为交叉验证结果；右图为系数估计结果，两根有颜色的线为与因变量有关的因子（signal variables），灰线表示无关变量（noise variables）。这里我们发现最佳的 λ 值比较大，而且最佳 λ 值只保留了两个显著变量，因此证明了交叉验证的结果很好地给出了一个合适的 λ 值，这很有挑战，因为这个例子中有 p=45 个因子，但是只有 n=50 个观测值。作为比较，最小二乘结果（右图横坐标为1）仅仅给一个显著变量分配了一个较大的系数（红线）。

Dimension Reduction Methods

上面提到的两种方法，subset selection 和 shrinkage methods ，都还是使用的原来的预测因子，X₁ , X₂ …… X_p 。 下面会提到如何 transform the predictors ，然后将转换后的预测因子用于最小二乘模型。我们一般将这些方法称为 dimension reduction methods。

我们假设 Z₁ , Z₂ …… Z_M 表示 M 个原来 p 个预测因子的线性组合 (linear combinations) ，其中 M < p 。即：

其中 Φ_jm 为常数。

把上面转换后的预测因子，带入到线性回归模型中，如下

这样就从估计 p+1 个系数降低至估计 M+1 个系数。根据上面的公式进行推导得到

因此，

通过降维的方法可能会增加系数估计的 bias，但是当相比对 n 而言，p 很大时，通过降维的方法可以显著降低 variance 。

所有的降维方法均分为两步。第一步，获得 Z₁ , Z₂ …… Z_M 这些新的预测因子；第二步，通过这些 M 个因子来拟合模型。但是如何获取新的预测因子，或者说如何得到这些 Φ_jm 的值，这存在很多方法。这里只考虑两个主要的方法：principal components (主成分) 和 partial least squares 。

Principal Components Regression

主成分分析 (principal components analysis, PCA) 是一个流行的降维方法。这里只讨论 PCA 在回归分析中的作用。第一主成分是指观测值在这个成分上的变异最丰富的成分。例如下图为两个因子的散点图，横坐标为城市的人口数量，纵坐标为某个特定公司的广告费用，这里总共是100个城市的结果。绿色实线表示第一主成分方向 (first principal component direction)，我们可以如果把所有的点均投影到这条线上，投影到第一主成分上得到的点的方差是最大的。如果我们把这些点投射到别的直线上，得到的变异均会更低。

第一主成分的公式如下，即 Φ₁₁ = 0.839 和 Φ₂₁ = 0.544 。( Φ₁₁ 等称为 principal component loadings ) （我不太清楚上面的直线和这个公式之间有什么关系）

我大致把 Z1 的三维图片画了出来，是一个平面，代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

install.packages("barsurf")
library("barsurf")

install.packages("plotrix")
library("plotrix")

install.packages("scatterplot3d")
library("scatterplot3d")

install.packages("plot3D")
library("plot3D")

install.packages("MBA")
library("MBA")

x <- seq(10, 70)
y <- seq(0,35)
f <- function(x, y) {
  r <- 0.839*(x-40) + 0.544*(y-15)
  return(r)
}
z <- outer(x, y, f)
op <- par(bg = "white")
nrz <- nrow(z)
ncz <- ncol(z)
jet.colors <- colorRampPalette(c("gray80", "gray10"))
nbcol <- 100
color <- jet.colors(nbcol)
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
facetcol <- cut(zfacet, nbcol)

persp(x, y, z,
  theta = 30, phi = 30,
  expand = 0.5, col = color[facetcol]
)

图片如下：

在所有的可能的线性组合中（满足条件， Φ₁₁² + Φ₂₁² = 1 ），只有上面的这种线性组合可以产生最高的方差（Var(Z₁) 最大）。

这里 pop 和 ad 都是长度为100的向量，那么 Z₁ 也是一个长度为100的向量。例如，对于第 i 个样本，公式如下

这些值 Z₁₁, …… , Z_n1 称为 principal component scores 。

我们对于 PCA 还有另外一种解释：第一主成分定义的直线是和数据越近越好。如下图的左图，数据与第一主成分的距离为直线距离，体现为黑色虚线，直线距离平方和最小的直线为第一主成分直线（这不就是最小二乘回归嘛）。右图为将左图向右下方旋转，旋转到第一主成分和横坐标方向一致。

第二主成分 Z₂ 需要和 Z₁ 无关，如图 6.14 所示（两条线垂直）。 第二主成分公式如下：

因为这里我们只有两个预测因子，因此前两个主成分就包含了全部信息。但是，这里第一主成分包含了最多的信息。图 6.15 右图看出，second principal component scores 几乎为0，表示第二主成分几乎没有什么信息了。

下图为第一主成分与第二主成分与两个性状的相关。

The Principal Components Regression Approach

主成分回归 (PCR) 包括创建前 M 个主成分 Z₁ , Z₂ …… Z_M ，然后用这些主成分作为预测因子加入到回归模型中。这么做的基础是，我们认为少数的几个主成分，可以体现预测因子的大部分变异，最终体现在与反应变量的关系上。

下图显示了两次模拟数据中主成分回归中使用的主成分数量的影响。模拟数据中 n=50, p=45，左图中的模拟数据设置所有预测因子都有影响，右图中的模拟数据只有2个预测因子有作用。我们看主成分回归确实有作用。

但是，同时比对 PCR，岭回归，lasso 三种方法。我们可以看到 PCR 方法相比于另外两种方法并没有优势。

一般来说，仅在当只有少数几个主成分可以捕获自变量的绝大部分变异（同时可以捕获因变量的绝大部分变异）时，PCR 分析效果才比较好。比如图 6.19 的情况，前5个主成分的 bias 就降低至0，这种情况下 PCR 分析更加合适。

我们需要注意地是，虽然 PCR 降低了估计参数的数目，但是它不是feature selection method ，主成分回归分析还是用到了所有的 p 个自变量。因此，PCR 分析从原理上更接近于岭回归，而不是 lasso。

在实际的主成分回归分析中，M 的值通常通过 CV 获得。Credit 数据的PCR分析结果如下图，我们可以看到最佳的M = 10，这就相当于没有进行什么降维，因为 p=11。

当我们进行 PCR 分析时，我们一般建议对预测因子进行 standardizing

如果不进行标准化，那么 high-variance 的变量会在主成分分析中倾向于占据更重要的地位，会对最终的模型结果产生影响。但是，如果所有的变量均采用相同的单位，那么你也可能不选择进行标准化这一步。

Partial Least Squares

PCR 分析中需要确认线性组合Z₁ , Z₂ …… Z_M ，来表示这些预测因子 。这个过程我们认为是 unsupervised ，因为这个过程因变量 Y 不会参与，或者说不会 supervise the identification of the principal components 。因此，PCR 分析有一个问题，可以很好地解释自变量变异的主成分不见得可以很好地预测因变量。

我们这里推出，偏最小二乘回归 , 一种替代 PCR 分析的 监督方法。PLS 方法创建新变量的过程时会利用因变量 Y 的值。粗略地说，PLS 方法倾向于找到同时可以解释自变量和因变量的新变量。

PLS 首先将 p 个自变量和因变量进行标准化。然后PLS 的第一主成分的 Φ_j1 设置为 Y ~ X_j 的简单线性回归的系数（公式如下）。因此，第一主成分的主要权重给了与因变量高度相关的变量。

下图显示了 p=2 的 PLS 与 PCR 的第一主成分的直线，绿色实线为 PLS 的第一主成分，虚线为 PCR 的第一主成分。PLS 的方向添加了 pop 的权重，降低了ad 的权重，表明 pop 与因变量的关系更强。

PLS 在 chemometrics (化学计量学) 上很流行。实际上，PLS 的效果不见得比 岭回归或 PCR 方法更好。虽然 PLS 可以减少 bias ，但是它可能会增加 variance ，

Considerations in High Dimensions

High-Dimensional Data

大部分传统的统计方法都是适用于 low_dimensional 的情况，即 n >> p 。 因为以前的计算任务都是这种情况的，因此发展起来的统计方法也是仅适用于这种情况的。

在过去的20年中，新技术改变了数据收集的方式，因此 p 非常大的情况也就出现了，但是由于收集数据是有成本的，因此 n 一般是有限的。ISL 书中举了两个例子

1. 我们可能会通过50万个 SNP 来预测血压性状，这里 n ≈ 200，p ≈ 500,000
2. 一个市场调研公司对客户的上网购物模式感兴趣，可能会把用户在搜索引擎上的所有搜索过的词汇当作特征 (features) 。对于某个特定用户，任何一个 search term 会被标记为 present (0) 和 absent (1) ，会创建一个非常大的二变量的 feature vector 。这里 n ≈ 1000，但是 p 非常大。

p>n (或 p≈n )的数据集一般称为 high-dimensional 。传统方法并不适合用于分析这些数据，例如最小二乘线性回归，原因来自于 bias-variance trade-off 和 overfitting 。

What Goes Wrong in High Dimensions?

这里我们用最小二乘回归举例。当 p 和 n 一样大时，最小二乘回归无法实现。理由很简单，无论因变量与自变量的真实关系如何，最小二乘法的结果都会完美拟合训练集，以至于残差均为0。

下图为 p=1 的两种情况，左边 n =20，右边 n=2，可以看到右边完美拟合了这两个点（无论这两个点是什么情况）。这就是 overfitting ，如果在 test data 里这个拟合结果就会很差。

下图进一步说明了p很大时最小二乘分析的风险，其中 n = 20， p 从 1-20 之间变化，所有的自变量均和因变量完全没有关系。如图，随着变量数目的增加，R² 逐渐趋于1，而 Training MSE 趋于0。由于所有的自变量都是噪音变量，随着自变量的提高， Test MSE 也是逐渐提高的。

如果分析人员很粗心，只看 training data 的情况，那么很容易得到包含所有变量的模型最好的结论。如果 p 很大，C_p , AIC , BIC 方法也不合适，因为估计的残差方差是有问题的，这里就是0。校正R² 也有问题，由于 RSS 是0，校正R² 就会是1。

Regression in High Dimensions

在 ISL 书中其实提到了很多应用于 p 很大 时的回归方法，例如 forward stepwise selection， 岭回归， lasso， 主成分回归。这些都可以用，都可以避免 overfitting 。

下图说明了 lasso 的表型。这里 p =20, 50, 2000 ，其中只有20个自变量是有效的， n=100 。然后通过一个无关的 test set 来评估效果。下面三个图，横坐标用的是自由度（系数个数，包括截距），不用 λ , 便于理解；纵坐标为 test MSE。

从下图我们可以得到三点：

1. lasso 有作用
2. 合适的λ值很重要
3. 随着 p 的增加， test error 倾向于增加，除非新增的因子和反应变量真的有关系。

第三点实际上是高维数据分析的一个重要的原则，称之为维度陷阱 。我们可能会本能地认为随着自变量数目的增加，模型拟合的效果会越好。 但是，比对下图地最左图与最右图结果，我们发现情况相反，p=2000的 test error 几乎是 p=20 的两倍。

一般来说，新增 signal features 会增加模型拟合效果；但是新增 noise features 会降低模型拟合效果。因为新增 noise features 会增加维度，使得 overfitting 的风险剧增。

因此，可以收集成千上万 features 的新技术是一把双刃剑：如果这些 features 都和关心的问题有关，那么就可以改善模型；但是如果这些 features 是噪音变量，那么把这些 features 加入到模型中就可能带来负优化。

Intercepting Results in High Dimensions

当我们对 high-dimensional setting 执行岭回归，lasso 等分析时，在解释结果时需要小心。在前面，我们提到了 multicollinearity ，意思是自变量之间的相关关系。在高维数据中， multicollinearity 的问题更加极端：所有的自变量都可以写成其他自变量的线性组合。这意味着，我们永远不知道哪些变量是真的与因变量有关，哪些仅仅是由于与signal features 相关导致与因变量相关。

例如，假设我们想要通过 50万个SNP来预测血压，forward stepwise selection 方法选择了 17个 SNP 加入到预测模型中。我们可能很容易错误地推断出这17个 SNP 预测血压地效果比其他 SNP 更有效。实际上，可能你可以选出很多个 17 SNPs sets，预测效果都和你选的模型一样。如果我们采用一个不同的 training data，我们可能会得到一个完全不同的 SNP 集合，这并不是说我们之前得到的模型是错的或无效的，事实上可能之前的模型非常有效，仅仅是我们不要过分解读我们获得的结果，要认识到我们仅仅是得到了 one of many possible models for predicting blood pressure , 而且需要后期更多数据的验证。

另外，在评价高维数据的模拟拟合效果的时候，需要注意之前的所有的训练集的统计量已经全部失效了。因此，很重要地是，我们需要采用 test set 的结果，或 cross-validation errors。

R代码

Subset Selection Methods

Best Subset Selection

查看数据，清除缺失值

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

> library(ISLR)
Warning message:
程辑包‘ISLR’是用R版本4.0.5 来建造的 
> fix(Hitters)
> names(Hitters)
 [1] "AtBat"     "Hits"      "HmRun"     "Runs"      "RBI"       "Walks"     "Years"    
 [8] "CAtBat"    "CHits"     "CHmRun"    "CRuns"     "CRBI"      "CWalks"    "League"   
[15] "Division"  "PutOuts"   "Assists"   "Errors"    "Salary"    "NewLeague"
> dim(Hitters)
[1] 322  20
> sum(is.na(Hitters$Salary))
[1] 59
> Hitters = na.omit(Hitters)
> dim(Hitters)
[1] 263  20
> sum(is.na(Hitters))
[1] 0

leaps 包 的 regsubsets() 函数可以执行 best subset selction，每个因子数目下，最好的模型通过 RSS 确定（regsubsets() 默认的最大因子数目为8，可以通过nvmax 设定）。最终结果可以通过 summary() 查看

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

> library(leaps)
Warning message:
程辑包‘leaps’是用R版本4.0.5 来建造的 
> regfit.full = regsubsets(Salary ~ ., Hitters)
> # 默认最大8个因子
> summary(regfit.full)
Subset selection object
Call: regsubsets.formula(Salary ~ ., Hitters)
19 Variables  (and intercept)
           Forced in Forced out
AtBat          FALSE      FALSE
Hits           FALSE      FALSE
HmRun          FALSE      FALSE
Runs           FALSE      FALSE
RBI            FALSE      FALSE
Walks          FALSE      FALSE
Years          FALSE      FALSE
CAtBat         FALSE      FALSE
CHits          FALSE      FALSE
CHmRun         FALSE      FALSE
CRuns          FALSE      FALSE
CRBI           FALSE      FALSE
CWalks         FALSE      FALSE
LeagueN        FALSE      FALSE
DivisionW      FALSE      FALSE
PutOuts        FALSE      FALSE
Assists        FALSE      FALSE
Errors         FALSE      FALSE
NewLeagueN     FALSE      FALSE
1 subsets of each size up to 8
Selection Algorithm: exhaustive
         AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns CRBI CWalks
1  ( 1 ) " "   " "  " "   " "  " " " "   " "   " "    " "   " "    " "   "*"  " "   
2  ( 1 ) " "   "*"  " "   " "  " " " "   " "   " "    " "   " "    " "   "*"  " "   
3  ( 1 ) " "   "*"  " "   " "  " " " "   " "   " "    " "   " "    " "   "*"  " "   
4  ( 1 ) " "   "*"  " "   " "  " " " "   " "   " "    " "   " "    " "   "*"  " "   
5  ( 1 ) "*"   "*"  " "   " "  " " " "   " "   " "    " "   " "    " "   "*"  " "   
6  ( 1 ) "*"   "*"  " "   " "  " " "*"   " "   " "    " "   " "    " "   "*"  " "   
7  ( 1 ) " "   "*"  " "   " "  " " "*"   " "   "*"    "*"   "*"    " "   " "  " "   
8  ( 1 ) "*"   "*"  " "   " "  " " "*"   " "   " "    " "   "*"    "*"   " "  "*"   
         LeagueN DivisionW PutOuts Assists Errors NewLeagueN
1  ( 1 ) " "     " "       " "     " "     " "    " "       
2  ( 1 ) " "     " "       " "     " "     " "    " "       
3  ( 1 ) " "     " "       "*"     " "     " "    " "       
4  ( 1 ) " "     "*"       "*"     " "     " "    " "       
5  ( 1 ) " "     "*"       "*"     " "     " "    " "       
6  ( 1 ) " "     "*"       "*"     " "     " "    " "       
7  ( 1 ) " "     "*"       "*"     " "     " "    " "       
8  ( 1 ) " "     "*"       "*"     " "     " "    " "       
> # 给定最大因子数
> regfit.full = regsubsets(Salary ~., data=Hitters, nvmax=19)
> reg.summary = summary(regfit.full)

summary() 也可以返回 R² , 校正R² ，C_p 等结果，我们可以根据这些统计量挑选最终的模型。

1
2
3
4
5
6
7

> names(reg.summary)
[1] "which"  "rsq"    "rss"    "adjr2"  "cp"     "bic"    "outmat" "obj"   
> #查看 R2 结果
> #可以看到随着因子数量增加，R2无脑增加
> reg.summary$rsq
 [1] 0.3214501 0.4252237 0.4514294 0.4754067 0.4908036 0.5087146 0.5141227 0.5285569 0.5346124 0.5404950 0.5426153 0.5436302
[13] 0.5444570 0.5452164 0.5454692 0.5457656 0.5459518 0.5460945 0.5461159

1
2
3
4
5
6
7
8
9
10

#将所有统计量画在一张图上
par(mfrow=c(2,2))
plot(1:19,reg.summary$rss, xlab = "NUmber of Variables",ylab = "RSS", type = "l")
plot(1:19,reg.summary$adjr2, xlab = "NUmber of Variables",ylab = "Adjusted RSq", type = "l")
which.max(reg.summary$adjr2)
points(11,reg.summary$adjr2[11], col="red", cex=2,pch=20)
plot(1:19,reg.summary$cp, xlab = "NUmber of Variables",ylab = "Cp", type = "l")
points(which.min(reg.summary$cp), reg.summary$cp[which.min(reg.summary$cp)], col="red", cex=2, pch=20)
plot(1:19,reg.summary$bic, xlab = "NUmber of Variables",ylab = "BIC", type = "l")
points(which.min(reg.summary$bic), reg.summary$bic[which.min(reg.summary$bic)], col="red", cex=2, pch=20)

手动画图结果如下：

regsubsets() 函数本身有一个内建函数 plot() ，也可以画图

1
2
3
4
5
6

# regsubsets() 内建函数
plot(regfit.full, scale="r2")
plot(regfit.full, scale="adjr2")
plot(regfit.full, scale="Cp")
plot(regfit.full, scale="bic")

但是画出的图看不懂

Forward and Backward Stepwise Selection

我们也可以用 regsubsets() 执行这两种算法，只要新增一个参数 method

1
2
3
4
5
6

# 向前/后
regfit.fwd = regsubsets(Salary~., data=Hitters, nvmax=19, method='forward')
summary(regfit.fwd)
regfit.bwd = regsubsets(Salary~., data=Hitters, nvmax=19, method='backward')
summary(regfit.bwd)

Choosing Among Models Using the Validation Set Approach and Cross-Validation

上面是通过 C_p , BIC, adjusted R² 这三个统计量来选择模型，这里是用验证集或交叉验证来实现。

这两种方法的思路都是，先找到最佳因子数目，然后用全部数据和 best subset 方法，找最佳因子数目下的模型。

为了适用验证集检验的方式，我们首先要把总数据拆分成训练集和验证集。首先根据训练集找到每个因子数量下的最好模型，再用不同因子数量的模型预测 test data， 计算 test MSE，找 test MSE 最小的因子数量。最后需要用全部数据使用 best subset 重新估计最佳因子数目的模型，因为用的数据越多，估计的参数越准，前面训练集估计参数只用了一部分数据（有可能最后一次找到的因子和前面训练集的因子不一样）。

10倍交叉验证的方法：首先随机分成10分，每一次随机挑出一份作为验证集，另外九份作为训练集。每一次计算出所有因子水平的 test MSE，最终结果是一个 10×19 的矩阵，其中元素 (i, j) 是第i 个 fold 的 j 个因子水平下最好模型的 test MSE。计算每一列的均值，查找最小值，即为最好的因子数目。最后需要用全部数据使用 best subset 重新估计最佳因子数目的模型，比如这里交叉验证得到的最佳因子数目为 11， 然后得到最终模型的代码如下：

Ridge Regression and the Lasso

这里我们用 glmnet 包来实现岭回归和 lasso 算法。

model.matrix() 创建因子矩阵特别好用，会自动将分类变量转化为哑变量。glmnet() 只接受哑变量和连续变量。

1
2
3
4
5

# 生成因子与因变量

x = model.matrix(Salary~., Hitters)[,-1]
y = Hitters$Salary

Ridge Regression

glmnet() 有一个 alpha 参数可以指定模型，如果 alpha=0 ，那么就是岭回归；如果 alpha=1，那么就是 lasso。

我们先看岭回归

1
2
3
4

# 岭回归
library(glmnet)
grid = 10^seq(10,-2,length=100)
ridge.mod = glmnet(x,y,alpha=0,lambda=grid)

glmnet() 函数一般会指定选择 λ 值的范围。这里是想指定λ的范围为 10¹⁰ 到 10^-2 ，覆盖了从只包含截距到普通的最小二乘回归的所有范围。glmnet() 默认会对所有变量进行标准化，如果想要关闭这一点，可以使用 standardize=FALSE 。

每一个λ值都会有相应的岭回归系数，可以通过coef() 查看，这里是20×100的矩阵，20 是因子数（包括截距），100 是 λ 值的数目。

1
2

> dim(coef(ridge.mod))
[1]  20 100

查看某一个λ值的结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14

> ridge.mod$lambda[50]
[1] 11497.57
> coef(ridge.mod)[,50]
  (Intercept)         AtBat          Hits         HmRun          Runs           RBI 
407.356050200   0.036957182   0.138180344   0.524629976   0.230701523   0.239841459 
        Walks         Years        CAtBat         CHits        CHmRun         CRuns 
  0.289618741   1.107702929   0.003131815   0.011653637   0.087545670   0.023379882 
         CRBI        CWalks       LeagueN     DivisionW       PutOuts       Assists 
  0.024138320   0.025015421   0.085028114  -6.215440973   0.016482577   0.002612988 
       Errors    NewLeagueN 
 -0.020502690   0.301433531 
> #计算岭回归的惩罚项
> sqrt(sum(coef(ridge.mod)[-1,50]^2))
[1] 6.360612

我们可以用交叉验证的方式来选择 λ 值，我们可以用 cv.glmnet() 函数，这个函数默认采取十倍的交叉验证，可以通过 nfolds 参数修改。这里我们通过 set.seed() 函数使我们的结果可重复，因为 CV 中参验群体的划分是随机的。通过交叉验证找到了最佳的的 λ 值后，我们需要使用全部数据，使用最佳的λ值来估计模型系数。

代码略

The Lasso

lassso 算法，只要改成 alpha = 1 即可。