Rozważania teoretyczne na temat znaczenia (szkodliwości) wspołiniowości zmiennych niezależnych (multicollinearity) zostały wyjaśnione w poprzedniej części.

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

Tym razem będziemy analizowali zagadnienie zanieczyszczenia miasta.
Źródło danych można znaleźć tutaj.

import pandas as pd
df = pd.read_csv('c:/TF/AirQ_filled.csv')
df.head(3)

Na początku sprawdzimy kompletność danych oraz ich format.

del df['Unnamed: 0']
df.isnull().sum()

Date             0
Time             0
CO(GT)           0
PT08.S1(CO)      0
C6H6(GT)         0
PT08.S2(NMHC)    0
NOx(GT)          0
PT08.S3(NOx)     0
NO2(GT)          0
PT08.S4(NO2)     0
PT08.S5(O3)      0
T                0
RH               0
AH               0
dtype: int64

df.dtypes

Date              object
Time              object
CO(GT)           float64
PT08.S1(CO)      float64
C6H6(GT)         float64
PT08.S2(NMHC)    float64
NOx(GT)          float64
PT08.S3(NOx)     float64
NO2(GT)          float64
PT08.S4(NO2)     float64
PT08.S5(O3)      float64
T                float64
RH               float64
AH               float64
dtype: object

Dane są kompletne i mają właściwy format do prowadzenia dalszej analizy.

Zadanie

Naszym zadaniem jest podział na trzy klasy zanieczyszczenia oznaczonego jako PT08.S5(03). Następnie musimy zbudować model klasyfikacji logistycznej oparty na pozostałych zmiennych w postaci poziomów zanieczyszczeń.

Korelacja zmiennych niezależnych ze zmienną zależną.

CORREL = df.corr()

import seaborn as sns
sns.heatmap(CORREL, annot=True, cbar=False, cmap="coolwarm")

Powyższa macierz wskazuje, że istnieje współliniowość, która negatywnie wpłynie na poziom klasyfikacji. Nie zważając na to kontynujemy realizację naszego zadania. Poniżej przeprowadziłem analizę korelacji zmiennej zależnej: 'PT08.S5(O3)’ ze zmiennymi niezależnymi. W wiekszości występuje wysoki poziom korelacji co jest zjawiskiem pozytywnym.

CORREL2 = df.corr().sort_values('PT08.S5(O3)')
CORREL2['PT08.S5(O3)'].plot(kind='bar')

C:ProgramDataAnaconda3libsite-packagesipykernel_launcher.py:1: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  """Entry point for launching an IPython kernel.

dtype('O')

CO_GT                   float64
PT08_S1_CO              float64
C6H6_GT                 float64
PT08_S2_NMHC            float64
NOx_GT                  float64
PT08_S3_NOx             float64
NO2_GT                  float64
PT08_S4_NO2             float64
T                       float64
RH                      float64
AH                      float64
Categores_PT08_S5_O3     object
dtype: object

Zbiór X treningowy:  (6269, 11)
Zbiór X testowy:     (3088, 11)
Zbiór y treningowy:  (6269,)
Zbiór y testowy:     (3088,)

C:ProgramDataAnaconda3libsite-packagessklearnlinear_modellogistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)

GridSearchCV(cv=2, error_score='raise-deprecating',
             estimator=LogisticRegression(C=1.0, class_weight=None, dual=False,
                                          fit_intercept=True,
                                          intercept_scaling=1, l1_ratio=None,
                                          max_iter=100, multi_class='warn',
                                          n_jobs=None, penalty='l2',
                                          random_state=None, solver='warn',
                                          tol=0.0001, verbose=0,
                                          warm_start=True),
             iid='warn', n_jobs=5,
             param_grid={'C': array([1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01, 1.e+02])},
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring='roc_auc', verbose=0)

array([[1414,  131],
       [ 176, 1367]], dtype=int64)

                  precision    recall  f1-score   support

(220.999, 964.0]       0.89      0.92      0.90      1545
 (964.0, 2523.0]       0.91      0.89      0.90      1543

        accuracy                           0.90      3088
       macro avg       0.90      0.90      0.90      3088
    weighted avg       0.90      0.90      0.90      3088

[648.230997747899, 11.43831319456862, 7.939972023928728, 46.44356640618427, 66.47134191824114, 5.220873142678388, 5.43267527445646, 3.714154563031008, 14.52279194616348, 14.359778937836424, 8.176237173633366, 10.619321015188527]

df['Categores_PT08.S5(O3)'] = pd.qcut(df['PT08.S5(O3)'],2)
df['Categores_PT08.S5(O3)'].value_counts().to_frame()

Badana kategoria PT08.S5(O3) została podzielona na równe przedziały. Mniejsza o to czy ten podział jest uzasadniony chemicznie. Najważniejsze że teraz możemy utworzyć model klasyfikacji oparty na regresji logistycznej.
Równy podział zbiorów oddala nas od negatywnego zjawiska niezbilansowania zbiorów i koniczności tworzenia oversampling. Dzisiaj nie to jest najważniejsze lecz przećwiczenie radzenia sobie ze zjawiskiej multicollinearity.

Zbudjemy model klasyfikacji tak jakby nie dotyczył nas problem multicollinearity.

Wskazujemy zmienne niezależne:¶

KOT = df[['CO(GT)', 'PT08.S1(CO)', 'C6H6(GT)', 'PT08.S2(NMHC)',
       'NOx(GT)', 'PT08.S3(NOx)', 'NO2(GT)', 'PT08.S4(NO2)', 
       'T', 'RH', 'AH','Categores_PT08.S5(O3)']]
KOT.columns = ['CO_GT', 'PT08_S1_CO', 'C6H6_GT', 'PT08_S2_NMHC',
       'NOx_GT', 'PT08_S3_NOx', 'NO2_GT', 'PT08_S4_NO2', 
       'T', 'RH', 'AH','Categores_PT08_S5_O3']

Przekształcamy zmienną zależną: Categores_PT08.S5(O3) na zmienną ciągłą.

KOT['Categores_PT08_S5_O3'] = KOT['Categores_PT08_S5_O3'].astype(str)
KOT['Categores_PT08_S5_O3'].dtypes

C:ProgramDataAnaconda3libsite-packagesipykernel_launcher.py:1: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  """Entry point for launching an IPython kernel.

dtype('O')

CO_GT                   float64
PT08_S1_CO              float64
C6H6_GT                 float64
PT08_S2_NMHC            float64
NOx_GT                  float64
PT08_S3_NOx             float64
NO2_GT                  float64
PT08_S4_NO2             float64
T                       float64
RH                      float64
AH                      float64
Categores_PT08_S5_O3     object
dtype: object

Zbiór X treningowy:  (6269, 11)
Zbiór X testowy:     (3088, 11)
Zbiór y treningowy:  (6269,)
Zbiór y testowy:     (3088,)

C:ProgramDataAnaconda3libsite-packagessklearnlinear_modellogistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)

GridSearchCV(cv=2, error_score='raise-deprecating',
             estimator=LogisticRegression(C=1.0, class_weight=None, dual=False,
                                          fit_intercept=True,
                                          intercept_scaling=1, l1_ratio=None,
                                          max_iter=100, multi_class='warn',
                                          n_jobs=None, penalty='l2',
                                          random_state=None, solver='warn',
                                          tol=0.0001, verbose=0,
                                          warm_start=True),
             iid='warn', n_jobs=5,
             param_grid={'C': array([1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01, 1.e+02])},
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring='roc_auc', verbose=0)

array([[1414,  131],
       [ 176, 1367]], dtype=int64)

                  precision    recall  f1-score   support

(220.999, 964.0]       0.89      0.92      0.90      1545
 (964.0, 2523.0]       0.91      0.89      0.90      1543

        accuracy                           0.90      3088
       macro avg       0.90      0.90      0.90      3088
    weighted avg       0.90      0.90      0.90      3088

[648.230997747899, 11.43831319456862, 7.939972023928728, 46.44356640618427, 66.47134191824114, 5.220873142678388, 5.43267527445646, 3.714154563031008, 14.52279194616348, 14.359778937836424, 8.176237173633366, 10.619321015188527]

[200.1579775895654, 5.895905610386816, 5.807909737315139, 3.1029777937073924, 2.7926706710126528, 2.4675790616762594]

from sklearn.model_selection import train_test_split 


y = KOT['Categores_PT08_S5_O3']
X = KOT.drop('Categores_PT08_S5_O3', axis=1)


Xtrain, Xtest, ytrain, ytest = train_test_split(X,y, test_size=0.33, stratify = y, random_state = 148)

KOT.dtypes

CO_GT                   float64
PT08_S1_CO              float64
C6H6_GT                 float64
PT08_S2_NMHC            float64
NOx_GT                  float64
PT08_S3_NOx             float64
NO2_GT                  float64
PT08_S4_NO2             float64
T                       float64
RH                      float64
AH                      float64
Categores_PT08_S5_O3     object
dtype: object

print ('Zbiór X treningowy: ',Xtrain.shape)
print ('Zbiór X testowy:    ', Xtest.shape)
print ('Zbiór y treningowy: ', ytrain.shape)
print ('Zbiór y testowy:    ', ytest.shape)

Zbiór X treningowy:  (6269, 11)
Zbiór X testowy:     (3088, 11)
Zbiór y treningowy:  (6269,)
Zbiór y testowy:     (3088,)

C:ProgramDataAnaconda3libsite-packagessklearnlinear_modellogistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)

GridSearchCV(cv=2, error_score='raise-deprecating',
             estimator=LogisticRegression(C=1.0, class_weight=None, dual=False,
                                          fit_intercept=True,
                                          intercept_scaling=1, l1_ratio=None,
                                          max_iter=100, multi_class='warn',
                                          n_jobs=None, penalty='l2',
                                          random_state=None, solver='warn',
                                          tol=0.0001, verbose=0,
                                          warm_start=True),
             iid='warn', n_jobs=5,
             param_grid={'C': array([1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01, 1.e+02])},
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring='roc_auc', verbose=0)

array([[1414,  131],
       [ 176, 1367]], dtype=int64)

                  precision    recall  f1-score   support

(220.999, 964.0]       0.89      0.92      0.90      1545
 (964.0, 2523.0]       0.91      0.89      0.90      1543

        accuracy                           0.90      3088
       macro avg       0.90      0.90      0.90      3088
    weighted avg       0.90      0.90      0.90      3088

[648.230997747899, 11.43831319456862, 7.939972023928728, 46.44356640618427, 66.47134191824114, 5.220873142678388, 5.43267527445646, 3.714154563031008, 14.52279194616348, 14.359778937836424, 8.176237173633366, 10.619321015188527]

[200.1579775895654, 5.895905610386816, 5.807909737315139, 3.1029777937073924, 2.7926706710126528, 2.4675790616762594]

C:ProgramDataAnaconda3libsite-packagessklearnlinear_modellogistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)

GridSearchCV(cv=2, error_score='raise-deprecating',
             estimator=LogisticRegression(C=1.0, class_weight=None, dual=False,
                                          fit_intercept=True,
                                          intercept_scaling=1, l1_ratio=None,
                                          max_iter=100, multi_class='warn',
                                          n_jobs=None, penalty='l2',
                                          random_state=None, solver='warn',
                                          tol=0.0001, verbose=0,
                                          warm_start=True),
             iid='warn', n_jobs=5,
             param_grid={'C': array([1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01, 1.e+02])},
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring='roc_auc', verbose=0)

array([[1387,  158],
       [ 182, 1361]], dtype=int64)

import numpy as np
from sklearn import model_selection
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV

Parameteres = {'C': np.power(10.0, np.arange(-3, 3))}
LR = LogisticRegression(warm_start = True)
LR_Grid = GridSearchCV(LR, param_grid = Parameteres, scoring = 'roc_auc', n_jobs = 5, cv=2)

LR_Grid.fit(Xtrain, ytrain)

C:ProgramDataAnaconda3libsite-packagessklearnlinear_modellogistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)

GridSearchCV(cv=2, error_score='raise-deprecating',
             estimator=LogisticRegression(C=1.0, class_weight=None, dual=False,
                                          fit_intercept=True,
                                          intercept_scaling=1, l1_ratio=None,
                                          max_iter=100, multi_class='warn',
                                          n_jobs=None, penalty='l2',
                                          random_state=None, solver='warn',
                                          tol=0.0001, verbose=0,
                                          warm_start=True),
             iid='warn', n_jobs=5,
             param_grid={'C': array([1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01, 1.e+02])},
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring='roc_auc', verbose=0)

array([[1414,  131],
       [ 176, 1367]], dtype=int64)

                  precision    recall  f1-score   support

(220.999, 964.0]       0.89      0.92      0.90      1545
 (964.0, 2523.0]       0.91      0.89      0.90      1543

        accuracy                           0.90      3088
       macro avg       0.90      0.90      0.90      3088
    weighted avg       0.90      0.90      0.90      3088

[648.230997747899, 11.43831319456862, 7.939972023928728, 46.44356640618427, 66.47134191824114, 5.220873142678388, 5.43267527445646, 3.714154563031008, 14.52279194616348, 14.359778937836424, 8.176237173633366, 10.619321015188527]

[200.1579775895654, 5.895905610386816, 5.807909737315139, 3.1029777937073924, 2.7926706710126528, 2.4675790616762594]

C:ProgramDataAnaconda3libsite-packagessklearnlinear_modellogistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)

GridSearchCV(cv=2, error_score='raise-deprecating',
             estimator=LogisticRegression(C=1.0, class_weight=None, dual=False,
                                          fit_intercept=True,
                                          intercept_scaling=1, l1_ratio=None,
                                          max_iter=100, multi_class='warn',
                                          n_jobs=None, penalty='l2',
                                          random_state=None, solver='warn',
                                          tol=0.0001, verbose=0,
                                          warm_start=True),
             iid='warn', n_jobs=5,
             param_grid={'C': array([1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01, 1.e+02])},
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring='roc_auc', verbose=0)

array([[1387,  158],
       [ 182, 1361]], dtype=int64)

                  precision    recall  f1-score   support

(220.999, 964.0]       0.88      0.90      0.89      1545
 (964.0, 2523.0]       0.90      0.88      0.89      1543

        accuracy                           0.89      3088
       macro avg       0.89      0.89      0.89      3088
    weighted avg       0.89      0.89      0.89      3088

Ocena modelu regresji logistycznej bez redukcji zjawiska multicollinearity¶

ypred = LR_Grid.predict(Xtest)

from sklearn.metrics import classification_report, confusion_matrix
from sklearn import metrics

co_matrix = metrics.confusion_matrix(ytest, ypred)
co_matrix

array([[1414,  131],
       [ 176, 1367]], dtype=int64)

print(classification_report(ytest, ypred))

                  precision    recall  f1-score   support

(220.999, 964.0]       0.89      0.92      0.90      1545
 (964.0, 2523.0]       0.91      0.89      0.90      1543

        accuracy                           0.90      3088
       macro avg       0.90      0.90      0.90      3088
    weighted avg       0.90      0.90      0.90      3088

[648.230997747899, 11.43831319456862, 7.939972023928728, 46.44356640618427, 66.47134191824114, 5.220873142678388, 5.43267527445646, 3.714154563031008, 14.52279194616348, 14.359778937836424, 8.176237173633366, 10.619321015188527]

[200.1579775895654, 5.895905610386816, 5.807909737315139, 3.1029777937073924, 2.7926706710126528, 2.4675790616762594]

C:ProgramDataAnaconda3libsite-packagessklearnlinear_modellogistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)

GridSearchCV(cv=2, error_score='raise-deprecating',
             estimator=LogisticRegression(C=1.0, class_weight=None, dual=False,
                                          fit_intercept=True,
                                          intercept_scaling=1, l1_ratio=None,
                                          max_iter=100, multi_class='warn',
                                          n_jobs=None, penalty='l2',
                                          random_state=None, solver='warn',
                                          tol=0.0001, verbose=0,
                                          warm_start=True),
             iid='warn', n_jobs=5,
             param_grid={'C': array([1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01, 1.e+02])},
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring='roc_auc', verbose=0)

array([[1387,  158],
       [ 182, 1361]], dtype=int64)

                  precision    recall  f1-score   support

(220.999, 964.0]       0.88      0.90      0.89      1545
 (964.0, 2523.0]       0.90      0.88      0.89      1543

        accuracy                           0.89      3088
       macro avg       0.89      0.89      0.89      3088
    weighted avg       0.89      0.89      0.89      3088

Model wykazuje się doskonałymi zdolnościmi prognostycznymi w klasyfikacji zanieczyszczenia. Nie oznacza to jednak, że wszystko jest doskonale. Wśród zmiennych opisujących występuje wysoka multicollinearity. Jest to niekorzystne zjawiski, które teraz wyelinimujemy.

Eliniminacja zjawiska multicollinearity¶

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.

Kroki wdrażania VIF

Uruchom regresję wieloraką.
Oblicz współczynniki VIF.
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.

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.

from patsy import dmatrices
from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.formula.api as smf

lm = smf.ols(formula = 'Categores_PT08_S5_O3 ~ CO_GT+PT08_S1_CO+C6H6_GT+PT08_S2_NMHC+NOx_GT+PT08_S3_NOx+NO2_GT+PT08_S4_NO2+T+RH+AH', data = KOT).fit()
y, X = dmatrices('Categores_PT08_S5_O3 ~ CO_GT+PT08_S1_CO+C6H6_GT+PT08_S2_NMHC+NOx_GT+PT08_S3_NOx+NO2_GT+PT08_S4_NO2+T+RH+AH', data = KOT, return_type = "dataframe")
vif = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
print(vif)

[648.230997747899, 11.43831319456862, 7.939972023928728, 46.44356640618427, 66.47134191824114, 5.220873142678388, 5.43267527445646, 3.714154563031008, 14.52279194616348, 14.359778937836424, 8.176237173633366, 10.619321015188527]

[200.1579775895654, 5.895905610386816, 5.807909737315139, 3.1029777937073924, 2.7926706710126528, 2.4675790616762594]

C:ProgramDataAnaconda3libsite-packagessklearnlinear_modellogistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)

GridSearchCV(cv=2, error_score='raise-deprecating',
             estimator=LogisticRegression(C=1.0, class_weight=None, dual=False,
                                          fit_intercept=True,
                                          intercept_scaling=1, l1_ratio=None,
                                          max_iter=100, multi_class='warn',
                                          n_jobs=None, penalty='l2',
                                          random_state=None, solver='warn',
                                          tol=0.0001, verbose=0,
                                          warm_start=True),
             iid='warn', n_jobs=5,
             param_grid={'C': array([1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01, 1.e+02])},
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring='roc_auc', verbose=0)

array([[1387,  158],
       [ 182, 1361]], dtype=int64)

                  precision    recall  f1-score   support

(220.999, 964.0]       0.88      0.90      0.89      1545
 (964.0, 2523.0]       0.90      0.88      0.89      1543

        accuracy                           0.89      3088
       macro avg       0.89      0.89      0.89      3088
    weighted avg       0.89      0.89      0.89      3088

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

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

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

RFE_list

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:¶

C6H6_GT
PT08_S2_NMHC
PT08_S4_NO2Eliminuje też zmienne gdzie korelacja ze zmienną zależną było bardzo małe a zmienne miały VIF powyżej 5:
T
RH
AH

from patsy import dmatrices
from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.formula.api as smf

lm = smf.ols(formula = 'Categores_PT08_S5_O3 ~ CO_GT+PT08_S1_CO+NOx_GT+PT08_S3_NOx+NO2_GT', data = KOT).fit()
y, X = dmatrices('Categores_PT08_S5_O3 ~ CO_GT+PT08_S1_CO+NOx_GT+PT08_S3_NOx+NO2_GT', data = KOT, return_type = "dataframe")
vif = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
print(vif)

[200.1579775895654, 5.895905610386816, 5.807909737315139, 3.1029777937073924, 2.7926706710126528, 2.4675790616762594]

C:ProgramDataAnaconda3libsite-packagessklearnlinear_modellogistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)

GridSearchCV(cv=2, error_score='raise-deprecating',
             estimator=LogisticRegression(C=1.0, class_weight=None, dual=False,
                                          fit_intercept=True,
                                          intercept_scaling=1, l1_ratio=None,
                                          max_iter=100, multi_class='warn',
                                          n_jobs=None, penalty='l2',
                                          random_state=None, solver='warn',
                                          tol=0.0001, verbose=0,
                                          warm_start=True),
             iid='warn', n_jobs=5,
             param_grid={'C': array([1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01, 1.e+02])},
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring='roc_auc', verbose=0)

array([[1387,  158],
       [ 182, 1361]], dtype=int64)

                  precision    recall  f1-score   support

(220.999, 964.0]       0.88      0.90      0.89      1545
 (964.0, 2523.0]       0.90      0.88      0.89      1543

        accuracy                           0.89      3088
       macro avg       0.89      0.89      0.89      3088
    weighted avg       0.89      0.89      0.89      3088

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

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

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

RFE_list

Teraz zmienne, pomimo wysokich stanów przekraczających 5, to są: CO_GT i PT08_S1_CO są akceptowane w zakresie poziomu multicollinearity. Tworzę nowy model klasyfikacji na okrojonym zbiorze zmienych opisujących.

KOT2 = KOT[['Categores_PT08_S5_O3', 'CO_GT', 'PT08_S1_CO','NOx_GT', 'PT08_S3_NOx', 'NO2_GT']]

KOT2.sample(3)

Budujemy model klasyfikacji po VIF¶

#from sklearn.model_selection import train_test_split 


y = KOT2['Categores_PT08_S5_O3']
X = KOT2.drop('Categores_PT08_S5_O3', axis=1)


Xtrain, Xtest, ytrain, ytest = train_test_split(X,y, test_size=0.33, stratify = y, random_state = 148)

#import numpy as np
#from sklearn import model_selection
#from sklearn.pipeline import make_pipeline
#from sklearn.linear_model import LogisticRegression
#from sklearn.model_selection import GridSearchCV

Parameteres = {'C': np.power(10.0, np.arange(-3, 3))}
LR = LogisticRegression(warm_start = True)
LR_VIF = GridSearchCV(LR, param_grid = Parameteres, scoring = 'roc_auc', n_jobs = 5, cv=2)

LR_VIF.fit(Xtrain, ytrain)

C:ProgramDataAnaconda3libsite-packagessklearnlinear_modellogistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.
  FutureWarning)

GridSearchCV(cv=2, error_score='raise-deprecating',
             estimator=LogisticRegression(C=1.0, class_weight=None, dual=False,
                                          fit_intercept=True,
                                          intercept_scaling=1, l1_ratio=None,
                                          max_iter=100, multi_class='warn',
                                          n_jobs=None, penalty='l2',
                                          random_state=None, solver='warn',
                                          tol=0.0001, verbose=0,
                                          warm_start=True),
             iid='warn', n_jobs=5,
             param_grid={'C': array([1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01, 1.e+02])},
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring='roc_auc', verbose=0)

array([[1387,  158],
       [ 182, 1361]], dtype=int64)

                  precision    recall  f1-score   support

(220.999, 964.0]       0.88      0.90      0.89      1545
 (964.0, 2523.0]       0.90      0.88      0.89      1543

        accuracy                           0.89      3088
       macro avg       0.89      0.89      0.89      3088
    weighted avg       0.89      0.89      0.89      3088

Ocena modelu regresji logistycznej bez redukcji zjawiska multicollinearity¶

ypred = LR_VIF.predict(Xtest)

from sklearn.metrics import classification_report, confusion_matrix
from sklearn import metrics

co_matrix = metrics.confusion_matrix(ytest, ypred)
co_matrix

array([[1387,  158],
       [ 182, 1361]], dtype=int64)

print(classification_report(ytest, ypred))

                  precision    recall  f1-score   support

(220.999, 964.0]       0.88      0.90      0.89      1545
 (964.0, 2523.0]       0.90      0.88      0.89      1543

        accuracy                           0.89      3088
       macro avg       0.89      0.89      0.89      3088
    weighted avg       0.89      0.89      0.89      3088

Nowy model zbudowany na pięciu niezwiązanych przez multicollinearity jest nieznacznie gorszy (na poziomie 1-2

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

Uruchom regresję wieloraką.
Oblicz współczynniki VIF.
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.

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)

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

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

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.

df.isnull().sum()

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

df.dtypes

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

df.columns

Index(['nr.', 'factorA', 'factorB', 'citric catoda', 'residual butanol',
       'caroton', 'stable nodinol', 'sulfur in nodinol', 'density', 'pH',
       'noracid', 'lacapon', 'quality class'],
      dtype='object')

df.columns=['nr.', 'factorA', 'factorB', 'citric_catoda', 'residual_butanol',
       'caroton', 'stable_nodinol', 'sulfur_in_nodinol', 'density', 'pH',
       'noracid', 'lacapon', 'quality_class']

#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.

df.columns

Index(['nr.', 'factorA', 'factorB', 'citric_catoda', 'residual_butanol',
       'caroton', 'stable_nodinol', 'sulfur_in_nodinol', 'density', 'pH',
       'noracid', 'lacapon', 'quality_class'],
      dtype='object')

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]

Index(['Intercept', 'factorA', 'factorB', 'citric_catoda', 'residual_butanol',
       'caroton', 'stable_nodinol', 'density', 'sulfur_in_nodinol', 'pH',
       'noracid', 'lacapon'],
      dtype='object')

[809.2230489081707, 1.6530213053949512, 2.1806661580635818, 1.093561268844909, 1.3702675584858117, 1.9475860748759193, 2.0181057450324356, 1.6077595870352759, 1.3354874384986997, 1.2793132834168612]

Ustawiamy macierz wyników.

X.columns

Index(['Intercept', 'factorA', 'factorB', 'citric_catoda', 'residual_butanol',
       'caroton', 'stable_nodinol', 'density', 'sulfur_in_nodinol', 'pH',
       'noracid', 'lacapon'],
      dtype='object')

vif  =  np.round(vif, decimals=2)

vif = list(map(float, vif))
name = list(X)

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

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

RFE_list

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.¶

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]

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

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

	Unnamed: 0	Date	Time	CO(GT)	PT08.S1(CO)	C6H6(GT)	PT08.S2(NMHC)	NOx(GT)	PT08.S3(NOx)	NO2(GT)	PT08.S4(NO2)	PT08.S5(O3)	T	RH	AH
0	0	10/03/2004	18.00.00	2.6	1360.0	11.9	1046.0	166.0	1056.0	113.0	1692.0	1268.0	13.6	48.9	0.7578
1	1	10/03/2004	19.00.00	2.0	1292.0	9.4	955.0	103.0	1174.0	92.0	1559.0	972.0	13.3	47.7	0.7255
2	2	10/03/2004	20.00.00	2.2	1402.0	9.0	939.0	131.0	1140.0	114.0	1555.0	1074.0	11.9	54.0	0.7502

	Categores_PT08_S5_O3	CO_GT	PT08_S1_CO	NOx_GT	PT08_S3_NOx	NO2_GT
790	(220.999, 964.0]	1.4	1035.0	96.0	1222.0	82.0
2484	(220.999, 964.0]	0.9	890.0	93.0	922.0	67.0
4306	(220.999, 964.0]	0.5	828.0	72.0	1281.0	48.0

	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

	name	vif
0	Intercept	648.23
1	CO_GT	11.44
2	PT08_S1_CO	7.94
3	C6H6_GT	46.44
4	PT08_S2_NMHC	66.47
5	NOx_GT	5.22
6	PT08_S3_NOx	5.43
7	NO2_GT	3.71
8	PT08_S4_NO2	14.52
9	T	14.36
10	RH	8.18
11	AH	10.62

	name	vif
0	Intercept	200.16
1	CO_GT	5.90
2	PT08_S1_CO	5.81
3	NOx_GT	3.10
4	PT08_S3_NOx	2.79
5	NO2_GT	2.47

	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

	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

vif - THE DATA SCIENCE LIBRARY

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

Zadanie

Wskazujemy zmienne niezależne:¶

Ocena modelu regresji logistycznej bez redukcji zjawiska multicollinearity¶

Eliniminacja zjawiska multicollinearity¶

Variance Inflation Factor (VIF)¶

Jeszcze raz tworzymy macierz korelacji vif i sprawdzamy zmienne po wyeliminowaniu:¶

Budujemy model klasyfikacji po VIF¶

Ocena modelu regresji logistycznej bez redukcji zjawiska multicollinearity¶

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

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.¶

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

Variance Inflation Factor (VIF)

Tworzymy model regresji liniowej

VIF (Variation Inflation Factor)¶

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