多重共线性py的相关实现

多重共线性

理论部分：多重共线性问题

以线性回归为例

#导入相关库
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score

#数据准备
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])

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)

# 系数的均方误差MSE
print(((lr1.coef_-coef0)**2).sum()/8)
# 测试集准确率（R2）
print(lr1.score(X1[~training],y[~training]))
# 平均测试集准确率
print(cross_val_score(lr1,X1,y,cv=5).mean())
'''
0.30205546186487053
0.9320526019314885
0.9546828999707184
'''

0.30205546186487053
0.9320526019314885
0.9546828999707184

#构建两份对比数据X2
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)])

利用新构建的数据再次拟合线性回归模型

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_val_score(lr2,X2,y,cv=5).mean())
'''
1.8934092769989603
0.929592629455834
0.9519696022183932
'''

1.8934092769989603
0.929592629455834
0.9519696022183932

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_val_score(lr3,X3,y,cv=5).mean())
'''
0.2677097747038732
0.9350851381708499
0.9547601223319211

'''

0.2677097747038732
0.9350851381708499
0.9547601223319211

可以发现，构建的第二份数据的均方误差明显上升，而第三份数据的均方误差以及测试集准确率与第一份数据无过度差异，表明，共线性并非会一定对数据集的表现产生影响。

#查看一下每份份数据的系数
print(lr1.coef_)
print(lr2.coef_)
print(lr3.coef_)

[ 6.05684499  5.52349054  6.42462281  8.40395492  8.96814584  9.93046512
 10.64546334 12.66833867]
[ 5.9317617   3.62119675  4.59353762  6.52605017  8.97919766  9.96109222
 10.59202948 12.69886599  0.07089803  1.80734528]
[ 5.94537999  5.47303855  6.53890852  8.30882708  9.01350639 10.14582669
 10.69893617 12.74174099 -1.18705515 -0.28565136]

共线性问题的检验

VIF的计算

import matplotlib.pyplot as plt

vif1=np.zeros((10,1))
for i in range(10):
    tmp=[k for k in range(10) if k!=i]
    lr1.fit(X2[:,tmp],X2[:,i])
    vifi=1/(1-lr1.score(X2[:,tmp],X2[:,i])) #VIF的公式
    vif2[i]=vifi

vif2=np.zeros((10,1))
for i in range(10):
    tmp=[k for k in range(10) if k!=i]
    lr2.fit(X2[:,tmp],X2[:,i])
    vifi=1/(1-lr2.score(X2[:,tmp],X2[:,i])) #VIF的公式
    vif2[i]=vifi

vif3=np.zeros((10,1))
for i in range(10):
    tmp=[k for k in range(10) if k!=i]
    lr3.fit(X3[:,tmp],X3[:,i])
    vifi=1/(1-lr3.score(X3[:,tmp],X3[:,i]))
    vif3[i]=vifi  
plt.figure()
ax = plt.gca()
ax.plot(vif1)
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明显过高，其中，第9个特征的VIF最高。
对比之下，发现第三份数据的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()

在岭回归中，变量出现波动表明，它们之间存在一定的共线性，由上图可以看到第1，8，9的波动较大，考虑剔除。

#alpha取0.1时，岭回归估计的系数如下：
print(coefs[0])

[ 4.09491447  0.7971758   3.63635471  5.28177248  8.95731907 10.12079948
 10.24244483 12.24172518  2.08787846  2.85162068]

Lasso回归情况类似，可以用sklearn中的linear_model.Lasso来学习；

？？对于逻辑回归任务，sklearn函数内部提供了L1或L2正则化方案，通过它们也可以去检测共线性问题。

总结

1、共线性问题的存在不一定会对模型的效果产生影响，是否要对共线性进行处理，视情况而定。
2、检验共线性问题建议用VIF，客观且清晰。

参考

讲讲共线性问题