Enunciado

Construir una función R que calcule el área bajo la Curva Operativa Característica (AUC). La función admitirá como argumentos de entrada, al menos los siguientes:

Un vector con las probabilidades estimadas de pertenencia a una determinada clase.
Un vector con la clase real a la que pertenece cada caso.

Se valorará que la función sea también capaz de dibujar la curva ROC.

Utilizar la función con algunos de los modelos vistos en el Tema 3, sobre datos de libre elección. Comparar los resultados obtenidos con los que se obtienen con alguna de las funciones R que calculan el criterio AUC.

Datos: Moneyball

A principios de la década de 2000, Billy Beane y Paul DePodesta trabajaron para los Oakland Athletics revolucionando el mundo del deporte y de la estadística deportiva. Este conjunto de datos contiene parte de la información que estaba disponible para Beane y DePodesta a principios de la década de 2000. Los datos son una adecuación de los que se encuentran en Kaggle:

https://www.kaggle.com/wduckett/moneyball-mlb-stats-19622012

Las variables son las siguientes:

  - Equipo
  - Liga
  - Año
  - Carreras anotadas (RS)
  - Ejecuciones permitidas (RA)
  - Victorias (W)
  - Porcentaje en base (OBP)
  - Porcentaje de slugging (SLG)
  - Promedio de bateo (BA)
  - Juegos jugados (G)
  - Playoffs (variable binaria a clasificar)

Team	League	Year	RS	RA	W	OBP	SLG	BA	G	Playoffs
ARI	NL	2012	734	688	81	0.328	0.418	0.259	162	0
ATL	NL	2012	700	600	94	0.320	0.389	0.247	162	1
BAL	AL	2012	712	705	93	0.311	0.417	0.247	162	1
BOS	AL	2012	734	806	69	0.315	0.415	0.260	162	0
CHC	NL	2012	613	759	61	0.302	0.378	0.240	162	0
CHW	AL	2012	748	676	85	0.318	0.422	0.255	162	0

summary(datos)

##      Team              League               Year            RS        
##  Length:1232        Length:1232        Min.   :1962   Min.   : 463.0  
##  Class :character   Class :character   1st Qu.:1977   1st Qu.: 652.0  
##  Mode  :character   Mode  :character   Median :1989   Median : 711.0  
##                                        Mean   :1989   Mean   : 715.1  
##                                        3rd Qu.:2002   3rd Qu.: 775.0  
##                                        Max.   :2012   Max.   :1009.0  
##        RA               W              OBP              SLG        
##  Min.   : 472.0   Min.   : 40.0   Min.   :0.2770   Min.   :0.3010  
##  1st Qu.: 649.8   1st Qu.: 73.0   1st Qu.:0.3170   1st Qu.:0.3750  
##  Median : 709.0   Median : 81.0   Median :0.3260   Median :0.3960  
##  Mean   : 715.1   Mean   : 80.9   Mean   :0.3263   Mean   :0.3973  
##  3rd Qu.: 774.2   3rd Qu.: 89.0   3rd Qu.:0.3370   3rd Qu.:0.4210  
##  Max.   :1103.0   Max.   :116.0   Max.   :0.3730   Max.   :0.4910  
##        BA               G         Playoffs
##  Min.   :0.2140   Min.   :158.0   0:988   
##  1st Qu.:0.2510   1st Qu.:162.0   1:244   
##  Median :0.2600   Median :162.0           
##  Mean   :0.2593   Mean   :161.9           
##  3rd Qu.:0.2680   3rd Qu.:162.0           
##  Max.   :0.2940   Max.   :165.0

Con eso en mente, nos plantearíamos la pregunta: ¿Llega un equipo a los playoffs?.

Modelos de clasificación

Partición entrenamiento-test

Particionaremos los datos en un conjunto entrenamiento (70%) y otro test (30%) sobre el que evaluaremos nuestros modelos.

set.seed("75895206")
n<- nrow(datos)
indin<- 1:n
nent<-ceiling(0.7*n)
ntest<- n-nent
indient<- sort(sample(indin,nent))
inditest<- setdiff(indin,indient)

#Posición de la variable dependiente en la tabla:
vd = 11

Para la clasificación de los datos emplearemos los siguientes modelos:

Naive Bayes (Gauss)

library(e1071) 
modeloNB_G    <- naiveBayes(Playoffs ~ ., data = datos[indient,])
preditest_NBG <- predict(modeloNB_G,datos[inditest,-vd])
probs_NBG     <- predict(modeloNB_G,datos[inditest,-vd],type="raw")[,2]

Naive Bayes (Poisson)

library(naivebayes)
modeloNB_P    <- naive_bayes(Playoffs ~ ., data = datos[indient,],usepoisson=TRUE)
preditest_NBP <- predict(modeloNB_P,datos[indient,-vd])
probs_NBP     <- predict(modeloNB_P,datos[inditest,-vd],type="prob")[,2]

Regresión logística (GLM)

modelo_GLM  <- glm(Playoffs~.,datos[indient,], family=binomial, maxit=10000)
probs_GLM   <- predict(modelo_GLM,datos[inditest,], type="response")

Para el desarrollo del ejercicio se tomará el modelo de regresión logística, al final se verá una comparación de los diferentes modelos.

Curvas ROC

Una Curva ROC (acrónimo de Receiver Operating Characteristic, o Característica Operativa del Receptor) es una representación gráfica de la sensibilidad frente a la especificidad para un sistema clasificador binario según se varía el umbral de discriminación.

Otra interpretación de este gráfico es la representación de la razón o ratio de verdaderos positivos (TVP = Tasa de Verdaderos Positivos) frente a la razón o ratio de falsos positivos (FPR = Razón de Falsos Positivos) también según se varía el umbral de discriminación (valor a partir del cual decidimos que un caso es un positivo).

AUC

El área bajo una curva ROC (AUC) se usa comúnmente en el aprendizaje automático para resumir el rendimiento de un modelo predictivo con un solo valor. Se utiliza en el análisis de clasificación con el fin de determinar específicamente de los modelos utilizados predice las clases mejor.

library(ROCR)

library(ROCR)

real <- datos[inditest,]$Playoffs

prediobj <-prediction(probs,real)
perf <-  performance(prediobj, "tpr","fpr")

plot(perf,
     main = "Curva ROC",
     xlab="Tasa de falsos positivos", 
     ylab="Tasa de verdaderos positivos")
abline(a=0,b=1,col="blue",lty=2)
grid()
auc <- as.numeric(performance(prediobj,"auc")@y.values)
legend("bottomright",legend=paste(" AUC =",round(auc,4)))

library(pROC)

library(pROC)

roc <- roc(real,probs)

plot(roc,col="red",lwd=2,main="ROC test")
legend("bottomright",legend=paste("AUC=",round(auc(roc),4)))

Curva ROC script

#MATRIZ DE CONFUSIÓN:

label <- datos[inditest,]$Playoffs

#Sensitividad vs Especificidad
TFP <- NULL #1-especificidad
TVP <- NULL #sensibilidad

epsilon <- 0.001
breaks  <-  seq(0 + epsilon, 1 - epsilon, epsilon)
#breaks  <-  seq(0,1, epsilon)
n <- length(breaks)

for (i in 1:n) {
  preditest <- probs>breaks[i]  
  confutest <- table(label,preditest)
  diagonal <- diag(prop.table(confutest, 1))
  TFP[i] = 1 - diagonal[1]
  TVP[i] = diagonal[2]
}

x = c(1,TFP,0)
y = c(1,TVP,0)

plot(x,y,type="l",
         main="Curva ROC",
         xlab="Tasa de falsos positivos", 
         ylab="Tasa de verdaderos positivos",
         lwd = 2,col="cornflowerblue")
abline(a=0, b=1, col="blue", lty=2)

Para valores de \(\varepsilon\) pequeños y número de observaciones grandes, la función tarda en computar. En AUC_script.R se ha comprobado los tiempos de computo con ayuda de la librería bench y la función system_time().

Podemos comparar nuestra aproximación con la del paquete ROCR:

plot(perf,
     main="Curva ROC",
     xlab="Tasa de falsos positivos", 
     ylab="Tasa de verdaderos positivos",
     lwd = 2)
lines(x,y,type="l",
      lty=2,lwd=2,col="red")
abline(a=0,b=1,col="blue",lty=2)
legend(0.7,0.2,
       legend = c("ROCR","Aproximación"),
       col=c("black", "red"), lty=1:2, cex=0.8)

Función AUC

library(ROCR)

library(ROCR)
prediobj  <- prediction(probs,real)
AUC_ROCR  <- as.numeric(performance(prediobj,"auc")@y.values)
AUC_ROCR

## [1] 0.9660934

library(pROC)

library(pROC)
roc      <- roc(real,probs)
AUC_pROC <- auc(roc)
AUC_pROC

## Area under the curve: 0.9661

Adecuación del factor para los siguientes apartados:

label = factor(label,labels = c(0,1))
label = as.numeric(as.character(label))

Fórmulas de integración numérica

TFP <- NULL 
TVP <- NULL 
epsilon <- 0.001
breaks  <-  seq(0 + epsilon, 1 - epsilon, epsilon)
n <- length(breaks)

for (i in 1:n) {
  preditest <- probs>breaks[i]  
  confutest <- table(label,preditest)
  diagonal <- diag(prop.table(confutest, 1))
  TFP[i] = 1 - diagonal[1]
  TVP[i] = diagonal[2]
}

TFP = sort(TFP)
TVP = sort(TVP)

x=c(0,TFP,1)
y=c(0,TVP,1)

Fórmula del trapecio

\[ \int_a^b f(t)dt \approx (b-a)\frac{f(a)+f(b)}{2} \]

sum(diff(x)*(y[1:(length(y)-1)] + y[2:length(y)])/2,na.rm = T)

## [1] 0.9656781

Fórmula del rectángulo inferior

\[ \int_a^b f(t)dt \approx (b-a)f(a) \]

sum(diff(x)*y[-length(y)],na.rm = T)

## [1] 0.9626734

Fórmula del rectángulo superior

\[ \int_a^b f(t)dt \approx (b-a)f(b) \]

sum(diff(x)*y[-1],na.rm = T)

## [1] 0.9686828

Estadístico U Mann-Whitney

De Wikipedia sabemos que: \[AUC = \frac{U}{n_1n_2}\]

dónde \(U\) es el estadístico U de Mann-Whitney, también conocida como estadística de prueba de suma de rango de Wilcoxon.

Usando `wilcox.test()`

auc_wmw <- function(labels, scores){
  labels <- as.logical(labels)
  pos <- scores[labels]
  neg <- scores[!labels]
  U <- as.numeric(wilcox.test(pos, neg)$statistic)
  U/(length(pos) * length(neg))
}
auc_wmw(label, probs)

## [1] 0.9660934

Calculando el estadístico

auc_wmw2 <- function(labels, scores){
  labels <- as.logical(labels)
  n1 <- sum(labels)
  n2 <- sum(!labels)
  R1 <- sum(rank(scores)[labels])
  U1 <- R1 - n1 * (n1 + 1)/2
  U1/(n1 * n2)
}
auc_wmw2(label,probs)

## [1] 0.9660934

En la librería library(auRoc) podemos encontrar funciónes más desarrolladas.

Visión probabilística

La interpretación probabilística es que si elige aleatoriamente un caso positivo y un caso negativo, la AUC da la probabilidad de que el caso positivo supere al caso negativo según el clasificador.

auc_probability <- function(labels, scores, N=1e7){
  labels <- as.logical(labels)
  pos <- sample(scores[labels], N, replace=TRUE)
  neg <- sample(scores[!labels], N, replace=TRUE)
  # sum( (1 + sign(pos - neg))/2)/N # does the same thing
  (sum(pos > neg) + sum(pos == neg)/2) / N # give partial credit for ties
}
auc_probability(label,probs)

## [1] 0.9661445

Comparación de rangos

AUC_rank_comparison <- function(labels, scores){
  score_order <- order(scores, decreasing=TRUE)
  labels <- as.logical(labels[score_order])
  scores <- scores[score_order]
  pos_scores <- scores[labels]
  neg_scores <- scores[!labels]
  n_pos <- sum(labels)
  n_neg <- sum(!labels)
  M <- outer(sum(labels):1, 1:sum(!labels), 
             function(i, j) (1 + sign(pos_scores[i] - neg_scores[j]))/2)
  AUC <- mean(M)
  AUC
}
AUC_rank_comparison(label,probs)

## [1] 0.9660934

Versión compacta

funcion_AUC = function(labels,scores,
                       method="ROCR",
                       plot.ROC=FALSE,epsilon=0.001,N=1e7){
  
  TFP <- NULL 
  TVP <- NULL 
  epsilon <- epsilon
  breaks  <-  seq(0 + epsilon, 1 - epsilon, epsilon)
  n <- length(breaks)
  
  for (i in 1:n) {
    preditest <- scores>breaks[i]  
    confutest <- table(labels,preditest)
    diagonal <- diag(prop.table(confutest, 1))
    TFP[i] = 1 - diagonal[1]
    TVP[i] = diagonal[2]
  }
  
  x=c(0,sort(TFP),1)
  y=c(0,sort(TVP),1)
  
  if (method=="ROCR") {
    if (!require(ROCR)) install.packages("ROCR")
    library(ROCR)
    prediobj  <- prediction(scores,labels)
    AUC       <- as.numeric(performance(prediobj,"auc")@y.values)
  }
  
  if (method=="pROC") {
    if (!require(ROCR)) install.packages("pROC")
    library(pROC)
    roc <- roc(labels,scores)
    AUC <- auc(roc)
  }
  
  if (method=="integracion") {
    AUC <- sum(diff(x)*(y[1:(length(y)-1)] + y[2:length(y)])/2,na.rm = T)
  }
  
  label = factor(label,labels = c(0,1))
  label = as.numeric(as.character(label))
  
  if (method=="UMW") {
    labels <- as.logical(labels)
    pos <- scores[labels]
    neg <- scores[!labels]
    U   <- as.numeric(wilcox.test(pos, neg)$statistic)
    AUC <- U/(length(pos) * length(neg))
  }
  
  if (method=="prob") {
    labels <- as.logical(labels)
    pos <- sample(scores[labels], N, replace=TRUE)
    neg <- sample(scores[!labels], N, replace=TRUE)
    AUC <- (sum(pos > neg) + sum(pos == neg)/2) / N
  }
  
  if (method=="rangos"){
    score_order <- order(scores, decreasing=TRUE)
    labels <- as.logical(labels[score_order])
    scores <- scores[score_order]
    pos_scores <- scores[labels]
    neg_scores <- scores[!labels]
    n_pos <- sum(labels)
    n_neg <- sum(!labels)
    M <- outer(sum(labels):1, 1:sum(!labels), 
               function(i, j) (1 + sign(pos_scores[i] - neg_scores[j]))/2)
    AUC <- mean(M)
  }
  
  if(plot.ROC){
    plot(x,y,type="l",
         main="Curva ROC",
         xlab="Tasa de falsos positivos", 
         ylab="Tasa de verdaderos positivos",
         lwd = 2,col="cornflowerblue")
    abline(a=0, b=1, col="blue", lty=2)
    legend(0.7, 0.1,
           legend = paste("AUC=",round(AUC,4)))
  }
  
  return(AUC)
  
}

Comparación de los distintos modelos

tabla <-  matrix(NA,nrow = 6,ncol = 3)
tabla <- data.frame(tabla)
colnames(tabla) <-  c("NaiveBayes_G","NaiveBayes_P","GLM")
row.names(tabla) <- c("ROCR","pROC","integracion","UMW","prob","rangos")
tabla$NaiveBayes_G <- c(
  funcion_AUC(label,probs_NBG,method = "ROCR"),
  funcion_AUC(label,probs_NBG,method = "pROC"),
  funcion_AUC(label,probs_NBG,method = "integracion"),
  funcion_AUC(label,probs_NBG,method = "UMW"),
  funcion_AUC(label,probs_NBG,method = "prob"),  
  funcion_AUC(label,probs_NBG,method = "rangos")
)

tabla$NaiveBayes_P <- c(
  funcion_AUC(label,probs_NBP,method = "ROCR"),
  funcion_AUC(label,probs_NBP,method = "pROC"),
  funcion_AUC(label,probs_NBP,method = "integracion"),
  funcion_AUC(label,probs_NBP,method = "UMW"),
  funcion_AUC(label,probs_NBP,method = "prob"),  
  funcion_AUC(label,probs_NBP,method = "rangos")
)

tabla$GLM <- c(
  funcion_AUC(label,probs_GLM,method = "ROCR"),
  funcion_AUC(label,probs_GLM,method = "pROC"),
  funcion_AUC(label,probs_GLM,method = "integracion"),
  funcion_AUC(label,probs_GLM,method = "UMW"),
  funcion_AUC(label,probs_GLM,method = "prob"),  
  funcion_AUC(label,probs_GLM,method = "rangos")
)

tabla

	NaiveBayes_G	NaiveBayes_P	GLM
ROCR	0.9346297	0.9306234	0.9660934
pROC	0.9346297	0.9306234	0.9660934
integracion	0.9346052	0.9180184	0.9656781
UMW	0.9346297	0.9306234	0.9660934
prob	0.9347422	0.9307716	0.9660372
rangos	0.9346297	0.9306234	0.9660934

En este caso, la regresión logística parece ser el mejor clasificador.

Bibliografía

Wikipedia. https://es.wikipedia.org/wiki/Curva_ROC
Kaggle https://www.kaggle.com/
Introducción a las curvas ROC. http://ares.inf.um.es/00Rteam/pub/clas/II_ROC.html#el-making-of-de-las-curvas-roc
Calculating AUC: the area under a ROC Curve. https://www.r-bloggers.com/calculating-auc-the-area-under-a-roc-curve/
AUC Meets the Wilcoxon-Mann-Whitney U-Statistic. https://blog.revolutionanalytics.com/2017/03/auc-meets-u-stat.html
CRAN: ROCR https://cran.r-project.org/web/packages/ROCR/ROCR.pdf
CRAN: pROC https://cran.r-project.org/web/packages/pROC/pROC.pdf

ECII

Área bajo la Curva Operativa Característica AUC

Alberto Torrejón Valenzuela
Doble Grado Matemáticas y Estadística

Enunciado

Datos: Moneyball

Modelos de clasificación

Partición entrenamiento-test

Curvas ROC

AUC

library(ROCR)

library(pROC)

Curva ROC script

Función AUC

library(ROCR)

library(pROC)

Fórmulas de integración numérica

Fórmula del trapecio

Fórmula del rectángulo inferior

Fórmula del rectángulo superior

Estadístico U Mann-Whitney

Usando `wilcox.test()`

Calculando el estadístico

Visión probabilística

Comparación de rangos

Versión compacta

Comparación de los distintos modelos

Bibliografía

ECII

Área bajo la Curva Operativa Característica AUC

Alberto Torrejón ValenzuelaDoble Grado Matemáticas y Estadística

Enunciado

Datos: Moneyball

Modelos de clasificación

Partición entrenamiento-test

Curvas ROC

AUC

library(ROCR)

library(pROC)

Curva ROC script

Función AUC

library(ROCR)

library(pROC)

Fórmulas de integración numérica

Fórmula del trapecio

Fórmula del rectángulo inferior

Fórmula del rectángulo superior

Estadístico U Mann-Whitney

Usando wilcox.test()

Calculando el estadístico

Visión probabilística

Comparación de rangos

Versión compacta

Comparación de los distintos modelos

Bibliografía

Alberto Torrejón Valenzuela
Doble Grado Matemáticas y Estadística

Usando `wilcox.test()`