Como implementar o algoritmo LASSO em Python e compará-lo ao OLS?

Question

Considere a questão 5.18 do livro Modern Multivariate Statistical Techniques de Alan Julian Izenman, que versa sobre o uso da base de dados "Boston House Prices" e de algoritmos de seleção de variáveis.
Interprete os resultados e documente cuidadosamente o seu material e Implementação.
Compare esses resultados com as estimativas de OLS usando os dados descritos acima.

Answer 1

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