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.
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.
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)
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()
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 |
A regressão linear será realizada utilizando o método dos mínimos quadrados. Apresentamos abaixo os coeficientes calculados.
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
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$ ...
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
OLS | |
---|---|
MSE | 21.9 |
R² | 0.741 |
... e o gráfico com a diferença entre o preço real e o calculado.
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))
É 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.
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
OLS | OLS - CV | |
---|---|---|
MSE | 21.9 | 37.2 |
R² | 0.741 | 0.56 |
É possível verificar que o $R^2$ foi baixo para os dois últimos grupos.
scatterPlot(y, p)
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$.
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
OLS | OLS - CV | LASSO | |
---|---|---|---|
MSE | 21.9 | 37.2 | 22.7 |
R² | 0.741 | 0.56 | 0.732 |
scatterPlot(y, las.predict(x))
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.
coef_df['LASSO'] = np.append([las.intercept_], las.coef_)
coef_df.T
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 |
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.
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
OLS | OLS - CV | LASSO | LASSO - CV | |
---|---|---|---|---|
MSE | 21.9 | 37.2 | 22.7 | 35.5 |
R² | 0.741 | 0.56 | 0.732 | 0.58 |
scatterPlot(y, p)
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.
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()
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.
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
OLS | OLS - CV | LASSO | LASSO - CV | LASSO - CV (α selection) | |
---|---|---|---|---|---|
MSE | 21.9 | 37.2 | 22.7 | 35.5 | 38.6 |
R² | 0.741 | 0.56 | 0.732 | 0.58 | 0.542 |
scatterPlot(y, p)
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.
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
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 |
B² | -0.318 |
LSTAT² | 1.15 |
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
OLS | OLS - CV | LASSO | LASSO - CV | LASSO - CV (α selection) | OLS X² | |
---|---|---|---|---|---|---|
MSE | 21.9 | 37.2 | 22.7 | 35.5 | 38.6 | 14.3 |
R² | 0.741 | 0.56 | 0.732 | 0.58 | 0.542 | 0.831 |
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.
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.
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
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 |
B² | -0.318 | -0.236 |
LSTAT² | 1.15 | 0.808 |
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
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 |
R² | 0.741 | 0.56 | 0.732 | 0.58 | 0.542 | 0.831 | 0.807 |
scatterPlot(y, las.predict(x2))
Faremos agora o teste de validação cruzada para o modelo com termos de segunda ordem.
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
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 |
R² | 0.741 | 0.56 | 0.732 | 0.58 | 0.542 | 0.831 | 0.807 | 0.638 | 0.746 |
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).