O código em Python é apresentado abaixo com os comentários.
Para uma melhor visualização do estudo, sugiro acessar o arquivo Housing.ipynb
# coding: utf-8
# # 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](http://scikit-learn.org/) 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)
# 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()
# ## 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
# 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
# ... e o gráfico com a diferença entre o preço real e o calculado.
# In[5]:
import matplotlib.pyplot as plt
#get_ipython().magic('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
# É 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
# 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
# ## 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
# 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
# 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
# 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
# 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
# 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
# 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
# ## 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).