Este tutorial é baseado nos exemplos do livro An Introduction to Statistical Learning, de Gareth James, Daniela Witten, Trevor Hastie e Rob Tibshirani (https://www.statlearning.com/) e no tutorial ROC Curves and AUC for Models Used for Binary Classification, de Jacob Goldstein-Greenwood, da University of Virginia Library.
Aqui, voltamos a trabalhar com a simple idéia de usar de um conjunto
de validação (ou teste) para estimar as taxas de erro de teste
resultantes do ajuste de vários modelos lineares ao conjunto de dados
Auto, que contém \(392\)
observações.
Antes de começarmos, usamos a função set.seed() para
definir uma semente para o gerador de números aleatórios do
R. Geralmente, é uma boa ideia usar set.seed()
ao realizar uma análise como a validação cruzada que contém um elemento
de aleatoriedade, para que os resultados obtidos possam ser reproduzidos
precisamente em um momento posterior.
Como já fizemos anteriormente, começamos usando a função
sample() para dividir o conjunto de observações em duas
metades, selecionando um subconjunto aleatório de \(196\) observações entre \(392\) (\(50\%\)), que servirá como nosso
conjunto de treinamento:
library(ISLR2)
head(Auto)
mpg cylinders displacement horsepower weight acceleration year origin
1 18 8 307 130 3504 12.0 70 1
2 15 8 350 165 3693 11.5 70 1
3 18 8 318 150 3436 11.0 70 1
4 16 8 304 150 3433 12.0 70 1
5 17 8 302 140 3449 10.5 70 1
6 15 8 429 198 4341 10.0 70 1
name
1 chevrolet chevelle malibu
2 buick skylark 320
3 plymouth satellite
4 amc rebel sst
5 ford torino
6 ford galaxie 500
set.seed(1)
train <- sample(392, 196)
No script acima, usamos um atalho no comando
sample(); consulte ?sample para entender o que
o script está fazendo.
Em seguida, usamos a opção subset em lm()
para ajustar uma regressão linear usando apenas as observações
correspondentes ao nosso conjunto de treinamento:
lm.fit <- lm(mpg ~ horsepower, data = Auto, subset = train)
Agora usamos a função predict() para estimar a resposta
para todas as \(392\) observações e
usamos a função mean() para calcular erro quadrático
médio do modelo baseando-se nas observações no conjunto de
validação ou teste. Note que o índice -train abaixo
seleciona apenas as observações que não estão no conjunto de
treinamento:
attach(Auto)
mean((mpg - predict(lm.fit, Auto))[-train]^2)
[1] 23.26601
Vemos que o erro de teste estimado para nosso modelo de regressão linear é \(23.27\).
Podemos usar a função poly() para estimar o erro de
teste para uma regressão quadrática e uma regressão
cúbica:
lm.fit2 <- lm(mpg ~ poly(horsepower, 2), data = Auto,
subset = train)
mean((mpg - predict(lm.fit2, Auto))[-train]^2)
[1] 18.71646
lm.fit3 <- lm(mpg ~ poly(horsepower, 3), data = Auto,
subset = train)
mean((mpg - predict(lm.fit3, Auto))[-train]^2)
[1] 18.79401
Essas taxas de erro são \(18.72\) e \(18.79\), respectivamente.
Claro que se escolhermos um conjunto de treinamento diferente (e consequentemente um conjunto de validação ou teste diferente), obteremos erros diferentes:
set.seed(2)
train <- sample(392, 196)
lm.fit <- lm(mpg ~ horsepower, subset = train)
mean((mpg - predict(lm.fit, Auto))[-train]^2)
[1] 25.72651
lm.fit2 <- lm(mpg ~ poly(horsepower, 2), data = Auto,
subset = train)
mean((mpg - predict(lm.fit2, Auto))[-train]^2)
[1] 20.43036
lm.fit3 <- lm(mpg ~ poly(horsepower, 3), data = Auto,
subset = train)
mean((mpg - predict(lm.fit3, Auto))[-train]^2)
[1] 20.38533
Usando esta divisão das observações em um conjunto de treinamento e um conjunto de validação/teste, encontramos que as taxas de erro do conjunto de validação/teste para os modelos com termos lineares, quadráticos e cúbicos são \(25.73\), \(20.43\) e \(20.39\), respectivamente.
Esses resultados são consistentes com nossas descobertas anteriores:
um modelo que prevê mpg usando uma função quadrática de
horsepower tem um desempenho melhor do que um modelo que
envolve apenas uma função linear de horsepower, e há pouca
evidência a favor de um modelo que usa uma função cúbica de
horsepower—simples, não? Até agora, nada de novo…
A métrica LOOCV pode ser calculada automaticamente para qualquer
modelo linear generalizado usando as funções glm() e
cv.glm(). Por exemplo, usamos a função glm()
em um tutorial anterior para realizar regressão logística passando o
argumento family = "binomial". Mas se usarmos
glm() para ajustar um modelo sem passar o argumento
family, ele realizará regressão linear, assim como a função
lm(). Por exemplo:
glm.fit <- glm(mpg ~ horsepower, data = Auto)
coef(glm.fit)
(Intercept) horsepower
39.9358610 -0.1578447
e…
lm.fit <- lm(mpg ~ horsepower, data = Auto)
coef(lm.fit)
(Intercept) horsepower
39.9358610 -0.1578447
…produzem modelos de regressão linear idênticos. Neste tutorial,
realizaremos a regressão linear usando a função glm() em
vez da função lm() porque a primeira pode ser usada em
conjunto com a função cv.glm(), do pacote boot
(e cv, por sinal, significa cross-validation):
library(boot)
glm.fit <- glm(mpg ~ horsepower, data = Auto)
cv.err <- cv.glm(Auto, glm.fit)
cv.err$delta
[1] 24.23151 24.23114
A função cv.glm() produz uma lista com vários
componentes. Os dois números no vetor delta contêm os
resultados da validação cruzada (cross-validation). Neste caso,
os números são idênticos (até duas casas decimais) e correspondem à
estatística LOOCV. Abaixo, discutimos uma situação em que os dois
números diferem. De qualquer forma, nossa estimativa de validação
cruzada para o erro de teste é aproximadamente \(24.23\).
Podemos repetir este procedimento para ajustes polinomiais cada vez
mais complexos. Para automatizar o processo, usamos a função
for() para iniciar um “for loop”
(for (i in 1:10) {...}) que ajusta iterativamente
regressões polinomiais para polinômios de ordem \(i=1\) a \(i=10\), calcula o erro de validação cruzada
associado a cada \(i\) e o armazena no
\(i\)-ésimo elemento do vetor
cv.error. Começamos inicializando o vetor:
cv.error <- rep(0, 10)
for (i in 1:10) {
glm.fit <- glm(mpg ~ poly(horsepower, i), data = Auto)
cv.error[i] <- cv.glm(Auto, glm.fit)$delta[1]
}
cv.error
[1] 24.23151 19.24821 19.33498 19.42443 19.03321 18.97864 18.83305 18.96115
[9] 19.06863 19.49093
plot(cv.error, type="b", col="red", lwd=2, pch=19)
Vemos uma queda acentuada no erro médio quádratico
(cv.error) de teste estimado entre os ajustes linear e
quadrático, mas depois nenhuma melhora clara ao usar polinômios de ordem
superior…
A função cv.glm() também pode ser usada para implementar
a validação cruzada \(k\)-fold. Abaixo, usamos \(k=10\), uma escolha comum para \(k\), no conjunto de dados
Auto. Mais uma vez, definimos uma semente aleatória e
inicializamos um vetor no qual armazenaremos os erros de validação
cruzada correspondentes aos ajustes polinomiais de ordens \(1\) a \(10\):
set.seed(17)
cv.error.10 <- rep(0, 10)
for (i in 1:10) {
glm.fit <- glm(mpg ~ poly(horsepower, i), data = Auto)
cv.error.10[i] <- cv.glm(Auto, glm.fit, K = 10)$delta[1]
}
cv.error.10
[1] 24.27207 19.26909 19.34805 19.29496 19.03198 18.89781 19.12061 19.14666
[9] 18.87013 20.95520
plot(cv.error.10, type="b", col="red", lwd=2, pch=19)
Observe que o tempo de computação do \(k\)-fold é menor do que o da LOOCV—faz sentido, certo? De qualquer forma, ainda vemos pouca evidência de que usar termos polinomiais cúbicos ou de ordem superior leve a um erro de teste menor do que simplesmente usar um ajuste quadrático…
Vamos agora nos aprofundar um pouco mais no processo de cálculo e visualização da Área Sob a Curva Característica de Operação do Receptor (AUC-ROC). A Curva ROC (e AUC) é uma ferramenta fundamental para avaliar o desempenho de modelos de classificação binária (e na verdade, existem extensões desse método que permitem utilizá-lo para validação de classificações multiclasse, vamos ver um exemplo disso mais abaixo).
Já que já utilizamos o pacote ROCR para criar curvas ROC
no tutorial anterior, vamos focar em um pacote diferente dessa vez:
pROC. (Note que o pacote pROC, por sua vez,
utiliza o pacote ggplot2 para plotar os resultados, então
preciamos dele também.)
Começamos carregando os pacotes e dados simulados que usaremos nesta parte do tutorial, produzidos pela Universidade de Virginia:
library(pROC)
Type 'citation("pROC")' for a citation.
Attaching package: 'pROC'
The following objects are masked from 'package:stats':
cov, smooth, var
library(ggplot2)
Attaching package: 'ggplot2'
The following object is masked from 'Auto':
mpg
sim_dat <- read.csv('http://static.lib.virginia.edu/statlab/materials/data/roc_sim_dat.csv', header = T)
Usamos pROC para calcular a taxa de verdadeiros
positivos dividida pela taxa de falsos positivos através
de uma variedade de limiares de classificação (note os comentários no
script abaixo):
sim_roc <- roc(response = sim_dat$actual_outcome,
predictor = sim_dat$predicted_prob_of_Yes,
levels = c('No', 'Yes'))
Setting direction: controls < cases
# response: vetor dos resultados reais
# predictor: vetor de probabilidades, uma para cada observação
# levels: valores da resposta a serem considerados como "No"/0/"Falha"/etc. e "Yes"/1/"Sucesso"/etc., respectivamente
# Plotar a curva ROC:
# `ggroc()` imprime a curva ROC usando `ggplot2` (você pode imprimir a curva usando gráficos de base do R com `plot.roc()` se preferir)
# `legacy.axes = TRUE` garante que o eixo x seja exibido como 'Taxa de Falsos Positivos (TFP)' (de 0 a 1) em vez de 'Especificidade' (1-TFP) de (1 a 0)
ggroc(sim_roc, legacy.axes = TRUE) +
labs(x = 'Taxa de Falsos Positivos', y = 'Taxa de Verdadeiros Positivos', title = 'Curva ROC Simulada')
Como sabemos, a AUC pode auxiliar na comparação do desempenho geral
de modelos utilizados para classificação binária. Por exemplo, considere
as duas regressões logísticas a seguir baseadas no conjunto de dados
titanic3, disponibilizado por Frank Harrell, da
Universidade de Vanderbilt. Esses dados contém informações sobre se os
passageiros do Titanic sobreviveram ao naufrágio ou morreram no mar,
assim como dados sobre a classe da cabine do passageiro (1ª, 2ª ou 3ª),
idade e sexo:
A primeira regressão logística prediz a sobrevivência
(survived: \(1\)=sobreviveu ou \(0\)=morreu) a partir da classe da cabine do
passageiro (pclass: 1ª, 2ª ou 3ª).
A segunda regressão logística prediz a sobrevivência a partir da classe da cabine do passageiro, idade do passageiro e sexo do passageiro.
Abaixo, selecionamos nosso subconjunto com as variáveis de nosso
interesse (note o uso do operador |> no script
abaixo—ele é basicamente equivalente ao %>%):
library(Hmisc)
Attaching package: 'Hmisc'
The following objects are masked from 'package:base':
format.pval, units
Hmisc::getHdata(titanic3)
titanic3 <- titanic3[, c('survived', 'pclass', 'sex', 'age')] |> na.omit()
knitr::kable(head(titanic3, n = 4), align = rep('l', ncol(titanic3)))
| survived | pclass | sex | age |
|---|---|---|---|
| 1 | 1st | female | 29.0000 |
| 1 | 1st | male | 0.9167 |
| 0 | 1st | female | 2.0000 |
| 0 | 1st | male | 30.0000 |
# Primeiro modelo
titanic_mod1 <- glm(survived ~ pclass, data = titanic3, family = binomial)
summary(titanic_mod1)
Call:
glm(formula = survived ~ pclass, family = binomial, data = titanic3)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.5638 0.1234 4.568 0.00000493 ***
pclass2nd -0.8024 0.1754 -4.574 0.00000479 ***
pclass3rd -1.6021 0.1599 -10.019 < 0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1414.6 on 1045 degrees of freedom
Residual deviance: 1305.9 on 1043 degrees of freedom
AIC: 1311.9
Number of Fisher Scoring iterations: 4
# Segundo modelo
titanic_mod2 <- glm(survived ~ pclass + sex + age, data = titanic3, family = binomial)
summary(titanic_mod2)
Call:
glm(formula = survived ~ pclass + sex + age, family = binomial,
data = titanic3)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.522074 0.326702 10.781 < 0.0000000000000002 ***
pclass2nd -1.280570 0.225538 -5.678 0.0000000136 ***
pclass3rd -2.289661 0.225802 -10.140 < 0.0000000000000002 ***
sexmale -2.497845 0.166037 -15.044 < 0.0000000000000002 ***
age -0.034393 0.006331 -5.433 0.0000000556 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1414.62 on 1045 degrees of freedom
Residual deviance: 982.45 on 1041 degrees of freedom
AIC: 992.45
Number of Fisher Scoring iterations: 4
Podemos usar cada modelo para gerar uma probabilidade de
sobrevivência para cada passageiro. Com as probabilidades de cada
modelo, assim como os resultados reais de sobrevivência, usaremos o
pacote pROC como antes para gerar as taxas de
verdadeiros positivos e as taxas de falsos positivos em
uma variedade de limiares de classificação. Aqui, utilizamos o pacote
gridExtrapara plotar dois gráficos do ggroc()
lado a lado:
library(gridExtra)
# Curva ROC para o primeiro modelo
probs_mod1 <- predict(titanic_mod1, type = 'response')
roc_mod1 <- roc(response = titanic3$survived, predictor = probs_mod1)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
mod1_plot <- ggroc(roc_mod1, legacy.axes = TRUE) +
labs(x = 'Taxa de Falsos Positivos', y = 'Taxa de Verdadeiros Positivos',
title = 'Curva ROC\nsurvived ~ pclass') +
annotate('text', x = .5, y = .5, label = paste0('AUC: ', round(auc(roc_mod1), digits = 2)))
# Curva ROC para o segundo modelo
probs_mod2 <- predict(titanic_mod2, type = 'response')
roc_mod2 <- roc(response = titanic3$survived, predictor = probs_mod2)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
mod2_plot <- ggroc(roc_mod2, legacy.axes = TRUE) +
labs(x = 'Taxa de Falsos Positivos', y = 'Taxa de Verdadeiros Positivos',
title = 'Curva ROC\nsurvived ~ pclass + age + sex') +
annotate('text', x = .5, y = .5, label = paste0('AUC: ', round(auc(roc_mod2), digits = 2)))
grid.arrange(mod1_plot, mod2_plot, nrow = 1)
Vemos que com a adição de idade e sexo como preditores, a AUC aumenta
em cerca de 25%. Neste caso, o segundo modelo, que inclui
age e sex, é definitivamente melhor (um
resultado dificilmente surpreendente, dado o que sabemos sobre o
naufrágio e seus sobreviventes).
Retornando mais uma vez à interpretação probabilística da AUC: se pegássemos aleatoriamente uma pessoa que sobreviveu e uma pessoa que morreu dos dados, o segundo modelo seria muito mais propenso do que o primeiro a classificar a probabilidade de sobrevivência da pessoa que viveu como maior do que a da pessoa que morreu.
Lembre-se que, em geral, quanto maior a AUC, melhor é o modelo de classificação. No entanto, buscar uma AUC cada vez maior pode levar ao sobreajuste (overfitting) do modelo. Os conjuntos de dados estão repletos de idiossincrasias (peculiaridades específicas e potencialmente irrelevantes presentes em um determinado conjunto de dados), e inserir preditor após preditor no seu modelo na tentativa de extrair um valor de AUC de \(0,87\) em vez de \(0,86\), por exemplo, pode degradar a capacidade de generalização e a utilidade prática do modelo—por isso a importância da validação de modelos!
Como vimos anteriormente, enquanto as curvas ROC podem ser utilizadas de forma diretas para classificação binária, estendê-las para classificação multiclasse apresenta desafios adicionais. A seguir, exploraremos como gerar e interpretar curvas ROC para classificação multiclasse.
Sem entrar em detalhes, existem diversas abordagens comuns para estender o uso das curvas ROC à classificação multiclasse (pesquise a respeito se for de seu interesse):
Aqui, vamos manipular o conjunto de dados titanic3 para
introduzir uma nova categoria na variável survived. De
forma aleatória, designaremos algumas observações como
"Almost_died". Com essa alteração, criaremos um problema de
classificação multiclasse com \(3\)
classes: "Died", "Almost_died" e
"Survived".
# Converter 'survived' para fator com rótulos descritivos: 0 = "Died", 1 = "Survived"
titanic3$survived <- factor(titanic3$survived, levels = c(0, 1), labels = c("Died", "Survived"))
# Criar a nova classe "almost_died":
# 1. Gerar índices para 1/3 dos dados
n_rows <- nrow(titanic3)
num_almost_died <- round(n_rows / 3)
set.seed(123) # Para reprodutibilidade
indices_almost_died <- sample(1:n_rows, num_almost_died)
# 2. Adicionar "Almost_died" como um novo nível ao fator 'survived'
# Isso é crucial para que você possa atribuir a nova categoria
levels(titanic3$survived) <- c(levels(titanic3$survived), "Almost_died")
# 3. Atribuir "almost_died" aos índices selecionados
titanic3$survived[indices_almost_died] <- "Almost_died"
# Verificar a contagem de cada categoria em 'survived'
print(table(titanic3$survived))
Died Survived Almost_died
422 275 349
Já que a variável survived agora possui \(3\) classes, a regressão logística binária
implementada com a função glm() (e
family = binomial) não é mais apropriada. Para lidar com
este problema de classificação multiclasse, utilizaremos uma máquina de
vetores de suporte, ou support vector machine (SVM), com kernel
radial, como fizemos no tutorial anterior:
library(caret) #instale se necessário para usar a função `createDataPartition()`
Loading required package: lattice
Attaching package: 'lattice'
The following object is masked from 'package:boot':
melanoma
# Criar índices para a partição dos dados de treinamento e teste:
set.seed(123)
train_indices <- createDataPartition(titanic3$survived, p = 0.8, list = FALSE)
# Criar os conjuntos de treinamento e teste:
titanic_train <- titanic3[train_indices, ]
titanic_test <- titanic3[-train_indices, ]
library(e1071)
Attaching package: 'e1071'
The following object is masked from 'package:Hmisc':
impute
svm_model_titanic <- svm(survived ~ pclass + sex + age, data = titanic_train,
kernel = "radial", cost = 1, gamma = 0.1,
probability = TRUE) # Necessário para a análise ROC a seguir
Agora, usamos nossa SVM para fazer previsões utilizando o conjunto de teste e aproveitamos para gerar uma matriz de confusão:
predictions_svm_test <- predict(svm_model_titanic, newdata = titanic_test, probability = TRUE)
# Avaliar o desempenho do modelo no conjunto de teste com uma matriz de confusão:
confusion_matrix_test <- table(predicted = predictions_svm_test, actual = titanic_test$survived)
print(confusion_matrix_test)
actual
predicted Died Survived Almost_died
Died 70 18 47
Survived 13 36 21
Almost_died 1 1 1
Vemos que o desempenho da nossa SVM não é dos melhores, o que já era
esperado, visto que a categoria Almost_died foi introduzida
de forma aleatória nos dados. No entanto, o objetivo aqui é puramente
demonstrativo: utilizá-la como exemplo para computar a AUC-ROC em um
cenário de classificação multiclasse. Realizamos essa computação com o
script a seguir (que é um pouco mais complexo que nosso exemplo
anterior de classificação binária). Para isso, adotamos a abordagem
Um-contra-o-Resto (One-vs-Rest):
library(pROC)
# Obter as probabilidades previstas para o conjunto de teste (necessárias para o calculo da AUC-ROC):
probabilities_svm_test <- attr(predictions_svm_test, "probabilities")
# Inicializar uma lista para armazenar os objetos ROC para cada classe:
roc_objects <- list()
# Inicializar um vetor para armazenar os valores AUC para cada classe:
auc_values <- numeric(length(levels(titanic_test$survived)))
names(auc_values) <- levels(titanic_test$survived)
for (i in seq_along(levels(titanic_test$survived))) {
class_name <- levels(titanic_test$survived)[i]
# Criar rótulos binários: 1 se a observação pertencer à classe atual, 0 caso contrário
binary_labels <- as.numeric(titanic_test$survived == class_name)
# A coluna de probabilidades para a classe atual
class_probabilities <- probabilities_svm_test[, class_name]
# Calcular a curva ROC para a classe atual
roc_obj <- roc(response = binary_labels, predictor = class_probabilities, levels = c(0, 1))
# Armazenar o objeto ROC e o valor AUC
roc_objects[[class_name]] <- roc_obj
auc_values[class_name] <- auc(roc_obj)
}
Setting direction: controls < cases
Setting direction: controls < cases
Setting direction: controls < cases
round(auc_values, digits = 2)
Died Survived Almost_died
0.69 0.80 0.49
Os valores de AUC obtidos são coerentes: o AUC de \(0,49\) para a classe
Almost_died indica que a SVM não se sai melhor que uma
classificação aleatória. Isso é esperado, visto que criamos essa classe
de forma intencionalmente aleatória.
A seguir plotamos as curvas ROC:
# Plotar as curvas ROC individuais (OvR)
plot(roc_objects[[1]], col = "red", main = "Curvas ROC por Classe (One-vs-Rest)",
xlab = "Taxa de Falsos Positivos (1 - Especificidade)",
ylab = "Taxa de Verdadeiros Positivos (Sensibilidade)")
lines(roc_objects[[2]], col = "blue")
lines(roc_objects[[3]], col = "green")
legend("bottomright",
legend = levels(titanic_test$survived),
col = c("red", "blue", "green"),
lwd = 2)
grid() # Adiciona uma grade ao gráfico para melhor visualização
Como mencionamos acima, existem diversas abordagens comuns para estender o uso das curvas ROC à classificação multiclasse. A seguir, calculamos as chamadas macro-média e micro-média associadas a nossa SVM. Note que essas métricas tendem a ser aplicadas exclusivamente à AUC:
# Macro-média AUC (média simples das AUCs por classe):
mean(auc_values)
[1] 0.6601806
# Micro-média AUC (requer a matriz de probabilidades e o vetor de classes reais):
multiclass.roc(response = titanic_test$survived, predictor = probabilities_svm_test)
Call:
multiclass.roc.default(response = titanic_test$survived, predictor = probabilities_svm_test)
Data: multivariate predictor probabilities_svm_test with 3 levels of titanic_test$survived: Survived, Died, Almost_died.
Multi-class area under the curve: 0.6591
Note que a micro-média AUC é calculada internamente pela função
multiclass.roc(), do pacote pROC, e difere (ao
menos um pouco) da macro-média AUC.
Agora é sua vez de validar alguns modelos. Você pode fazer isso
utilizando seus próprios dados, dados disponíveis na internet ou em
pacotes do R. Outra opção é utilizar os dados sobre
desmatamento municipal em 2004, disponíveis em https://thaleswest.wixsite.com/home/tutorials (note que
o site contém a descrição das variáveis deste conjunto de dados).
Já que esse é nosso último tutorial, queremos que você utilize o que aprendeu ao longo desse curso. Siga os seguintes passos: