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, multicolinearidade, heterocedasticidade, 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).

[sourcecode language=”r”] install.packages(“caret”) install.packages(“earth”) [/sourcecode]

Pacotes instalados, vamos carregar as bibliotecas:

[sourcecode language=”r”] library(caret) library(earth)[/sourcecode]

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.

[sourcecode language=”r”] # Get data tambore_parnaiba_2013_data <- read.csv(‘https://raw.githubusercontent.com/fclesio/learning-space/master/Datasets/imoveis-tambore-parnaiba-2013.csv’); [/sourcecode]

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

[sourcecode language=”r”] # Check data tambore_parnaiba_2013_data [/sourcecode]

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:

[sourcecode language=”r”] # Eliminate Bairro variable tambore_parnaiba_2013_data$Bairro <- NULL [/sourcecode]

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:

[sourcecode language=”r”] # 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.) [/sourcecode]

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.

[sourcecode language=”r”] # Summarize the fit summary(fit,digit=3) [/sourcecode]

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

[code language=”text”] 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 [/code]

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.

[sourcecode language=”text”] 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 [/sourcecode]

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:

[sourcecode language=”r”] # Plot plotmo(fit) [/sourcecode]

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.

[sourcecode language=”r”] # Variable importance evimp(fit) [/sourcecode]

Esses são os resultados:

[sourcecode language=”text”] 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 [/sourcecode]

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.

[sourcecode language=”r”] # 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 [/sourcecode]

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.

[sourcecode language=”r”] # 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) [/sourcecode]

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:

[sourcecode language=”r”] # Final Predictions FinalPredictions$resuduals <- FinalPredictions$ValorVenda - FinalPredictions$PredictValorVenda # Final table with results FinalPredictions [/sourcecode]

E esse é o output do objeto FinalPredictions:

[sourcecode language=”text”] 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 [/sourcecode]

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:

Colocando essa representação no R, temos:

[sourcecode language=”r”] # Root Mean Square Error rmse <- sqrt(mean(((FinalPredictions$ValorVenda - FinalPredictions$PredictValorVenda)^2))) # Show results rmse [/sourcecode]

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!