作者 | 宋老師

來源 | JSong的數(shù)據(jù)科學(xué)小站

多重共線性是使用線性回歸算法時經(jīng)常要面對的一個問題。在其他算法中，例如決策樹和貝葉斯，前者的建模過程是逐步遞進，每次拆分只有一個變量參與，這種建模機制含有抗多重共線性干擾的功能；后者干脆假定變量之間是相互獨立的，因此從表面上看，也沒有多重共線性的問題。但是對于回歸算法，不論是一般回歸，邏輯回歸，或存活分析，都要同時考慮多個預(yù)測因子，因此多重共線性是不可避免需要面對的，在很多時候，多重共線性是一個普遍的現(xiàn)象。在構(gòu)造預(yù)測模型時如何處理多重共線性是一個比較微妙的議題。既不能不加控制，又不能一刀切，認為凡是多重共線性就應(yīng)該消除。

1、共線性的原理

假設(shè)有k個自變量的多元線性回歸模型：

其中誤差項是一個期望值為0且服從正態(tài)分布的隨機變量：

則利用最小二乘法可得參數(shù)的估計值為：

該求解公式唯一的條件是矩陣X是列滿秩的，不然會有無窮多解：

當各變量之間存在共線性問題，即各變量之間存在部分線性相關(guān)時，例如：

易知此時X近乎是不滿秩的（實際情況很難完全共線性），X^TX近乎是奇異的，X的最小奇異值會非常小，那它的影響到底有多大呢？我們先從矩陣計算的角度來看。

1.1 擾動分析

對于一個方程或者系統(tǒng)而言，當輸入有一個非常微小的擾動時，我們希望方程或系統(tǒng)的輸出變化也非常微小，如果輸出的變化非常大，且不能被控制，那這個系統(tǒng)的預(yù)測就無效了，蝴蝶效應(yīng)講的就是這個。在矩陣計算中，這叫做擾動分析。

【擾動分析定理】設(shè)非奇異方陣A滿足方程

它的精確解為x* ，當A存在一個小擾動時，假設(shè) $\hat{x}$ 是新方程的解：

可以證明x* 的擾動滿足：

可以看到矩陣的條件數(shù)越大，擾動就越大，即x的求解值會變得非常不準確?；氐缴厦嬷v的線性回歸問題，容易證明最小二乘法的解滿足下面的正定方程：

此時

當方程有共線性問題時，X的最小特征值非常小，相應(yīng)的，上述的條件數(shù)會非常大。也就是說機器學(xué)習(xí)中的共線性問題實際上就是矩陣計算中的條件數(shù)問題。

從實際應(yīng)用的角度，一般若K<100，則認為多重共線性的程度很小，若是100<=K<=1000，則認為存在一般程度上的多重共線性，若是K>1000，則就認為存在嚴重的多重共線性。

1.2 方差分析

再從統(tǒng)計學(xué)的角度來看共線性?？梢宰C明參數(shù)的協(xié)方差矩陣為

又對任意的常數(shù)矩陣A和隨機變量x有

代入上式即可得

具體到每個參數(shù)，有：

其中是將第i個變量作為因變量，其他k-1個變量作為自變量進行線性回歸獲得的，且令

為方差膨脹因子(variance inflation factor，VIF)。當

時，即當?shù)趇個變量和其他變量之間存在線性關(guān)系時，VIF趨于無窮大。所以 VIF 的大小反應(yīng)了變量的共線性程度。一般地，當VIF大于5或10時，認為模型存在嚴重的共線性問題。

同時考慮參數(shù)顯著性檢驗的 t 統(tǒng)計量：

當存在共線性時，參數(shù)的標準差偏大，相應(yīng)的 t 統(tǒng)計量 會偏小，這樣容易淘汰一些不應(yīng)淘汰的解釋變量，使統(tǒng)計檢驗的結(jié)果失去可靠性。

另外考慮線性回歸的殘差

其中M是一個投影矩陣，且滿足

易證明

而矩陣M的范數(shù)與X的條件數(shù)毫無關(guān)系，于是可以得出共線性并不影響模型的訓(xùn)練精度。但是對于泛化精度，由于參數(shù)的估計已經(jīng)不準確啦，所以泛化誤差肯定要差些，具體差多少，我還很難用公式表示出來。

總結(jié)一下，共線性問題對線性回歸模型有如下影響：

• 參數(shù)的方差增大；
• 難以區(qū)分每個解釋變量的單獨影響；
• 變量的顯著性檢驗失去意義；
• 回歸模型缺乏穩(wěn)定性。樣本的微小擾動都可能帶來參數(shù)很大的變化；
• 影響模型的泛化誤差。

2、共線性問題的解決方法

根據(jù)上一節(jié)的描述，共線性問題有如下幾種檢驗方法：

• 相關(guān)性分析。檢驗變量之間的相關(guān)系數(shù)；
• 方差膨脹因子VIF。當VIF大于5或10時，代表模型存在嚴重的共線性問題；
• 條件數(shù)檢驗。當條件數(shù)大于100、1000時，代表模型存在嚴重的共線性問題。

當變量數(shù)不多，樣本數(shù)不是很大時，上述的方法是沒問題的，檢驗?zāi)硞€變量有共線性問題時，可以結(jié)合實際業(yè)務(wù)考慮直接剔除該變量。但是有的時候變量數(shù)大到有上千個，VIF的計算需要建立上千個回歸模型（條件數(shù)僅能判定是否存在共線性，但不能找到對應(yīng)的變量），這將耗費很長時間。

事實上我們可以從模型角度來直接規(guī)避共線性問題。

2.1 PCA等降維法

主成分分析法作為多元統(tǒng)計分析的一種常用方法在處理多變量問題時具有其一定的優(yōu)越性，其降維的優(yōu)勢是明顯的，主成分回歸方法對于一般的多重共線性問題還是適用的，尤其是對共線性較強的變量之間。當采取主成分提取了新的變量后，往往這些變量間的組內(nèi)差異小而組間差異大，起到了消除共線性的問題。

2.2 逐步回歸法

逐步回歸（Stepwise Regression）是一種常用的消除多重共線性、選取“最優(yōu)”回歸方程的方法。其做法是將逐個引入自變量，引入的條件是該自變量經(jīng)F檢驗是顯著的，每引入一個自變量后，對已選入的變量進行逐個檢驗，如果原來引入的變量由于后面變量的引入而變得不再顯著，那么就將其剔除。引入一個變量或從回歸方程中剔除一個變量，為逐步回歸的一步，每一步都要進行F 檢驗，以確保每次引入新變量之前回歸方程中只包含顯著的變量。這個過程反復(fù)進行，直到既沒有不顯著的自變量選入回歸方程，也沒有顯著自變量從回歸方程中剔除為止。

• 第一：建立全部變量的回歸方程
• 第二：分別建立單獨的回歸方程，依照t檢驗和擬合度依次加入各變量來構(gòu)建回歸方程
• 第三：判斷新引入的變量，對于之前的系數(shù)影響是否顯著，是否符合實際以及對于擬合度的變量，來選擇是否將變量引入模型中。

2.3 嶺回歸、L2正則化(ridge regression)

嶺回歸是一種可用于共線性數(shù)據(jù)分析的有偏估計回歸方法，它是一種改良的最小二乘估計法，通過放棄最小二乘法的無偏性，以損失部分信息、降低精度為代價獲得回歸系數(shù)更為符合實際、更可靠的回歸方法，對條件數(shù)很大（病態(tài)數(shù)據(jù)）的擬合要強于最小二乘法。

在線性回歸問題中，最小二乘法實際上是最小化問題：

而嶺回歸則是加入了L2懲罰項：

這樣參數(shù)的方差不會過大，且隨著懲罰項系數(shù)C的增大，共線性的影響將越來也小。在這個過程中，可以記錄 (嶺跡)的變化情況，通過對嶺跡的波動來判斷我們是否要剔除該變量。

那為什么說嶺回歸能解決共線性問題呢？從矩陣計算的角度來看，L2正則化下方程的解為：

在上一節(jié)我們講到共線性代表正定矩陣X^T^X的條件數(shù)很大:

而當條件數(shù)很大時，矩陣的逆的數(shù)值計算也是非常不準確的，但是當我們給矩陣加上一個單位矩陣時，奇異性（不可逆）問題就完全沒有啦。

進一步考慮對懲罰項對奇異值的影響，假設(shè)X的奇異值（SVD）分解為：

則容易證明

其中D是對角矩陣，且滿足

其反應(yīng)了懲罰項是如何影響到條件數(shù)的。

2.4 LASSO回歸

LASSO回歸和嶺回歸類似，只不過將懲罰項由L2范數(shù)改為了L1范數(shù)

L1范數(shù)沒有L2范數(shù)那么圓潤，畢竟存在不可導(dǎo)點，而且在L1范數(shù)下LASSO回歸也給不出解析解啦，但是相對于嶺回歸，LASSO估計的參數(shù)能更容易收斂到0

2.5 ElasticNet回歸等

ElasticNet回歸同時兼顧了L1和L2懲罰項：

當許多變量是相關(guān)的時候，Elastic-net是有用的。Lasso一般會隨機選擇其中一個，而Elastic-net則會選在兩個。

除此之外，還有L0范數(shù)（非零元的個數(shù)）、L1/2范數(shù)等。

3、Python實踐

首先捏造一份好的數(shù)據(jù)，樣本量為100，特征數(shù)為8，且滿足方程：

其中誤差項是期望為0，標準差為1.5的正態(tài)分布隨機變量。

import numpy as npfrom sklearn.linear_model import LinearRegressionfrom sklearn import cross_validation

coef0=np.array([5,6,7,8,9,10,11,12])
X1=np.random.rand(100,8)
y=np.dot(X1,coef0)+np.random.normal(0,1.5,size=100)
training=np.random.choice([True,False],p=[0.8,0.2],size=100)
lr1=LinearRegression()
lr1.fit(X1[training],y[training])# 系數(shù)的均方誤差MSEprint(((lr1.coef_-coef0)**2).sum()/8)# 測試集準確率（R2）print(lr1.score(X1[~training],y[~training]))# 平均測試集準確率print(cross_validation.cross_val_score(lr1,X1,y,cv=5).mean())

此時平均準確率為0.934955，擬合的系數(shù)MSE為0.203657

然后我們基于這份數(shù)據(jù)另外構(gòu)造出兩份數(shù)據(jù)，第二份數(shù)據(jù)增加兩個隨機的特征用作對比，第一份數(shù)據(jù)則增加兩個共線性特征：

X2=np.column_stack([X1,np.dot(X1[:,[0,1]],np.array([1,1]))+np.random.normal(0,0.05,size=100)])
X2=np.column_stack([X2,np.dot(X2[:,[1,2,3]],np.array([1,1,1]))+np.random.normal(0,0.05,size=100)])
X3=np.column_stack([X1,np.random.rand(100,2)])

先來看下它們的條件數(shù)

>>>print(np.linalg.cond(X1))
>>>print(np.linalg.cond(X2))
>>>print(np.linalg.cond(X3))6.29077685383110.9306124087.25066276479

可以看到X2的條件數(shù)很搭，最小奇異值為0.213，此時還不至于完全共線性。

拿這兩份數(shù)據(jù)重新用線性回歸擬合模型。

lr2=LinearRegression()
lr2.fit(X2[training],y[training])# 系數(shù)的均方誤差MSEprint(((lr2.coef_[:8]-coef0)**2).sum()/8)# 測試集準確率（R2）print(lr2.score(X2[~training],y[~training]))# 平均測試集準確率print(cross_validation.cross_val_score(lr2,X2,y,cv=5).mean())

lr3=LinearRegression()
lr3.fit(X3[training],y[training])# 系數(shù)的均方誤差MSEprint(((lr3.coef_[:8]-coef0)**2).sum()/8)# 測試集準確率（R2）print(lr3.score(X3[~training],y[~training]))# 平均測試集準確率print(cross_validation.cross_val_score(lr3,X3,y,cv=5).mean())

對于第二份共線性構(gòu)造數(shù)據(jù)X2，有平均測試集準確率為0.932070，擬合的參數(shù)MSE為7.697837?？梢钥吹組SE增加了很多，準確率也下降了0.2%，測試擬合的系數(shù)為：

>>>print(lr2.coef_)
[ 10.506618 11.467777 6.35562175 7.56698262 9.44509206
 9.81032939 11.66187822 12.29728702 -5.07439399 0.02649089]

在來看對比用的數(shù)據(jù)X3，其平均測試集準確率為0.934952，參數(shù)MSE為0.171651，與X1無異。

以上是直接的結(jié)果，我們再來看VIF

import matplotlib.pyplot as plt
vif2=np.zeros((10,1))for i in range(10):
 tmp=[k for k in range(10) if k!=i]
 clf.fit(X2[:,tmp],X2[:,i])
 vifi=1/(1-clf.score(X2[:,tmp],X2[:,i]))
 vif2[i]=vifi

vif3=np.zeros((10,1))for i in range(10):
 tmp=[k for k in range(10) if k!=i]
 clf.fit(X3[:,tmp],X3[:,i])
 vifi=1/(1-clf.score(X3[:,tmp],X3[:,i]))
 vif3[i]=vifi 
plt.figure()
ax = plt.gca()
ax.plot(vif2)
ax.plot(vif3)
plt.xlabel('feature')
plt.ylabel('VIF')
plt.title('VIF coefficients of the features')
plt.axis('tight')
plt.show()

可以看到第0、1、2、3、8、9個特征的VIF都過高。且可以看出第1個特征相對第0、2、3個特征的VIF較高。

最后我們試著用模型的方法來檢測共線性問題

from sklearn.linear_model import Ridge
plt.figure()
n_alphas = 20alphas = np.logspace(-1,4,num=n_alphas)
coefs = []for a in alphas:
 ridge = Ridge(alpha=a, fit_intercept=False)
 ridge.fit(X2, y)
 coefs.append(ridge.coef_)
ax = plt.gca()
ax.plot(alphas, coefs)
ax.set_xscale('log')
handles, labels = ax.get_legend_handles_labels()
plt.legend(labels=[0,1,2,3,4,5,6,7,8,9])
plt.xlabel('alpha')
plt.ylabel('weights')
plt.title('Ridge coefficients as a function of the regularization')
plt.axis('tight')
plt.show()

其中當alpha取0.1時，嶺回歸估計的系數(shù)分別為

>>>print(coefs[0])
[ 2.70748655 0.95748918 3.53687372 5.2073456 8.70186695
 9.84484102 10.67351759 11.74614246 2.46502016 3.19919212]

可以看到第0、1、2、3、8、9個變量都出現(xiàn)了波動，代表它們之間存在一定的共線性。觀察嶺跡，我們可以考慮剔除其中波動比較大的第1、8、9個變量。

另外Lasso回歸類似，可以用sklearn中的linear_model.Lasso來學(xué)習(xí)，這里就不展示了。最后對于邏輯回歸任務(wù)，sklearn函數(shù)內(nèi)部提供了L1或L2正則化方案，通過它們也可以去檢測共線性問題。