Curvas ROC y áreas bajo la curva

library("pROC")
library("randomForest")

set.seed(420)
num_samples <- 100

weight <- sort(rnorm(n = num_samples, mean = 172, sd = 29))

obese <- ifelse(test = (runif(n = num_samples) < (rank(weight)/100)),
                yes = 1, no = 0)


plot(x = weight, y = obese)

glm_fit <- glm(obese ~ weight, family = "binomial")
plot(x = weight, y = obese)
lines(weight, glm_fit$fitted.values)

El primer resultado nos dice cuántos resultados fueron clasificados como obsesos y cuantos como no obesos.

Tambien reporta el area bajo la curva

par(pty = "s") # square
roc(obese, glm_fit$fitted.values, plot = TRUE, legacy.axes = TRUE,
    percent = TRUE, xlab = "Porcentaje Falsos positivos",
    ylab = "Porcentaje verdaderos postivios", col = "#377eb8", lwd = 2,
    print.auc = TRUE)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## 
## Call:
## roc.default(response = obese, predictor = glm_fit$fitted.values,     percent = TRUE, plot = TRUE, legacy.axes = TRUE, xlab = "Porcentaje Falsos positivos",     ylab = "Porcentaje verdaderos postivios", col = "#377eb8",     lwd = 2, print.auc = TRUE)
## 
## Data: glm_fit$fitted.values in 45 controls (obese 0) < 55 cases (obese 1).
## Area under the curve: 82.91%

roc_info <- roc(obese, glm_fit$fitted.values, legacy.axes = TRUE)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

roc_df <- data.frame(tpp =roc_info$sensitivities * 100,
                     fpp = (1- roc_info$specificities) * 100,
                     thresholds = roc_info$thresholds)

# Auc parcial
# Util cuando queremos enfocarnos  en una parte de la curva ROC que unicamente permite
# pequeños numeros de falsos positivos

par(pty = "s") # square
roc(obese, glm_fit$fitted.values, plot = TRUE, legacy.axes = TRUE,
    percent = TRUE, xlab = "Porcentaje Falsos positivos",
    ylab = "Porcentaje verdaderos postivios", col = "#377eb8", lwd = 2,
    print.auc = TRUE, print.auc.x =45, partial.auc = c(100, 90), # en terminos de especificidad
    auc.polygon = TRUE, auc.polygon.col = "#377eb850")

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## 
## Call:
## roc.default(response = obese, predictor = glm_fit$fitted.values,     percent = TRUE, plot = TRUE, legacy.axes = TRUE, xlab = "Porcentaje Falsos positivos",     ylab = "Porcentaje verdaderos postivios", col = "#377eb8",     lwd = 2, print.auc = TRUE, print.auc.x = 45, partial.auc = c(100,         90), auc.polygon = TRUE, auc.polygon.col = "#377eb850")
## 
## Data: glm_fit$fitted.values in 45 controls (obese 0) < 55 cases (obese 1).
## Partial area under the curve (specificity 100%-90%): 4.727%

rfmodel <- randomForest(factor(obese) ~ weight)

par(pty = "s") # square
roc(obese, glm_fit$fitted.values, plot = TRUE, legacy.axes = TRUE,
    percent = TRUE, xlab = "Porcentaje Falsos positivos",
    ylab = "Porcentaje verdaderos postivios", col = "#377eb8", lwd = 2,
    print.auc = TRUE)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## 
## Call:
## roc.default(response = obese, predictor = glm_fit$fitted.values,     percent = TRUE, plot = TRUE, legacy.axes = TRUE, xlab = "Porcentaje Falsos positivos",     ylab = "Porcentaje verdaderos postivios", col = "#377eb8",     lwd = 2, print.auc = TRUE)
## 
## Data: glm_fit$fitted.values in 45 controls (obese 0) < 55 cases (obese 1).
## Area under the curve: 82.91%

plot.roc(obese, rfmodel$votes[,1], percent=TRUE, col="#4daf4a",lwd= 2,
         print.auc =TRUE, add=TRUE,print.auc.y = 40)

## Setting levels: control = 0, case = 1

## Setting direction: controls > cases