Photo by Colin Watts on Unsplash

Random Forests 101 com R

Uma introdução sobre o time

O texto de hoje é sobre um dos modelos preditivos mais relevantes atualmente capaz de rivalizar com redes neurais complexas e que deixa muitos modelos de regressão linear no chinelo! Hoje vamos de Floretas Aleatórias! (Mais um exemplo de como a tradução direta do inglês para o português é complicado!).

Este tipo de algoritmo de machine learning se mostra excepcionalmente relevante em problemas de classificação, quando o resultado esperado não é um valor numérico, mas uma informação categórica, vulgo qualitativa.

O nome (esquisito em português) vem de sua base, as árvores de decisão, e utilizando uma série de regras de separação dos dados estas “florestas” são capazes de predições altamente acuradas.

Falando rapidamente das árvores de decisão, elas apresentam uma forma simples de classificação, separando as informações em galhos até chegar em um nível onde a informação não pode mais ser segmentada. Um exemplo lúdico, pode ser aquela brincadeira da adivinhação, onde você começa perguntando: “Sou um animal ou uma pessoa”, “Sou homem ou mulher”, “Sou uma criança ou um adulto?” e através destas perguntas você começa a estreitar as opções até chegar no melhor resultado possível, se fosse uma brinca com meu filho este resultado geralmente é Homem Aranha!

Cada pergunta desta pode ser comparada aos galhos de uma árvore de decisão e as folhas, seu nível mais baixo são as respostas possíveis. Abaixo uma representação de uma árvore de decisão olhando para os dados criada com base nos dados que exploraremos daqui a pouco!

Caso queira se aprofundar um pouco mais no tema, deixo o texto do Gabriel Stankevix aqui do blog Medium, que possui uma infinidade de conteúdos interessantes sobre data science e machine learning.

As random forests se apoiam nas árvores de decisão criando centenas de versões diferentes destas árvores com base nos dados e nas variáveis informadas. Para isso elas utilizam:

Seleção aleatória dos dados: de forma geral 2/3 dos dados informados serão selecionados de forma aleatória para serem utilizados na etapa de aprendizagem e o outro 1/3 será utilizado para validar o resultado de cada uma das árvores de decisão criadas.

Grande número de árvores: no mundo real, quanto maior uma floresta melhor para nós e no caso das random forests geralmente se trabalha com centenas arvores de decisão, não há um número padrão embora a maior parte dos exemplos encontrados trabalham com 500 árvores por modelo. Cada uma receberá uma amostra aleatória dos dados e fará predições baseadas nestas informações

Seleção aleatória de variáveis: à cada execução de uma árvore um número aleatório de variáveis de sua base de dados é selecionado para ser utilizado. Este número é calculado de acordo com a quantidade total de colunas que forem informadas na criação do modelo.

Resultado baseado na média de acerto de cada árvore: ao final da execução do algoritmo é feita uma média da taxa de acerto de cada uma das centenas de árvores de decisão e das variáveis utilizadas em cada uma, com este resultado as melhores árvores serão utilizadas para os cálculos de predição bem como as variáveis classificadas de acordo com sua utilização taxa de assertividade.

Para ficar mais claro, nada melhor que um exemplo escrito em R :)

pacotes <- c("stringr","tidyr","ggplot2", "plotly", "dplyr", "randomForest", "caTools", "data.table", "kableExtra", "caret", "rpart", "rpart.plot", "funModeling", "janitor")

if(sum(as.numeric(!pacotes %in% installed.packages())) != 0){
  instalador <- pacotes[!pacotes %in% installed.packages()]
  for(i in 1:length(instalador)) {
    install.packages(instalador, dependencies = T)
    break()}
  sapply(pacotes, require, character = T) 
} else {
  sapply(pacotes, require, character = T) 
}

Para este exemplo, utilizarei uma base de dados (disponível no Kaggle) com 11 informações diferentes de 4909 pacientes indicando aqueles que tiveram e não tiveram um ataque do coração.

Para sua utilização foi necessário transformar as variáveis categóricas em fatores e remover a coluna ID, que não será utilizada em nossa análise. Também foram removidas todas as linhas que possuíam alguma informação em branco.

dados = 
  fread("dataset-stroke-data.csv", header = T) %>%
  select(everything(), -id) %>%
  mutate(
    hypertension = as.factor(hypertension),
    heart_disease = as.factor(heart_disease),
    smoking_status = as.factor(smoking_status),
    Residence_type = as.factor(Residence_type),
    work_type = as.factor(work_type),
    ever_married = as.factor(ever_married),
    bmi = as.numeric(ifelse(bmi == "N/A", NA, bmi)),
    stroke = as.factor(ifelse(stroke==0, "N", "Y"))) %>%
  drop_na

summary(dados)
##     gender               age        hypertension heart_disease ever_married
##  Length:4909        Min.   : 0.08   0:4458       0:4666        No :1705    
##  Class :character   1st Qu.:25.00   1: 451       1: 243        Yes:3204    
##  Mode  :character   Median :44.00                                          
##                     Mean   :42.87                                          
##                     3rd Qu.:60.00                                          
##                     Max.   :82.00                                          
##          work_type    Residence_type avg_glucose_level      bmi       
##  children     : 671   Rural:2419     Min.   : 55.12    Min.   :10.30  
##  Govt_job     : 630   Urban:2490     1st Qu.: 77.07    1st Qu.:23.50  
##  Never_worked :  22                  Median : 91.68    Median :28.10  
##  Private      :2811                  Mean   :105.31    Mean   :28.89  
##  Self-employed: 775                  3rd Qu.:113.57    3rd Qu.:33.10  
##                                      Max.   :271.74    Max.   :97.60  
##          smoking_status stroke  
##  formerly smoked: 837   N:4700  
##  never smoked   :1852   Y: 209  
##  smokes         : 737           
##  Unknown        :1483           
##                                 
##

Vamos fazer uma análise exploratória dos dados e verificar se temos alguma informação inconsistente em nossas variáveis categóricas

#Seleção do nome das variáveis categóricas
dados %>% 
  select_if(is.factor) %>% 
  names()
## [1] "hypertension"   "heart_disease"  "ever_married"   "work_type"     
## [5] "Residence_type" "smoking_status" "stroke"
dados %>% tabyl(hypertension) %>% adorn_pct_formatting() %>% arrange(desc(n))
##  hypertension    n percent
##             0 4458   90.8%
##             1  451    9.2%
dados %>% tabyl(heart_disease) %>% adorn_pct_formatting() %>% arrange(desc(n))
##  heart_disease    n percent
##              0 4666   95.0%
##              1  243    5.0%
dados %>% tabyl(ever_married) %>% adorn_pct_formatting() %>% arrange(desc(n))
##  ever_married    n percent
##           Yes 3204   65.3%
##            No 1705   34.7%
dados %>% tabyl(Residence_type) %>% adorn_pct_formatting() %>% arrange(desc(n))
##  Residence_type    n percent
##           Urban 2490   50.7%
##           Rural 2419   49.3%
dados %>% tabyl(smoking_status) %>% adorn_pct_formatting() %>% arrange(desc(n))
##   smoking_status    n percent
##     never smoked 1852   37.7%
##          Unknown 1483   30.2%
##  formerly smoked  837   17.1%
##           smokes  737   15.0%
dados %>% tabyl(stroke) %>% adorn_pct_formatting() %>% arrange(desc(n))
##  stroke    n percent
##       N 4700   95.7%
##       Y  209    4.3%
dados %>% tabyl(work_type) %>% adorn_pct_formatting() %>% arrange(desc(n))
##      work_type    n percent
##        Private 2811   57.3%
##  Self-employed  775   15.8%
##       children  671   13.7%
##       Govt_job  630   12.8%
##   Never_worked   22    0.4%

A variável work_typepossui valores relativamente bem distribuídos, com exceção das pessoas que nunca trabalharam (Never_worked), com apenas 22 registros ou 0.04% de toda a base. Vamos optar pela remoção destes registros para reduzir a quantidade de categorias de nossa base, dado que representam um percentual muito pequeno de dados.

dados = 
  dados %>%
  filter(work_type != "Never_worked")

Uma vez que temos nossa base pronta, seguiremos para a separação de dados para treino de nosso modelo e para testarmos sua acurácia. Mesmo que o algoritmo faça esta verificação interna utilizarei 10% da base para mostrar como ficou resultado.

amostra = sample.split(dados$stroke, SplitRatio = .90, )

treino = subset(dados, amostra == T)
teste  = subset(dados, amostra == F)

A random forest

Embora o modelo faça o cálculo para definir a quantidade variáveis de nossa base que será usada em cada árvore, vamos fazer o cálculo do % de erro que cada combinação gera ao final de uma execução. Para isso vou utilizar o método tuneRF.

mtry = 
  tuneRF(
    #lista das variáveis que serão utilizadas para a especificação do modelo
    treino %>% select(-stroke), 
    # coluna que desejamos prever
    treino$stroke, 
    # quantidade de árvores que devem ser geradas no modelo
    ntreeTry=500,
    stepFactor=1.5,
    improve=0.01)
print(mtry)
##       mtry   OOBError
## 2.OOB    2 0.04297408
## 3.OOB    3 0.04297408
## 4.OOB    4 0.04365621

Algo importante das random foreres e de outros modelos de machine learning é que nem sempre todas as variáveis disponíveis precisam ser utilizadas, neste exemplo poderemos trabalhar com duas variáveis em cada árvore para aumentarmos a chance de um modelo eficaz.

modelo = 
  randomForest(
    # a fórmula da árvore, a informação de ataque cardíaco (stroke) em função de todas as variáveis de nossa base de dados
    formula = stroke ~ .,
    # a base de treino com 4418 registros
    data = treino, 
    # quantidade de colunas que serão utilizadas em cada árvore
    mtry = 2,
    # quantidade de árvores que devem ser criadas pelo algorítmo 
    ntree=500)

Vamos verificar, de acordo com o modelo final quais variáveis foram mais importantes em todas as 500 combinações realizadas:

varImpPlot(modelo, sort = T)

As variáveis avg_glucose_level(média dos níveis de açúcar no sangue), age(idade) e bmi(índice de massa corporal) foram as variáveis mais relevantes enquanto a ever_married (se foi casado alguma vez a menos relevante).

Com o modelo criado podemos realizar a predição de resultado com base nos valores da base de teste e comparar o resultado predito com os valores reais.

pred = predict(modelo, teste, type = "class")

Para verificar o resultado utilizaremos uma matriz de confusão, que faz a comparação entre o real e o predito

confusionMatrix(pred, teste$stroke)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   N   Y
##          N 468  21
##          Y   0   0
##                                           
##                Accuracy : 0.9571          
##                  95% CI : (0.9351, 0.9732)
##     No Information Rate : 0.9571          
##     P-Value [Acc > NIR] : 0.5577          
##                                           
##                   Kappa : 0               
##                                           
##  Mcnemar's Test P-Value : 1.275e-05       
##                                           
##             Sensitivity : 1.0000          
##             Specificity : 0.0000          
##          Pos Pred Value : 0.9571          
##          Neg Pred Value :    NaN          
##              Prevalence : 0.9571          
##          Detection Rate : 0.9571          
##    Detection Prevalence : 1.0000          
##       Balanced Accuracy : 0.5000          
##                                           
##        'Positive' Class : N               
##

Este exemplo possui uma acurácia de 95,7%, nada mal! Mas nossa matriz de confusão não traz bons resultados!

Como podemos observar na imagem abaixo, vemos que na linha N e coluna N (valores preditos e reais, respectivamente) que indica pacientes que não tiveram um ataque cardíaco o modelo foi capaz de acertar os 468 registros, ou seja 100% dos casos.

Por outro lado, temos 21 pacientes que tiveram um ataque do coração e o modelo indicou outro resultado, por isso na coluna Y e linha Y da tabela vemos um número 0, mostrando que o modelo errou feio!

##           Reference
## Prediction   N   Y
##          N 468  21
##          Y   0   0

Para resolver este problema é temos alguns caminhos:

1. rever a parametrização do algoritmo, existe um universo de opções que não foram exploradas neste texto e que podem aumentar sua acurácia
2. Aumentar a quantidade dados de treino para que o modelo possa aprender melhor com os dados históricos
3. Avaliar a substituição de variáveis, como o caso da ever_married ou inclusão de novas informações que possam ajudar na previsão de resultados

Por não ser este o objetivo do texto não vamos entrar no detalhe, mas fica a sugestão de contribuição sobre o que pode ser feito para aumentar a acurácia nos casos de em que o modelo falhou e se estiver interessado podemos escrever uma versão 2.0 deste texto :)

Algumas Limitações

Como você viu as random forestes não são a bala de prata dos modelos preditivos e por isso é preciso considera alguns pontos de atenção ao utilizá-las:

1. Elas não muito boas em modelos de regressão, onde o resultado esperado é um valor numérico como o valor de uma ação ou o preço da gasolina
2. Elas podem ser tendenciosas especialmente em casos em que variáveis categóricas possuem muitos valores (categorias), por exemplo, um campo de tipo de cliente que possua 10 ou mais opções, o tipo que possuir mais registros em sua base pode contaminar a predição dos demais.

Mesmo com estes pontos ainda sim as random forests representam um avanço nas técnicas de machine learning e uma ferramenta incrível para a criação de modelos capazes de auxiliar na tomada de decisão.

Se você encontrou algum erro neste texto de conceito ou nos scripts por favor me avise através do e-mail marcosmhs@live.com e terei prazer em fazer a correção e incluí-lo(a) nos créditos.