covid-19: uma primeira abordagem

Esse post é uma análise em conjunto dos seguintes colaboradores:

• Paulo C. Marques F.
• Hedibert F. Lopes (site)
• Tiago Mendonça dos Santos

Nesse post apresentamos um classificador baseado nos dados que o Hospital Israelita Albert Einstein publicou no Kaggle. Após uma longa análise exploratória, chegamos a um conjuto de dados com apenas 6 preditoras (faça o download dos dados clicando aqui). A partir dessa seleção, utilizamos os modelos de Floresta Aleatória com o pacote ranger e boosting com o pacote xgboost. Posteriormente, avaliamos a performance preditiva desses modelos com método de reamostragem.

A performance preditiva será avaliada com base nas seguintes métricas:

• sensibilidade (taxa de verdadeiro positivo): probabilidade de classificar um indivíduo como positivo para o vírus dado que esse indivíduo é, de fato, positivo.

• especificidade (taxa verdadeiro negativo): probabilidade de classificar um indivíduo como negativo para o vírus dado que esse indivíduo é, de fato, negativo.

• AUC: área sob a curva ROC. A curva ROC é um gráfico dado pela \(\text{sensibilidade}\) no eixo y e \(1 - \text{especificidade}\) no eixo x. Queremos que esse número seja o mais próximo possível de 1.

O controle dessas métricas é de extrema importância. Note que a consequência de ter uma alta taxa de falso negativo (ou seja, 1 - sensibilidade) seria dispensar muitas pessoas doentes (aumentando o espalhamento do vírus) e deixá-las sem tratamento. Já um valor alto da taxa de falso positivo (ou seja, 1 - especificidade) levaria a uma sobrecarga desnecessária do hospital.

O critério para classificar um indivíduo como positivo será dado pelo valor de corte na probabilidade que maximiza a soma da sensibilidade e especificidade. Embora o banco de dados processado apresente um tamanho reduzido (n = 501 pacientes) e seja desbalanceado (apenas 13.8% dos pacientes testaram positivo para o vírus), os classificadores apresentam um desempenho preditivo regular. No entanto, para esse contexto, é um desempenho abaixo do desejável.

A seguir faremos a leitura dos dados, a função para o ajuste dos modelos (Floresta Aleatória e boosting) que será utilizada em cada um dos bancos de dados obtidos por reamostragem (ver esse post para mais detalhes) e a replicação do procedimento para \(\text{N} = 10^3\) amostras.

set.seed(1234)

library(tidyverse)
library(tidymodels)
library(pROC)
library(doParallel)

db <- read.csv("covid-19.csv") %>%  # the data after some feature engineering
  mutate(result = factor(result))

registerDoParallel()

# ajustar hiperparâmetros do boosting
splits_cv <- vfold_cv(db, v = 10, strata = result)

fit_bst <- boost_tree(tree_depth = tune(), learn_rate = tune()) %>% 
            set_engine("xgboost") %>% 
            set_mode("classification")

bst_grid <- tune_grid(fit_bst, 
                      result ~ .,
                      splits_cv, 
                      grid = 30)

best <- bst_grid %>% 
         select_best("roc_auc")

# finaliza modelo boosting
fit_bst <- finalize_model(fit_bst, parameters = best)


# função para ajustar os modelos
ajustes <- function(df){
  
  tr   <- analysis(df)   
  test <- assessment(df) 
  
  # random forest
  rf <- rand_forest() %>% 
          set_engine("ranger") %>% 
          set_mode("classification") %>% 
          fit(result ~ ., data = tr)
  
  roc_rf <- bind_cols(obs = test$result,
                      predict(rf, new_data = test, type = "prob") %>% 
                      select(.pred_POS)) %>% 
              roc(obs, .pred_POS)

  best_rf <- which.max(roc_rf$sensitivities + roc_rf$specificities)
  
  # boosting
  bst <- fit(fit_bst,
             result ~ ., 
             data = tr)
  
  roc_bst <- bind_cols(obs = test$result,
                       predict(bst, new_data = test, type = "prob") %>% 
                       select(.pred_POS)) %>% 
              roc(obs, .pred_POS)
  
  best_bst <- which.max(roc_bst$sensitivities + roc_bst$specificities)
  
  # medidas de desempenho
  tibble(modelo = "floresta aleatória", 
         auc = roc_rf$auc[[1]], 
         best = roc_rf$thresholds[best_rf],
         sensibilidade = roc_rf$sensitivities[best_rf], 
         especificidade = roc_rf$specificities[best_rf]) %>% 
  bind_rows(tibble(modelo = "boosting", 
                   auc = roc_bst$auc[[1]], 
                   best = roc_bst$thresholds[best_bst],
                   sensibilidade = roc_bst$sensitivities[best_bst], 
                   especificidade = roc_bst$specificities[best_bst]))

}


# método de reamostragem para estimar medidas
splits <- bootstraps(db, times = 10^3, strata = result)

resultados <- splits %>% 
                mutate(medidas = map(splits, ajustes))

As medidas de desempenho dos dois modelos são apresentadas na tabela abaixo.

resultados %>%
    unnest(medidas) %>%
    select(modelo:especificidade, -best) %>%
    pivot_longer(-modelo, names_to = "medidas") %>%
    group_by(modelo, medidas) %>%
    summarise(media = 100*mean(value),
              lim_inferior = 100*quantile(value, .025),
              lim_superior = 100*quantile(value, .975))  

modelo	medidas	media	lim_inferior	lim_superior
boosting	auc	86.5	79.3	92.0
boosting	especificidade	77.8	61.9	91.4
boosting	sensibilidade	85.0	66.7	100.0
floresta aleatória	auc	88.1	82.7	93.2
floresta aleatória	especificidade	80.3	64.6	92.3
floresta aleatória	sensibilidade	86.1	70.0	100.0

Caso tenha alguma crítica, sugestão ou comentário, me envie uma mensagem!