推薦 ：為你介紹7種流行的線性回歸收縮與選擇方法（附程式碼）

作者：Michał Oleszak 翻譯：張恬鈺 校對：吳金笛

本文5000字，建議閱讀12分鐘。

本文討論了幾種子集和收縮方法：最佳子集回歸, 嶺回歸, LASSO, 彈性網, 最小角度回歸, 主成分回歸和偏最小二乘。

本文討論了七種流行的收縮和選擇方法的數學屬性和實際的Python應用。

在本文中，我們將介紹七種流行的子集選擇和線性回歸收縮方法。在介紹了證明需要這些方法的主題之後，我們將逐一研究每種方法，包括數學屬性和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

譯者簡介：張恬鈺，上海交通大學本科物理專業，Emory University生物統計碩士在讀。以後想繼續在生物統計方向深造。希望能留在美國學習和工作。希望能和廣大的資料愛好者做朋友！分享