多元线性回归是我们在数据分析中经常用到的一个方法,很多人在遇到多维数据时基本上无脑使用该方法,而在用多元线性回归之后所得到的结果又并不总是完美的,其问题实际上并不出在方法上,而是出在数据上。当数据涉及的维度过多时,我们就很难保证维度之间互不相关,而这些维度又都对结果产生一定影响,当一组维度或者变量之间有较强的相关性时,就认为是一种违背多元线性回归模型基本假设的情形。今天我们就讲解一下如何用VIF方法消除多维数据中多重共线性的问题。
首先介绍一下多重共线性。
多元回归模型有一个基本假设,就是要求设计矩阵X
的秩rank(X)=p+1
,其中p
是维度数,即要求X中的列向量之间线性无关。如果存在不全为零的p+1
个数c0、c1、c2、...、cp
,使得c0 + c1xi1 + c2xi2 + ... + cpxip = 0
,i=1, 2, 3, ..., n
,则自变量x1、x2、...、xp
之间存在多重共线性(multi-collinearity),因为实际问题中,完全多重共线性不太常见,所以上式中的等号经常用约等号。多重共线性到底会带来什么问题呢,笔者下面就用一个实际例子来告诉大家。
这个例子来自1994年统计摘要,是一个中国民航客运量的回归模型,统计了1978至1993年的各年数据。该模型以民航客运量为因变量y
,以国民收入、消费额、民航航线里程、来华旅游入境人数作为影响客运量的主要因素,其中y
的单位是万人,x1
表示国民收入(亿元),x2
表示消费额(亿元),x3
表示铁路客运量(万人),x4
表示民航航线里程(万公里),x5
表示来华旅游入境人数(万人)。该数据集如图1所示,一共16行、7列(包括前面的year
,虽然这一列用不到)。
图1. 数据集截图
我们用该数据集来做一个多元线性回归模型,主要使用statsmodels,代码如下。首先是导入各种库。
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
接下来是读取数据,并生成自变量X
和因变量y
。
file = r'C:\Users\Desktop\data.xlsx'
data = pd.read_excel(file) #读取数据
y = data['y'] #因变量数据
X = data.loc[:, 'x1':] #自变量数据
然后是生成多元回归模型,并输出结果,结果如图2所示。
X = sm.add_constant(X) #加上一列常数1,这是回归模型中的常数项
reg = sm.OLS(y, X) #生成回归模型
model = reg.fit() #拟合数据
model.summary() #生成结果
图2. 原数据的多元回归模型结果
图2中的参数较多,如果大家对这个结果不太明白,可以参考笔者之前给公众号写的文章《详解用statsmodels进行回归分析》。从图2中可以得出,我们的模型的回归方程为y = 450.9 + 0.354x1 - 0.561x2 - 0.0073x3 + 21.578x4 + 0.435x5
,看到这里,估计很多人就看出一些问题了。这个回归模型从表面上看没有什么太大的问题,但是仔细分析一下其实际意义,就能看出一些问题,x2
是消费额,从经济学角度分析,消费额与民航客运量之间的关系应该是正相关的,也就是x2
前面的系数应该是正的,但这里却是负值,问题出在哪?这就是变量之间的多重共线性造成的。
多重共线性的影响就在于此,我们的模型结果中每一个参数都能通过检验,而且模型整体的线性显著性也很好(比如该例中R-squared
值为0.998,效果非常好),但其部分参数的实际意义却和我们的常识是相违背的,而这种情况我们往往很难察觉,很多人看到自己的模型在数学角度上没有任何问题,就直接拿去用了,结果总是得到错误的结论。
那么如何来诊断多重共线性呢?笔者今天就介绍一下VIF
方法。VIF
全称是Variance Inflation Factor
,即方差扩大因子,我们对自变量X
作中心标准化,则X
变为Xs
,然后可以得到Xs’ Xs = (rij)
,这个就是自变量的相关阵。如图3所示,式(1)中C的主对角线元素VIFj=cjj
,就是自变量xj
的方差扩大因子,式(2)中的Rj^2
是自变量xj
对其余p-1
个自变量的复决定系数,式(2)也可以作为方差扩大因子VIFj
的定义,可知VIFj
是大于等于1的。
图3. VIF方法部分公式
Rj^2
度量了自变量xj
与其余p-1
个自变量的线性相关程度,这种相关程度越强,说明自变量之间的多重共线性就越严重,Rj^2
越接近1,VIFj
就越大。反之,xj
与其余p-1
个自变量之间的线性相关程度越弱,多重共线性就越弱,Rj^2
就越接近于0,VIFj
就约接近于1。由此可见,VIFj
的大小反映了自变量之间是否存在多重共线性,可由它来度量多重共线性的严重程度,那么VIFj
多大才算是有严重的多重共线性呢?根据统计学中的使用经验,当VIFj
大于等于10的时候,就说明自变量xj
与其余自变量之间存在严重的多重共线性,且这种多重共线性会过度地影响最小二乘估计值。
在了解了VIF
的概念之后,我们就用代码来展示一下如何诊断并消除多重共线性。这里笔者依然使用前面的数据,但加入了VIF
检测,同时给出消除多重共线性后的结果,全部代码如下。
file = r'C:\Users\Desktop\data.xlsx'
data = pd.read_excel(file)
y = data['y']
X = data.loc[:, 'x1':]
X = sm.add_constant(X)
def process(data, col):
data = data.loc[:, col] #读取对应列标数据
vif = [variance_inflation_factor(data.values, i) for i in range(data.shape[1])][1:]
if max(vif) >= 10:
index = np.argmax(vif)+1 #得到最大值的标号
del col[index] #删除vif值最大的一项
return process(data, col) #递归过程
else:
vif = [variance_inflation_factor(data.values, i) for i in range(data.shape[1])][1:]
return col, vif
cols = ['const', 'x1', 'x2', 'x3', 'x4', 'x5']
cols, vif = process(X, cols)
reg = sm.OLS(y, X[cols])
model = reg.fit()
model.summary()
这里我们从process
这个函数开始讲起,process
需要两个参数,一个是data
,就是要输入的数据,另一个是col
,就是数据的columns
(即数据的列标题),我们这里默认使用的数据集是pandas.DataFrame
格式的,所以数据都是有columns
的。在process
函数中,data = data.loc[:, col]
就是读取只含有col
列标的那些数据, vif = [variance_inflation_factor(data.values, i) for i in range(data.shape[1])][1:]
这行代码就是计算vif
的过程,variance_inflation_factor
函数需要输入两个参数,分别是数据和每列数据的标号,这个标号也是从0开始的。而最终我们取得的vif
结果是去掉第一项的,因为第一项对应数据集中const
那一列,这一列因为都是1,所以在vif
结果中要去掉,但在计算时要保留。而得到vif
之后,我们要找出vif
中数据最大的一项,判断其是否大于等于10。如果是,就找到其对应的标号,利用np.argmax
即可,然后删除col
中这一项,再把所得的结果带入到process
函数中,形成递归;如果不是,则直接返回col
和vif
这两个结果。
最终我们得到的cols
是['const', 'x3', 'x4', 'x5']
,const
就是前面X = sm.add_constant(X)
中加入的常数项一列,这个const
列标是自动添加的,我们在这里仍沿用这个叫法,这列数据在VIF
方法中只参与计算,但其值不用于比较大小。我们可以看到这里的结果去掉了x1
和x2
这两列数据,消除多重共线性最好的方式就是把那些造成多重共线性的维度(自变量)直接去掉,vif
是[1.9836946236748652, 6.6499090855830225, 8.513876170172715]
,vif
中所有数值都在10以内,说明目前已经消除了多重共线性。
然后用剩下的这些数据进行建模,得到多元回归模型,其结果如图4所示。该模型为y = 591.9 - 0.0104x3 + 26.4358x4 + 0.3174x5
,该模型无论从数学上还是经济意义上,都能合理有效地进行解释。
图4. 用VIF法处理后的模型结果
判断数据是否具有多重共线性实际上有多种方法,比如特征根判定法、直接判断法等,本文主要讲解如何用VIF
法来诊断多重共线性,有兴趣的读者也可以把此方法和其他方法进行一下对比学习。
作者简介:Mort,数据分析爱好者,擅长数据可视化,比较关注机器学习领域,希望能和业内朋友多学习交流。