Problem występowania współliniowości zmiennych niezależnych w regresji logistycznej. Przykład użycia testu VIF.

PL202001081905

Jednym z najważniejszych warunków w modelu regresji logistycznej jest brak współliniowości (multicollinearity) zmiennych niezależnych. Oznacza to, że model powinien mieć niewielką współliniowość zmiennych opisujących lub jej brak.

Multicollinearity zmiennych niezależnych prowadzi do niewiarygodnych i niestabilnych oszacowań. Model nie jest stabilny, niewielkie zmiany w zmiennych niezależnych mogą diametralnie zmieniać wyniki klasyfikacji modelu.

Mimo że regresja logistyczna nie jest regresją lecz klasyfikacją, obowiązują tu podobne zasady. Ważna jest zmiana zmiennej zależnej dla każdej 1 zmiany jednostki w zmiennej niezależnej, gdy wszystkie pozostałe zmienne niezależne są stałe.

Chodzi o to, żeby zmienna jednej zmiennej nie wiązała się ze zmianą innych zmiennych niezależnych. Gdy zmienne niezależne są skorelowane, oznacza to, że zmiany w jednej zmiennej są powiązane ze zmianami w innej zmiennej. Im silniejsza korelacja, tym trudniej jest zmienić jedną zmienną bez zmiany drugiej. Modelowi trudno jest oszacować związek między każdą zmienną niezależną a zmienną zależną niezależnie, ponieważ wtedy gdy zmienne niezależne są skorelowane zmieniają się zgodnie.

Dowodem na istnienie współliniowości (multicollinearity) zmiennych niezależnych w modelu regresji logistycznej jest ich skorelowanie. Ta korelacja stanowi problem, ponieważ zmienne niezależne powinny być niezależne. Jeśli stopień korelacji między zmiennymi jest wystarczająco wysoki, może to powodować problemy w dopasowaniu modelu i interpretacji wyników.

Jak radzić sobie z problemem multicollinearity w modelu regresji logistycznej?

Multicolinearity można sprawdzić robiąc prostą macierz regresji dla zmiennych niezależnych.
Korfelacja zmiennych niezależnych to nie to samo co multicollinearity ponieważ multicollinearity może się pojawić w modelu, nawet gdy izolowane pary zmiennych nie są współliniowe (skorelowane), pojawia się interakcja 3 zmiennych itd.
Dlatego multicollinearity jest czasem bardziej trudne do wykrycia.

Variance Inflation Factor (VIF)

Variance Inflation Factor (VIF) to miara multicollinearity między zmiennymi predykcyjnymi w regresji wielokrotnej. Określa ilościowo nasilenie multicollinearity w zwykłej analizie regresji metodą najmniejszych kwadratów . Zapewnia wskaźnik, który mierzy, o ile wariancja (kwadrat odchylenia standardowego oszacowania ) szacowanego współczynnika regresji jest zwiększona z powodu kolinearności.

Pierwiastek kwadratowy współczynnika inflacji wariancji wskazuje, o ile większy błąd standardowy wzrasta w porównaniu z tym, czy zmienna ta ma korelację 0 z innymi zmiennymi predykcyjnymi w modelu.

*Przykład
Jeśli współczynnik inflacji wariancji zmiennej predykcyjnej wynosił 5,27 (√5,27 = 2,3), oznacza to, że błąd standardowy współczynnika tej zmiennej predykcyjnej jest 2,3 razy większy niż w przypadku, gdy zmienna predyktorowa ma 0 korelacji z innymi zmiennymi predykcyjnymi.*

(źródło:https://en.wikipedia.org/wiki/Variance_inflation_factor)

Kroki wdrażania VIF

  1. Uruchom regresję wieloraką.
  2. Oblicz współczynniki VIF.
  3. Sprawdź współczynniki dla każdej zmiennej predykcyjnej, jeśli VIF wynosi między 5-10, prawdopodobnie występuje multicollinearity i powinieneś rozważyć usunięcie tej zmiennej.

źródło: https://etav.github.io/python/vif_factor_python.html

Czyli przed rozpoczęciem pracy nad modelem regresji logistycznej należy utworzyć model regresji wielorakiej i przez VIF wyselekcjonować, wybrać dane do modelu regresji logistycznej.

Jak wspominają tutaj: https://www.ibm.com/support/pages/multicollinearity-diagnostics-logistic-regression-nomreg-or-plum

Procedury regresji dla zmiennych zależnych kategorycznie nie mają diagnostyki kolinearności. W tym celu można jednak użyć procedury regresji liniowej. Statystyka kolinearności w regresji dotyczy relacji między predyktorami, ignorując zmienną zależną. Możesz więc uruchomić REGRESSION z tą samą listą predyktorów i zmiennej zależnej, jakiej chcesz użyć w REGRESJI LOGISTYCZNEJ (na przykład) i zażądać diagnostyki kolinearności. Uruchom regresję logistyczną, aby uzyskać właściwe współczynniki, przewidywane prawdopodobieństwa itp. Po podjęciu niezbędnych decyzji (porzucenie predyktorów itp.) Wynikających z analizy kolinearności.

In [1]:
import pandas as pd
import numpy as np

df = pd.read_csv('c:/2/poliaxid.csv', index_col=0)
#del df['Unnamed: 0.1']
df.head(5)
Out[1]:
nr. factorA factorB citric catoda residual butanol caroton stable nodinol sulfur in nodinol density pH noracid lacapon quality class
0 0 4.933333 0.466667 0.000000 1.266667 0.050667 7.333333 22.666667 0.665200 2.340000 0.373333 6.266667 1
1 1 5.200000 0.586667 0.000000 1.733333 0.065333 16.666667 44.666667 0.664533 2.133333 0.453333 6.533333 1
2 2 5.200000 0.506667 0.026667 1.533333 0.061333 10.000000 36.000000 0.664667 2.173333 0.433333 6.533333 1
3 3 7.466667 0.186667 0.373333 1.266667 0.050000 11.333333 40.000000 0.665333 2.106667 0.386667 6.533333 2
4 4 4.933333 0.466667 0.000000 1.266667 0.050667 7.333333 22.666667 0.665200 2.340000 0.373333 6.266667 1

Z ciekawości sprawdzamy jaki jest poziom korelacji pomiędzy zmiennymi niezależnymi.

In [2]:
import seaborn as sns
import matplotlib.pyplot as plt
plt.figure(figsize=(10,10))


CORREL = df.corr()
sns.heatmap(CORREL, cmap="YlGnBu", annot=True, cbar=False)
plt.show
Out[2]:
<function matplotlib.pyplot.show(*args, **kw)>

Widać wyraźnie, że niektóre zmienne są ze sobą silnie skorelowane. Nie możemy ich ślepo eliminować z modelu. Następnym krokiem jest stworzenie prostego modelu regresji wielorakiej OLS.

Rutynowo sprawdzamy kompletność i format danych.

In [3]:
df.isnull().sum()
Out[3]:
nr.                  0
factorA              0
factorB              0
citric catoda        0
residual butanol     0
caroton              0
stable nodinol       0
sulfur in nodinol    0
density              0
pH                   0
noracid              0
lacapon              0
quality class        0
dtype: int64
In [4]:
df.dtypes
Out[4]:
nr.                    int64
factorA              float64
factorB              float64
citric catoda        float64
residual butanol     float64
caroton              float64
stable nodinol       float64
sulfur in nodinol    float64
density              float64
pH                   float64
noracid              float64
lacapon              float64
quality class          int64
dtype: object

Budujemy regresję liniową w bibliotece statmodel.

Tworzymy model regresji liniowej

In [5]:
df.columns
Out[5]:
Index(['nr.', 'factorA', 'factorB', 'citric catoda', 'residual butanol',
       'caroton', 'stable nodinol', 'sulfur in nodinol', 'density', 'pH',
       'noracid', 'lacapon', 'quality class'],
      dtype='object')
In [6]:
df.columns=['nr.', 'factorA', 'factorB', 'citric_catoda', 'residual_butanol',
       'caroton', 'stable_nodinol', 'sulfur_in_nodinol', 'density', 'pH',
       'noracid', 'lacapon', 'quality_class']
In [7]:
#import statsmodels.api as sm
#from statsmodels.formula.api import ols

VIF (Variation Inflation Factor)

Dzięki temu czynnikowi jesteśmy w stanie wskazać, która zmienna powinna zostać wyeliminowana z modelu.

In [8]:
df.columns
Out[8]:
Index(['nr.', 'factorA', 'factorB', 'citric_catoda', 'residual_butanol',
       'caroton', 'stable_nodinol', 'sulfur_in_nodinol', 'density', 'pH',
       'noracid', 'lacapon', 'quality_class'],
      dtype='object')
In [9]:
from patsy import dmatrices
from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.formula.api as smf

lm = smf.ols(formula = 'quality_class ~ factorA+factorB+citric_catoda+residual_butanol+caroton+stable_nodinol+density+sulfur_in_nodinol+pH+noracid+lacapon', data = df).fit()
y, X = dmatrices('quality_class ~ factorA+factorB+citric_catoda+residual_butanol+caroton+stable_nodinol+density+sulfur_in_nodinol+pH+noracid+lacapon', data = df, return_type = "dataframe")
vif = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
print(vif)
[1709157.8270185539, 7.777537721773227, 1.7885509640263009, 3.1329923025205733, 1.7030403470293642, 1.4834333245817521, 1.9734391365941484, 6.3474498618969335, 2.20192092692558, 3.323580122673945, 1.429921652882204, 3.034070351762669]

Ustawiamy macierz wyników.

In [10]:
X.columns
Out[10]:
Index(['Intercept', 'factorA', 'factorB', 'citric_catoda', 'residual_butanol',
       'caroton', 'stable_nodinol', 'density', 'sulfur_in_nodinol', 'pH',
       'noracid', 'lacapon'],
      dtype='object')
In [11]:
vif  =  np.round(vif, decimals=2) 
In [12]:
vif = list(map(float, vif))
name = list(X)
In [13]:
s1=pd.Series(name,name='name')
s2=pd.Series( vif,name='vif')

RFE_list = pd.concat([s1,s2], axis=1)

RFE_list
Out[13]:
name vif
0 Intercept 1709157.83
1 factorA 7.78
2 factorB 1.79
3 citric_catoda 3.13
4 residual_butanol 1.70
5 caroton 1.48
6 stable_nodinol 1.97
7 density 6.35
8 sulfur_in_nodinol 2.20
9 pH 3.32
10 noracid 1.43
11 lacapon 3.03

Interpretacja: wynik w postaci wektora reprezentuje zmienną w określonej kolejności jak w modelu. W zaleceniu VIF wskazano, że jeśli współczynnik przypisany do zmiennej jest większy niż 5, zmienna ta jest wysoce skorelowana z innymi zmiennymi i powinna zostać wyeliminowana z modelu.

Test wykazał, że zmienne: density oraz factorA powinny zostać usunięta z modelu.
jeszcze raz oglądamy macierz korelacji

Jeszcze raz tworzymy macierz korelacji vif i sprawdzamy zmienne po wyeliminowaniu: factorA i density.

In [14]:
lm = smf.ols(formula = 'quality_class ~ factorB+citric_catoda+residual_butanol+caroton+stable_nodinol+sulfur_in_nodinol+pH+noracid+lacapon', data = df).fit()
y, X = dmatrices('quality_class ~ factorB+citric_catoda+residual_butanol+caroton+stable_nodinol+sulfur_in_nodinol+pH+noracid+lacapon', data = df, return_type = "dataframe")
vif2 = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
print(vif2)
[809.2230489081707, 1.6530213053949512, 2.1806661580635818, 1.093561268844909, 1.3702675584858117, 1.9475860748759193, 2.0181057450324356, 1.6077595870352759, 1.3354874384986997, 1.2793132834168612]
In [15]:
X2 = 'Intercept', 'factorB', 'citric_catoda', 'residual_butanol','caroton', 'stable_nodinol', 'sulfur_in_nodinol', 'pH','noracid', 'lacapon'


vif  =  np.round(vif2, decimals=2) 
vif = list(map(float, vif2))
name = list(X2)

s1=pd.Series(name,name='name')
s2=pd.Series( vif2,name='vif')

RFE_list = pd.concat([s1,s2], axis=1)

RFE_list
Out[15]:
name vif
0 Intercept 809.223049
1 factorB 1.653021
2 citric_catoda 2.180666
3 residual_butanol 1.093561
4 caroton 1.370268
5 stable_nodinol 1.947586
6 sulfur_in_nodinol 2.018106
7 pH 1.607760
8 noracid 1.335487
9 lacapon 1.279313

Teraz wyraźnie zmniejszył się poziom vif, oznacza to że zmniejszył się stopień multicollinearity.

ciąg dalszy tutaj:

Przykład klasyfikacji wykonanej za pomocą regresji logistycznej. Eliminacja współliniowości zmiennych niezależnych za pomocą VIF