为你介绍7种流行的线性回归收缩与选择方法（附代码）

在本文中，我们将介绍七种流行的子集选择和线性回归收缩方法。在介绍了证明需要这些方法的主题之后，我们将逐一研究每种方法，包括数学属性和Python应用程序。

为什么收缩或子集，这是什么意思？

在线性回归上下文中，子集意味着从可用变量中选择要包含在模型中的子集，从而减少其维数。另一方面，收缩意味着减小系数估计的大小（将它们缩小到零）。请注意，如果系数缩小到恰好为零，则相应的变量将退出模型。因此，这种情况也可以看作是一种子集。

收缩和选择旨在改进简单的线性回归。关于为什么需要改进，这里有两个原因：

• 预测准确性：线性回归估计倾向于具有低偏差和高方差。降低模型复杂性（需要估计的参数数量）导致减少差异，但代价是引入更多偏差。如果我们能够找到总误差的最佳位置，那么偏差导致的误差加上方差的误差被最小化，这样我们就可以改进模型的预测。

• 模型的可解释性：由于预测变量太多，人类很难掌握变量之间的所有关系。在某些情况下，我们愿意确定影响最大的一小部分变量，为全局着想而牺牲一些细节。

设置和数据加载

在直接跳到方法本身之前，让我们先看看我们将要分析的数据集。它来自Stamey等人的一项研究（1989）。Stamey研究了不同临床测量对前列腺特异性抗原（PSA）水平的影响。任务是根据一组临床和人口统计学变量确定前列腺癌的风险因素。数据以及变量的一些描述可以在Hastie等人的网站以及 “统计学习的要素”教科书的数据部分找到。 

Hastie等人的网站

http://web.stanford.edu/~hastie/ElemStatLearn/

我们将首先导入本文中使用的模块，加载数据并将其拆分为训练和测试集，分别保留目标和特征。然后，我们将讨论每种收缩和选择方法，使其适合训练数据，并使用测试集检查它预测新数据的PSA水平的效果如何。

# Import necessary modules and set options

import pandas as pd

import numpy as np

import itertools

from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV, ElasticNetCV, LarsCV

from sklearn.cross_decomposition import PLSRegression

from sklearn.decomposition import PCA

from sklearn.pipeline import Pipeline

from sklearn.model_selection import GridSearchCV

import warnings

warnings.filterwarnings("ignore") 

# Load data

data = pd.read_csv("prostate_data", sep = "\t")

print(data.head())

# Train-test split

y_train = np.array(data[data.train == "T"]['lpsa'])

y_test = np.array(data[data.train == "F"]['lpsa'])

X_train = np.array(data[data.train == "T"].drop(['lpsa', 'train'], axis=1))

X_test = np.array(data[data.train == "F"].drop(['lpsa', 'train'], axis=1))

线性回归

让我们从简单的线性回归开始，这将构成我们的基准。它将目标变量y建模为p个预测变量或特征X的线性组合：

该模型具有必须从训练数据估计的p + 2个参数：

• p个特征系数β，每个变量一个，表示它们对目标的影响;

• 一个截距参数，表示为上面的β0，它是在所有X都为零的情况下的预测。不一定要将它包含在模型中，实际上在某些情况下应该将其删除（例如，如果想要包含表示分类变量级别的全套虚拟对象），但通常它会给模型更大的灵活性，正如你将在下一段中看到的;

• 高斯误差项的一个方差参数。

通常使用普通最小二乘法（OLS）估计这些参数。 OLS最小化残差平方和，由下式给出

以图形方式考虑这种最小化标准是有帮助的。只有一个预测变量X时，我们处于由预测变量和目标形成的2D空间中。在此设置中，模型在最接近所有数据点的X-Y空间中拟合这样的线：接近度测量为所有数据点的垂直距离的平方和 - 请参见下面的左图。如果有两个预测变量X1和X2，则空间增长到3D，现在模型拟合最接近3D空间中所有点的平面 - 请参见下面的右图。有了两个以上的特征，平面变成了有点抽象的超平面，但这个思想仍然是一样的。这些可视化也有助于了解截距如何为模型提供更大的灵活性：如果包含它，它允许线或平面不跨越空间的原点。

上述最小化问题证明具有解析解，并且β参数可以被计算为

在X矩阵中包括一列1可以表达上述公式中的β帽矢量的截距部分。 “β”上方的“帽子”表示它是基于训练数据的估计值。

偏差-方差权衡

在统计学中，要考虑估计量的两个关键特征：偏差和方差。偏差是真实总体参数和预期估计量之间的差异。它衡量估计数的不准确性。另一方面，方差衡量它们之间的差异。

显然，如果模型太大，偏差和方差都会损害模型的预测性能。然而，线性回归更受到方差的影响，同时具有低偏差。如果模型中存在许多预测特征或者它们彼此高度相关，则尤其如此。这就是用到子集化和正则化来修正的地方。它们允许以引入一些偏差为代价来减少方差，最终减少模型的总误差。

在详细讨论这些方法之前，让我们将线性回归拟合到前列腺数据中并检查其样本外的平均预测误差（MAE）。

linreg_model = LinearRegression(normalize=True).fit(X_train, y_train)

linreg_prediction = linreg_model.predict(X_test)

linreg_mae = np.mean(np.abs(y_test - linreg_prediction))

linreg_coefs = dict(

    zip(['Intercept'] + data.columns.tolist()[:-1],

        np.round(np.concatenate((linreg_model.intercept_, linreg_model.coef_),

        axis=None), 3))

)

print('Linear Regression MAE: {}'.format(np.round(linreg_mae, 3)))

print('Linear Regression coefficients:')

linreg_coefs

最佳子集回归

选择线性回归变量子集的直接方法是尝试所有可能的组合，并选择一个最小化某些标准的组合。这就是最佳子集回归的目标。对于每个k∈{1,2，...，p}，其中p是可用特征的总数，它选择大小为k的子集，其给出最小的残差平方和。然而，平方和不能用作确定k本身的标准，因为它必然随k减小：模型中包含的变量越多，其残差越小。但这并不能保证更好的预测性能。这就是为什么应该使用另一个标准来选择最终模型的原因。对于专注于预测的模型，测试数据上的（可能是交叉验证的）错误是常见的选择。

由于最佳子集回归没有在任何Python包中实现，我们必须手动循环k和k大小的所有子集。以下代码块完成了这项工作。

results = pd.DataFrame(columns=['num_features', 'features', 'MAE'])

# Loop over all possible numbers of features to be included

for k in range(1, X_train.shape[1] + 1):

    # Loop over all possible subsets of size k

    for subset in itertools.combinations(range(X_train.shape[1]), k):

        subset = list(subset)

        linreg_model = LinearRegression(normalize=True).fit(X_train[:, subset], y_train)

        linreg_prediction = linreg_model.predict(X_test[:, subset])

        linreg_mae = np.mean(np.abs(y_test - linreg_prediction))

        results = results.append(pd.DataFrame([{'num_features': k,

                                                'features': subset,

                                                'MAE': linreg_mae}]))

# Inspect best combinations

results = results.sort_values('MAE').reset_index()

print(results.head())

# Fit best model

best_subset_model = LinearRegression(normalize=True).fit(X_train[:, results['features'][0]], y_train)

best_subset_coefs = dict(

    zip(['Intercept'] + data.columns.tolist()[:-1],

        np.round(np.concatenate((best_subset_model.intercept_, best_subset_model.coef_), axis=None), 3))

)

print('Best Subset Regression MAE: {}'.format(np.round(results['MAE'][0], 3)))

print('Best Subset Regression coefficients:')

best_subset_coefs

岭回归

最佳子集回归的一个缺点是它没有告诉我们关于从模型中排除的变量对响应变量的影响。岭回归提供了这种难以选择的变量的替代方案，这些变量将它们分解为模型中包括的和不包括的。相反，它惩罚系数以将它们缩小到零。不完全为零，因为这意味着从模型中移除，但是在零方向上，这可以被视为以连续方式降低模型的复杂性，同时将所有变量保持在模型中。

在岭回归中，线性回归损失函数以这样的方式增强，不仅可以最小化残差平方和，还可以惩罚参数估计的大小：

解决这个最小化问题可得到βs的分析公式：

其中I表示单位矩阵。惩罚项λ是要选择的超参数：其值越大，系数越向零收缩。从上面的公式可以看出，当λ变为零时，加性罚分消失，β-ridge与线性回归中的β-OLS相同。另一方面，随着λ增长到无穷大，β-ridge接近零：在足够高的惩罚下，系数可以任意地收缩到接近零。

但这种收缩是否真的会导致减少模型的方差并以承诺的方式引入一些偏差呢？是的，确实如此，从岭回归估计的偏差和方差的公式中可以清楚地看出：随着λ的增加，偏差也随之增加，而方差则下降！

现在，如何选择λ的最佳值？进行交叉验证尝试一组不同的值，并选择一个最小化测试数据上交叉验证错误的值。幸运的是，Python的scikit-learn可以为我们做到这一点。

ridge_cv = RidgeCV(normalize=True, alphas=np.logspace(-10, 1, 400))

ridge_model = ridge_cv.fit(X_train, y_train)

ridge_prediction = ridge_model.predict(X_test)

ridge_mae = np.mean(np.abs(y_test - ridge_prediction))

ridge_coefs = dict(

    zip(['Intercept'] + data.columns.tolist()[:-1],

        np.round(np.concatenate((ridge_model.intercept_, ridge_model.coef_),

                                axis=None), 3))

)

print('Ridge Regression MAE: {}'.format(np.round(ridge_mae, 3)))

print('Ridge Regression coefficients:')

ridge_coefs

LASSO

Lasso，或最小绝对收缩和选择算子，在本质上与岭回归非常相似。它也为损失函数的非零系数增加了一个惩罚，但与惩罚平方系数之和（所谓的L2惩罚）的岭回归不同，LASSO惩罚它们的绝对值之和（L1惩罚）。因此，对于λ的高值，许多系数在LASSO下完全归零，在岭回归中从未如此。

它们之间的另一个重要区别是它们如何解决这些特征之间的多重共线性问题。在岭回归中，相关变量的系数趋于相似，而在LASSO中，其中一个通常为零，另一个则赋值为整个影响。因此，如果存在相同值的许多大参数，即当大多数预测器真正影响响应时，预期岭回归将有更好的效果。另一方面，当存在少量重要参数且其他参数接近零时，即当只有少数预测因子实际影响响应时，LASSO将占据首位。

然而，在实践中，人们不知道参数的真实值。因此，岭回归和LASSO之间的选择可以基于样本外预测误差。另一种选择是将这两种方法合二为一 - 见下一节！

LASSO的损失函数如下：

与岭回归不同，这种最小化问题无法通过分析解决。幸运的是，有数值算法可以处理它。

lasso_cv = LassoCV(normalize=True, alphas=np.logspace(-10, 1, 400))

lasso_model = lasso_cv.fit(X_train, y_train)

lasso_prediction = lasso_model.predict(X_test)

lasso_mae = np.mean(np.abs(y_test - lasso_prediction))

lasso_coefs = dict(

    zip(['Intercept'] + data.columns.tolist()[:-1],

        np.round(np.concatenate((lasso_model.intercept_, lasso_model.coef_), axis=None), 3))

)

print('LASSO MAE: {}'.format(np.round(lasso_mae, 3)))

print('LASSO coefficients:')

lasso_coefs

弹性网

弹性网首先是对LASSO的批评而产生的，LASSO的变量选择过于依赖数据，因而不稳定。它的解决方案是将岭回归和LASSO的惩罚结合起来，以获得两全其美的效果。弹性网旨在最大限度地减少包括L1和L2惩罚的损失函数：

其中α是岭回归（当它为零时）和LASSO（当它为1时）之间的混合参数。可以使用基于scikit-learn的基于交叉验证的超左侧调整来选择最佳α。

elastic_net_cv = ElasticNetCV(normalize=True, alphas=np.logspace(-10, 1, 400),

                              l1_ratio=np.linspace(0, 1, 100))

elastic_net_model = elastic_net_cv.fit(X_train, y_train)

elastic_net_prediction = elastic_net_model.predict(X_test)

elastic_net_mae = np.mean(np.abs(y_test - elastic_net_prediction))

elastic_net_coefs = dict(

    zip(['Intercept'] + data.columns.tolist()[:-1],

        np.round(np.concatenate((elastic_net_model.intercept_,

                                 elastic_net_model.coef_), axis=None), 3))

)

print('Elastic Net MAE: {}'.format(np.round(elastic_net_mae, 3)))

print('Elastic Net coefficients:')

elastic_net_coefs

最小角度回归

到目前为止，我们已经讨论了一种子集化方法，最佳子集回归和三种收缩方法：岭回归，LASSO及其组合，弹性网络。本节专门介绍位于子集和收缩之间的方法：最小角度回归（LAR）。该算法以空模型开始，所有系数等于零，然后迭代地工作，在每个步骤将一个变量的系数移向其最小二乘值。

更具体地说，LAR从识别与响应最相关的变量开始。然后，它将该变量的系数连续地移向其最小平方值，从而降低其与演化残差的相关性。一旦另一个变量在与残差的相关性方面“赶上”，该过程就会暂停。然后，第二个变量加入有效集，即具有非零系数的变量集，并且它们的系数以保持它们的相关性连接和减少的方式一起移动。继续该过程直到所有变量都在模型中，并以完全最小二乘拟合结束。名称“最小角度回归”来自算法的几何解释，其中给定步骤处的新拟合方向与已经具有非零系数的每个特征形成最小角度。

下面的代码块将LAR应用于前列腺数据。

LAR_cv = LarsCV(normalize=True)

LAR_model = LAR_cv.fit(X_train, y_train)

LAR_prediction = LAR_model.predict(X_test)

LAR_mae = np.mean(np.abs(y_test - LAR_prediction))

LAR_coefs = dict(

    zip(['Intercept'] + data.columns.tolist()[:-1],

        np.round(np.concatenate((LAR_model.intercept_, LAR_model.coef_), axis=None), 3))

)

print('Least Angle Regression MAE: {}'.format(np.round(LAR_mae, 3)))

print('Least Angle Regression coefficients:')

LAR_coefs

主成分回归

我们已经讨论了选择变量（子集）和降低系数（收缩率）的方法。本文中介绍的最后两种方法采用了稍微不同的方法：它们将原始要素的输入空间挤压到较低维度的空间中。主要是，他们使用X创建一小组新特征Z，它们是X的线性组合，然后在回归模型中使用它们。

这两种方法中的第一种是主成分回归。它应用主成分分析，这种方法允许获得一组新特征，彼此不相关且具有高方差（以便它们可以解释目标的方差），然后将它们用作简单线性回归中的特征。这使得它类似于岭回归，因为它们都在原始特征的主成分空间上运行（对于基于PCA的岭回归推导，参见本文底部的Sources中的[1]）。不同之处在于PCR丢弃具有最少信息功能的组件，而岭回归只是将它们收缩得更强。

要重新获得的组件数量可以视为超参数，并通过交叉验证进行调整，如下面的代码块中的情况。

regression_model = LinearRegression(normalize=True)

pca_model = PCA()

pipe = Pipeline(steps=[('pca', pca_model), ('least_squares', regression_model)])

param_grid = {'pca__n_components': range(1, 9)}

search = GridSearchCV(pipe, param_grid)

pcareg_model = search.fit(X_train, y_train)

pcareg_prediction = pcareg_model.predict(X_test)

pcareg_mae = np.mean(np.abs(y_test - pcareg_prediction))

n_comp = list(pcareg_model.best_params_.values())[0]

pcareg_coefs = dict(

   zip(['Intercept'] + ['PCA_comp_' + str(x) for x in range(1, n_comp + 1)],

       np.round(np.concatenate((pcareg_model.best_estimator_.steps[1][1].intercept_, 

                                pcareg_model.best_estimator_.steps[1][1].coef_), axis=None), 3))

)

print('Principal Components Regression MAE: {}'.format(np.round(pcareg_mae, 3)))

print('Principal Components Regression coefficients:')

pcareg_coefs

偏最小二乘法

本文讨论的最终方法是偏最小二乘法（PLS）。与主成分回归类似，它也使用原始要素的一小组线性组合。不同之处在于如何构建这些组合。虽然主成分回归仅使用X自身来创建派生特征Z，但偏最小二乘另外使用目标y。因此，在构建Z时，PLS寻找具有高方差的方向（因为这些可以解释目标中的方差）以及与目标的高相关性。与主成分分析形成对比，主成分分析仅关注高差异。

在算法的驱动下，第一个新特征z1被创建为所有特征X的线性组合，其中每个X由其内积与目标y加权。然后，y在z1上回归，给出PLSβ系数。最后，所有X都相对于z1正交化。然后，该过程重新开始z2并继续，直到获得Z中所需的组件数量。像往常一样，这个数字可以通过交叉验证来选择。

可以证明，尽管PLS根据需要缩小了Z中的低方差分量，但它有时会使高方差分量膨胀，这可能导致在某些情况下更高的预测误差。这似乎是我们的前列腺数据的情况：PLS在所有讨论的方法中表现最差。

pls_model_setup = PLSRegression(scale=True)

param_grid = {'n_components': range(1, 9)}

search = GridSearchCV(pls_model_setup, param_grid)

pls_model = search.fit(X_train, y_train)

pls_prediction = pls_model.predict(X_test)

pls_mae = np.mean(np.abs(y_test - pls_prediction))

pls_coefs = dict(

  zip(data.columns.tolist()[:-1],

      np.round(np.concatenate((pls_model.best_estimator_.coef_), axis=None), 3))

)

print('Partial Least Squares Regression MAE: {}'.format(np.round(pls_mae, 3)))

print('Partial Least Squares Regression coefficients:')

pls_coefs

回顾与结论

由于模型参数存在很大的方差，线性模型具有许多可能相关的特征，导致预测精度和模型的可解释性方面较差。这可以通过减少方差来缓解，这种方差只会以引入一些偏差为代价。然而，找到最佳的偏差 - 方差权衡可以优化模型的性能。

允许实现此目的的两大类方法是子集和收缩。前者选择变量的子集，而后者将模型的系数缩小为零。这两种方法都会降低模型的复杂性，从而导致参数方差的减少。

本文讨论了几种子集和收缩方法：

• 最佳子集回归迭代所有可能的特征组合以选择最佳特征组合;

• 岭回归惩罚平方系数值（L2惩罚），强制它们很小;

• LASSO惩罚系数的绝对值（L1惩罚），这可以迫使它们中的一些精确为零;

• 弹性网结合了L1和L2的惩罚，享受了Ridge和Lasso的精华;

• 最小角度回归适用于子集和收缩之间：它迭代地工作，在每个步骤中添加其中一个特征的“某个部分”;

• 主成分回归执行PCA将原始特征压缩为一小部分新特征，然后将其用作预测变量;

• 偏最小二乘也将原始特征概括为较小的新特征子集，但与PCR不同，它也利用目标构建它们。

正如您在上面运行代码块时从应用程序看到的前列腺数据，大多数这些方法在预测准确性方面表现相似。前5种方法的误差范围在0.467和0.517之间，击败最小二乘误差为0.523。最后两个，PCR和PLS，表现更差，可能是由于数据中没有那么多特征，因此降维的收益是有限的。

谢谢阅读！我希望你学到了新东西！

来源

• Hastie, T., Tibshirani, R., & Friedman, J. H. (2009). The elements of statistical learning: data mining, inference, and prediction. 2nd ed. New York: Springer.

• https://www.datacamp.com/community/tutorials/tutorial-ridge-lasso-elastic-net

原文标题：

A Comparison of Shrinkage and Selection Methods for Linear Regression

原文链接：

https://towardsdatascience.com/a-comparison-of-shrinkage-and-selection-methods-for-linear-regression-ee4dd3a71f16