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

1绍哎、共線性的原理

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

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

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

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

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

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

1.1 擾動分析

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

---

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

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

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

其中

是非奇異方陣的條件數(shù)，且此時矩陣范數(shù)等價于矩陣最大的奇異值龄糊，即矩陣的條件數(shù)等價于 最大奇異值/最小奇異值

---

可以看到矩陣的條件數(shù)越大逆粹，擾動就越大，即x的求解值會變得非常不準確炫惩∑У回到上面講的線性回歸問題，容易證明最小二乘法的解滿足下面的正定方程：

此時

當(dāng)方程有共線性問題時他嚷，X的最小特征值非常小蹋绽，相應(yīng)的搔扁，上述的條件數(shù)會非常大。也就是說機器學(xué)習(xí)中的共線性問題實際上就是矩陣計算中的條件數(shù)問題蟋字。從實際應(yīng)用的角度，一般若K<100扭勉，則認為多重共線性的程度很小鹊奖，若是100<=K<=1000，則認為存在一般程度上的多重共線性涂炎，若是K>1000忠聚，則就認為存在嚴重的多重共線性。

1.2 方差分析

再從統(tǒng)計學(xué)的角度來看共線性唱捣×襟埃可以證明參數(shù)$\theta$的協(xié)方差矩陣為

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

代入上式即可得

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

其中$R_i²$是將第i個變量$x_i$作為因變量震缭，其他k-1個變量作為自變量進行線性回歸獲得的$R²$赂毯，且令

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

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

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

當(dāng)存在共線性時，參數(shù)的標(biāo)準差偏大朝群，相應(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朴肺。當(dāng)VIF大于5或10時，代表模型存在嚴重的共線性問題坚洽；
• 條件數(shù)檢驗戈稿。當(dāng)條件數(shù)大于100、1000時讶舰，代表模型存在嚴重的共線性問題鞍盗。

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

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

2.1 PCA等降維法

主成分分析法作為多元統(tǒng)計分析的一種常用方法在處理多變量問題時具有其一定的優(yōu)越性杠娱，其降維的優(yōu)勢是明顯的挽牢，主成分回歸方法對于一般的多重共線性問題還是適用的，尤其是對共線性較強的變量之間摊求。當(dāng)采取主成分提取了新的變量后禽拔，往往這些變量間的組內(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的增大碗誉，共線性的影響將越來也小。在這個過程中父晶，可以記錄$\theta(k)$(嶺跡)的變化情況哮缺，通過對嶺跡的波動來判斷我們是否要剔除該變量。

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

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

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

進一步考慮對懲罰項對奇異值的影響直撤，假設(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懲罰項：

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

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

3婉陷、Python實踐

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

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

import numpy as np
from sklearn.linear_model import LinearRegression
from 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ù)的均方誤差MSE
print(((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.29077685383
110.930612408
7.25066276479

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

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

lr2=LinearRegression()
lr2.fit(X2[training],y[training])
# 系數(shù)的均方誤差MSE
print(((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ù)的均方誤差MSE
print(((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流礁。可以看到MSE增加了很多罗丰，準確率也下降了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較高。

10個特征的VIF

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

from sklearn.linear_model import Ridge
plt.figure()
n_alphas = 20
alphas = 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()

嶺回歸各個系數(shù)的嶺跡

其中當(dāng)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正則化方案，通過它們也可以去檢測共線性問題刃宵。

參考文獻

[1]. variance inflation factor
[2]. 多重共線性的解決方法之——嶺回歸與LASSO
[3]. ridge regression