講講共線性問題

作者：JSong，時間：2017.12.18，公眾號：JSong老師

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

1、共線性的原理

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

\[y=\theta_0+\theta_1x_1+\cdots+\theta_kx_k+\varepsilon=X\theta+\varepsilon \]

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

\[\varepsilon\sim\mathcal{N}(0,\sigma^2) \]

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

\[\hat{\theta}=X^{\dagger}y=(X^TX)^{-1}X^Ty \]

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

\[\hat{\theta}=X^{\dagger}y+(I-X^{\dagger}X)\xi \]

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

\[x_3=x_2+x_1+\varepsilon \]

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

1.1 擾動分析

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

---

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

\[Ax=y \]

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

\[(A+\delta\,A)\hat{x}=y \]

可以證明 \(x^*\) 的擾動滿足：

\[\frac{||\delta\,x||}{||\hat{x}||}\leq\kappa(A)\frac{||\delta\,A||}{||A||} \]

可以看到矩陣的條件數越大，擾動就越大，即x的求解值會變得非常不准確。回到上面講的線性回歸問題，容易證明最小二乘法的解滿足下面的正定方程：

\[X^{T}X\hat{\theta}=X^{T}y \]

此時

\[\kappa(X^{T}X)=\frac{\lambda_{max}(X^{T}X)}{\lambda_{min}(X^{T}X)}=\frac{\sigma^2_{max}(X)}{\sigma^2_{min}(X)} \]

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

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

1.2 方差分析

再從統計學的角度來看共線性。可以證明參數\(\theta\)的協方差矩陣為

\[Var(\hat{\theta})=Var(\hat{\theta}-\theta)=Var[(X^TX)^{-1}X^T\varepsilon] \]

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

\[Var(Ax)=A\cdot\,Var(x)\cdot\,A^T \]

代入上式即可得

\[Var(\hat{\theta})=\sigma^2(X^{T}X)^{-1} \]

具體到每個參數，有：

\[Var(\hat{\theta_i})=\frac{\sigma^2}{(n-1)Var(x_j)}\cdot\frac{1}{1-R_i^2} \]

其中\(R_i^2\)是將第i個變量\(x_i\)作為因變量，其他k-1個變量作為自變量進行線性回歸獲得的\(R^2\)，且令

\[\text{VIF}_i=\frac{1}{1-R_i^2} \]

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

\[R_i^2\longrightarrow\,1 \]

時，即當第i個變量和其他變量之間存在線性關系時，VIF趨於無窮大。所以 VIF 的大小反應了變量的共線性程度。一般地，當VIF大於5或10時，認為模型存在嚴重的共線性問題。

同時考慮參數顯著性檢驗的 t 統計量：

\[t=\frac{\hat{\theta_i}}{std(\hat{\theta_i})}\sim\,t(n-k-1) \]

當存在共線性時，參數的標准差偏大，相應的 t 統計量 會偏小，這樣容易淘汰一些不應淘汰的解釋變量，使統計檢驗的結果失去可靠性。

另外考慮線性回歸的殘差

\[\hat{\varepsilon}=y-X\hat{\theta}=M\varepsilon \]

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

\[M=I-X(X^{T}X)^{-1}X^{T} \]

易證明

\[||\hat{\varepsilon}||_{2}^{2}=\varepsilon^{T}M\varepsilon\leq||M||_F^2\cdot||\varepsilon||_2^2=(n-k)||\varepsilon||_2^2 \]

而矩陣M的范數與X的條件數毫無關系，於是可以得出共線性並不影響模型的訓練精度。但是對於泛化精度，由於參數的估計已經不准確啦，所以泛化誤差肯定要差些，具體差多少，我還很難用公式表示出來。

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

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

---

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

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

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

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

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

2.1 PCA等降維法

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

2.2 逐步回歸法

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

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

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

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

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

\[\hat{\theta}=\arg_{\theta}\min||X\theta-y||^2_2 \]

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

\[\hat{\theta}=\arg_{\theta}\min||X\theta-y||^2_2+C||\theta||^2_2 \]

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

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

\[\hat{\theta}=(X^{T}X+CI)^{-1}X^{T}y \]

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

\[\kappa(X^{T}X)=\frac{\lambda_{max}(X^{T}X)}{\lambda_{min}(X^{T}X)} \]

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

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

\[X=U\Sigma\,V^{T} \]

則容易證明

\[\hat{\theta}=(X^{T}X+CI)^{-1}X^{T}y=VDU^{T}y \]

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

\[D_{ii}=\frac{\sigma_i}{\sigma_i^2+C} \]

其反應了懲罰項是如何影響到條件數的。

2.4 LASSO回歸

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

\[\hat{\theta}=\arg_{\theta}\min||X\theta-y||^2_2+C||\theta||^2_1 \]

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

2.5 ElasticNet回歸等

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

\[\hat{\theta}=\arg_{\theta}\min||X\theta-y||^2_2+C_1||\theta||^2_1+C_2||\theta||^2_2 \]

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

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

3、Python實踐

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

\[y=5x_0+6x_1+7x_2+8x_3+9x_4+10x_5+11x_6+12x_7+\varepsilon \]

其中誤差項是期望為0，標准差為1.5的正態分布隨機變量。

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])
# 系數的均方誤差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，擬合的系數MSE為0.203657

然后我們基於這份數據另外構造出兩份數據，第二份數據增加兩個隨機的特征用作對比，第一份數據則增加兩個共線性特征：

\[x_8=x_0+x_1+\varepsilon \]

\[x_9=x_1+x_2+x_3+\varepsilon \]

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.29077685383
110.930612408
7.25066276479

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

拿這兩份數據重新用線性回歸擬合模型。

lr2=LinearRegression()
lr2.fit(X2[training],y[training])
# 系數的均方誤差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])
# 系數的均方誤差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())

對於第二份共線性構造數據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
clf=LinearRegression()
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 = 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()

其中當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、1、2、3、8、9個變量都出現了波動，代表它們之間存在一定的共線性。觀察嶺跡，我們可以考慮剔除其中波動比較大的第1、8、9個變量。

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

參考文獻

[1]. variance inflation factor
[2]. 多重共線性的解決方法之——嶺回歸與LASSO
[3]. [ridge regression](https://en.wikipedia.org/wiki/Tikhonov_regularization)

---