作者 | 宋老師

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

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

1、共線性的原理

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

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

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

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

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

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

1.1 擾動(dòng)分析

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

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

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

可以證明x* 的擾動(dòng)滿足：

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

此時(shí)

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

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

1.2 方差分析

再?gòu)慕y(tǒng)計(jì)學(xué)的角度來(lái)看共線性。可以證明參數(shù)的協(xié)方差矩陣為

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

代入上式即可得

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

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

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

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

同時(shí)考慮參數(shù)顯著性檢驗(yàn)的 t 統(tǒng)計(jì)量：

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

另外考慮線性回歸的殘差

其中M是一個(gè)投影矩陣，且滿足

易證明

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

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

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

2、共線性問(wèn)題的解決方法

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

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

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

事實(shí)上我們可以從模型角度來(lái)直接規(guī)避共線性問(wèn)題。

2.1 PCA等降維法

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

2.2 逐步回歸法

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

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

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

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

在線性回歸問(wèn)題中，最小二乘法實(shí)際上是最小化問(wèn)題：

而嶺回歸則是加入了L2懲罰項(xiàng)：

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

那為什么說(shuō)嶺回歸能解決共線性問(wèn)題呢？從矩陣計(jì)算的角度來(lái)看，L2正則化下方程的解為：

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

而當(dāng)條件數(shù)很大時(shí)，矩陣的逆的數(shù)值計(jì)算也是非常不準(zhǔn)確的，但是當(dāng)我們給矩陣加上一個(gè)單位矩陣時(shí)，奇異性（不可逆）問(wèn)題就完全沒(méi)有啦。

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

則容易證明

其中D是對(duì)角矩陣，且滿足

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

2.4 LASSO回歸

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

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

2.5 ElasticNet回歸等

ElasticNet回歸同時(shí)兼顧了L1和L2懲罰項(xiàng)：

當(dāng)許多變量是相關(guān)的時(shí)候，Elastic-net是有用的。Lasso一般會(huì)隨機(jī)選擇其中一個(gè)，而Elastic-net則會(huì)選在兩個(gè)。

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

3、Python實(shí)踐

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

其中誤差項(xiàng)是期望為0，標(biāo)準(zhǔn)差為1.5的正態(tài)分布隨機(jī)變量。

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)# 測(cè)試集準(zhǔn)確率（R2）print(lr1.score(X1[~training],y[~training]))# 平均測(cè)試集準(zhǔn)確率print(cross_validation.cross_val_score(lr1,X1,y,cv=5).mean())

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

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

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)])

先來(lái)看下它們的條件數(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í)還不至于完全共線性。

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

lr2=LinearRegression()
lr2.fit(X2[training],y[training])# 系數(shù)的均方誤差MSEprint(((lr2.coef_[:8]-coef0)**2).sum()/8)# 測(cè)試集準(zhǔn)確率（R2）print(lr2.score(X2[~training],y[~training]))# 平均測(cè)試集準(zhǔn)確率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)# 測(cè)試集準(zhǔn)確率（R2）print(lr3.score(X3[~training],y[~training]))# 平均測(cè)試集準(zhǔn)確率print(cross_validation.cross_val_score(lr3,X3,y,cv=5).mean())

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

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

在來(lái)看對(duì)比用的數(shù)據(jù)X3，其平均測(cè)試集準(zhǔn)確率為0.934952，參數(shù)MSE為0.171651，與X1無(wú)異。

以上是直接的結(jié)果，我們?cè)賮?lái)看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個(gè)特征的VIF都過(guò)高。且可以看出第1個(gè)特征相對(duì)第0、2、3個(gè)特征的VIF較高。

最后我們?cè)囍媚Ｐ偷姆椒▉?lái)檢測(cè)共線性問(wèn)題

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()

其中當(dāng)alpha取0.1時(shí)，嶺回歸估計(jì)的系數(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個(gè)變量都出現(xiàn)了波動(dòng)，代表它們之間存在一定的共線性。觀察嶺跡，我們可以考慮剔除其中波動(dòng)比較大的第1、8、9個(gè)變量。

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