Um estudo sobre modelos de aprendizagem baseados em árvores com desafio do Kaggle
Um estudo aplicado de modelos de aprendizagem baseados em árvores utilizando a base de dados do Kaggle para prever o preço final de casas residenciais em Ames, Iowa, utilizando uma variedade de aspectos
Post publicado originalmente em: https://gomesfellipe.github.io/post/2018-08-31-modelos-em-arvore/modelos-em-arvore/
KAGGLE
Segundo o Wikipédia: “Kaggle é a maior comunidade mundial de cientistas de dados e machine learning.” Aprendo muito estudando as resoluções de alguns competidores pois lá é possível conferir tanto as metodologias utilizadas pelos competidores quando os códigos e é notável o cuidado dos participantes para que seja possível a reprodutibilidade dos resultados, o que pode impulsionar o aprendizado.
O Kaggle trabalha com a ideia de gamificação, que é um assunto do qual já escrevi em um post sobre gamificação e porque aprender R é tão divertido e gosto deste conceito de se criar jogos para motivar e engajar as pessoas em atividades profissionais e a ideia de se estar em um jogo possibilita doses de motivação especialmente a quem gosta de competir.
A plataforma é focada em competições que envolvem modelagem preditiva, que julgam apenas o seu desempenho preditivo, embora a inteligibilidade não deixe de ser importante. Neste post farei também a modelagem descritiva com modelos de aprendizagem baseados em árvores, na qual o principal objetivo será obter informações sobre os dados para o ajuste dos modelos preditivos que iremos submeter à competição do Kaggle House Prices: Advanced Regression Techniques.
A diferença entre modelos preditivos e descritivos não é tão rigorosa assim pois algumas das técnicas podem ser utilizadas para ambos e geralmente um modelo pode servir para ambos os propósitos (mesmo que de de forma insuficiente).
Além dos modelos de machine learning baseados em árvores, também será ajustado um modelo de regressão linear multivariado para compararmos os resultados dos ajustes e submeter nossas previsões no site do kaggle.
Os pacotes que serão utilizados serão os seguintes:
library(purrr) # Programacao funciona
library(broom) # Arrumar outputs
library(dplyr) # Manipulacao de dados
library(magrittr) # pipes
library(funModeling) # df_status()
library(plyr) # revalue()
library(gridExtra) # Juntar ggplots
library(reshape) # funcao melt()
library(rpart) # Arvore de Decisoes
library(rpart.plot) # Plot da Arvore de Decisoes
library(data.table) # aux na manipulacao do heatmap
library(readr) # Leitura da base de dados
library(stringr) # Manipulacao de strings
library(ggplot2) # Graficos elegantes
library(caret) # Machine Learning
library(GGally) # up ggplot
library(ggfortify) # autoplot()BASE DE DADOS
A base de dados deste post vem de uma competição ótima para estudantes de ciência de dados de dados com alguma experiência com R ou Python e noções básicas de machine learning e estatística.
Pode ser útil para aqueles que desejam expandir seu conjunto de habilidades em uma tarefa de regressão, quando a variável y que desejamos estimar é do tipo numérico (contínuo ou discreto).
Trata-se do conjunto de dados Ames Housing que foi compilado por Dean De Cock para uso em educação de ciência de dados.
train <- read_csv("train.csv")
test <- read_csv("test.csv")
full <- bind_rows(train, test) id <- test$Id
full %<>% select(-Id)
DESCRIÇÃO DA COMPETIÇÃO
Traduzido do site oficial do kaggle:
“Peça a um comprador que descreva a casa dos seus sonhos, e eles provavelmente não começarão com a altura do teto do porão ou a proximidade de uma ferrovia leste-oeste. Mas o conjunto de dados desta competição de playground prova que muito mais influencia as negociações de preço do que o número de quartos ou uma cerca branca.
Com 79 variáveis explicativas descrevendo (quase) todos os aspectos de casas residenciais em Ames, Iowa, esta competição desafia você a prever o preço final de cada casa.”
Portanto, primeiramente vamos entender o comportamento da variável resposta, depois buscar quais dessas 79 variáveis explicativas são mais importantes para representar a variação do preço de venda das casas através dos métodos baseados em árvores e por fim ajustar os modelos propostos e submeter nossas estimativas no site!
ANÁLISE EXPLORATÓRIA DOS DADOS
Antes de pensar em ajustar algum modelo é extremamente necessário entender como se comportam os dados, portanto, tanto a variável resposta quanto as variáveis explicativas serão avaliadas.
Para que o post não fique tão grande apresentarei apenas uma introdução da análise exploratória dos dados e já foi apresentar os resultados da parte dos ajustes dos modelos. Convido quem gosta de análise de dados a conferir a análise exploratória com mais detalhes no link original do post:
VARIÁVEL RESPOSTA:
SalePrice - o preço de venda da propriedade em dólares. Essa é a variável de destino que estamos tentando prever.
Note que a distribuição dos dados referentes ao preço de venda se distribui de maneira assimétrica e não possuem evidências de normalidade dos dados. Apesar dos métodos baseados em árvore se tratarem de técnicas não paramétricas essa transformação será feita pois ao final deste post desejo comparar os resultados com um modelo de regressão linear múltipla.
ÁRVORE DE DECISÃO
Uma técnica muito popular que é mais comumente usada para resolver tarefas de classificação de dados porém a árvore conhecida como CART (Classification and Regression Trees)(Breiman, 1986) lida com todos os tipos de atributos (incluindo atributos numéricos que são tratados a partir da criação de intervalos). Para seu ajuste é possível realizar podas e produzir árvores binárias.
A construção da árvore é realizada por meio do algoritmo que iterativamente analisa os atributos descritivos de um conjunto de dados previamente rotulado. Sua popularidade como apoio para a tomada de decisão se deve principalmente ao fato da fácil visualização do conhecimento gerado e o fácil entendimento.
Outra característica legal da árvore de decisões é que ela permite ajustar um modelo sem um pré-processamento detalhado, pois é fácil de ajustar, aceita valores faltantes e é de fácil interpretação, veja:
library(rpart) # minsplit: o número mínimo de observações em um n
# cp: parametro de complexidade q controla o tamanho da arvor
control <-rpart.control(minsplit =10,cp = 0.006)rpartFit <- rpart(exp(SalePrice) ~ . , train, method = "anova", control = control) rpart.plot::rpart.plot(rpartFit,cex = 0.6)
No topo, vemos o primeiro nó com 100% das observações, que representa o total da base (100%). Em seguida, vemos que a primeira variável que determina o preço de venda das casas SalePrice é a variável OverallQual. As casas que apresentaram OverallQual < 7.5 ocorrem em maior proporção do que as que tiveram OverallQual>7.5. A interpretação pode continuar dessa forma recursivamente.
É possível notar que as variáveis OverallQual,Neighborhood,1stFlrSF,2ndFlrSF,GrLivArea, BsmtFinSF1foram as que melhor representaram os dados de acordo com os parâmetros que determinamos para ajustar esta árvore, vejamos com mais detalhes se existe relação linear e intensidade e direção dessa relação com o coeficiente de correlação de Pearson entre estas variáveis dois a dois e em relação à variável resposta:
devtools::source_url("https://raw.githubusercontent.com/gomesfellipe/functions/master/correlations_for_ggpairs.R")
train %>%
select(SalePrice, OverallQual, `1stFlrSF`, `2ndFlrSF`, GrLivArea, BsmtFinSF1) %>%
ggpairs(lower = list(continuous = my_fn))+
theme_bw()
Com esta figura temos muitas informações, destaca-se que todas essas variáveis possuem algum tipo de relação linear com a variável resposta, a menor correlação observada foi com o BsmtFinSF1 e a variável que apresentou a maior correlação foi a OverallQual. Atenção para a correlação entre SalePricee OverallQual, pois Overallqual parece ser uma variável ordinal e uma outra medida de correlação que melhor representaria esta relação é o coeficiente de correlação de Spearman, veja:
cor(full$SalePrice, full$OverallQual, method = "spearman", use = "complete.obs")## [1] 0.8098286
Um pouco diferente do resultado da correlação de Pearson pois avalia relações lineares, já a correlação de Spearman avalia relações monótonas, sejam elas lineares ou não.
Após a análise exploratória e inclusão dos valores faltantes, seguiremos com o ajuste dos modelos
MACHINE LEARNING COM ALGORITMOS DE APRENDIZAGEM BASEADOS EM ÁRVORES
Os métodos baseados em árvores fornecem modelos preditivos de alta precisão, estabilidade e facilidade de interpretação. Ao contrário dos modelos lineares, eles são capazes de lidar bem com relações não-lineares além de poderem ser adaptados para resolver tanto problemas de classificação quanto problemas de regressão.
Algoritmos como árvores de decisão, random forest e “gradient boosting” estão sendo muito usados em todos os tipos de problemas de data science e é notável o uso desses algorítimos para resolver os desafios do Kaggle. Para resolver este problema utilizaremos estes três algoritmos e ao final, pegando carona na seleção de variáveis para os algoritmos de árvore, será ajustado um modelo de regressão linear para compararmos e conferirmos a significância estatística de cada uma das variáveis.
VARIMP COM RANDOM FOREST
Um dos benefícios da floresta aleatória é o poder de lidar com grande conjunto de dados com maior dimensionalidade e identificar as variáveis a importância das variáveis, que pode ser uma característica muito útil porém deve ser feita com cautela.
Veja uma reflexão (traduzida) da nota de Leo Breiman (Universidade da Califórnia em Berkeley)
O ajuste da árvore será feito com o pacote caret e o estudo de estimativas de erro foi definido como o Out of bag que remove a necessidade de um conjunto de teste pois é o erro médio de previsão em cada amostra de treinamento xi , usando apenas as árvores que não tinham xi em sua amostra de bootstrap.
set.seed(1)
control <- trainControl(method = "oob",verboseIter = F)
rfFit1 <- train(SalePrice ~. ,
data=train,
method="rf",
metric = "Rsquared",
trControl = control
) randomForest::varImpPlot(rfFit1$finalModel)
rfFit1$finalModel$importance %>%
as.data.frame %>%
mutate(row = rownames(.)) %>%
arrange(desc(IncNodePurity)) %>%
as_tibble()## # A tibble: 218 x 2
## IncNodePurity row
## <dbl> <chr>
## 1 77.6 OverallQual
## 2 33.7 GrLivArea
## 3 16.3 YearBuilt
## 4 10.5 KitchenQual
## 5 9.57 GarageCars
## 6 8.88 TotalBsmtSF
## 7 7.48 GarageArea
## 8 7.42 `1stFlrSF`
## 9 5.84 ExterQualTA
## 10 4.41 FireplaceQu
## # ... with 208 more rows
Após inspecionar a importância das variáveis vamos selecionar as seguintes variáveis:
full %<>%
select(
SalePrice , Neighborhood, OverallQual , GrLivArea ,
YearBuilt , KitchenQual, GarageCars , GarageArea ,
`1stFlrSF` , ExterQual , BsmtFinSF1 , FireplaceQu ,
BsmtQual , `2ndFlrSF` , CentralAir , GarageFinish,
YearRemodAdd, FullBath , GarageYrBlt , Fireplaces ,
LotFrontage , BsmtUnfSF , TotalBsmtSF , BsmtFinType1,
OpenPorchSF , GarageType , BsmtExposure, OverallCond ,
TotalBsmtSF , LotArea
)Portanto, vamos definir novamente o conjunto de dados de treino e de teste:
train <- full[1:nrow(train),] %>% as.data.frame()
test <- full[(nrow(train)+1):nrow(full),-1] %>% as.data.frame()AJUSTANDO MODELOS
ÁRVORE DE DECISÃO
O modelo de árvore de decisão já foi comentado e deixei algumas referências ao final do post portanto vejamos a seguir o ajusto no R. Segundo a documentação:
cp: parâmetro de complexidade. No nosso caso isso significa que o R2 total deve aumentar em cp em cada etapa. O principal papel desse parâmetro é economizar tempo de computação removendo as divisões que obviamente não valem a pena. Essencialmente, informamos ao programa que qualquer divisão que não melhore o ajuste por cp provavelmente será eliminada por validação cruzada, e que, portanto, o programa não precisa buscá-la.
Para pesquisa de grade existem duas maneiras de ajustar um algoritmo no pacote caret: permitir que o sistema faça isso automaticamente ou especificar o tuneGride manualmente onde cada parâmetro do algoritmo pode ser especificado como um vetor de valores possíveis. Confira o ajuste manual em R:
set.seed(1)
control <- trainControl(method = "cv", number = 5,verboseIter = F)
tunegrid <- expand.grid(cp=seq(0.001, 0.01, 0.001)) rpartFit2 <- train(y=train$SalePrice, x=train[,-1],
method="rpart",
trControl=control,
tuneGrid=tunegrid,
metric = "Rsquared"
)
rpartFit2## CART
##
## 1460 samples
## 28 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 1169, 1166, 1169, 1167, 1169
## Resampling results across tuning parameters:
##
## cp RMSE Rsquared MAE
## 0.001 0.1851270 0.7850311 0.1335729
## 0.002 0.1882550 0.7774668 0.1360632
## 0.003 0.1937508 0.7642278 0.1406901
## 0.004 0.1986244 0.7518272 0.1451432
## 0.005 0.1983841 0.7524198 0.1466253
## 0.006 0.2018946 0.7436975 0.1489733
## 0.007 0.2029032 0.7412966 0.1498322
## 0.008 0.2036704 0.7386382 0.1504721
## 0.009 0.2063662 0.7321483 0.1516738
## 0.010 0.2088097 0.7261248 0.1541896
##
## Rsquared was used to select the optimal model using the largest value.
## The final value used for the model was cp = 0.001.
Podemos conferir os resultados novamente de maneira visual:
rpart.plot(rpartFit2$finalModel, cex = 0.5)
Gerando arquivo para submissão no kaggle:
id %>% cbind(predict(rpartFit2, test) %>% exp) %>%
`colnames<-`(c("Id", "SalePrice")) %>
write.csv("rpartFit2.csv",row.names = F)BAGGING
“Bagging” é usado quando desejamos reduzir a variação de uma árvore de decisão. Ela combina o resultado de vários modelos onde todas as variáveis são considerados para divisão um nó. Em R:
set.seed(1)
control <- trainControl(method = "cv", number = 5,verboseIter = F) treebagFit <- train(y=train$SalePrice,
x=train[,-1],
method = "treebag",
metric = "Rsquared",
trControl=control
)
treebagFit## Bagged CART
##
## 1460 samples
## 28 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 1169, 1166, 1169, 1167, 1169
## Resampling results:
##
## RMSE Rsquared MAE
## 0.180257 0.800383 0.1276895
Note que o R2 aumentou e o RMSE diminuiu após o uso desta técnica.
Resultados para enviar para o Kaggle:
id %>% cbind(predict(treebagFit, test)%>% exp) %>%
`colnames<-`(c("Id", "SalePrice")) %>%
write.csv("treebagFit.csv",row.names = F)RANDOM FOREST
A principal diferença entre “bagging” e o algoritmo Random Forest é que em randomForest, apenas um subconjunto de características é selecionado aleatoriamente em cada divisão em uma árvore de decisão enquanto que no bagging todos os recursos são usados.
Para pesquisa de grade especificaremos um vetor com os possíveis valores, pois o default adotado para o parâmetro mtry é mtry = p/3 (Número de variáveis amostradas aleatoriamente como candidatos em cada divisão), onde p é o número de variáveis e pode ser que o modelo se ajuste melhor aos dados ao utilizar outro valor.
Veja:
set.seed(1)
tunegrid <- expand.grid(mtry = seq(4, ncol(train) * 0.8, 2))
control <- trainControl(method = "cv", number = 5,verboseIter = F)
rfFit <- train(SalePrice ~. ,
data=train,
method="rf",
metric = "Rsquared",
tuneGrid=tunegrid,
trControl=control
)
rfFit## Random Forest
##
## 1460 samples
## 28 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 1169, 1166, 1169, 1167, 1169
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared MAE
## 4 0.1455418 0.8777812 0.09715090
## 6 0.1414106 0.8818890 0.09430431
## 8 0.1390004 0.8845912 0.09260696
## 10 0.1384602 0.8847086 0.09192925
## 12 0.1375739 0.8856181 0.09146023
## 14 0.1379341 0.8845199 0.09191786
## 16 0.1375007 0.8848814 0.09173468
## 18 0.1378254 0.8841412 0.09181452
## 20 0.1370225 0.8852069 0.09167087
## 22 0.1377031 0.8836794 0.09185160
##
## Rsquared was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 12.
Resultados para enviar para o Kaggle:
id %>% cbind(predict(rfFit, test) %>% exp) %>%
`colnames<-`(c("Id", "SalePrice")) %>%
write.csv("rfFit.csv",row.names = F)GBM
Diferentemente do “bagging”, o “boosting” é uma técnica de ensemble (conjunto) na qual os preditores não são feitos independentemente, mas sequencialmente. Na imagem a seguir é possível ver uma representação visual dessa diferença:
A imagem foi obtida neste artigo: Gradient Boosting from scratch, recomendo a leitura pois da uma boa intuição de como o algoritmo funciona.
Para a pesquisa de grade vamos permitir que o sistema faça isso automaticamente configurando apenas o tuneLength para indicar o número de valores diferentes para cada parâmetro do algoritmo.
set.seed(1)
control <- trainControl(method = "cv", number = 5,verboseIter = F)
gbmFit <- train(SalePrice~.,
data=train,
method = "gbm",
trControl=control,
tuneLength=5,
metric = "Rsquared",
verbose = FALSE
)
gbmFit## Stochastic Gradient Boosting
##
## 1460 samples
## 28 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 1169, 1166, 1169, 1167, 1169
## Resampling results across tuning parameters:
##
## interaction.depth n.trees RMSE Rsquared MAE
## 1 50 0.1730397 0.8351648 0.12135999
## 1 100 0.1444704 0.8699373 0.10144622
## 1 150 0.1369571 0.8808203 0.09661622
## 1 200 0.1343402 0.8848707 0.09422986
## 1 250 0.1337297 0.8860971 0.09358998
## 2 50 0.1483094 0.8669459 0.10310990
## 2 100 0.1332364 0.8868509 0.09277565
## 2 150 0.1308947 0.8906276 0.09104835
## 2 200 0.1300717 0.8921879 0.09024389
## 2 250 0.1300771 0.8921040 0.09007936
## 3 50 0.1419465 0.8747933 0.09804802
## 3 100 0.1319456 0.8889471 0.09055141
## 3 150 0.1299294 0.8925420 0.08939164
## 3 200 0.1294453 0.8936116 0.08892228
## 3 250 0.1286343 0.8949347 0.08860049
## 4 50 0.1371505 0.8816338 0.09454808
## 4 100 0.1305124 0.8909805 0.09050102
## 4 150 0.1294474 0.8932055 0.08960923
## 4 200 0.1284459 0.8950245 0.08918657
## 4 250 0.1280869 0.8959699 0.08897876
## 5 50 0.1356299 0.8840098 0.09265884
## 5 100 0.1302069 0.8918284 0.08938808
## 5 150 0.1297269 0.8929259 0.08900940
## 5 200 0.1288664 0.8946107 0.08839819
## 5 250 0.1294697 0.8939757 0.08878455
##
## Tuning parameter 'shrinkage' was held constant at a value of 0.1
##
## Tuning parameter 'n.minobsinnode' was held constant at a value of 10
## Rsquared was used to select the optimal model using the largest value.
## The final values used for the model were n.trees = 250,
## interaction.depth = 4, shrinkage = 0.1 and n.minobsinnode = 10.
Note que este foi o modelo que apresentou os melhores resultados quanto só R2 e ao RMSE em comparação com os outros modelos.
Submissão para Kaggle:
id %>% cbind(predict(gbmFit, test) %>% exp) %>%
`colnames<-`(c("Id", "SalePrice")) %>%
write.csv("gbmFit.csv", row.names = F)REGRESSÃO LINEAR
Por fim faremos o ajuste de um modelo de regressão linear multivariado utilizando o pacote caret.
Utilizaremos validação cruzada separando nossa amostra em 5 e utilizaremos o método lmStepAICque realiza a seleção do modelo escalonado pelo critério de informação de Akaike - AIC.
set.seed(1)
control <- trainControl(method = "cv", number = 5,verboseIter = F)
lmFit <- train(SalePrice~.,data=train,
method = "lmStepAIC",
trControl=control,
metric = "Rsquared",
trace=F
)
lmFit## Linear Regression with Stepwise Selection
##
## 1460 samples
## 28 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 1169, 1166, 1169, 1167, 1169
## Resampling results:
##
## RMSE Rsquared MAE
## 0.1476205 0.8625503 0.09569549
Note que o ajuste do modelo se apresenta de maneira satisfatória com R2 e RMSE semelhantes aos modelos de bagging e boosting e além disso, diferente dos modelos baseados em árvore, com este ajuste é possível notar a significância estatística de cada parâmetro ajustado, o que possibilita tanto o uso tanto como modelo preditivo quanto como modelo descritivo. Veja:
ggcoef(
lmFit$finalModel,
vline_color = "red",
errorbar_color = "blue",
errorbar_height = .25, shape = 18,
size=2,
color="black",
exclude_intercept = TRUE,
mapping = aes(x = estimate, y = term, size = p.value))+
scale_size_continuous(trans = "reverse")+
theme_bw()
Note que o intercepto β0 foi retirado da imagem pois é muito superior aos demais coeficientes. Note também que βi informa quão sensível é y, no caso log(SalePrice) às variações de cara umas das xi,j variáveis explicativas. Mais concretamente, se xi,j aumenta em uma unidade, o valor de y varia em β1unidades.
Uma rápida Análise dos Resíduos:
lmFit$finalModel %>%
autoplot(which = 1:2) +
theme_bw()
É possível notar que parece haver alguns outliers em ambas as figuras. Na primeira é possível notar uma nuvem de pontos aleatórios em torno de zero porém na segunda figura nota-se que alguns valores não estão de acordo com os quantils teóricos de uma distribuição normal, o que pode prejudicar nossa interpretação dos coeficientes do modelo. Vamos encerrar o modelo por aqui mesmo e ver como ele se sai na competição do Kaggle, preparando a submissão:
id %>% cbind(predict(lmFit, test) %>% exp ) %>%
`colnames<-`(c("Id", "SalePrice")) %>%
write.csv("lmFit.csv",row.names = F)O score obtido com esta submissão no Kaggle foi muito próximo dos modelos baseados e árvore e o tempo computacional para este ajuste foi bem menor.
COMPARANDO AJUSTES
Vejamos a seguir uma comparação entre estes modelos com as funções fornecidas pelo pacote caret.
resamps <- resamples(list(rpart = rpartFit2,
treebag = treebagFit,
rf = rfFit,
gbm = gbmFit,
lm = lmFit
))
bwplot(resamps)
Com este gráfico é possível notar que o modelo de regressão linear múltipla apresentou resultados semelhantes aos de bagging e boosting.
É importante frisar que a maneira como as variáveis foram selecionadas para o modelo de regressão linear múltipla através da importância das variáveis obtida com o modelo randomForest não é um padrão e existem diversos outros modos estatísticos de se de determinar a significância e a relação das variáveis para o modelo.
Um possível problema neste método é que não detecta a multicolinearidade, que ocorre quando as variáveis explicativas estão fortemente correlacionadas entre si e a análise de regressão linear pode ficar confusa e desprovida de significado, pois há dificuldade em distinguir o efeito de uma ou outra variável explicativa sobre a variável resposta Y devido à variâncias muito elevadas ou sinais inconsistentes.
Essa proposta de aprender se divertindo e de maneira produtiva me deixa muito empolgado, espero que tenham se divertido como eu me diverti fazendo este post!
REFERÊNCIAS:
- DataCamp Course:Machine Learning with Tree-Based Models in R
- Data Science for Business
- Learn ML Algorithms by coding: Decision Trees
- DataCamp Tutorials: Decision Trees in R
- The caret Package — Max Kuhn
- Um tutorial completo sobre modelagem baseada em árvores de decisão (códigos R e Python)
- Tuning Machine Learning Models Using the Caret R Package
- Gradient Boosting from scratch
- Tune Machine Learning Algorithms in R (random forest case study)
- Random Forests — Leo Breiman and Adele Cutler
- An Introduction to Recursive Partitioning Using the RPART Routines — CRAN
