Machine-Learning

Objetivo

O objetivo deste trabalho é avaliar a base de dados "Boston House Prices" utilizando o algoritmo de seleção de variáveis LASSO e comparar o resultado com o modelo tradicional de regressão linear.

Carga inicial

Será utilizada a biblioteca scikit-learn para o trabalho.

A partir da biblioteca, é possível carregar diretamente os dados que serão utilizados. Esses dados serão carregados em um DataFrame para facilitar a visualização e a manipulação.

In [1]:
from sklearn.datasets import load_boston

boston = load_boston()
print(boston.DESCR)
Boston House Prices dataset

Notes
------
Data Set Characteristics:  

    :Number of Instances: 506 

    :Number of Attributes: 13 numeric/categorical predictive
    
    :Median Value (attribute 14) is usually the target

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

    :Missing Attribute Values: None

    :Creator: Harrison, D. and Rubinfeld, D.L.

This is a copy of UCI ML housing dataset.
http://archive.ics.uci.edu/ml/datasets/Housing


This dataset was taken from the StatLib library which is maintained at Carnegie Mellon University.

The Boston house-price data of Harrison, D. and Rubinfeld, D.L. 'Hedonic
prices and the demand for clean air', J. Environ. Economics & Management,
vol.5, 81-102, 1978.   Used in Belsley, Kuh & Welsch, 'Regression diagnostics
...', Wiley, 1980.   N.B. Various transformations are used in the table on
pages 244-261 of the latter.

The Boston house-price data has been used in many machine learning papers that address regression
problems.   
     
**References**

   - Belsley, Kuh & Welsch, 'Regression diagnostics: Identifying Influential Data and Sources of Collinearity', Wiley, 1980. 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.
   - many more! (see http://archive.ics.uci.edu/ml/datasets/Housing)

In [2]:
import pandas as pd
boston_df = pd.DataFrame(boston.data)
boston_df.columns = boston.feature_names
boston_df['Price'] = boston.target

pd.options.display.float_format = '{:,.3g}'.format # Just Format!
boston_df.head()
Out[2]:
CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT Price
0 0.00632 18 2.31 0 0.538 6.58 65.2 4.09 1 296 15.3 397 4.98 24
1 0.0273 0 7.07 0 0.469 6.42 78.9 4.97 2 242 17.8 397 9.14 21.6
2 0.0273 0 7.07 0 0.469 7.18 61.1 4.97 2 242 17.8 393 4.03 34.7
3 0.0324 0 2.18 0 0.458 7 45.8 6.06 3 222 18.7 395 2.94 33.4
4 0.0691 0 2.18 0 0.458 7.15 54.2 6.06 3 222 18.7 397 5.33 36.2

Regressão Linear

A regressão linear será realizada utilizando o método dos mínimos quadrados. Apresentamos abaixo os coeficientes calculados.

In [3]:
from sklearn.linear_model import LinearRegression
import numpy as np

lr = LinearRegression()

_x = boston_df.drop('Price', axis=1)
y = boston_df['Price']

# Normalization
mu = np.mean(_x, axis=0)
sigma = np.std(_x, axis=0)
x = (_x - mu) / sigma

lr.fit(x, y)

coef_df = pd.DataFrame(index=np.append(['Intercept'], boston.feature_names))
coef_df['OLS'] = np.append([lr.intercept_], lr.coef_) 
coef_df.T
Out[3]:
Intercept CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT
OLS 22.5 -0.92 1.08 0.143 0.682 -2.06 2.67 0.0211 -3.1 2.66 -2.08 -2.06 0.857 -3.75

Apresentamos abaixo a média do quadrado dos resídudos, o coeficiente de determinação $R^2$ ...

In [4]:
from sklearn.metrics import mean_squared_error, r2_score

mse_ls = mean_squared_error(y, lr.predict(x))
r2_ls = r2_score(y, lr.predict(x))

result_df = pd.DataFrame(index=['MSE', 'R\u00b2'])
result_df['OLS'] = [mse_ls, r2_ls]
result_df
Out[4]:
OLS
MSE 21.9
0.741

... e o gráfico com a diferença entre o preço real e o calculado.

In [5]:
import matplotlib.pyplot as plt
%matplotlib inline

def scatterPlot(actual, predicted):
    plt.scatter(actual, predicted)
    range = [actual.min(), actual.max()]
    plt.plot(range, range, 'black')
    plt.xlabel("Preco Real")
    plt.ylabel("Preco Calculado")
    plt.show()
    
scatterPlot(y, lr.predict(x))

Validação Cruzada

É possível verificar no gráfico acima que, para os preços mais altos, o modelo não representou bem os dados.

Faremos uma validação cruzada para verificar esse comportamento. Os dados serão divididos em 5 grupos.

In [6]:
from sklearn.cross_validation import KFold

kf = KFold(len(x), n_folds=5)
p = np.zeros_like(y)
for k, (train, test) in enumerate(kf):
    lr.fit(x.ix[train], y[train])
    p[test] = lr.predict(x.ix[test])
    print("[Fold {0}] R\u00b2: {1:.2f}".
          format(k+1, lr.score(x.ix[test], y[test])))

mse_ls_cv = mean_squared_error(y, p)
r2_ls_cv = r2_score(y, p)

result_df['OLS - CV'] = [mse_ls_cv, r2_ls_cv]
result_df
[Fold 1] R²: 0.64
[Fold 2] R²: 0.71
[Fold 3] R²: 0.59
[Fold 4] R²: 0.08
[Fold 5] R²: -0.26
Out[6]:
OLS OLS - CV
MSE 21.9 37.2
0.741 0.56

É possível verificar que o $R^2$ foi baixo para os dois últimos grupos.

In [7]:
scatterPlot(y, p)

Regressão LASSO

Utilizaremos a técnica Least Absolute Shrinkage and Selection Operator - LASSO para selecionar as variáveis explicativas utilizando um parâmetro de regularização $\alpha = 0.13$.

In [8]:
from sklearn.linear_model import Lasso

alpha = 0.13

las = Lasso(alpha=alpha)
las.fit(x, y)

mse_lasso = mean_squared_error(y, las.predict(x))
r2_lasso = r2_score(y, las.predict(x))

result_df['LASSO'] = [mse_lasso, r2_lasso]
result_df
Out[8]:
OLS OLS - CV LASSO
MSE 21.9 37.2 22.7
0.741 0.56 0.732
In [9]:
scatterPlot(y, las.predict(x))

Comparação dos coeficientes

Abaixo podemos verificar que, ao penalizar os estimadores com um $\alpha > 0$, quase todos os coeficientes ficaram com valores menores em módulo. O coeficiente de AGE ficou igual a 0, ou seja, o algoritmo eliminou um parâmetro da regressão.

In [10]:
coef_df['LASSO'] = np.append([las.intercept_], las.coef_) 
coef_df.T
Out[10]:
Intercept CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT
OLS 22.5 -0.92 1.08 0.143 0.682 -2.06 2.67 0.0211 -3.1 2.66 -2.08 -2.06 0.857 -3.75
LASSO 22.5 -0.537 0.596 -0.0366 0.65 -1.43 2.87 -0 -2.21 0.755 -0.479 -1.88 0.74 -3.73

Validação Cruzada - LASSO

Faremos também um teste de validação cruzada, dividindo os dados em 5 grupos.

É possível verificar que o $R^2$ continua baixo para os dois últimos grupos, porém o $R^2$ geral da validação cruzada foi maior que o da validação cruzada no método de mínimos quadrados ordinários.

Concluímos que, apesar de termos um ajuste pior da amostra de treino utilizando LASSO, o modelo se comporta melhor nos dados fora da amostra.

In [11]:
kf = KFold(len(x), n_folds=5)
p = np.zeros_like(y)
for k, (train, test) in enumerate(kf):
    las.fit(x.ix[train], y[train])
    p[test] = las.predict(x.ix[test])
    print("[Fold {0}] \u03b1: {1:.2f}, R\u00b2: {2:.2f}".
          format(k+1, alpha, las.score(x.ix[test], y[test])))

mse_lasso_cv = mean_squared_error(y, p)
r2_lasso_cv = r2_score(y, p)

result_df['LASSO - CV'] = [mse_lasso_cv, r2_lasso_cv]
result_df
[Fold 1] α: 0.13, R²: 0.69
[Fold 2] α: 0.13, R²: 0.74
[Fold 3] α: 0.13, R²: 0.58
[Fold 4] α: 0.13, R²: 0.06
[Fold 5] α: 0.13, R²: -0.04
Out[11]:
OLS OLS - CV LASSO LASSO - CV
MSE 21.9 37.2 22.7 35.5
0.741 0.56 0.732 0.58
In [12]:
scatterPlot(y, p)

Caminho LASSO

Com a visualização do LASSO Path, podemos ver o comportamento de cada coeficiente da regressão de acordo com a variação do parâmetro de regularização $\alpha$.

Quando $\alpha$ é muito grande, a penalização é significativa e todos os coeficientes ficam muito próximos de 0. A medida que o $\alpha$ tende a 0, os coeficientes se estabilizam e tendem aos valores calculados pelo método dos mínimos quadrados ordinários.

In [13]:
las = Lasso()
alphas = np.logspace(-3, 1.5, 1000)
alphas, coefs, _ = las.path(x, y, alphas=alphas)

fig, ax = plt.subplots()
ax.plot(alphas, coefs.T)
ax.set_xscale('log')
ax.set_xlim(alphas.max(), alphas.min())
plt.vlines(alpha, coefs.min(), coefs.max())
plt.xlabel("Alpha")
plt.ylabel("Coeficientes")

plt.show()

LASSO com seleção automática do $\alpha$

Faremos uma implementação de validação cruzada para encontrar o melhor parâmetro de regularização $\alpha$ de acordo com os dados da amostra.

In [14]:
from sklearn.linear_model import LassoCV

las = LassoCV()
kf = KFold(len(x), n_folds=5)
p = np.zeros_like(y)

for k, (train, test) in enumerate(kf):
    las.fit(x.ix[train], y[train])
    p[test] = las.predict(x.ix[test])
    print("[Fold {0}] \u03b1: {1:.2f}, R\u00b2: {2:.2f}".
          format(k+1, las.alpha_, las.score(x.ix[test], y[test])))

mse_lasso_cv_2 = mean_squared_error(y, p)
r2_lasso_cv_2 = r2_score(y, p)

result_df['LASSO - CV (\u03b1 selection)'] = [mse_lasso_cv_2, r2_lasso_cv_2]
result_df
[Fold 1] α: 0.76, R²: 0.73
[Fold 2] α: 0.66, R²: 0.70
[Fold 3] α: 0.67, R²: 0.40
[Fold 4] α: 0.02, R²: 0.08
[Fold 5] α: 0.16, R²: -0.03
Out[14]:
OLS OLS - CV LASSO LASSO - CV LASSO - CV (α selection)
MSE 21.9 37.2 22.7 35.5 38.6
0.741 0.56 0.732 0.58 0.542
In [15]:
scatterPlot(y, p)

Regressão com termos de segunda ordem

Repetiremos as regressões adicionando o quadrado das variáveis explicativas ao modelo. Vamos verificar se o modelo de segunda ordem pode predizer melhor o preço dos imóveis.

In [16]:
col2 = ['{0}\u00b2'.format(c) for c in list(x.columns)]
df2 = x.copy()**2
df2.columns = col2

x2 = pd.concat([x, df2], axis=1)

lr = LinearRegression()
lr.fit(x2, y)

coef_df_2 = pd.DataFrame(index=np.append(['Intercept'], x2.columns))
coef_df_2['OLS X\u00b2'] = np.append([lr.intercept_], lr.coef_) 
coef_df_2
Out[16]:
OLS X²
Intercept 18.7
CRIM -2.96
ZN -0.988
INDUS 0.0311
CHAS 0.0536
NOX -2.62
RM 1.89
AGE 0.0627
DIS -3.24
RAD 3.83
TAX -2.23
PTRATIO -1.32
B -0.36
LSTAT -5.32
CRIM² 0.213
ZN² 0.397
INDUS² 0.331
CHAS² 0.182
NOX² -0.0441
RM² 0.859
AGE² 0.13
DIS² 0.599
RAD² -0.667
TAX² 0.427
PTRATIO² 0.58
-0.318
LSTAT² 1.15
In [17]:
mse_ls = mean_squared_error(y, lr.predict(x2))
r2_ls = r2_score(y, lr.predict(x2))

result_df['OLS X\u00b2'] = [mse_ls, r2_ls]
result_df
Out[17]:
OLS OLS - CV LASSO LASSO - CV LASSO - CV (α selection) OLS X²
MSE 21.9 37.2 22.7 35.5 38.6 14.3
0.741 0.56 0.732 0.58 0.542 0.831
In [18]:
scatterPlot(y, lr.predict(x2))

Conforme esperado, vemos que o modelo de segunda ordem se encaixa melhor. Porém, faremos mais a frente a validação cruzada para verificar se não está ocorrendo overfitting.

Regressão LASSO com termos de segunda ordem

Dessa vez escolhemos um valor maior para $\alpha$, visto que há agora um maior número de variáveis a serem eliminadas pelo modelo.

Observamos um maior número de coeficientes iguais a 0.

In [19]:
alpha = 0.20

las = Lasso(alpha=alpha)
las.fit(x2, y)

coef_df_2['LASSO X\u00b2'] = np.append([las.intercept_], las.coef_) 
coef_df_2
Out[19]:
OLS X² LASSO X²
Intercept 18.7 21.2
CRIM -2.96 -0
ZN -0.988 0
INDUS 0.0311 -0
CHAS 0.0536 0
NOX -2.62 -0
RM 1.89 2.23
AGE 0.0627 -0
DIS -3.24 -0.74
RAD 3.83 0
TAX -2.23 -0.216
PTRATIO -1.32 -1.39
B -0.36 0
LSTAT -5.32 -4.97
CRIM² 0.213 -0.108
ZN² 0.397 0.145
INDUS² 0.331 0
CHAS² 0.182 0.215
NOX² -0.0441 -0.619
RM² 0.859 1.03
AGE² 0.13 0
DIS² 0.599 -0
RAD² -0.667 -0
TAX² 0.427 0.053
PTRATIO² 0.58 0.0311
-0.318 -0.236
LSTAT² 1.15 0.808
In [20]:
mse_lasso = mean_squared_error(y, las.predict(x2))
r2_lasso = r2_score(y, las.predict(x2))

result_df['LASSO X\u00b2'] = [mse_lasso, r2_lasso]
result_df
Out[20]:
OLS OLS - CV LASSO LASSO - CV LASSO - CV (α selection) OLS X² LASSO X²
MSE 21.9 37.2 22.7 35.5 38.6 14.3 16.3
0.741 0.56 0.732 0.58 0.542 0.831 0.807
In [21]:
scatterPlot(y, las.predict(x2))

Validação Cruzada

Faremos agora o teste de validação cruzada para o modelo com termos de segunda ordem.

In [22]:
kf = KFold(len(x2), n_folds=5)
p = np.zeros_like(y)
for k, (train, test) in enumerate(kf):
    lr.fit(x2.ix[train], y[train])
    p[test] = lr.predict(x2.ix[test])
    print("[Fold {0}] R\u00b2: {1:.2f}".
          format(k+1, lr.score(x2.ix[test], y[test])))

mse_ls_cv = mean_squared_error(y, p)
r2_ls_cv = r2_score(y, p)

result_df['OLS X\u00b2 - CV'] = [mse_ls_cv, r2_ls_cv]

print('-------------------')

p = np.zeros_like(y)
for k, (train, test) in enumerate(kf):
    las.fit(x2.ix[train], y[train])
    p[test] = las.predict(x2.ix[test])
    print("[Fold {0}] \u03b1: {1:.2f}, R\u00b2: {2:.2f}".
          format(k+1, alpha, las.score(x2.ix[test], y[test])))

mse_lasso_cv = mean_squared_error(y, p)
r2_lasso_cv = r2_score(y, p)

result_df['LASSO X\u00b2 - CV'] = [mse_lasso_cv, r2_lasso_cv]

result_df
[Fold 1] R²: 0.68
[Fold 2] R²: 0.53
[Fold 3] R²: 0.82
[Fold 4] R²: 0.42
[Fold 5] R²: -0.25
-------------------
[Fold 1] α: 0.20, R²: 0.80
[Fold 2] α: 0.20, R²: 0.83
[Fold 3] α: 0.20, R²: 0.81
[Fold 4] α: 0.20, R²: 0.36
[Fold 5] α: 0.20, R²: 0.50
Out[22]:
OLS OLS - CV LASSO LASSO - CV LASSO - CV (α selection) OLS X² LASSO X² OLS X² - CV LASSO X² - CV
MSE 21.9 37.2 22.7 35.5 38.6 14.3 16.3 30.6 21.4
0.741 0.56 0.732 0.58 0.542 0.831 0.807 0.638 0.746

Conclusão

Observamos na tabela anterior que nos testes de validação cruzada a regressão LASSO teve um desempenho superior em relação ao método dos mínimos quadrados ordinários. Podemos concluir, portanto, que a seleção de variáveis feita pela regularização melhorou o modelo para previsões (dados fora da amostra).