Factor de inflación varianza en Python

Estoy tratando de calcular el factor de inflación de varianza (VIF) para cada columna en un conjunto de datos simple en python:

abcd 1 2 4 4 1 2 6 3 2 3 7 4 3 2 8 5 4 1 9 4 

Ya hice esto en R usando la función vif de la biblioteca usdm que da los siguientes resultados:

 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 

Sin embargo, cuando hago lo mismo en python usando la función statsmodel vif , mis resultados son:

 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 

Los resultados son muy diferentes, aunque las entradas son las mismas. En general, los resultados de la función VIF de statsmodel parecen estar equivocados, pero no estoy seguro de si esto se debe a la forma en que lo llamo o si es un problema con la función en sí.

Esperaba que alguien pudiera ayudarme a averiguar si estaba llamando incorrectamente a la función statsmodel o explicar las discrepancias en los resultados. Si es un problema con la función, ¿existen alternativas VIF en python?

Creo que la razón de esto se debe a una diferencia en el OLS de Python. OLS, que se utiliza en el cálculo del factor de inflación de la varianza de Python, no agrega una intercepción de forma predeterminada. Sin embargo, definitivamente quieres una intercepción allí.

Lo que querría hacer es agregar una columna más a su matriz, ck, llena de unas para representar una constante. Este será el término de intercepción de la ecuación. Una vez hecho esto, sus valores deben coincidir correctamente.

Editado: sustituye los ceros por unos.

Como lo mencionaron otros y en este post de Josef Perktold, el autor de la función, variance_inflation_factor espera la presencia de una constante en la matriz de variables explicativas. Se puede usar add_constant desde statsmodels para agregar la constante requerida al dataframe antes de pasar sus valores a la función.

 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 

Creo que también podría agregar la constante a la columna más a la derecha del dataframe usando 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 

El código fuente en sí es bastante conciso:

 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 

También es bastante simple modificar el código para devolver todos los VIF como una serie:

 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 

Para futuros usuarios de este hilo (como yo):

 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() 

Este código da

 array([22.95, 3. , 12.95, 3. ]) 

[EDITAR]

En respuesta a un comentario, traté de usar DataFrame tanto como sea posible (se requiere numpy para invertir una matriz).

 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) 

El código da

  abcd 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 

Los elementos diagonales dan VIF.

Ejemplo para los datos de Boston :

El VIF se calcula por regresión auxiliar, por lo que no depende del ajuste real.

Vea abajo:

 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() 

En caso de que no quieras lidiar con variance_inflation_factor y add_constant y solo quieras usar la fórmula. Por favor considere la siguiente función.

 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 

Escribí esta función basada en algunos otros mensajes que vi en Stack y CrossValidated. Muestra las características que están por encima del umbral y devuelve un nuevo dataframe con las características eliminadas.

 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