# LARS vs. OLS (MMST - Ex. 5.18)

+1 voto
116 visitas

Considere a questão 5.18 do livro Modern multivariate
statistical techniques - Alan Julian Izenman que versa sobre o uso
da base de dados "Boston housing data"e variable selection
o seu material e implementacão.

website lib.stat.cmu.edu/datasets/boston corrected.txt. There are
506 observations on census tracts in the Boston Standard Metropolitan
Statistical Area (SMSA) in 1970. The response variable is the logarithm of
the median value of owner-occupied homes in thousands of dollars; there
are 13 input variables (plus information on location of each observation).
Compute the OLS estimates and compare them with those obtained from
the following variable-selection algorithms: Forwards Selection (stepwise),
Cp, the Lasso, LARS, and Forwards Stagewise.

(c) Como implementar o algoritmo de selecão de variáveis em
Python conhecido como LARS em python? Compare esses
acima.

## 2 Respostas

+1 voto
respondida Jun 21, 2018 por (46 pontos)

Melhor resposta

Lista de Variáveis

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


from sklearn.datasets import load_boston
import pandas as pd
import numpy as np

print(boston.DESCR)

bh = pd.DataFrame(boston.data)
bh.columns = boston.feature_names
bh['PRICE'] = boston.target

CRIM    ZN  INDUS  CHAS    NOX     RM   AGE     DIS  RAD    TAX  \
0  0.00632  18.0   2.31   0.0  0.538  6.575  65.2  4.0900  1.0  296.0
1  0.02731   0.0   7.07   0.0  0.469  6.421  78.9  4.9671  2.0  242.0
2  0.02729   0.0   7.07   0.0  0.469  7.185  61.1  4.9671  2.0  242.0
3  0.03237   0.0   2.18   0.0  0.458  6.998  45.8  6.0622  3.0  222.0
4  0.06905   0.0   2.18   0.0  0.458  7.147  54.2  6.0622  3.0  222.0

PTRATIO       B  LSTAT  PRICE
0     15.3  396.90   4.98   24.0
1     17.8  396.90   9.14   21.6
2     17.8  392.83   4.03   34.7
3     18.7  394.63   2.94   33.4
4     18.7  396.90   5.33   36.2


OLS:

   if __name__ == '__main__':

import statsmodels.api as sm
import patsy as ps

y,X = ps.dmatrices('PRICE ~ CRIM + ZN + INDUS + CHAS + NOX + RM + AGE + DIS + RAD + TAX + PTRATIO + B + LSTAT',data=bh, return_type='dataframe')
model = sm.OLS(y,X)
results = model.fit()
print (results.summary())


                            OLS Regression Results
==============================================================================
Dep. Variable:                  PRICE   R-squared:                       0.741
Method:                 Least Squares   F-statistic:                     108.1
Date:                Mon, 25 Jun 2018   Prob (F-statistic):          6.95e-135
Time:                        20:41:52   Log-Likelihood:                -1498.8
No. Observations:                 506   AIC:                             3026.
Df Residuals:                     492   BIC:                             3085.
Df Model:                          13
Covariance Type:            nonrobust
==============================================================================
coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     36.4911      5.104      7.149      0.000      26.462      46.520
CRIM          -0.1072      0.033     -3.276      0.001      -0.171      -0.043
ZN             0.0464      0.014      3.380      0.001       0.019       0.073
INDUS          0.0209      0.061      0.339      0.735      -0.100       0.142
CHAS           2.6886      0.862      3.120      0.002       0.996       4.381
NOX          -17.7958      3.821     -4.658      0.000     -25.302     -10.289
RM             3.8048      0.418      9.102      0.000       2.983       4.626
AGE            0.0008      0.013      0.057      0.955      -0.025       0.027
DIS           -1.4758      0.199     -7.398      0.000      -1.868      -1.084
RAD            0.3057      0.066      4.608      0.000       0.175       0.436
TAX           -0.0123      0.004     -3.278      0.001      -0.020      -0.005
PTRATIO       -0.9535      0.131     -7.287      0.000      -1.211      -0.696
B              0.0094      0.003      3.500      0.001       0.004       0.015
LSTAT         -0.5255      0.051    -10.366      0.000      -0.625      -0.426
==============================================================================
Omnibus:                      178.029   Durbin-Watson:                   1.078
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              782.015
Skew:                           1.521   Prob(JB):                    1.54e-170
Kurtosis:                       8.276   Cond. No.                     1.51e+04
==============================================================================


Calcular o MSE:

from sklearn.linear_model import LinearRegression
lr = LinearRegression()
lr.fit(X, y)

mseFull = np.mean((bh['PRICE'] - lr.predict(X))**2)
print(mseFull)


21.897779217687486


LARS:

#1) Lars normal
reg = linear_model.Lars(n_nonzero_coefs=500)
reg.fit(X, y)
print(reg.coef_)

#2) Lasso Lars
reg2 = linear_model.LassoLars(alpha=0.01)
reg2.fit(X,y)
print(reg2.coef_)

#3) Lasso Lars com Cross Validation
reg3 = linear_model.LassoLarsCV(cv=10)
reg3.fit(X,y)
print(reg3.coef_)


 #1) Lasso normal

[ 0.00000000e+00 -1.09250387e-01  4.78372644e-02  3.57286797e-02
2.67472193e+00 -1.83209467e+01  3.79660038e+00  1.06365664e-03
-1.48741200e+00  3.21253220e-01 -1.31885495e-02 -9.61744609e-01
9.47486584e-03 -5.26838797e-01]

#2) Lasso Lars

[ 0.00000000e+00 -3.58135035e-02  1.30928584e-02  0.00000000e+00
2.35445155e+00 -8.56012099e+00  4.23361219e+00  0.00000000e+00
-7.43136732e-01  0.00000000e+00  0.00000000e+00 -8.19078277e-01
7.27593433e-03 -5.21059152e-01]

#3) Lasso Lars com Cross Validation

[ 0.00000000e+00 -3.98460117e-02  1.66050711e-02 -4.45454004e-03
2.44200356e+00 -9.47813285e+00  4.21500719e+00  0.00000000e+00
-8.33093415e-01  2.70507003e-03  0.00000000e+00 -8.22840000e-01
7.42443460e-03 -5.21319026e-01]


Calcular o MSE:

#1) Lars normal
mseLARS = np.mean((bh['PRICE'] - reg.predict(X))**2)
print(mseLARS)
print(reg.score(X,y)) #R2

#2) Lasso Lars
mseLARS2 = np.mean((bh['PRICE'] - reg2.predict(X))**2)
print(mseLARS2)
print(reg2.score(X,y))  #R2

#3) Lasso Lars com Cross Validation
mseLARS3 = np.mean((bh['PRICE'] - reg3.predict(X))**2)
print(mseLARS3)
print(reg3.score(X,y)) #R2


#1) Lasso normal
21.90246113588432
0.7405522827510782

#2)   Lasso Lars
23.678624594592755
0.7195125670787672

#3) Lasso Lars com Cross Validation
23.469597857245414
0.7219886134697318


Entre os modelos de Lars, o Lars normal apresentou uma pequena vantagem nos resultados em relação ao Lasso Lars e ao Lasso Lars com Cross Validations. O modelo em questão apresentou um menor EQM (21.90) e um maior R2 (0.74) em relação aos outros modelos.

Porém, o Lars normal e o OLS apresentaram resultados parecidos (assim como os outros modelos), visto que o R2 do OLS foi de 0.74 e seu EQM de 21.89.

Ainda, como o Felipe observa abaixo, os resultados do OLS para um processo de validação "out-of-sample" foram os piores (pior R2 e maior EQM) em relação aos modelos Lars. O que reforça o poder de previsão dos modelos Lars em relação ao OLS.

comentou Jun 26, 2018 por (76 pontos)
Pedro, boa implementação.
Uma dúvida: X2 seria a matriz com todas as colunas da base de dados menos a coluna PRICE?

Um passo adicional seria realizar um processo de validação cruzada, ou também comparando a performance do modelo em dados "fora da amostra".
Postarei a implementação em forma de resposta, a fim de proporcionar maior clareza.
comentou Jun 26, 2018 por (46 pontos)
X2 era a matriz X. Já arrumei o código.

0 votos
respondida Jun 26, 2018 por (76 pontos)

Conforme discutido no comentário da resposta do Pedro, segue implementação da validação "out-of-sample" dos dois modelos - note que os comentários estão contidos junto ao código:

    ### 1) IMPORTAÇÃO DOS DADOS ###

print(boston.DESCR)

import pandas as pd

bh = pd.DataFrame(boston.data)
bh.columns = boston.feature_names
bh['PRICE'] = boston.target

### 2) OLS ###

from sklearn.linear_model import LinearRegression
import numpy as np

lr = LinearRegression()

x = bh.drop('PRICE', axis=1)
y = bh['PRICE']

lr.fit(x, y)

print('Coefficients: \n', lr.coef_)

mseOLS = np.mean((bh['PRICE'] - lr.predict(x))**2)
R2OLS = lr.score(x,y)
print(mseOLS) ## MSE do modelo OLS ##
print(R2OLS)  ## R² do modelo OLS ##

### 3) LARS ###

import time
import matplotlib.pyplot as plt
from sklearn.linear_model import LassoLarsCV
from sklearn import linear_model

## Computing regularization path using the Lars lasso... ##
t1 = time.time()
model = LassoLarsCV(cv=10).fit(x, y)
t_lasso_lars_cv = time.time() - t1

# Display results
m_log_alphas = -np.log10(model.cv_alphas_)

plt.figure()
plt.plot(m_log_alphas, model.mse_path_, ':')
plt.plot(m_log_alphas, model.mse_path_.mean(axis=-1), 'k',
label='Average across the folds', linewidth=2)
plt.axvline(-np.log10(model.alpha_), linestyle='--', color='k',
label='alpha CV')
plt.legend()

plt.xlabel('-log(alpha)')
plt.ylabel('Mean square error')
plt.title('Mean square error on each fold: Lars (train time: %.2fs)'
% t_lasso_lars_cv)
plt.axis('tight')
plt.ylim(ymin, ymax)

plt.show()

## Computing LARS w/ alpha selected via CV (alpha = 0,01)... ##

reg = linear_model.LassoLars(alpha=0.01)
reg.fit(x, y)
print('Coefficients: \n', reg.coef_)

mseLARS = np.mean((bh['PRICE'] - reg.predict(x))**2)
R2LARS = reg.score(x,y)

print(mseLARS) ## MSE do modelo LARS ##
print(R2LARS) ## R² do modelo LARS ##

### 4) CROSS-VALIDATION OLS ###

from sklearn.model_selection import train_test_split

X_train, X_test, Y_train, Y_test = train_test_split(x, bh.PRICE, test_size=0.33,
random_state=5)

print (X_train.shape)
print (X_test.shape)
print (Y_train.shape)
print (Y_test.shape)

lm = LinearRegression()
lm.fit(X_train, Y_train)

print("Fit a model X_train, and calculate MSE with Y_train:", np.mean((Y_train -
lm.predict(X_train))**2))
print("Fit a model X_train, and calculate MSE with X_test, Y_test:", np.mean((Y_test -
lm.predict(X_test)) ** 2) )

print(lm.score(X_test, Y_test)) ### R² do OLS a partir dos dados 'fora da amostra' ###

import matplotlib.pyplot as plt

plt.scatter(lm.predict(X_train), lm.predict(X_train) - Y_train, c = 'b', s = 40, alpha = 0.5)
plt.scatter(lm.predict(X_test), lm.predict(X_test) - Y_test, c = 'g', s = 40)
plt.hlines(y = 0, xmin=0,xmax=50)
plt.title('Plotagem dos resíduos a partir dos dados de treinamento (azul) e teste
(verde)')
plt.ylabel('Resíduos')
plt.show()

### 5) CROSS VALIDATION LARS ###

x2 = bh.drop(['INDUS', 'AGE', 'RAD', 'TAX'], axis = 1)

X_train2, X_test2, Y_train2, Y_test2 = train_test_split(x2, bh.PRICE, test_size=0.33,
random_state=5)

print (X_train2.shape)
print (X_test2.shape)
print (Y_train2.shape)
print (Y_test2.shape)

lm2 = LinearRegression()
lm2.fit(X_train2, Y_train2)

print("Fit a model X_train2, and calculate MSE with Y_train2:", np.mean((Y_train2 -
lm2.predict(X_train2))**2))
print("Fit a model X_train2, and calculate MSE with X_test2, Y_test2:",
np.mean((Y_test2 - lm2.predict(X_test2)) ** 2) )

print(lm2.score(X_test2, Y_test2)) ### R² do LARS a partir dos dados 'fora da amostra'
###
plt.scatter(lm2.predict(X_train2), lm2.predict(X_train2) - Y_train2, c = 'b', s = 40, alpha =
0.5)
plt.scatter(lm2.predict(X_test2), lm2.predict(X_test2) - Y_test2, c = 'g', s = 40)
plt.hlines(y = 0, xmin=0,xmax=50)
plt.title('Plotagem dos resíduos a partir dos dados de treinamento (azul) e teste
(verde)')
plt.ylabel('Resíduos')
plt.show()


Após as observações da Júlia, corrigi meu código.
Note que o LARS e o OLS produziram resultados bem parecidos sobre os dados da amostra: MSE (LARS) = 23.69, MSE (OLS) = 21.9; R² (LARS) = 0.72, R² (OLS) = 0.74.

Observando os resultados num processo de validação "out-of-sample", notamos que o modelo OLS apresenta menor R² (0.69) e maior MSE (28.54), enquanto o modelo LARS apresenta um MSE menor (quase nulo) e um R² maior (1). Evidenciando o maior poder preditivo para este segundo modelo.

Por meio da plotagem dos resíduos a partir dos dados de treinamento e de teste dos dois modelos, notamos no caso do LARS que as variáveis excluídas no processo de regularização melhoraram seu potencial preditivo - observe que os resíduos se concentram em torno do zero e não há nenhuma aparente estrutura.

OLS:

LARS:

LARS CV 10-fold Path:

comentou Jul 4, 2018 por (41 pontos)
Por que forçar o LARS a escolher apenas um coeficiente?

(em reg = linear_model.Lars(n_nonzero_coefs=1)

acho que tinha que deixar esse parâmetro no default (que é 500) para que o LARS selecione tantos coeficientes quanto necessários.

Além disso, acho que seria legal usar o LassoLars (reg = linear_model.LassoLars(alpha=.1)) porque daí dá para escolher o alpha. Ou melhor ainda, usar o LarsLassoCV porque daí ele faz o gráfico de MSE por alpha e escolhe o alpha que gera o menor MSE.
comentou Jul 5, 2018 por (76 pontos)