Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

LARS vs. OLS (MMST - Ex. 5.18)

+1 voto
62 visitas
perguntada Jun 21 em Ciência da Computação por Pedro Campelo (41 pontos)  
editado Jul 6 por Pedro Campelo

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
algorithms. Interprete os resultados e documente cuidadosamente
o seu material e implementacão.

5.18 The Boston housing data can be downloaded from the StatLib
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
resultados com as estimativas de OLS usando os dados descritos
acima.

Compartilhe

2 Respostas

+1 voto
respondida Jun 21 por Pedro Campelo (41 pontos)  
selecionada Jul 5 por Pedro Campelo
 
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
              - 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

Importar a base de dados:

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

boston = load_boston()
print(boston.DESCR)

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

print(bh.head())

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

Resultado:

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  PRICE   R-squared:                       0.741
Model:                            OLS   Adj. R-squared:                  0.734
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)

Resultado:

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

Resultados:

 #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

Resultados:

#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 por Felipe Carneiro (26 pontos)  
editado Jun 26 por Felipe Carneiro
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 por Pedro Campelo (41 pontos)  
X2 era a matriz X. Já arrumei o código.

Obrigado, Felipe.
0 votos
respondida Jun 26 por Felipe Carneiro (26 pontos)  
editado Jul 5 por Felipe Carneiro

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

    from sklearn.datasets import load_boston

    boston = load_boston()
    print(boston.DESCR)

    import pandas as pd

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

    print(bh.head())

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

A imagem será apresentada aqui.

LARS:

A imagem será apresentada aqui.

LARS CV 10-fold Path:

A imagem será apresentada aqui.

comentou Jul 4 por Julia Regina Scotti (31 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 por Felipe Carneiro (26 pontos)  
Júlia, bem notado.
Modifiquei meu código e, de fato, os resultados agora parecem fazer mais sentido!
comentou Jul 5 por Pedro Campelo (41 pontos)  
Obrigado pelo comentário, Julia. Não tinha notado esse erro e já arrumei a questão para os 3 modelos citados do Lars.
...