做者 | 宋老師
來源 | JSong的數據科學小站
多重共線性是使用線性迴歸算法時常常要面對的一個問題。在其餘算法中,例如決策樹和貝葉斯,前者的建模過程是逐步遞進,每次拆分只有一個變量參與,這種建模機制含有抗多重共線性干擾的功能;後者乾脆假定變量之間是相互獨立的,所以從表面上看,也沒有多重共線性的問題。可是對於迴歸算法,不管是通常回歸,邏輯迴歸,或存活分析,都要同時考慮多個預測因子,所以多重共線性是不可避免須要面對的,在不少時候,多重共線性是一個廣泛的現象。在構造預測模型時如何處理多重共線性是一個比較微妙的議題。既不能不加控制,又不能一刀切,認爲凡是多重共線性就應該消除。算法
假設有k個自變量的多元線性迴歸模型:小程序
其中偏差項是一個指望值爲0且服從正態分佈的隨機變量:app
則利用最小二乘法可得參數的估計值爲:dom
該求解公式惟一的條件是矩陣X是列滿秩的,否則會有無窮多解:機器學習
當各變量之間存在共線性問題,即各變量之間存在部分線性相關時,例如:函數
易知此時X近乎是不滿秩的(實際狀況很難徹底共線性),X^TX近乎是奇異的,X的最小奇異值會很是小,那它的影響到底有多大呢?咱們先從矩陣計算的角度來看。學習
對於一個方程或者系統而言,當輸入有一個很是微小的擾動時,咱們但願方程或系統的輸出變化也很是微小,若是輸出的變化很是大,且不能被控制,那這個系統的預測就無效了,蝴蝶效應講的就是這個。在矩陣計算中,這叫作擾動分析。測試
【擾動分析定理】設非奇異方陣A知足方程spa
它的精確解爲x* ,當A存在一個小擾動時,假設 $\hat{x}$ 是新方程的解:3d
能夠證實x* 的擾動知足:
能夠看到矩陣的條件數越大,擾動就越大,即x的求解值會變得很是不許確。回到上面講的線性迴歸問題,容易證實最小二乘法的解知足下面的正定方程:
此時
當方程有共線性問題時,X的最小特徵值很是小,相應的,上述的條件數會很是大。也就是說機器學習中的共線性問題實際上就是矩陣計算中的條件數問題。
從實際應用的角度,通常若K<100,則認爲多重共線性的程度很小,如果100<=K<=1000,則認爲存在通常程度上的多重共線性,如果K>1000,則就認爲存在嚴重的多重共線性。
再從統計學的角度來看共線性。能夠證實參數的協方差矩陣爲
又對任意的常數矩陣A和隨機變量x有
代入上式便可得
具體到每一個參數,有:
其中是將第i個變量做爲因變量,其餘k-1個變量做爲自變量進行線性迴歸得到的,且令
爲方差膨脹因子(variance inflation factor,VIF)。當
時,即當第i個變量和其餘變量之間存在線性關係時,VIF趨於無窮大。因此 VIF 的大小反應了變量的共線性程度。通常地,當VIF大於5或10時,認爲模型存在嚴重的共線性問題。
同時考慮參數顯著性檢驗的 t 統計量:
當存在共線性時,參數的標準差偏大,相應的 t 統計量 會偏小,這樣容易淘汰一些不該淘汰的解釋變量,使統計檢驗的結果失去可靠性。
另外考慮線性迴歸的殘差
其中M是一個投影矩陣,且知足
易證實
而矩陣M的範數與X的條件數毫無關係,因而能夠得出共線性並不影響模型的訓練精度。可是對於泛化精度,因爲參數的估計已經不許確啦,因此泛化偏差確定要差些,具體差多少,我還很難用公式表示出來。
總結一下,共線性問題對線性迴歸模型有以下影響:
根據上一節的描述,共線性問題有以下幾種檢驗方法:
當變量數很少,樣本數不是很大時,上述的方法是沒問題的,檢驗某個變量有共線性問題時,能夠結合實際業務考慮直接剔除該變量。可是有的時候變量數大到有上千個,VIF的計算須要創建上千個迴歸模型(條件數僅能斷定是否存在共線性,但不能找到對應的變量),這將耗費很長時間。
事實上咱們能夠從模型角度來直接規避共線性問題。
主成分分析法做爲多元統計分析的一種經常使用方法在處理多變量問題時具備其必定的優越性,其降維的優點是明顯的,主成分迴歸方法對於通常的多重共線性問題仍是適用的,尤爲是對共線性較強的變量之間。當採起主成分提取了新的變量後,每每這些變量間的組內差別小而組間差別大,起到了消除共線性的問題。
逐步迴歸(Stepwise Regression)是一種經常使用的消除多重共線性、選取「最優」迴歸方程的方法。其作法是將逐個引入自變量,引入的條件是該自變量經F檢驗是顯著的,每引入一個自變量後,對已選入的變量進行逐個檢驗,若是原來引入的變量因爲後面變量的引入而變得再也不顯著,那麼就將其剔除。引入一個變量或從迴歸方程中剔除一個變量,爲逐步迴歸的一步,每一步都要進行F 檢驗,以確保每次引入新變量以前回歸方程中只包含顯著的變量。這個過程反覆進行,直到既沒有不顯著的自變量選入迴歸方程,也沒有顯著自變量從迴歸方程中剔除爲止。
嶺迴歸是一種可用於共線性數據分析的有偏估計迴歸方法,它是一種改良的最小二乘估計法,經過放棄最小二乘法的無偏性,以損失部分信息、下降精度爲代價得到迴歸係數更爲符合實際、更可靠的迴歸方法,對條件數很大(病態數據)的擬合要強於最小二乘法。
在線性迴歸問題中,最小二乘法其實是最小化問題:
而嶺迴歸則是加入了L2懲罰項:
這樣參數的方差不會過大,且隨着懲罰項係數C的增大,共線性的影響將越來也小。在這個過程當中,能夠記錄 (嶺跡)的變化狀況,經過對嶺跡的波動來判斷咱們是否要剔除該變量。
那爲何說嶺迴歸能解決共線性問題呢?從矩陣計算的角度來看,L2正則化下方程的解爲:
在上一節咱們講到共線性表明正定矩陣X^T^X的條件數很大:
而當條件數很大時,矩陣的逆的數值計算也是很是不許確的,可是當咱們給矩陣加上一個單位矩陣時,奇異性(不可逆)問題就徹底沒有啦。
進一步考慮對懲罰項對奇異值的影響,假設X的奇異值(SVD)分解爲:
則容易證實
其中D是對角矩陣,且知足
其反應了懲罰項是如何影響到條件數的。
LASSO迴歸和嶺迴歸相似,只不過將懲罰項由L2範數改成了L1範數
L1範數沒有L2範數那麼圓潤,畢竟存在不可導點,並且在L1範數下LASSO迴歸也給不出解析解啦,可是相對於嶺迴歸,LASSO估計的參數能更容易收斂到0
ElasticNet迴歸同時兼顧了L1和L2懲罰項:
當許多變量是相關的時候,Elastic-net是有用的。Lasso通常會隨機選擇其中一個,而Elastic-net則會選在兩個。
除此以外,還有L0範數(非零元的個數)、L1/2範數等。
首先捏造一份好的數據,樣本量爲100,特徵數爲8,且知足方程:
其中偏差項是指望爲0,標準差爲1.5的正態分佈隨機變量。
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])# 係數的均方偏差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,擬合的係數MSE爲0.203657
而後咱們基於這份數據另外構造出兩份數據,第二份數據增長兩個隨機的特徵用做對比,第一份數據則增長兩個共線性特徵:
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)])
先來看下它們的條件數
>>>print(np.linalg.cond(X1)) >>>print(np.linalg.cond(X2)) >>>print(np.linalg.cond(X3))6.29077685383110.9306124087.25066276479
能夠看到X2的條件數很搭,最小奇異值爲0.213,此時還不至於徹底共線性。
拿這兩份數據從新用線性迴歸擬合模型。
lr2=LinearRegression() lr2.fit(X2[training],y[training])# 係數的均方偏差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])# 係數的均方偏差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())
對於第二份共線性構造數據X2,有平均測試集準確率爲0.932070,擬合的參數MSE爲7.697837。能夠看到MSE增長了不少,準確率也降低了0.2%,測試擬合的係數爲:
>>>print(lr2.coef_) [ 10.506618 11.467777 6.35562175 7.56698262 9.44509206 9.81032939 11.66187822 12.29728702 -5.07439399 0.02649089]
在來看對比用的數據X3,其平均測試集準確率爲0.934952,參數MSE爲0.171651,與X1無異。
以上是直接的結果,咱們再來看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、一、二、三、八、9個特徵的VIF都太高。且能夠看出第1個特徵相對第0、二、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時,嶺迴歸估計的係數分別爲
>>>print(coefs[0]) [ 2.70748655 0.95748918 3.53687372 5.2073456 8.70186695 9.84484102 10.67351759 11.74614246 2.46502016 3.19919212]
能夠看到第0、一、二、三、八、9個變量都出現了波動,表明它們之間存在必定的共線性。觀察嶺跡,咱們能夠考慮剔除其中波動比較大的第一、八、9個變量。
另外Lasso迴歸相似,能夠用sklearn中的linear_model.Lasso來學習,這裏就不展現了。最後對於邏輯迴歸任務,sklearn函數內部提供了L1或L2正則化方案,經過它們也能夠去檢測共線性問題。
掃碼進入CDA官方小程序,解鎖更多新鮮資訊和優質內容,還有免費試聽課程,不要錯過喲!