Árvores de regressão básicas particionam um conjunto de dados em grupos menores e então ajustam um modelo simples (constante) para cada subgrupo. Infelizmente, um modelo de árvore única tende a ser altamente instável e um preditor pobre. No entanto, por bootstrap agregando (ensacando) árvores de regressão, essa técnica pode se tornar bastante poderosa e eficaz. Além disso, isso fornece a base fundamental de modelos baseados em árvore mais complexos, como florestas aleatórias e máquinas de aumento de gradiente.

Exemplo

Este tutorial aproveita os seguintes pacotes. A maioria desses pacotes está desempenhando um papel de suporte, enquanto a ênfase principal estará no pacote rpart. 

# Pacotes necessários

library(rsample)      
library(dplyr)       
library(rpart)       
library(rpart.plot)  
library(ipred)       
library(caret)

Para ilustrar vários conceitos usaremos os dados do Ames Housing que foram incluídos no pacote AmesHousing. 

library(AmesHousing)
ls.str(AmesHousing::make_ames())
## Alley :  Factor w/ 3 levels "Gravel","No_Alley_Access",..: 2 2 2 2 2 2 2 2 2 2 ...
## Bedroom_AbvGr :  int [1:2930] 3 2 3 3 3 3 2 2 2 3 ...
## Bldg_Type :  Factor w/ 5 levels "OneFam","TwoFmCon",..: 1 1 1 1 1 1 5 5 5 1 ...
## Bsmt_Cond :  Factor w/ 6 levels "Excellent","Fair",..: 3 6 6 6 6 6 6 6 6 6 ...
## Bsmt_Exposure :  Factor w/ 5 levels "Av","Gd","Mn",..: 2 4 4 4 4 4 3 4 4 4 ...
## Bsmt_Full_Bath :  num [1:2930] 1 0 0 1 0 0 1 0 1 0 ...
## Bsmt_Half_Bath :  num [1:2930] 0 0 0 0 0 0 0 0 0 0 ...
## Bsmt_Qual :  Factor w/ 6 levels "Excellent","Fair",..: 6 6 6 6 3 6 3 3 3 6 ...
## Bsmt_Unf_SF :  num [1:2930] 441 270 406 1045 137 ...
## BsmtFin_SF_1 :  num [1:2930] 2 6 1 1 3 3 3 1 3 7 ...
## BsmtFin_SF_2 :  num [1:2930] 0 144 0 0 0 0 0 0 0 0 ...
## BsmtFin_Type_1 :  Factor w/ 7 levels "ALQ","BLQ","GLQ",..: 2 6 1 1 3 3 3 1 3 7 ...
## BsmtFin_Type_2 :  Factor w/ 7 levels "ALQ","BLQ","GLQ",..: 7 4 7 7 7 7 7 7 7 7 ...
## Central_Air :  Factor w/ 2 levels "N","Y": 2 2 2 2 2 2 2 2 2 2 ...
## Condition_1 :  Factor w/ 9 levels "Artery","Feedr",..: 3 2 3 3 3 3 3 3 3 3 ...
## Condition_2 :  Factor w/ 8 levels "Artery","Feedr",..: 3 3 3 3 3 3 3 3 3 3 ...
## Electrical :  Factor w/ 6 levels "FuseA","FuseF",..: 5 5 5 5 5 5 5 5 5 5 ...
## Enclosed_Porch :  int [1:2930] 0 0 0 0 0 0 170 0 0 0 ...
## Exter_Cond :  Factor w/ 5 levels "Excellent","Fair",..: 5 5 5 5 5 5 5 5 5 5 ...
## Exter_Qual :  Factor w/ 4 levels "Excellent","Fair",..: 4 4 4 3 4 4 3 3 3 4 ...
## Exterior_1st :  Factor w/ 16 levels "AsbShng","AsphShn",..: 4 14 15 4 14 14 6 7 6 14 ...
## Exterior_2nd :  Factor w/ 17 levels "AsbShng","AsphShn",..: 11 15 16 4 15 15 6 7 6 15 ...
## Fence :  Factor w/ 5 levels "Good_Privacy",..: 5 3 5 5 3 5 5 5 5 5 ...
## Fireplace_Qu :  Factor w/ 6 levels "Excellent","Fair",..: 3 4 4 6 6 3 4 4 6 6 ...
## Fireplaces :  int [1:2930] 2 0 0 2 1 1 0 0 1 1 ...
## First_Flr_SF :  int [1:2930] 1656 896 1329 2110 928 926 1338 1280 1616 1028 ...
## Foundation :  Factor w/ 6 levels "BrkTil","CBlock",..: 2 2 2 2 3 3 3 3 3 3 ...
## Full_Bath :  int [1:2930] 1 1 1 2 2 2 2 2 2 2 ...
## Functional :  Factor w/ 8 levels "Maj1","Maj2",..: 8 8 8 8 8 8 8 8 8 8 ...
## Garage_Area :  num [1:2930] 528 730 312 522 482 470 582 506 608 442 ...
## Garage_Cars :  num [1:2930] 2 1 1 2 2 2 2 2 2 2 ...
## Garage_Cond :  Factor w/ 6 levels "Excellent","Fair",..: 6 6 6 6 6 6 6 6 6 6 ...
## Garage_Finish :  Factor w/ 4 levels "Fin","No_Garage",..: 1 4 4 1 1 1 1 3 3 1 ...
## Garage_Qual :  Factor w/ 6 levels "Excellent","Fair",..: 6 6 6 6 6 6 6 6 6 6 ...
## Garage_Type :  Factor w/ 7 levels "Attchd","Basment",..: 1 1 1 1 1 1 1 1 1 1 ...
## Gr_Liv_Area :  int [1:2930] 1656 896 1329 2110 1629 1604 1338 1280 1616 1804 ...
## Half_Bath :  int [1:2930] 0 0 1 1 1 1 0 0 0 1 ...
## Heating :  Factor w/ 6 levels "Floor","GasA",..: 2 2 2 2 2 2 2 2 2 2 ...
## Heating_QC :  Factor w/ 5 levels "Excellent","Fair",..: 2 5 5 1 3 1 1 1 1 3 ...
## House_Style :  Factor w/ 8 levels "One_and_Half_Fin",..: 3 3 3 3 8 8 3 3 3 8 ...
## Kitchen_AbvGr :  int [1:2930] 1 1 1 1 1 1 1 1 1 1 ...
## Kitchen_Qual :  Factor w/ 5 levels "Excellent","Fair",..: 5 5 3 1 5 3 3 3 3 3 ...
## Land_Contour :  Factor w/ 4 levels "Bnk","HLS","Low",..: 4 4 4 4 4 4 4 2 4 4 ...
## Land_Slope :  Factor w/ 3 levels "Gtl","Mod","Sev": 1 1 1 1 1 1 1 1 1 1 ...
## Latitude :  num [1:2930] 42.1 42.1 42.1 42.1 42.1 ...
## Longitude :  num [1:2930] -93.6 -93.6 -93.6 -93.6 -93.6 ...
## Lot_Area :  int [1:2930] 31770 11622 14267 11160 13830 9978 4920 5005 5389 7500 ...
## Lot_Config :  Factor w/ 5 levels "Corner","CulDSac",..: 1 5 1 1 5 5 5 5 5 5 ...
## Lot_Frontage :  num [1:2930] 141 80 81 93 74 78 41 43 39 60 ...
## Lot_Shape :  Factor w/ 4 levels "Regular","Slightly_Irregular",..: 2 1 2 1 2 2 1 2 2 1 ...
## Low_Qual_Fin_SF :  int [1:2930] 0 0 0 0 0 0 0 0 0 0 ...
## Mas_Vnr_Area :  num [1:2930] 112 0 108 0 0 20 0 0 0 0 ...
## Mas_Vnr_Type :  Factor w/ 5 levels "BrkCmn","BrkFace",..: 5 4 2 4 4 2 4 4 4 4 ...
## Misc_Feature :  Factor w/ 6 levels "Elev","Gar2",..: 3 3 2 3 3 3 3 3 3 3 ...
## Misc_Val :  int [1:2930] 0 0 12500 0 0 0 0 0 0 0 ...
## Mo_Sold :  int [1:2930] 5 6 6 4 3 6 4 1 3 6 ...
## MS_SubClass :  Factor w/ 16 levels "One_Story_1946_and_Newer_All_Styles",..: 1 1 1 1 6 6 12 12 12 6 ...
## MS_Zoning :  Factor w/ 7 levels "Floating_Village_Residential",..: 3 2 3 3 3 3 3 3 3 3 ...
## Neighborhood :  Factor w/ 29 levels "North_Ames","College_Creek",..: 1 1 1 1 7 7 17 17 17 7 ...
## Open_Porch_SF :  int [1:2930] 62 0 36 0 34 36 0 82 152 60 ...
## Overall_Cond :  Factor w/ 10 levels "Very_Poor","Poor",..: 5 6 6 5 5 6 5 5 5 5 ...
## Overall_Qual :  Factor w/ 10 levels "Very_Poor","Poor",..: 6 5 6 7 5 6 8 8 8 7 ...
## Paved_Drive :  Factor w/ 3 levels "Dirt_Gravel",..: 2 3 3 3 3 3 3 3 3 3 ...
## Pool_Area :  int [1:2930] 0 0 0 0 0 0 0 0 0 0 ...
## Pool_QC :  Factor w/ 5 levels "Excellent","Fair",..: 4 4 4 4 4 4 4 4 4 4 ...
## Roof_Matl :  Factor w/ 8 levels "ClyTile","CompShg",..: 2 2 2 2 2 2 2 2 2 2 ...
## Roof_Style :  Factor w/ 6 levels "Flat","Gable",..: 4 2 4 4 2 2 2 2 2 2 ...
## Sale_Condition :  Factor w/ 6 levels "Abnorml","AdjLand",..: 5 5 5 5 5 5 5 5 5 5 ...
## Sale_Price :  int [1:2930] 215000 105000 172000 244000 189900 195500 213500 191500 236500 189000 ...
## Sale_Type :  Factor w/ 10 levels "COD","Con","ConLD",..: 10 10 10 10 10 10 10 10 10 10 ...
## Screen_Porch :  int [1:2930] 0 120 0 0 0 0 0 144 0 0 ...
## Second_Flr_SF :  int [1:2930] 0 0 0 0 701 678 0 0 0 776 ...
## Street :  Factor w/ 2 levels "Grvl","Pave": 2 2 2 2 2 2 2 2 2 2 ...
## Three_season_porch :  int [1:2930] 0 0 0 0 0 0 0 0 0 0 ...
## Total_Bsmt_SF :  num [1:2930] 1080 882 1329 2110 928 ...
## TotRms_AbvGrd :  int [1:2930] 7 5 6 8 6 7 6 5 5 7 ...
## Utilities :  Factor w/ 3 levels "AllPub","NoSeWa",..: 1 1 1 1 1 1 1 1 1 1 ...
## Wood_Deck_SF :  int [1:2930] 210 140 393 0 212 360 0 0 237 140 ...
## Year_Built :  int [1:2930] 1960 1961 1958 1968 1997 1998 2001 1992 1995 1999 ...
## Year_Remod_Add :  int [1:2930] 1960 1961 1958 1968 1998 1998 2001 1992 1996 1999 ...
## Year_Sold :  int [1:2930] 2010 2010 2010 2010 2010 2010 2010 2010 2010 2010 ...

A ideia

Existem muitas metodologias para construir árvores de regressão, mas uma das mais antigas é conhecida como a abordagem de árvore de classificação e regressão (CART) desenvolvida por Breiman et al. (1984). Árvores de regressão básicas particionam um conjunto de dados em subgrupos menores e então ajustam uma constante simples para cada observação no subgrupo. O particionamento é obtido por partições binárias sucessivas (também conhecidas como particionamento recursivo) com base nos diferentes preditores. A constante a prever é baseada nos valores médios de resposta para todas as observações que se enquadram nesse subgrupo.

Existem várias vantagens nas árvores de regressão:

Mas também existem alguns pontos fracos significativos:

# Criando conjuntos de treinamento (70%) e teste (30%) para os dados AmesHousing::make_ames().
set.seed(123) # Use set.seed para reprodutibilidade 
ames_split <- initial_split(AmesHousing::make_ames(), prop = .7)
ames_train <- training(ames_split)
ames_test  <- testing(ames_split)

Implementação Básica

Podemos ajustar uma árvore de regressão usando rpart e depois visualizá-la usando rpart.plot. O processo de ajuste e a saída visual das árvores de regressão e árvores de classificação são muito semelhantes. Ambos usam o método da fórmula para expressar o modelo (semelhante ao lm). No entanto, ao ajustar uma árvore de regressão, precisamos definir method = “anova”. Por padrão, rpart fará uma estimativa inteligente de qual valor do método deve ser baseado no tipo de dados de sua coluna de resposta, mas é recomendado que você defina explicitamente o método por motivos de reprodutibilidade (uma vez que o adivinhador automático pode mudar no futuro). 

m1 <- rpart(formula = Sale_Price ~ ., data    = ames_train, method  = "anova")

Depois de ajustar nosso modelo, podemos obter um pico na saída de m1. Isso apenas explica as etapas das divisões.

Por exemplo, começamos com 2051 observações no nó raiz (bem no início) e a primeira variável que dividimos (a primeira variável que otimiza uma redução na SSE) é Overall_Qual.

Vemos que no primeiro nó todas as observações com:
Overall_Qual = Very_Poor,Poor,Fair,Below_Average,Average,Above_Average,Good
vão para o 2º (2)) ramal.

O número total de observações que seguem este ramo (1725), seu preço médio de venda (155077.10) e SSE (4.089851e+12) são listados. Se você olhar para o terceiro ramo (3)), verá que 326 observações com Overall_Qual = Very_Good, Excellent, Very_Excellent seguem este ramo e seus preços médios de venda são 305506.30 e o SEE nesta região é 2.902384e+12. Basicamente, isso está nos dizendo que a variável mais importante que tem a maior redução em SEE inicialmente é Overall_Qual com as casas na extremidade superior do espectro de qualidade tendo quase o dobro do preço médio de venda. 

m1
## n= 2051 
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##  1) root 2051 1.319672e+13 178987.30  
##    2) Overall_Qual=Very_Poor,Poor,Fair,Below_Average,Average,Above_Average,Good 1725 4.089851e+12 155077.10  
##      4) Neighborhood=North_Ames,Old_Town,Edwards,Sawyer,Mitchell,Brookside,Iowa_DOT_and_Rail_Road,South_and_West_of_Iowa_State_University,Meadow_Village,Briardale,Northpark_Villa,Blueste,Landmark 1039 1.414812e+12 131469.60  
##        8) Overall_Qual=Very_Poor,Poor,Fair,Below_Average 202 1.682257e+11  96476.75 *
##        9) Overall_Qual=Average,Above_Average,Good 837 9.395439e+11 139914.60  
##         18) First_Flr_SF< 1240 647 4.048057e+11 131629.20 *
##         19) First_Flr_SF>=1240 190 3.390731e+11 168128.90 *
##      5) Neighborhood=College_Creek,Somerset,Northridge_Heights,Gilbert,Northwest_Ames,Sawyer_West,Crawford,Timberland,Northridge,Stone_Brook,Clear_Creek,Bloomington_Heights,Veenker,Green_Hills 686 1.218974e+12 190832.40  
##       10) Gr_Liv_Area< 1725 497 5.267647e+11 176577.80  
##         20) Total_Bsmt_SF< 1217.5 317 2.021716e+11 163165.00 *
##         21) Total_Bsmt_SF>=1217.5 180 1.671280e+11 200199.30 *
##       11) Gr_Liv_Area>=1725 189 3.256624e+11 228316.80 *
##    3) Overall_Qual=Very_Good,Excellent,Very_Excellent 326 2.902384e+12 305506.30  
##      6) Overall_Qual=Very_Good 235 9.229383e+11 270486.30  
##       12) Gr_Liv_Area< 1956.5 146 3.185953e+11 241512.80 *
##       13) Gr_Liv_Area>=1956.5 89 2.807253e+11 318015.80 *
##      7) Overall_Qual=Excellent,Very_Excellent 91 9.469774e+11 395942.60  
##       14) Gr_Liv_Area< 1958 38 6.559161e+10 337195.90 *
##       15) Gr_Liv_Area>=1958 53 6.562129e+11 438062.80  
##         30) Neighborhood=Edwards,Somerset,Veenker 7 9.951422e+10 290836.40 *
##         31) Neighborhood=College_Creek,Old_Town,Northridge_Heights,Timberland,Northridge,Stone_Brook 46 3.818802e+11 460466.80  
##           62) Total_Bsmt_SF< 2394 39 1.996090e+11 437638.20 *
##           63) Total_Bsmt_SF>=2394 7 4.870803e+10 587655.30 *

Podemos visualizar nosso modelo com rpart.plot. rpart.plot tem muitas opções de plotagem, que deixaremos para o leitor explorar. No entanto, na impressão padrão, ele mostrará a porcentagem de dados que caem para aquele nó e o preço médio de venda para aquela filial.

Uma coisa que você pode notar é que esta árvore contém 11 nós internos, resultando em 12 nós terminais. Basicamente, esta árvore está particionando em 11 variáveis para produzir seu modelo. No entanto, existem 80 variáveis em ames_train. Então o que aconteceu? 

rpart.plot(m1)

Nos bastidores, rpart está aplicando automaticamente uma gama complexa de valores \(\alpha\) para podar a árvore. Para comparar o erro de cada valor \(\alpha\), rpart executa uma validação cruzada de 10 vezes para que o erro associado a um determinado valor \(\alpha\) é calculado nos dados de validação. Neste exemplo, encontramos retornos decrescentes após 12 nós terminais, conforme ilustrado abaixo (o eixo \(y\) é um erro de validação cruzada, o eixo \(x\) inferior é a complexidade de custo, valor \(\alpha\), o eixo \(x\) superior é o número de nós terminais (tamanho da árvore = |T|).

Você também pode notar a linha tracejada que atravessa o ponto |T|=9. Breiman et al. (1984) sugeriram que, na prática real, é comum usar a menor árvore dentro de 1 desvio padrão do erro de validação cruzada mínimo (também conhecida como regra 1-SE). Assim, poderíamos usar uma árvore com 9 nós terminais e esperar obter resultados semelhantes com uma pequena margem de erro. 

plotcp(m1)

Para ilustrar o ponto de selecionar uma árvore com 12 nós terminais (ou 9 se você seguir a regra 1-SE), podemos forçar rpart a gerar uma árvore completa usando cp = 0 (nenhuma penalidade resulta em uma árvore totalmente crescida). Podemos ver que depois de 12 nós terminais, vemos retornos decrescentes na redução de erros à medida que a árvore se torna mais profunda. Assim, podemos podar significativamente nossa árvore e ainda assim atingir o erro mínimo esperado. xval é o número de réplicas de validação cruzada. 

m2 <- rpart(formula = Sale_Price ~ ., data = ames_train, method  = "anova", control = list(cp = 0, xval = 10))
plotcp(m2)
abline(v = 12, lty = "dashed")

Portanto, por padrão, rpart está realizando algum ajuste automatizado, com uma subárvore ideal de 11 divisões, 12 nós terminais e um erro de validação cruzada de 0.241 (observe que esse erro é equivalente à estatística PRESS, mas não ao MSE). No entanto, podemos realizar ajustes adicionais para tentar melhorar o desempenho do modelo. 

m1$cptable
##            CP nsplit rel error    xerror       xstd
## 1  0.47015337      0 1.0000000 1.0013200 0.05921212
## 2  0.11033540      1 0.5298466 0.5311496 0.02925211
## 3  0.07823674      2 0.4195112 0.4213588 0.02860305
## 4  0.02777564      3 0.3412745 0.3435680 0.02028168
## 5  0.02452259      4 0.3134988 0.3165230 0.02007149
## 6  0.02326660      5 0.2889763 0.2997117 0.01974815
## 7  0.01706279      6 0.2657097 0.2764728 0.01969300
## 8  0.01482680      7 0.2486469 0.2764824 0.01969404
## 9  0.01324712      8 0.2338201 0.2574236 0.01802444
## 10 0.01193214      9 0.2205729 0.2532587 0.01762624
## 11 0.01012094     10 0.2086408 0.2444527 0.01761806
## 12 0.01000000     11 0.1985199 0.2411865 0.01669430

Afinação (Tuning)

Além do parâmetro de complexidade de custo (\(\alpha\)), também é comum ajustar:

  • minsplit: o número mínimo de pontos de dados necessários para tentar uma divisão antes de ser forçado a criar um nó terminal. O padrão é 20. Tornar isso menor permite que nós terminais que podem conter apenas um punhado de observações criem o valor previsto.
  • maxdepth: o número máximo de nós internos entre o nó raiz e os nós terminais. O padrão é 30, que é bastante liberal e permite a construção de árvores bastante grandes.

rpart usa um argumento de controle especial onde fornecemos uma lista de valores de hiperparâmetros. Por exemplo, se quiséssemos avaliar um modelo com minsplit = 10 e maxdepth = 12, poderíamos executar o seguinte: 

m3 <- rpart(formula = Sale_Price ~ ., data = ames_train, method  = "anova", 
    control = list(minsplit = 10, maxdepth = 12, xval = 10))

Embora útil, essa abordagem requer que você avalie manualmente vários modelos. Em vez disso, podemos realizar uma pesquisa de grade para pesquisar automaticamente em uma gama de modelos ajustados de forma diferente para identificar a configuração ideal do hiperparâmetro.

Para realizar uma pesquisa de grade, primeiro criamos nossa grade de hiperparâmetros. Neste exemplo, procuramos um intervalo de minsplit de 5-20 e variamos maxdepth de 8-15 (já que nosso modelo original encontrou uma profundidade ideal de 12). O que resulta são 128 combinações diferentes, o que requer 128 modelos diferentes. 

hyper_grid <- expand.grid(minsplit = seq(5, 20, 1), maxdepth = seq(8, 15, 1) ) 
head(hyper_grid)
##   minsplit maxdepth
## 1        5        8
## 2        6        8
## 3        7        8
## 4        8        8
## 5        9        8
## 6       10        8
# número total de combinações
nrow(hyper_grid)
## [1] 128

Para automatizar a modelagem, simplesmente configuramos um loop for e iteramos através de cada combinação de minsplit e maxdepth. Salvamos cada modelo em seu próprio item de lista. 

models <- list()
for (i in 1:nrow(hyper_grid)) {
  # obtemos valores de minsplit, maxdepth na linha i
  minsplit <- hyper_grid$minsplit[i]
  maxdepth <- hyper_grid$maxdepth[i]
  # treinamos um modelo e o armazenamos na lista
  models[[i]] <- rpart(
    formula = Sale_Price ~ .,
    data    = ames_train,
    method  = "anova",
    control = list(minsplit = minsplit, maxdepth = maxdepth) )
}

Podemos agora criar uma função para extrair o erro mínimo associado ao valor α da complexidade de custo ideal para cada modelo. Após uma pequena discussão de dados para extrair o valor α ideal e seu respectivo erro, adicionando-o de volta à nossa grade, e filtrar pelos 5 principais valores de erro mínimo, vemos que o modelo ideal faz uma pequena melhoria em relação ao nosso modelo anterior (xerror de 0,242 versus 0,272). 

# função para obter o cp ideal
get_cp <- function(x) {
  min    <- which.min(x$cptable[, "xerror"])
  cp <- x$cptable[min, "CP"] 
}

# função para obter erro mínimo
get_min_error <- function(x) {
  min    <- which.min(x$cptable[, "xerror"])
  xerror <- x$cptable[min, "xerror"] 
}

hyper_grid %>%
  mutate(
    cp    = purrr::map_dbl(models, get_cp),
    error = purrr::map_dbl(models, get_min_error)
    ) %>%
  arrange(error) %>%
  top_n(-5, wt = error)
##   minsplit maxdepth         cp     error
## 1        6       10 0.01000000 0.2331321
## 2        7       12 0.01012094 0.2335002
## 3       15       14 0.01012094 0.2352524
## 4        6        9 0.01000000 0.2356500
## 5       17       14 0.01000000 0.2357567

Se estivéssemos satisfeitos com esses resultados, poderíamos aplicar este modelo final ótimo e fazer uma previsão em nosso conjunto de teste. O RMSE final é 38595.29, o que sugere que, em média, nossos preços de venda previstos estão cerca de $ 38595,29 abaixo do preço de venda real. 

optimal_tree <- rpart(
    formula = Sale_Price ~ .,
    data    = ames_train,
    method  = "anova",
    control = list(minsplit = 11, maxdepth = 8, cp = 0.01)
    )

pred <- predict(optimal_tree, newdata = ames_test)
RMSE(pred = pred, obs = ames_test$Sale_Price)
## [1] 38595.29

Ensacamento

A ideia

Como mencionado anteriormente, os modelos de árvore única sofrem de alta variação. Embora podar a árvore ajude a reduzir essa variação, existem métodos alternativos que realmente exploram a variabilidade de árvores individuais de uma maneira que pode melhorar significativamente o desempenho além do das árvores individuais. A agregação bootstrap (bootstrap aggregating bagging) é uma dessas abordagens (originalmente proposta por Breiman, 1996).

O ensacamento combina e calcula a média de vários modelos. A média em várias árvores reduz a variabilidade de qualquer árvore e reduz o sobreajuste, o que melhora o desempenho preditivo. O ensacamento segue três etapas simples:

  • Crie \(m\) amostras bootstrap a partir dos dados de treinamento. As amostras inicializadas nos permitem criar muitos conjuntos de dados ligeiramente diferentes, mas com a mesma distribuição do conjunto de treinamento geral.
  • Para cada amostra bootstrap, treine uma única árvore de regressão não ajustada.
  • Predições individuais médias de cada árvore para criar um valor predito médio geral.

Este processo pode realmente ser aplicado a qualquer modelo de regressão ou classificação; no entanto, ele fornece o maior aprimoramento para modelos com alta variância. Por exemplo, modelos paramétricos mais estáveis, como regressão linear, tendem a apresentar menos melhorias no desempenho preditivo.

Um benefício do bagging é que, em média, uma amostra de bootstrap conterá 63% (2/3) dos dados de treinamento. Isso deixa cerca de 33% (1/3) dos dados fora da amostra inicializada. Chamamos isso de amostra out-of-bag (OOB). Podemos usar as observações OOB para estimar a precisão do modelo, criando um processo natural de validação cruzada.

Bagging com ipred

Ajustar um modelo de árvore ensacado é bastante simples. Em vez de usar rpart, usamos ipred::bagging. Usamos coob = TRUE para usar a amostra OOB para estimar o erro de teste. Vemos que nosso erro de estimativa inicial está perto de $2K menos do que o erro de teste que alcançamos com nossa única árvore ótima (36389.23 vs. 38595.29) 

# fazendo o bootstrap reproducivel
set.seed(123)

# treinando o modelo ensacado (bagged model)
bagged_m1 <- bagging(
  formula = Sale_Price ~ .,
  data    = ames_train,
  coob    = TRUE
)

bagged_m1
## 
## Bagging regression trees with 25 bootstrap replications 
## 
## Call: bagging.data.frame(formula = Sale_Price ~ ., data = ames_train, 
##     coob = TRUE)
## 
## Out-of-bag estimate of root mean squared error:  36389.23

Uma coisa a se notar é que, normalmente, quanto mais árvores, melhor. À medida que adicionamos mais árvores, estamos calculando a média de mais árvores únicas de alta variação. O resultado é que, logo no início, vemos uma redução dramática na variância (e, portanto, nosso erro) e, eventualmente, a redução no erro irá nivelar sinalizando um número apropriado de árvores para criar um modelo estável. Raramente você precisará de mais de 50 árvores para estabilizar o erro.

Por padrão, o bagging executa 25 amostras e árvores de bootstrap, mas podemos exigir mais. Podemos avaliar o erro em relação ao número de árvores conforme abaixo. Vemos que o erro está se estabilizando em cerca de 19 árvores, então provavelmente não obteremos muitas melhorias simplesmente ensacando mais árvores. 

# avaliando 10-50 árvores ensacadas
ntree <- 10:50

# criando um vetor vazio para armazenar valores OOB RMSE
rmse <- vector(mode = "numeric", length = length(ntree))

for (i in seq_along(ntree)) {
  # reprodutibilidade
  set.seed(123)
  
  # realizando o modelo ensacado
  model <- bagging(
  formula = Sale_Price ~ .,
  data    = ames_train,
  coob    = TRUE,
  nbagg   = ntree[i]
)
  # obtendo o erro OOB 
  rmse[i] <- model$err
}

plot(ntree, rmse, type = 'l', lwd = 2)
abline(v = 19, col = "red", lty = "dashed")
grid()

Ensacando com caret

O ensacamento com ipred é bastante simples; no entanto, existem alguns benefícios adicionais de ensacamento com cursor (caret):

  • É mais fácil realizar a validação cruzada. Embora possamos usar o erro OOB, executar a validação cruzada também fornecerá um entendimento mais robusto do verdadeiro erro de teste esperado.
  • Podemos avaliar a importância variável entre as árvores ensacadas.

Aqui, realizamos um modelo de validação cruzada de 10 vezes. Vemos que o RMSE com validação cruzada é de $ 35589.21. Também avaliamos as 20 principais variáveis de nosso modelo. A importância da variável para árvores de regressão é medida avaliando-se a quantidade total de SSE diminuída por divisões em um determinado preditor, calculada a média de todas as \(m\) árvores. Os preditores com o maior impacto médio para ESS são considerados os mais importantes. O valor de importância é simplesmente a diminuição média relativa em SSE em comparação com a variável mais importante (fornece uma escala de 0-100). 

# Especificando validação cruzada de 10 vezes
ctrl <- trainControl(method = "cv",  number = 10) 

# CV modelo ensacado (bagged)
bagged_cv <- train(
  Sale_Price ~ .,
  data = ames_train,
  method = "treebag",
  trControl = ctrl,
  importance = TRUE
  )

# avaliando os resultados
bagged_cv
## Bagged CART 
## 
## 2051 samples
##   80 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 1844, 1846, 1846, 1846, 1846, 1846, ... 
## Resampling results:
## 
##   RMSE      Rsquared   MAE     
##   35589.21  0.8110434  23998.56
# mostrando as variáveis mais importantes
plot(varImp(bagged_cv), 20)

Se compararmos isso com o teste estabelecido fora da amostra, vemos que nossa estimativa de erro com validação cruzada estava muito próxima. Reduzimos com sucesso nosso erro para cerca de $ 37.000. 

pred <- predict(bagged_cv, ames_test)
RMSE(pred, ames_test$Sale_Price)
## [1] 37333.42

Aprendendo mais

As árvores de decisão fornecem uma abordagem de modelagem muito intuitiva com vários benefícios flexíveis. Infelizmente, eles sofrem de alta variação; no entanto, ao combiná-los com o ensacamento, você pode minimizar essa desvantagem. Além disso, as árvores ensacadas fornecem a base fundamental para algoritmos mais complexos e poderosos, como florestas aleatórias e máquinas de aumento de gradiente.