Python中的方差膨胀因子
Posted
技术标签:
【中文标题】Python中的方差膨胀因子【英文标题】：Variance Inflation Factor in Python 【发布时间】：2017-07-28 06:09:01 【问题描述】：我正在尝试计算 python 中一个简单数据集中每一列的方差膨胀因子 (VIF)：
a b c d
1 2 4 4
1 2 6 3
2 3 7 4
3 2 8 5
4 1 9 4
我已经在 R 中使用 usdm library 中的 vif 函数完成了此操作，结果如下：
a <- c(1, 1, 2, 3, 4)
b <- c(2, 2, 3, 2, 1)
c <- c(4, 6, 7, 8, 9)
d <- c(4, 3, 4, 5, 4)
df <- data.frame(a, b, c, d)
vif_df <- vif(df)
print(vif_df)
Variables VIF
a 22.95
b 3.00
c 12.95
d 3.00
但是，当我在 python 中使用statsmodel vif function 执行相同操作时，我的结果是：
a = [1, 1, 2, 3, 4]
b = [2, 2, 3, 2, 1]
c = [4, 6, 7, 8, 9]
d = [4, 3, 4, 5, 4]
ck = np.column_stack([a, b, c, d])
vif = [variance_inflation_factor(ck, i) for i in range(ck.shape[1])]
print(vif)
Variables VIF
a 47.136986301369774
b 28.931506849315081
c 80.31506849315096
d 40.438356164383549
即使输入相同，结果也大不相同。一般来说，statsmodel VIF 函数的结果似乎是错误的，但我不确定这是因为我调用它的方式还是函数本身的问题。
我希望有人可以帮助我弄清楚我是否错误地调用了 statsmodel 函数或解释结果中的差异。如果是函数的问题，那么 python 中是否有任何 VIF 替代方案？
【问题讨论】：
【参考方案1】：另一种解决方案。以下代码给出了与 R car 包完全相同的 VIF 结果。
def calc_reg_return_vif(X, y):
"""
Utility function to calculate the VIF. This section calculates the linear
regression inverse R squared.
Parameters
----------
X : DataFrame
Input data.
y : Series
Target.
Returns
-------
vif : float
Calculated VIF value.
"""
X = X.values
y = y.values
if X.shape[1] == 1:
print("Note, there is only one predictor here")
X = X.reshape(-1, 1)
reg = LinearRegression().fit(X, y)
vif = 1 / (1 - reg.score(X, y))
return vif
def calc_vif_from_scratch(df):
"""
Calculating VIF using function from scratch
Parameters
----------
df : DataFrame
without target variable.
Returns
-------
vif : DataFrame
giving the feature - VIF value pair.
"""
vif = pd.DataFrame()
vif_list = []
for feature in list(df.columns):
y = df[feature]
X = df.drop(feature, axis="columns")
vif_list.append(calc_reg_return_vif(X, y))
vif["feature"] = df.columns
vif["VIF"] = vif_list
return vif
我已经在 Titanic 数据集上对其进行了测试。您可以在此处获取完整示例：https://github.com/tulicsgabriel/Variance-Inflation-Factor-VIF-
【讨论】：
【参考方案2】：正如其他人以及函数作者 Josef Perktold 在this post 中提到的，variance_inflation_factor 预计解释变量矩阵中存在一个常数。可以使用 statsmodels 中的add_constant 将所需的常量添加到数据帧，然后再将其值传递给函数。
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.tools.tools import add_constant
df = pd.DataFrame(
'a': [1, 1, 2, 3, 4],
'b': [2, 2, 3, 2, 1],
'c': [4, 6, 7, 8, 9],
'd': [4, 3, 4, 5, 4]
)
X = add_constant(df)
>>> pd.Series([variance_inflation_factor(X.values, i)
for i in range(X.shape[1])],
index=X.columns)
const 136.875
a 22.950
b 3.000
c 12.950
d 3.000
dtype: float64
我相信您也可以使用 assign 将常量添加到数据框的最右侧列：
X = df.assign(const=1)
>>> pd.Series([variance_inflation_factor(X.values, i)
for i in range(X.shape[1])],
index=X.columns)
a 22.950
b 3.000
c 12.950
d 3.000
const 136.875
dtype: float64
源代码本身比较简洁：
def variance_inflation_factor(exog, exog_idx):
"""
exog : ndarray, (nobs, k_vars)
design matrix with all explanatory variables, as for example used in
regression
exog_idx : int
index of the exogenous variable in the columns of exog
"""
k_vars = exog.shape[1]
x_i = exog[:, exog_idx]
mask = np.arange(k_vars) != exog_idx
x_noti = exog[:, mask]
r_squared_i = OLS(x_i, x_noti).fit().rsquared
vif = 1. / (1. - r_squared_i)
return vif
修改代码以将所有 VIF 作为一个系列返回也相当简单：
from statsmodels.regression.linear_model import OLS
from statsmodels.tools.tools import add_constant
def variance_inflation_factors(exog_df):
'''
Parameters
----------
exog_df : dataframe, (nobs, k_vars)
design matrix with all explanatory variables, as for example used in
regression.
Returns
-------
vif : Series
variance inflation factors
'''
exog_df = add_constant(exog_df)
vifs = pd.Series(
[1 / (1. - OLS(exog_df[col].values,
exog_df.loc[:, exog_df.columns != col].values).fit().rsquared)
for col in exog_df],
index=exog_df.columns,
name='VIF'
)
return vifs
>>> variance_inflation_factors(df)
const 136.875
a 22.950
b 3.000
c 12.950
Name: VIF, dtype: float64
根据@T_T的解决方案，也可以简单地做以下事情：
vifs = pd.Series(np.linalg.inv(df.corr().to_numpy()).diagonal(),
index=df.columns,
name='VIF')
【讨论】：
我认为在缺少值的情况下添加X = add_constant(df.dropna()) 是安全的。
感谢您的解决方案。我非常困惑为什么我的模型自变量的 VIF 如此之高，这就是我最终发表这篇文章的原因。尽管我很讨厌这样做，但我几乎很想在 R 中完成我的分析。
对于一些数据，不可能创建导致LinAlgError("Singular matrix")numpy.linalg.LinAlgError: Singular matrix 的逆矩阵。在这种情况下，将 inv() 替换为 pinv()。 pinv() 计算矩阵的 (Moore-Penrose) 伪逆。 pd.Series(np.linalg.pinv(X.corr().to_numpy()).diagonal(), index=X.columns, name='VIF') Out[13]: a 22.95 b 3.00 c 12.95 d 3.00【参考方案3】：
我认为这是由于 Python 的 OLS 不同所致。在python方差膨胀因子计算中使用的OLS默认不添加截距。但是，您肯定想要在那里拦截。
您想要做的是在矩阵 ck 中再添加一列，用一列填充以表示一个常数。这将是方程的截距项。完成此操作后，您的值应该正确匹配。
已编辑：用一替换零
【讨论】：
从所有变量中减去平均值是相似的。 错字：常量列应该用 1（而不是 0）填充。 对我的错字很好。使用修复程序编辑了我的原始帖子。 这是有道理的。添加一列 1 就可以了。谢谢！【参考方案4】：虽然已经很晚了，但我正在从给定的答案中添加一些修改。如果我们使用@Chef1075 解决方案，为了在消除多重共线性后获得最佳集合，那么我们将丢失相关的变量。我们只需要删除其中一个。为此，我使用@steve answer提供了以下解决方案：
import pandas as pd
from sklearn.linear_model import LinearRegression
def sklearn_vif(exogs, data):
'''
This function calculates variance inflation function in sklearn way.
It is a comparatively faster process.
'''
# initialize dictionaries
vif_dict, tolerance_dict = ,
# form input data for each exogenous variable
for exog in exogs:
not_exog = [i for i in exogs if i != exog]
X, y = data[not_exog], data[exog]
# extract r-squared from the fit
r_squared = LinearRegression().fit(X, y).score(X, y)
# calculate VIF
vif = 1/(1 - r_squared)
vif_dict[exog] = vif
# calculate tolerance
tolerance = 1 - r_squared
tolerance_dict[exog] = tolerance
# return VIF DataFrame
df_vif = pd.DataFrame('VIF': vif_dict, 'Tolerance': tolerance_dict)
return df_vif
df = pd.DataFrame(
'a': [1, 1, 2, 3, 4,1],
'b': [2, 2, 3, 2, 1,3],
'c': [4, 6, 7, 8, 9,5],
'd': [4, 3, 4, 5, 4,6],
'e': [8,8,14,15,17,20]
)
df_vif= sklearn_vif(exogs=df.columns, data=df).sort_values(by='VIF',ascending=False)
while (df_vif.VIF>5).any() ==True:
red_df_vif= df_vif.drop(df_vif.index[0])
df= df[red_df_vif.index]
df_vif=sklearn_vif(exogs=df.columns,data=df).sort_values(by='VIF',ascending=False)
print(df)
d c b
0 4 4 2
1 3 6 2
2 4 7 3
3 5 8 2
4 4 9 1
5 6 5 3
【讨论】：
那么，在这种情况下，d、c 和b 列不会导致多重共线性，对吧？
@AlvaroMartinez。对
@MdAsrafulKabir 我能问一下你为什么要这样做red_df_vif= df_vif.drop(df_vif.index[0])吗？因此，您计算 VIF，将它们从高到低排序；如果最高值大于 5 则将其删除并重新计算整个过程？
@MdAsrafulKabir 我能问一下你为什么要这样做red_df_vif= df_vif.drop(df_vif.index[0])吗？因此，您计算 VIF，将它们从高到低排序；如果最高值大于 5 则将其删除并重新计算整个过程？【参考方案5】：
我根据我在 Stack 和 CrossValidated 上看到的其他一些帖子编写了这个函数。它显示超过阈值的特征并返回一个新的数据框，其中删除了特征。
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.tools.tools import add_constant
def calculate_vif_(df, thresh=5):
'''
Calculates VIF each feature in a pandas dataframe
A constant must be added to variance_inflation_factor or the results will be incorrect
:param df: the pandas dataframe containing only the predictor features, not the response variable
:param thresh: the max VIF value before the feature is removed from the dataframe
:return: dataframe with features removed
'''
const = add_constant(df)
cols = const.columns
variables = np.arange(const.shape[1])
vif_df = pd.Series([variance_inflation_factor(const.values, i)
for i in range(const.shape[1])],
index=const.columns).to_frame()
vif_df = vif_df.sort_values(by=0, ascending=False).rename(columns=0: 'VIF')
vif_df = vif_df.drop('const')
vif_df = vif_df[vif_df['VIF'] > thresh]
print 'Features above VIF threshold:\n'
print vif_df[vif_df['VIF'] > thresh]
col_to_drop = list(vif_df.index)
for i in col_to_drop:
print 'Dropping: '.format(i)
df = df.drop(columns=i)
return df
【讨论】：
删除所有 VIF 高于阈值的变量是不正确的。正确的做法是去掉 VIF 最高的变量，然后对剩余变量重新计算 VIF，重复此步骤，直到没有剩余变量的 VIF 大于 thresh。例如，假设 x3=x2+x1，如果我们简单地删除所有具有高 VIF 的变量，则 x1/x2/x3 将被删除并且没有保留，我们可能会丢失一个重要的变量。 是的，同意环发。 @chef 和其他人 - 如果您只是从初始运行中删除高于 VIF 阈值的所有列，您将停止更多的变量。正如环发所说，这需要迭代完成。【参考方案6】：这里使用数据框python的代码：
创建数据
import numpy as npimport scipy as sp
a = [1, 1, 2, 3, 4]b = [2, 2, 3, 2, 1]c = [4, 6, 7, 8, 9]d = [4, 3, 4, 5, 4]
创建数据框
import pandas as pddata = pd.DataFrame()data["a"] = adata["b"] = bdata["c"] = cdata["d"] = d
计算 VIF
cc = np.corrcoef(data, rowvar=False)VIF = np.linalg.inv(cc)VIF.diagonal()
结果
array([22.95, 3. , 12.95, 3. ])
【讨论】：
【参考方案7】：如果您不想处理 variance_inflation_factor 和 add_constant。请考虑以下两个函数。
1.在 statasmodels 中使用公式：
import pandas as pd
import statsmodels.formula.api as smf
def get_vif(exogs, data):
'''Return VIF (variance inflation factor) DataFrame
Args:
exogs (list): list of exogenous/independent variables
data (DataFrame): the df storing all variables
Returns:
VIF and Tolerance DataFrame for each exogenous variable
Notes:
Assume we have a list of exogenous variable [X1, X2, X3, X4].
To calculate the VIF and Tolerance for each variable, we regress
each of them against other exogenous variables. For instance, the
regression model for X3 is defined as:
X3 ~ X1 + X2 + X4
And then we extract the R-squared from the model to calculate:
VIF = 1 / (1 - R-squared)
Tolerance = 1 - R-squared
The cutoff to detect multicollinearity:
VIF > 10 or Tolerance < 0.1
'''
# initialize dictionaries
vif_dict, tolerance_dict = ,
# create formula for each exogenous variable
for exog in exogs:
not_exog = [i for i in exogs if i != exog]
formula = f"exog ~ ' + '.join(not_exog)"
# extract r-squared from the fit
r_squared = smf.ols(formula, data=data).fit().rsquared
# calculate VIF
vif = 1/(1 - r_squared)
vif_dict[exog] = vif
# calculate tolerance
tolerance = 1 - r_squared
tolerance_dict[exog] = tolerance
# return VIF DataFrame
df_vif = pd.DataFrame('VIF': vif_dict, 'Tolerance': tolerance_dict)
return df_vif
2。在 sklearn 中使用LinearRegression：
# import warnings
# warnings.simplefilter(action='ignore', category=FutureWarning)
import pandas as pd
from sklearn.linear_model import LinearRegression
def sklearn_vif(exogs, data):
# initialize dictionaries
vif_dict, tolerance_dict = ,
# form input data for each exogenous variable
for exog in exogs:
not_exog = [i for i in exogs if i != exog]
X, y = data[not_exog], data[exog]
# extract r-squared from the fit
r_squared = LinearRegression().fit(X, y).score(X, y)
# calculate VIF
vif = 1/(1 - r_squared)
vif_dict[exog] = vif
# calculate tolerance
tolerance = 1 - r_squared
tolerance_dict[exog] = tolerance
# return VIF DataFrame
df_vif = pd.DataFrame('VIF': vif_dict, 'Tolerance': tolerance_dict)
return df_vif
示例：
import seaborn as sns
df = sns.load_dataset('car_crashes')
exogs = ['alcohol', 'speeding', 'no_previous', 'not_distracted']
[In] %%timeit -n 100
get_vif(exogs=exogs, data=df)
[Out]
VIF Tolerance
alcohol 3.436072 0.291030
no_previous 3.113984 0.321132
not_distracted 2.668456 0.374749
speeding 1.884340 0.530690
69.6 ms ± 8.96 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)
[In] %%timeit -n 100
sklearn_vif(exogs=exogs, data=df)
[Out]
VIF Tolerance
alcohol 3.436072 0.291030
no_previous 3.113984 0.321132
not_distracted 2.668456 0.374749
speeding 1.884340 0.530690
15.7 ms ± 1.4 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)
【讨论】：
检测多重共线性的截止值：VIF > 10 或 Tolerance @SandraGuerrero 确实是一个错字。 感谢您为解释这一点所做的努力。非常感谢！ 非常感谢你【参考方案8】：对于这个帖子的未来来者（比如我）：
import numpy as np
import scipy as sp
a = [1, 1, 2, 3, 4]
b = [2, 2, 3, 2, 1]
c = [4, 6, 7, 8, 9]
d = [4, 3, 4, 5, 4]
ck = np.column_stack([a, b, c, d])
cc = sp.corrcoef(ck, rowvar=False)
VIF = np.linalg.inv(cc)
VIF.diagonal()
这段代码给出了
array([22.95, 3. , 12.95, 3. ])
[编辑]
在回应评论时，我尝试尽可能多地使用DataFrame（需要numpy 来反转矩阵）。
import pandas as pd
import numpy as np
a = [1, 1, 2, 3, 4]
b = [2, 2, 3, 2, 1]
c = [4, 6, 7, 8, 9]
d = [4, 3, 4, 5, 4]
df = pd.DataFrame('a':a,'b':b,'c':c,'d':d)
df_cor = df.corr()
pd.DataFrame(np.linalg.inv(df.corr().values), index = df_cor.index, columns=df_cor.columns)
代码给出
a b c d
a 22.950000 6.453681 -16.301917 -6.453681
b 6.453681 3.000000 -4.080441 -2.000000
c -16.301917 -4.080441 12.950000 4.080441
d -6.453681 -2.000000 4.080441 3.000000
对角线元素给出 VIF。
【讨论】：
能否请您为数据框输入而不是 numpy 数组添加一个解决方案？ 看起来不错。将 VIF 作为一个系列获取：vifs = pd.Series(np.linalg.inv(df.corr().values).diagonal(), index=df_cor.index)
vif是逆corr矩阵的对角元素？编辑：是的，检查链接：documentation.statsoft.com/STATISTICAHelp.aspx?path=glossary/…【参考方案9】：
波士顿数据示例：
VIF 是通过辅助回归计算的，因此不依赖于实际拟合。
见下文：
from patsy import dmatrices
from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.api as sm
# Break into left and right hand side; y and X
y, X = dmatrices(formula="medv ~ crim + zn + nox + ptratio + black + rm ", data=boston, return_type="dataframe")
# For each Xi, calculate VIF
vif = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
# Fit X to y
result = sm.OLS(y, X).fit()
【讨论】：
以上是关于Python中的方差膨胀因子的主要内容，如果未能解决你的问题，请参考以下文章
R计算方差膨胀因子（VIF，Variance Inflation Factor）计算并解读VID与共线性（Multicollinearity）的关系实战