Prawidłowo utworzony model regresji liniowej jest odwzorowaniem rzeczywistości. Może być on podstawą systemów diagnostycznych lub automatyki przemysłowej. Stały monitoring właściwości procesu przez model statystyczny skutecznie wskazywać na pojawiające się anomalie i w ten sposób przepowiadać zbliżającą się awarię instalację przemysłowej.

W poprzedniej części tego artykułu omówiliśmy proces tworzenia modelu regresji jednoczynnikowej i regresji wielorakiej. Dzisiaj rozwiniemy ten temat pod kontem analizy modelu i analizy statystycznej zmiennych niezależnych.

Na początek przypomnimy warunki które muszą być spełnione aby utworzyć prawidłowy modelu regresji.

1. Między predyktorami i zmienną wynikową musi być zachowana relacja liniowa.
2. Nie może istnieć silna korelacja pomiędzy zmiennymi niezależnymi.
3. Wartości resztkowe powinny mieć rozkład normalność. Indywidualne reszty modelu powinny być podobne w swej zbiorowości oraz mieć rozkład normalny.
4. Homogeniczność wariancji - wariancja błędu powinna być stała.
5. Niezależność - błędy związane z jedną obserwacją nie mogą być skorelowane z błędami innych obserwacji.

Parametry modelu regresji wielorakiej

Analizę w języku Python zaczynamy od wgrania potrzebnych bibliotek oraz pobraniu źródła danych.

import pandas as pd
import numpy as np
import itertools
from itertools import chain, combinations
import statsmodels.formula.api as smf
import scipy.stats as scipystats
import statsmodels.api as sm
import statsmodels.stats.stattools as stools
import statsmodels.stats as stats 
from statsmodels.graphics.regressionplots import *
import matplotlib.pyplot as plt
import seaborn as sns
import copy
#from sklearn.cross_validation import train_test_split
import math
#18:00
## https://archive.ics.uci.edu/ml/datasets/Combined+Cycle+Power+Plant

df = pd.read_excel('c:/1/Folds5x2_pp.xlsx')

df.columns = ['Temperature', 'Exhaust_Vacuum', 'Ambient_Pressure', 'Relative_Humidity', 'Energy_output']

Tworzymy model wieloczynnikowej regresji liniowej.

lm = smf.ols(formula = 'Energy_output ~ Temperature + Exhaust_Vacuum + Relative_Humidity + Ambient_Pressure', data = df).fit()

Wyświetlenie podsumowania modelu.

lm.summary()

 Tworzymy słownik współczynników 

coef = lm.params.to_dict()

Teraz możemy sobie przywołać dowolny współczynnik przy x.

coef['Temperature']

-1.977513106635395

coef['Intercept']

454.6092743153136

 Wyświetlenie reszt modelu 

Reszty modelu (inaczej błędy modelu) to różnice pomiędzy wartościami empirycznymi a wartościami wynikającącymi z estymacji modelu.

lm.resid.sample(6)

 Wyświetlenia wyników predykcji modelu 

Generujemy wartości predykcyjne modelu

lm.predict()

Tworzymy dodatkowe kolumny w tabeli danych zawierającą predykcje modelu oraz wartości reszt, będących różnicą pomiędzy wartościami a ich estymacją.

df['Predict'] = pd.Series(lm.predict())
df['Resid'] = pd.Series(lm.resid)
df.sample(8)[['Energy_output','Predict','Resid']]

Punkt 1. Analiza normalności rozkładu reszt modelu regresji

sns.kdeplot(np.array(df['Resid']), bw=10)

Wykres rozkładu normalnego reszt modelu regresji. Niepokojące jest wydłużenie lewego ogona wykresu (zaznaczone na czerwono).

sns.distplot(np.array(df['Resid']))

 Test normalności rozkładu Anderson-Darling 

Test Anderson-Darling został opracowany w roku 1952 przez Teodora Andersona i Donalda Darlinga.

Jest on najczęściej wykorzystywanym testem wśród testów badających dopasowanie zestawu danych do rozkładu normalnego. Rozkład prawdopodobieństwa zdarzeń jest normalny jeżeli pokrywa się z linią przekątną.

Można skorzystać z prostej formy wykresu.

sm.qqplot(df['Resid'],color='r')

import pylab
scipystats.probplot(df['Resid'], dist="norm", plot=pylab)

 Wartość wskaźnika Anderson-Darling 

import scipy
scipy.stats.anderson(df['Resid'], dist='norm' ) 

Jeśli zwrócona statystyka (parametr: statistic) jest większa niż wartości krytyczne (parametr: critical_values, to dla odpowiedniego poziomu istotności (parametr: significance_level ) można odrzucić hipotezę zerową, mówiąca, że dane posiadają rozkład normalny.

Test Anderson-Darling odrzucił hipotezę zerową mówiącą, że rozkład reszt modelu ma jest normalny.

 Test Kołmogorowa-Smirnowa 

from scipy import stats
stats.kstest(df['Resid'], 'norm')

Wartość statystyki jest większa od pvalue co wskazuje na odrzucenie hipotezy zerowej mówiącej o normalności rozkładu reszt modelu.

Porównanie rozkładów wykonania i predykcji.

sns.kdeplot(np.array(df['Energy_output']), bw=10)
sns.distplot(np.array(df['Predict']))

Punkt 2. Testy na jednorodność wariancji reszt

Drugim głównych założeń dla zwykłej regresji jest jednorodność wariancji reszt. Jeśli model jest dobrze dopasowany, reszty powinny być jednorodne.

resid = lm.resid
plt.scatter(lm.predict(), resid)

Wykres wskazuje na częściową niejednorodność wariancji reszt modelu. Powód może być ten sam co w przypadku braku normalności rozkładu dla reszt modelu.

 

Test na autokorelacje reszt modelu

Jednym z założeń dotyczących modelu regresji jest niezależność błędów obserwacji, czyli fakt, czy występujące reszty w predykcji zmiennej zależnej są ze sobą skorelowane. Dobrze dopasowane modele regresji zakładają, że otrzymywane reszty - są niezależne od siebie. Innymi słowy, ich rozkład jest losowy, przypadkowy, bez stale występującego wzorca.

Jednym ze sposobów określenia niezależności błędów obserwacji jest wyznaczenie autokorelacji składnika resztowego, czyli korelacji r-Pearsona pomiędzy kolejnymi resztami, powstałymi z nieidealnego dopasowania modelu. Test Durbina-Watsona (statystyka) służy do oceny występowania korelacji pomiędzy resztami (błędami, składnikami resztowymi).

Źródło: http://www.naukowiec.org/wiedza/statystyka/test-durbina-watsona-niezaleznosc-bledow-obserwacji_423.html

 

import statsmodels
statsmodels.stats.stattools.durbin_watson(df['Resid'], axis=0)

2.032935807327476

Statystyka testu jest równa 2, co wskazujące na brak korelacji szeregowej pomiędzy resztami w predykcji zmiennej zależnej. Reszty więc nie są ze sobą skorelowane. Statystyka ta zawsze będzie zawierać się w przedziale od 0 do 4. Im statystyka bliższa 0, tym więcej dowodów na dodatnią korelację szeregową. Im bliżej 4, tym więcej dowodów na ujemną korelację szeregową.

Test współliniowości zmiennych niezależnych

Jednym z założeń analizy regresji jest brak występowania współliniowości predyktorów (zmiennych objaśniających). Wprowadzając do modelu regresji silnie skorelowane ze sobą zmienne wprowadzamy do modelu (przy każdej zmiennej) małą bądź zerową unikalną "część wyjaśnienia" zmiennej zależnej. W zależności od sposobu liczenia jeden z predyktorów silnie powiązanych straci swoją "moc" przewidywania zmiennej zależnej na rzecz innego predyktora, z racji występowania silnego związku pomiędzy nimi.

Współliniowość oznacza sytuację, w której predyktory wprowadzane do modelu regresji są ze sobą silnie skorelowane i taka sytuacja powoduje „pogorszenie się” parametrów modelu.

Najczęściej współliniowość szacowana jest poprzez dwa parametry: wskaźnik tolerancji oraz parametr VIF (CIW – czynnik inflacji wariancji). Należy zaznaczyć, że parametry te obliczane są dla każdego z predyktorów wprowadzonych do modelu. Dzięki temu możliwe jest oszacowanie, który z nich wprowadza do modelu współliniowość.

Źródło: http://www.naukowiec.org/wiedza/statystyka/wspolliniowosc-zmiennych-objasniajacych_954.html

 

from patsy import dmatrices
from statsmodels.stats.outliers_influence import variance_inflation_factor

lm = smf.ols(formula = 'Energy_output ~ Temperature + Exhaust_Vacuum + Relative_Humidity + Ambient_Pressure', data = df).fit()
y, X = dmatrices('Energy_output ~ Temperature + Exhaust_Vacuum + Relative_Humidity + Ambient_Pressure', data = df, return_type = "dataframe")
vif = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
print(vif)

Jedną z rekomendacji jest to, że jeśli VIF jest większy niż 5, to zmienna objaśniająca podana przez parametr: exog_idx jest wysoce współliniowa z innymi zmiennymi objaśniającymi, a estymatory parametrów będą miały z tego powodu duże błędy standardowe.

Źródło: https://www.statsmodels.org/stable/generated/statsmodels.stats.outliers_influence.variance_inflation_factor.html

Test wykazał wysoką współliniowość parametru Temperatura. Analiza ta potwierdza test korelacji zmiennych egzogenicznych wykonany w pierwszej części tego artykułu.

Wnioski końcowe

Aby model regresji wielorakiej mógł być wykorzystany produkcyjnie musi zmienna niezależna Temperatura musi zostać wyłączona ze zmiennych opisujących. Po wyłączeniu tej zmiennej wszystkie czynności opisane w tej publikacji powinny zostać powtórzone pod kontem analizy zgodności warunków regresji liniowej.

Możliwe jest również utworzenie modelu regresji jednoczynnikowej w oparciu o zmienne egzogeniczne Temperatura lub Próżnia spalin.