Varianz-Inflationsfaktor in Python
Ich versuche, den Varianz-Inflationsfaktor (VIF) für jede Spalte eines einfachen Datensatzes in Python zu berechnen:
a b c d
1 2 4 4
1 2 6 3
2 3 7 4
3 2 8 5
4 1 9 4
Ich habe dies in R bereits mit der vif-Funktion aus der usdm-Bibliothek durchgeführt , die die folgenden Ergebnisse liefert:
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
Wenn ich dasselbe in Python mit der vif-Funktion statsmodel mache, sind meine Ergebnisse jedoch:
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
Die Ergebnisse sind sehr unterschiedlich, obwohl die Eingaben gleich sind. Im Allgemeinen scheinen die Ergebnisse der statistischen VIF-Funktion falsch zu sein, aber ich bin nicht sicher, ob dies auf die Art und Weise zurückzuführen ist, wie ich sie anrufe, oder ob es ein Problem mit der Funktion selbst ist.
Ich hatte gehofft, jemand könnte mir helfen herauszufinden, ob ich die Statistikmodel-Funktion falsch aufgerufen habe, oder die Diskrepanzen in den Ergebnissen erklären. Gibt es bei der Funktion ein Problem mit der VIF-Alternative in Python?
Ich glaube, der Grund dafür liegt in einem Unterschied in Pythons OLS. OLS, das bei der Berechnung des Python-Varianz-Inflationsfaktors verwendet wird, fügt standardmäßig keinen Intercept hinzu. Sie wollen aber auf jeden Fall einen Abschnitt drin haben.
Sie möchten Ihrer Matrix eine weitere Spalte hinzufügen, ck, gefüllt mit einer, um eine Konstante darzustellen. Dies ist der Intercept-Term der Gleichung. Sobald dies geschehen ist, sollten Ihre Werte richtig übereinstimmen.
Editiert: Nullen durch Einsen ersetzt
Wie von anderen und in diesem Beitrag von Josef Perktold, dem Autor der Funktion, erwähnt, erwartet variance_inflation_factor das Vorhandensein einer Konstante in der Matrix der erklärenden Variablen. Sie können add_constant von statsmodels verwenden, um dem Datenrahmen die erforderliche Konstante hinzuzufügen, bevor Sie dessen Werte an die Funktion übergeben.
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
Ich glaube, Sie könnten die Konstante auch in die rechte rechte Spalte des Datenrahmens einfügen, indem Sie assign verwenden:
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
Der Quellcode selbst ist ziemlich knapp:
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
Es ist auch ziemlich einfach, den Code so zu ändern, dass alle VIFs als Serie zurückgegeben werden:
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
Für zukünftige Besucher dieses Threads (wie ich):
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()
Dieser Code gibt an
array([22.95, 3. , 12.95, 3. ])
[BEARBEITEN]
Als Antwort auf einen Kommentar habe ich versucht, so oft wie möglich DataFrame zu verwenden (numpy ist erforderlich, um eine Matrix zu invertieren).
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)
Der Code gibt an
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
Die diagonalen Elemente geben VIF an.
Beispiel für Boston Data :
VIFwird durch Hilfsregression berechnet und ist daher nicht von der tatsächlichen Anpassung abhängig.
Siehe unten:
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()
Ich habe diese Funktion basierend auf einigen anderen Beiträgen geschrieben, die ich auf Stack und CrossValidated gesehen habe. Es zeigt die Features, die über dem Schwellenwert liegen, und gibt ein neues Datenframe mit entfernten Features zurück.
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 X: the pandas dataframe
: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
Falls Sie sich nicht mit variance_inflation_factor und add_constant auseinandersetzen möchten und einfach die Formel verwenden möchten Bitte beachten Sie die folgende Funktion.
import pandas as pd
import seaborn as sns
import numpy as np
import statsmodels.formula.api as smf
# define function
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.2
'''
# initialize arrays
vif_array = np.array([])
tolerance_array = np.array([])
# 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_array = np.append(vif_array, vif).round(2)
# calculate tolerance
tolerance = 1-r_squared
tolerance_array = np.append(tolerance_array, tolerance)
# return VIF DataFrame
df_vif = pd.DataFrame({'VIF': vif_array, 'Tolerance': tolerance_array},
index=exogs)
return df_vif
[In]
df = sns.load_dataset('car_crashes')
exogs = ['alcohol', 'speeding', 'no_previous', 'not_distracted']
get_vif(exogs=exogs, data=df)
[Out]
VIF Tolerance
alcohol 3.44 0.291030
speeding 1.88 0.530690
no_previous 3.11 0.321132
not_distracted 2.67 0.374749