Um walkthrough com o MARS

Postado originalmente no How Movile Works.

Esse post tem sido adiado a no mínimo a uns 2 anos, desde que eu estava em alguns estudos para finalizar a minha dissertação de mestrado, mas como acredito que o conteúdo dele é atemporal, e por isso vou falar um pouquinho sobre essa técnica MARS aqui no How Movile Works.

Contudo, mesmo com todo o material que está na Web, há algumas coisas a serem ditas sobre  Multivariate Adaptive Regression Splines, ou simplesmente MARS, dado que essa é uma das técnicas mais subestimadas e desconhecidas de todo universo de algoritmos regressores de Machine Learning.

Uma das primeiras interações que eu tive com MARS foi em um dos testes com o software da Salford Systems, que detém a patente sobre essa técnica. MARS é uma forma de regressão criada em meados da década de 90 por Jerome Friedman.

Eu poderia fazer várias descrições, mas o sumário do artigo original já é matador nesse sentido:

A new method is presented for flexible regression modeling of high dimensional data. The model takes the form of an expansion in product spline basis functions, where the number of basis functions as well as the parameters associated with each one (product degree and knot locations) are automatically determined by the data. This procedure is motivated by the recursive partitioning approach to regression and shares its attractive properties. Unlike recursive partitioning, however, this method produces continuous models with continuous derivatives. It has more power and flexibility to model relationships that are nearly additive or involve interactions in at most a few variables. In addition, the model can be represented in a form that separately identifies the additive contributions and those associated with the different multivariable interactions

A primeira coisa que me chamou atenção nessa técnica de Machine Learning na época era devido ao fato de que a MARS foi desenhada para atacar dois dos maiores problemas que eu enfrentei, na época trabalhando no mercado financeiro, que eram: (i) bases de dados com alta dimensionalidade e (ii) não linearidades nas dinâmicas de recuperação de Non-Performing Loans.

Exemplificando o primeiro caso, eu trabalhava com uma base de aproximadamente 8 milhões de registros, que tinha mais de 400 colunas, e aí tem-se o que é de praxe: inúmeras variáveis com missing data, multicolinearidadeheterocedasticidade, junk data, e afins. Em resumo uma tremenda bagunça.

Já o segundo caso tratar de dados de dinâmicas de Non-Performing Loans (NPL) envolve diversas características comportamentais que podem ser caracterizadas através dos dados como a influência de um telefone na recuperação o impacto financeiro no envio de uma carta para o devedor, a influência do parcelamento em relação ao fluxo de recebimento a influência do valor de face da dívida na recuperação, et cetera.

Em termos computationais/algoritmicos o MARS é uma extensão das Classifications and Regresson Trees em que o split ao invés de ser no atributo ocorre no próprio conjunto de dados. Em outras palavras, as divisões ocorrem dentro do range de dados do atributo.

Essa divisão ocorre através do fitting de uma regressão linear em que essa divisão (split) é feito usando funções de base radial (RBF) em que cada uma dessas funções ocorrem em pares f(x) = {x − t SE x > t, SENAO 0} e g(x)= {t − x SE x < t, SENAO 0} em que o valor de t é chamado de knot.

Questões de negócios e teóricas colocadas, vamos para o código.

Em primeiro lugar vamos instalar os pacotes caret e earth (sim, devido a questões de copyright o MARS no R é conhecido como Earth).

install.packages("caret")
install.packages("earth")

Pacotes instalados, vamos carregar as bibliotecas:

library(caret)
library(earth)

Para esse pequeno experimento, vamos usar uma base de dados de alguns imóveis da região de Tamboré aqui no estado de São Paulo de 2013. Essa base de dados contém alguns anúncios do Zap Imóveis do ano de 2013 de alguns imóveis como apartamentos e casas em uma região nobre de São Paulo.

Fiz o scrapping dessa base na época para saber quais eram os fatores mais determinantes nos preços de casas nessa região da cidade.

A conclusão usando estatística descritiva simples e um pouquinho de Data Mining foi que a dinâmica de preços desses imóveis segue uma proposta quase que irracional em alguns momentos, mas esse não é o objetivo desse post, em que vamos usar essa base apenas para fins didáticos.

Primeiramente vamos realizar o upload dos dados.

# Get data
tambore_parnaiba_2013_data <- read.csv('https://raw.githubusercontent.com/fclesio/learning-space/master/Datasets/imoveis-tambore-parnaiba-2013.csv');

Para ver se deu tudo certo, vamos dar uma olhada no resultset gerado:

# Check data
tambore_parnaiba_2013_data

Antes do próximo comando duas observações com MARS: a) esse algoritmo não trabalha com variáveis categóricas (i.e. o ideal é criar variáveis dummy), e b) esse algoritmo não trabalha com missing values (i.e. por motivos óbvios não tem como realizar nenhuma regressão com algum valor que não existe.)

Tendo isso em mente, vamos retirar a variável bairro da base de dados:

# Eliminate Bairro variable
tambore_parnaiba_2013_data$Bairro <- NULL

Com a base de dados higienizada por assim dizer, vamos para a parte de implementação em si em que vamos submeter os dados para o ajuste do modelo. Em outras palavras: Vamos realizar o fitting do modelo.

Nesse próximo snippet vamos chamar a função que vai realizar o ajuste do modelo:

# Fit model(MARS)
fit <- earth(ValorVenda~.
             ,tambore_parnaiba_2013_data # Base de dados
             ,trace = 1                  # Trace: Detalhes do runtime do algoritmo
             ,ncross= 50                 # Número de validações cruzadas
             ,nfold = 10                 # Número de folds para cross validation
             ,pmethod="backward"         # Nesse caso o prunning será da última variável do modelo para a primeira
             ,nprune=15                  # Número de termos máximos após o prunning
             ,varmod.method="lm")        #Variance Model: Opção para a escolha do modelo de variancia (Nesse caso será o modelo linear, mas o glm pode ser escolhido.)

Algumas explicações antes de prosseguirmos:
– O número de validações cruzadas (Cross Validation e número de n folds é mais alquimia do que química: Isso é, vai de acordo com a heurística de trabalho de quem está parametrizando o modelo e tentativa e erro;
– Vamos usar backward devido ao fato de que é mais confortável saber o resultado com todas as variáveis dispostas no modelo e depois realizar o prunning do que o contrário. E além do mais como o conjunto de dados não é grande não haverá complicações em termos computacionais; e
– A aproximação via modelo linear usando MARS será escolhido apenas para fins de exemplo.

Agora vamos ver um resumo do modelo.

# Summarize the fit
summary(fit,digit=3)

Vamos quebrar cada parte do output descrevendo alguns dos parâmetros da saída do console.

               coefficients
(Intercept)           596993
Playground             55849
h(ValorM2-5248)          236
h(3-Suite)            -28711
h(132-AreaUtil)        -4118
h(AreaUtil-132)         2802
h(AreaUtil-158)         4014

Logo de cara podemos ver que temos quase 596,993 de valor de intercepto, que “significa” que se todos os outros atributos forem zero o valor previsto para um imóvel será de R$ 596.993. Tentar ao menos interpretar o intercepto é algo que não faz sentido algum mas coloquei essa exemplificação apenas para ficar mais didática a análise.

Os demais valores são as funções de base radial que são geradas pelo algoritmo MARS, mais os seus respectivos coeficientes.

Selected 7 of 9 terms, and 4 of 17 predictors
Termination condition: RSq changed by less than 0.001 at 9 terms
Importance: AreaUtil, ValorM2, Playground, Suite, Entrada-unused, Condominio-unused, ...
Number of terms at each degree of interaction: 1 6 (additive model)
GCV 4074691053  RSS 1.11565e+11  GRSq 0.9885132  RSq 0.9934504  CVRSq 0.9055663

Note: the cross-validation sd's below are standard deviations across folds

Cross validation:   nterms 5.73 sd 0.97    nvars 3.06 sd 0.77

     CVRSq    sd     MaxErr     sd
     0.906 0.157     650028 216758

varmod: method "lm"    min.sd 7.46e+03    iter.rsq 0.148

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)    46481.318      9082.7           20
ValorVenda         0.028        0.01           35

                              mean   smallest    largest      ratio
95% prediction interval   292420.1   223271.9   481074.8   2.154659

                                         68%     80%     90%     95%
response values in prediction interval   92     100     100     100

Em termos gerais o que esse bloco de código acima diz é que o modelo regressor selecionou 7 termos (conjunto de constantes, coeficientes e intercepto) e das 17 variáveis independentes foram escolhidas 4 (i.e. AreaUtil, ValorM2, Playground, Suite).

Na sequência ele diz que o critério de parada foi que o coeficiente de determinação mudou menos de 0.001 ao longo da cadeia de processamento do algoritmo.

Depois ele diz que a ordem da importância de cada variável é AreaUtil, ValorM2, Playground, Suite, Entrada-unused, e Condominio-unused.

Na sequência ele apresenta algumas estatísticas de Generalized Cross Validation (GCV), Residual sum-of-squares (RSS – Soma dos resíduos quadráticos), e o mais importante que é o Coeficiente de Determinação (RSq). O RSq nesse modelo está em 99,35%; o que significa em termos gerais que 99,35% do Valor de Venda está explicado pelas variáveis constantes no modelo.

É um bom indicativo de fitting (ajuste) do modelo, mas vamos ver posteriormente que nem sempre ter um coeficiente de determinação alto significa bons resultados práticos.

Dito isso, agora vamos ver alguns dos plots do earth.

Em primeiro lugar vamos a função de plot chamada “plotmo”. Essa função tem como principal objetivo mostrar a modificação da variável preditora se todas as outras variáveis permanecerem constantes. Vejamos como é o comportamento da variável dependente de acordo com a modificação de cada uma das variáveis independentes:

# Plot
plotmo(fit)

Screen Shot 2016-05-22 at 9.23.46 AM

Como podemos ver a variável área útil tem o maior impacto sobre a o preço de venda (variável dependente), seguido da variável valorM2. Quando pensamos na dinâmica do mercado imobiliário isso faz total sentido, dado que o valor do metro quadrado já incorpora financeiramente uma série de atributos do local em que o imóvel está sendo vendido como infraestrutura, transporte, saneamento básico, serviços como supermercados, postos de saúde, e amenidades em geral.

Outro fato interessante que podemos ver nesse gráfico é que as variáveis playground e suite tem uma baixa relevância no modelo como um todo, i.e. mesmo se o imóvel tiver 1 playground e 4 suites, o que vai determinar o valor final será majoritariamente a metragem e o valor do metro na região.

Para comprovar essas conjecturas citadas anteriormente de forma objetiva, o MARS tem uma funcionalidade de apresentar a importância de cada variável no modelo. Essa funcionalidade pode ser apresentada através da funcão evimp.

# Variable importance
evimp(fit)

Esses são os resultados:

          nsubsets   gcv    rss
AreaUtil          6 100.0  100.0
ValorM2           5  31.9   31.2
Playground        3   4.0    5.9
Suite             1   0.7    2.8

Como podemos ver, as variáveis AreaUtil e ValorM2 são as mais importantes para determinar o preço de venda de um apartamento.

Agora que o modelo está ajustado, e entendemos o que ele faz, vamos testar o poder preditivo de fato. Para isso vamos usar a função predict. Para fins de evitar overfitting, vamos tirar a variável ValorVenda do modelo para ter uma visão mais realista da aderência do modelo.

# New dataset called "data"
data <- tambore_parnaiba_2013_data

# Eliminate variables
data$ValorVenda <- NULL
data$Bairro <- NULL

# Predictions (Fitting data over the new dataset)
predictions <- predict(fit, data)

# Change the name of variable for binding
colnames(predictions) <- c("PredictValorVenda")

# See results
predictions

O que fizemos aqui foi ajustar uma os dados do dataset data e após isso eliminamos as variáveis descenessárias (nesse caso bairro que não tem nenhum poder perditivo, e ValorVenda para não causar overfitting). Na sequência usamos a propriedade predict para aplicar o modelo fit na base de dados data e um pequeno ajuste de colunas para colocar o nome da variável.

Agora vamos ver os resultados vis-a-vis do resultado da previsão com os dados originais.

# Binding predictions with original data
predit <- cbind(tambore_parnaiba_2013_data,predictions)

# Only the two variables
predit <- cbind(predit$ValorVenda, predit$PredictValorVenda) 

# Name ajustments
colnames(predit) <- c("ValorVenda","PredictValorVenda")

# Convert to a dataframe
FinalPredictions <-as.data.frame(predit)

Explicando o snippet acima: Usamos a função cbind para concatenar dois resultsets diferentes, e após isso reutilizamos o objeto predict que recebeu apenas as colunas de interesse. Realizamos um pequeno ajuste nos nomes das variáveis e transformamos o objeto predict em um dataframe para cálculos posteriores.

Finalmente vamos ver de forma objetiva a diferença entre o valor de venda e o valor previsto pelo modelo:

# Final Predictions
FinalPredictions$resuduals <- FinalPredictions$ValorVenda - FinalPredictions$PredictValorVenda

# Final table with results
FinalPredictions

E esse é o output do objeto FinalPredictions:

   ValorVenda PredictValorVenda    resuduals
1      400000          377165.5   22834.5395
2      412000          381283.5   30716.5381
3      424000          499672.9  -75672.8656
4      490000          438935.5   51064.5182
5      540000          519845.6   20154.3738
6      545000          540044.9    4955.1494
7      550000          533534.8   16465.2286
8      550000          533534.8   16465.2286
9      555000          595419.5  -40419.5357
10     560000          624130.8  -64130.8026
11     570000          539570.5   30429.4529
12     575500          595419.5  -19919.5357
13     585000          569056.5   15943.5033
14     595000          595419.5    -419.5357
15     599000          558401.5   40598.4771
16     619000          600189.4   18810.6230
17     650000          616607.3   33392.7016
18     660000          616607.3   43392.7016
19     664000          725694.9  -61694.8772
20     680000          640252.0   39747.9673
21     690000          719742.4  -29742.3607
22     720000          725694.9   -5694.8772
23     740000          748934.3   -8934.3359
24     740000          804783.3  -64783.3245
25     750000          789256.9  -39256.9344
26     750000          789256.9  -39256.9344
27     750000          768293.7  -18293.6808
28     750000          707397.5   42602.4565
29     779000          724713.0   54287.0472
30     780000          729681.6   50318.3761
31     800000          836143.0  -36143.0484
32     960000         1025613.3  -65613.2683
33    1100000         1043647.7   56352.2635
34    1150000         1074648.6   75351.4097
35    1200000         1279797.1  -79797.1005
36    1200000         1223948.1  -23948.1119
37    1380000         1346151.6   33848.4030
38    1400000         1424162.7  -24162.7348
39    1400000         1367103.0   32896.9721
40    1450000         1427100.0   22899.9506
41    1590000         1528703.1   61296.9004
42    1690000         1787723.5  -97723.5283
43    1691000         1695382.3   -4382.3368
44    1700000         1702937.2   -2937.2031
45    1700000         1736685.9  -36685.8630
46    1850000         1960909.0 -110909.0338
47    1900000         1941890.0  -41890.0183
48    2277000         2238520.2   38479.7998
49    2650000         2616112.0   33888.0410
50    2850000         2744780.8  105219.2248

Bem, alguns de vocês estão pensando “Nossa mas esse modelo não teve o valor de venda explicado em 99,35% pelos dados constantes no modelo na fase de treinamento?”. Pois bem, aqui vai uma lição sobre modelagem: Ajuste nem sempre significa aderência.

Colocando em outros termos, o ajuste do modelo e o uso de métricas com o coeficiente de determinação deve ser feito muito mais para balizar a modelagem (e.g. inclusão ou exclusão de variáveis, ajuste de escala de variáveis, normalização de dados, et cetera) do que como uma potencial expectativa de sucesso do modelo.

Dito isso, vamos calcular o Root Mean Square Error (RMSE, ou Raiz da média do erro quadrático) que é a raiz quadrada da média da raiz de todos os erros (i.e. diferença entre a variável a ser prevista (ValorVenda) entre o que foi previsto de fato(predicted)). Ou para quem se sente mais confortável com menos palavras e mais demonstrações matemáticas é essa formula abaixo:

RMSE

Colocando essa representação no R, temos:

# Root Mean Square Error
rmse <- sqrt(mean(((FinalPredictions$ValorVenda - FinalPredictions$PredictValorVenda)^2)))

# Show results
rmse

Como resultado final temos o RMSE de 47236.65. O que isso significa isoladamente? Absolutamente nada. Mesmo que essa medida de avaliação de modelos seja do tipo “menor melhor” em termos práticos essa é uma medida que serve muito mais para termos de comparação com outros modelos do que um número definitivo, até porque vendo os resultados da variável residuals temos uma dispersão muito grande nos erros.

Bem pessoal é isso, espero que tenha ficado claro a aplicação dessa técnica e que o MARS entre aí no cinto de utilidades de analistas de BI ou Data Scientists.

Forte abraço!

Um walkthrough com o MARS

Deixe o seu comentário inteligente e educado! :o)

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s