Rozmawiająć z osobami, które zawodowo wykorzystują metody analizy danych byłem zaskoczony jak wiele modeli predykcyjnych opartych jest na regresji liniowej. Jest to jedna z podstawowych technik w arsenale analityka, stosukowo prosta w implementacji oraz zrozumieniu, a jednak niezwykle efektywna i użyteczna. Dzięki bibliotece scikit-learn jesteśmy w stanie w kilku wierszach kodu python’a zaimplementować jej podstawowe rodzaje.

Regresja liniowa wprowadzenie

Celem regresji liniowej jest dopasowanie do danych doświadczalnych/historycznych prostej, która oddaje charakter tych danych. Zadanie polega na znalezieniu funkcji liniowej \(f(x) \rightarrow \mathbf{R}\) postaci:
\begin{equation}
f(x) = w_0 + w_1x_1+w_2x_2+…w_nx_n
\label{eq1}
\end{equation}

gdzie \(x \in R^N\) a \(w \in R^N\) jest wektorem zawierającym współczynniki prostej. I to właśnie na znalezieniu wektora \(w\) będziemy skupiać swoją uwagę.

Metoda najmniejszych kwadratów

Mając dane historyczne  \(X = \{x^1,x^2,…,x^k \} , x^i \in R^N\)  oraz odpowiadające im wartości \(Y=\{y^1, y^2,…,y^k\}\) będziemy poszukiwać najlepiej dopasowanego wektora \(w\). Co to dla nas znaczy, najlepiej dopasowanego?

Mianowicie takiego, dla którego wartości zwracane przez funkcję \(f\) i odpowiadająca wartości rzeczywiste najmniej się różnią. Załóżmy że mamy już wybrany wektor np. \(w=[1,0.5,0,0,2,0.3,…,1]\) możemy obliczyć:

\begin{align}
f(x^1)= & w_0 + w_1x^1_1+w_2x^1_2+…w_nx^1_n = w_0 + (x^1)^T*w\\
f(x^2)= & w_0 + w_1x^2_1+w_2x^2_2+…w_nx^2_n = w_0 + (x^2)^T*w \\
\cdots  &\\
f(x^k)= & w_0 + w_1x^k_1+w_2x^k_2+…w_nx^k_n =  w_0 + (x^k)^T*w\\
\end{align}

Wartości te możemy porównać z rzeczywistymi odpowiadającymi wartościami \(y\), stąd mamy

\begin{align}
e^1 = &f(x^1)- y^1 \\
e^2 = &f(x^2)- y^2\\
\cdots  &\\
e^k= &f(x^k)- y^k \\
\end{align}

Następnie sumując kwadraty tych błędów \(e^1, …, e^k\) otrzymamy wyrażenie:

$$E(w) =  \sum_{i=1}^{k}\left(f(x^i)-y^i\right)^2 $$

Oczywiście chcielibyśmy aby suma kwadratów błędów była jak najmniejsza, a jedyne co możemy zmieniać to wartości współczynników wektora \(w\). W ten sposób zdefiniowaliśmy zadanie minimalizacji \(E\) względem \(w\). Aby skrócić i ułatwić sobie zapis powyższe można zapisać w postaci macierzowej, jeżeli \(X\) zapiszemy jako macierz, w której poszczególne egzemplarze \(x^i\) zostaną ułożone w wierszach, wektor \(w\) jest wektorem kolumnowym, a \(Y\) potraktujemy także jako wektor kolumnowy całość możemy zapisać w postaci:

$$ \underset{w}{min} ||X*w-Y||_2^2 $$

 Rodzaje regresji liniowej

W zależności od definicji zadania minimalizacji możemy wyróżnić kilka typów regresji:

• regresja liniowa $$ \underset{w}{min} ||X*w-Y||_2^2 $$
• regresja grzbietowa (Ridge) $$ \underset{w}{min} ||X*w-Y||_2^2 + \alpha ||w||_2^2$$
• regresja Lasso $$ \underset{w}{min} ||X*w-Y||_2^2 + \alpha ||w||_1$$
• Elastic Net $$\underset{w}{min} { \frac{1}{2n_{samples}} ||X w – y||_2 ^ 2 + \alpha \rho ||w||_1 + \frac{\alpha(1-\rho)}{2} ||w||_2 ^ 2} $$

Kolejne trzy typy regresji dokładają do zadania optymalizacji dodatkowy czynnik (penalty, regularization) wpływający na cechy wektora \(w\). W regresji Ridge dodatkowo dodajemy warunek aby norma \(l_2\) wektora \(w\) była minimalizowana, dzięki czemu wymuszamy aby poszczególne współczynniki nie były za duże. Natomiast w regresji Lasso użycie normy \(l_1\) preferuje rozwiązania, w których wektor jest rzadki (ma dużo zer). Regresja Elastic Net łączy w sobie cechy Ridge oraz Lasso, poprzez zastosowanie dwóch typów norm, których istotność względem siebie kontroluje parametr \(\rho\)

Regresja liniowa przykłady w scikit-learn

No dobra, dość teorii! Czas zacząć kodować. Dzięki wykorzystaniu biblioteki scikit-learn wraz z zaimplementowanymi algorytmami regresji jesteśmy w stanie bardzo szybko nauczyć nasz model.

Implementacja dowolnego typu regresji jest bardzo prosta, tworzymy obiekt odpowiedniej klasy (LinearRegression, Ridge, Lasso, ElasticNet), na którym wywołujemy metodę fit podając jej za argumenty zbiór treningowy, docelowe wartości oraz ewentualne parametry (\(\alpha, \rho\)). Gdy model się nauczy to możemy wywołać metodę predict, która zwróci przewidywaną wartość. Ogólny szablon przedstawia się następująco

import numpy as np
from sklearn import datasets, linear_model

#########
# Load some data !
#...
#########
# Create linear regression object
regr = linear_model.LinearRegression()

# Train the model using the training sets
regr.fit(X_train, Y_train)
# Predict values
Y_predicted = regr.predict(X_test)

# The coefficients
print('Coefficients: \n', regr.coef_)
# The mean square error
error = np.mean((regr.predict(X_test) - Y_test) ** 2)
print("Residual sum of squares: {}".format(error))

Dodatkowo warto zobaczyć jak wyglądają parametry prostej (regr.coef) oraz obliczyć rzeczywisty popełniony błąd średnio-kwadratowy.

Przykład 1. Regresja liniowa dla zbioru Diabetes

Zróbmy wizualizację prostych regresji obliczonych dla poszczególnych zmiennych w zbiorze Diabetes.

import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model

# Load the diabetes dataset
diabetes = datasets.load_diabetes()

f, axarr = plt.subplots(5,2, sharex=True, sharey=True,figsize=(12,12))
for i in range(0,5):
    for j in range(0,2):
        # Use only one feature
        diabetes_X = diabetes.data[:, np.newaxis, i*2+j]
        
        # Split the data into training/testing sets
        diabetes_X_train = diabetes_X[:-20]
        diabetes_X_test = diabetes_X[-20:]
        
        # Split the targets into training/testing sets
        diabetes_y_train = diabetes.target[:-20]
        diabetes_y_test = diabetes.target[-20:]
        
        # Create linear regression object
        regr = linear_model.LinearRegression()
        
        # Train the model using the training sets
        regr.fit(diabetes_X_train, diabetes_y_train)
        
        # The coefficients
        print('Coefficients: \n', regr.coef_)
        # The mean square error
        print("Residual sum of squares: %.2f"
              % np.mean((regr.predict(diabetes_X_test) - diabetes_y_test) ** 2))
        # Explained variance score: 1 is perfect prediction
        print('Variance score: %.2f' % regr.score(diabetes_X_test, diabetes_y_test))
        
        # Plot outputs
        axarr[i,j].scatter(diabetes_X_test, diabetes_y_test,  color='red')
        axarr[i,j].plot(diabetes_X_test, regr.predict(diabetes_X_test), color='blue', linewidth=1)
        
plt.show()

W wyniku otrzymamy następującą wizualizację:

Przykład 2. Porównanie modeli regresji liniowej

W drugim przykładzie dokonamy porównania modeli regresji pod względem błędu dopasowania (wartości błędu średnio kwadratowego). Ponadto zobaczmy jak wygląda wektor kierunkowy poszczególnych prostych. Posłużymy się zbiorem Boston House Pricing, ma on 14 atrybutów opisujących domy (ilość łazienek itp.). Naszym zadaniem jest przewidzenie ostatniego 14, który jest ceną domu. Zbiór nie jest duży ma tylko 506 instancji, ale dla nas w zupełności wystarczy. Poniżej opis z pliku informacyjnego:

1. Title: Boston Housing Data

2. Sources:
   (a) Origin:  This dataset was taken from the StatLib library which is
                maintained at Carnegie Mellon University.
   (b) Creator:  Harrison, D. and Rubinfeld, D.L. 'Hedonic prices and the 
                 demand for clean air', J. Environ. Economics & Management,
                 vol.5, 81-102, 1978.
   (c) Date: July 7, 1993

3. Past Usage:
   -   Used in Belsley, Kuh & Welsch, 'Regression diagnostics ...', Wiley, 
       1980.   N.B. Various transformations are used in the table on
       pages 244-261.
    -  Quinlan,R. (1993). Combining Instance-Based and Model-Based Learning.
       In Proceedings on the Tenth International Conference of Machine 
       Learning, 236-243, University of Massachusetts, Amherst. Morgan
       Kaufmann.

4. Relevant Information:

   Concerns housing values in suburbs of Boston.

5. Number of Instances: 506

6. Number of Attributes: 13 continuous attributes (including "class"
                         attribute "MEDV"), 1 binary-valued attribute.

7. Attribute Information:

    1. CRIM   -   per capita crime rate by town
    2. ZN     -  proportion of residential land zoned for lots over 
                 25,000 sq.ft.
    3. INDUS  -   proportion of non-retail business acres per town
    4. CHAS   -   Charles River dummy variable (= 1 if tract bounds 
                 river; 0 otherwise)
    5. NOX    -   nitric oxides concentration (parts per 10 million)
    6. RM     -   average number of rooms per dwelling
    7. AGE    -   proportion of owner-occupied units built prior to 1940
    8. DIS    -   weighted distances to five Boston employment centres
    9. RAD    -   index of accessibility to radial highways
    10. TAX   -   full-value property-tax rate per $10,000
    11. PTRATIO-  pupil-teacher ratio by town
    12. B      -  1000(Bk - 0.63)^2 where Bk is the proportion of blacks  by town
    13. LSTAT -   % lower status of the population
    14. MEDV  -   Median value of owner-occupied homes in $1000's

8. Missing Attribute Values:  None.

W naszym przypadku dane załadujemy automatycznie przy wykorzystaniu scikit-learn. Klasa datasets ma szereg metod do załadowania danych w tym load_boston() ładującą właśnie nasz zbiór.

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets, linear_model as linm

# Reggression models
# http://scikit-learn.org/stable/modules/linear_model.html

# Load the diabetes dataset
boston = datasets.load_boston()
# print description
print(boston.DESCR)
# get the data
boston_X = boston.data
boston_Y = boston.target
# Split the data into training/testing sets
boston_X_train = boston_X[:-50]
boston_X_test = boston_X[-50:]

# Split the targets into training/testing sets
boston_y_train = boston_Y[:-50]
boston_y_test = boston_Y[-50:]

Po załadowaniu danych wydzielamy do oddzielnych zmiennych zmienne zależne i niezależne, boston_X zawiera dane, na podstawie których będziemy dokonywać przewidywania docelowej wartości boston_Y. Następnie dokonujemy podziału na dane treningowe i testowe, gdzie ostatnie 50 instancji to dane testowe.

Przypominam, że chcemy dokonać porównania różnych typów regresji. W tym celu deklarujemy słownik zawierający obiekty rozpatrywanych typów regresji: LinearRegression, Ridge, Lasso, ElasticNet. Następnie w pętli dla każdego obliczamy błąd i zapamiętujemy w zmiennej fit_results.

regressors = {}
regressors['LinReg'] =linm.LinearRegression()
regressors['Ridge'] = linm.Ridge(alpha = .5)
regressors['Lasso'] = linm.Lasso(alpha = 5.1)
regressors['ElNet'] =linm.ElasticNet(alpha = .5, l1_ratio=0.5)

fit_results={}


for key in regressors:
    # Train the model using the training sets
    regr = regressors[key]
    regr.fit(boston_X_train, boston_y_train)
    # mean square error
    mse = np.mean((regr.predict(boston_X_test) - boston_y_test) ** 2)
    w = regr.coef_
    # l1 norm
    wl1 = np.sum(np.abs(w))
    # l2 norm
    wl2 = np.sqrt(np.sum(w**2))
    fit_results[key]= { 'mse': mse, 'wl2': wl2, 'wl1': wl1, 'w': w}
    print("{}\n----------\n  mse={}\n  wl1={}\n  wl2={}\n  w={}\n ".format(key,mse,wl1,wl2,w))
    
    

Na końcu wyświetlmy obliczone błędy regresji na wykresie, tak aby łatwo można było dokonać porównania. W pętli pobieramy uprzednio obliczone wartości i rysujemy w postaci wykresu słupkowego. Dla każdego typu regresji mamy trzy parametry (mse, wl1, wl2), które odkładamy na w postaci trzech słupków. Każdy typ regresji otrzymał inny kolor.

groups = 3
index = np.arange(groups)
bar_width = .2
opacity = 0.4

fig, ax = plt.subplots(figsize=(15,9))

t=0
for key in regressors:
    results = fit_results[key]
    res_val = (results['mse'],results['wl1'],results['wl2'])
    plt.bar(index+ bar_width*t, res_val, bar_width,
                 alpha=opacity,
                 color=np.random.rand(3,1),
                 label=key)
    t+=1

#plt.xlabel('Modele regresji')
plt.title('Porownanie modeli regresji: MSE, wl1, wl2')
plt.xticks(index + (t-2)*bar_width, ('MSE', 'wl1', 'wl2'))
plt.legend()

plt.tight_layout()
plt.show()

W wyniku otrzymamy diagram dokonujący porównania trzech wartości:

• MSE – mean square error, czyli błędu dopasowania,
• wl1 – wartość normy l1 wektora w,
• wl2 – wartość normy l2 wektora w

Na powyższym wykresie mamy zobrazowany trzy rozpatrywane parametry MSE, normę l1 i l2 dla każdego typu regresji: ElasticNet(zielony), Linear regression (czerwony), Ridge regression (szary), Lasso (żółty).

Pokazuje on kilka istotnych informacji, które mogą posłużyć do wyboru modelu. Po pierwsze, zależy nam na tym, aby MSE był jak najmniejszy. W naszym przykładzie najmniejszą wartość osiągnęła Liniowa Regresja (kolor czerwony), tuż za nią Ridge (szary).

Warto także zwrócić uwagą na dwa pozostałe policzone parametry, normy l1 i l2 charakteryzujące znaleziony wektor 'w’. Ze względu na stabilność numeryczną zależy nam, aby wartości tego wektora nie były duże (norma l2 mała). Pod względem tego kryterium lepiej wypada ElasticNet i Lasso.

Natomiast w przypadku gdy mamy do czynienia z danymi wielowymiarowymi chcielibyśmy, aby wektor 'w’ był rzadki (norma l1 mała). W tym przypadku Lasso (kolor żółty) i ElasticNet (zielony) promują rozwiązania rzadkie.

Polecam poczytać o zaletach i wadach poszczególnych modeli:

• https://www.analyticsvidhya.com/blog/2015/08/comprehensive-guide-regression/
• http://www.datasciencecentral.com/profiles/blogs/10-types-of-regressions-which-one-to-use