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 ( tiagoms1@insper.edu.br)

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(doParallel)

db <- read.csv("covid-19.csv") # the data after some feature engineering

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))  

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