#http://ares.inf.um.es/00Rteam/pub/clas/II_ROC.html
#https://rpubs.com/JairoAyala/592802

#Making-of
# Crear datos
altE <- rnorm( 100, 173, 2 )
altS <- rnorm( 100, 167, 2 )
dE <- density( altE, from = min( altS ), to = max( altE ), n = 1000 )
dS <- density( altS, from = min( altS ), to = max( altE ), n = 1000 )

#Para los enfermos
muestra <- which( dE$y > 0.005 )
muestra <- sample( muestra, 50, replace = FALSE )
y <- vector()

for ( i in muestra ){
  y <- append( y, runif( 1, min( dE$y ), dE$y[i] )  ) 
}

#Para los sanos
muestra2 <- which( dS$y > 0.005 )
muestra2 <- sample( muestra2, 50, replace = FALSE )
ym <- vector()

for ( i in muestra2 ){
  ym <- append(ym, runif(1, min( dS$y ), dS$y[i] ) )
}

incr <- 0
intersX <- which( diff( dE$y < dS$y ) != 0 )

curvas <- function() {
  plot( dE$x, dE$y, type = "l", xlim = c( min( min( dE$x ), min( dS$x ) ),
                                          max( max( dE$x ), max( dS$x ) ) ), ylim=c( 0, 0.25 ) ) 
  lines( dS$x, dS$y )
  points( dE$x[muestra], y, pch = "+", col = 4 )
  text( 175, 0.25, labels = "Enfermos" )
  points( dS$x[muestra2], ym, pch = "-", col = 3 )
  text( 165, 0.25, labels = "Sanos" )
  text( c( 167, 167, 172, 171 ), c( 0, 0.10, 0.10, 0 ), labels = c("FN", "TN" , "TP", "FP" ) ) }

curvas()
abline( v = dE$x[intersX], col = 2 ) 

# número de FN
FN <- length( which( dE$x[muestra] < dE$x[intersX] ) )
# Número de TP
TP <- length( which( dE$x[muestra] >= dE$x[intersX] ) )
# Número de FP
FP <- length( which( dS$x[muestra2] >= dS$x[intersX] ) )
# Número de TN
TN <- length (which( dS$x[muestra2] < dS$x[intersX] ) )

( sensibilidad <- TP / ( FN + TP ) )
( especificidad <- TN / ( FP + TN ) )

curvas()
incr <- c( -2 )
intersX <- which( diff( dE$y < dS$y ) != 0 )
abline( v = dE$x[intersX] + incr, col = 3 )


incr <- c( -2 )
FN <- length( which( dE$x[muestra] < dE$x[intersX]+incr ) )
TP <- length( which( dE$x[muestra] >= dE$x[intersX]+incr ) )
FP <- length( which( dS$x[muestra2] >= dS$x[intersX]+incr ) )
TN <- length( which( dS$x[muestra2] < dS$x[intersX]+incr ) )

( sensibilidad2 <- TP / ( FN + TP ) )

( especificidad2 <- TN / ( FP + TN ) )

# Aumentando la sensibilidad
plot( seq( 0, 1, 0.1 ), seq( 0, 1, 0.1 ), type = "l", xlab = "1-especificidad", ylab = "Sensibilidad" ) #out
sn <- c( sensibilidad, sensibilidad2 )
es <- c( especificidad, especificidad2 )
points( 1-es, sn, pch = c( 1, 2 ), col = c(2,3) )
lines( 1-es, sn )

plot( dE$x, dE$y, type = "l", xlim = c( min( min( dE$x ), min( dS$x ) ),
                                        max( max( dE$x ), max( dS$x ) ) ),
      ylim=c( 0, 0.25 ) )

points( dE$x[muestra], y, pch = "+", col = 4 )
text( 175, 0.25, labels = "Enfermos" )
lines( dS$x, dS$y )
points( dS$x[muestra2], ym, pch = "-", col = 3 )
text(165, 0.25, labels = "Sanos" )

for (i in seq(-5,5,0.3)){
  abline( v = 170+i, col = 1 )
}

x <- vector()
y <- vector()

for ( i in seq( -10, 10, 0.1 ) ){
  incr <- i
  FN <- length( which( dE$x[muestra] < dE$x[intersX] + incr ) )
  TP <- length( which( dE$x[muestra] >= dE$x[intersX] + incr ) )
  FP <- length( which( dS$x[muestra2] >= dS$x[intersX] + incr ) )
  TN <- length( which( dS$x[muestra2] < dS$x[intersX] + incr ) )
  sensi <- TP / ( FN + TP )
  espec <- TN /( FP + TN )
  x <- append( x, espec ) 
  y <- append( y, sensi )
  
}

plot( seq(0, 1, 0.1 ) , seq(0, 1, 0.1 ), type = "n", xlab = "1-especificidad", ylab = "Sensibilidad" ) #out
abline( a = 0, b = 1 )
lines( 1-x, y) 
points(1-x, y, pch = "*" )


#con paquete pROC
datos <- data.frame( c( altS,altE ), as.factor( rep( c("S", "E" ), times = 1, each = 100 ) ) )
colnames( datos ) <- c( "val", "attr" )

rocobj <- roc(datos$attr, datos$val, auc = TRUE, ci = TRUE  )
print(rocobj)
plot.roc(rocobj,print.auc=T,print.thres = "best",
         col="blue",xlab="1-ESpecificidad",ylab="Sensibilidad",legacy.axes=T)

SensEspec <-  rocobj$sensitivities + rocobj$specificities 
maximo <- max( SensEspec )
numOrdenCutoff <- which( SensEspec == maximo )
rocobj$thresholds[numOrdenCutoff]

#con paquete ROCR

library('ROCR')
datos <- data.frame( c( altS,altE ), as.factor( rep( c("E", "S" ), times = 1, each = 100 ) ) )
colnames( datos ) <- c( "val", "attr" )
##ACa se necesita cambiar las etiquetas

pred <- prediction(predictions=datos$val,labels=datos$attr)
perf <- performance(pred, measure = "tpr", x.measure = "fpr" )

plot(perf, colorize = TRUE, type = "l" ) #out
abline(a = 0, b = 1 )

# Calcular el área bajo la curva
AUC <- performance(pred, measure = "auc")
AUCaltuca <- AUC@y.values

# Calcular el punto de corte óptimo
cost.perf <- performance( pred, measure = "cost" )
opt.cut <- pred@cutoffs[[1]][which.min(cost.perf@y.values[[1]])]
#coordenadas del punto de corte óptimo
x<-perf@x.values[[1]][which.min( cost.perf@y.values[[1]] ) ]
y<-perf@y.values[[1]][which.min( cost.perf@y.values[[1]] ) ]
points(x,y, pch=20, col="red")



#Ejercicio

imc <- read.table( "http://gauss.inf.um.es/datos/IMC.dat", sep = " ", header = TRUE )

rocimc <- roc( imc$attr, imc$val, auc = TRUE, ci = TRUE  )
print( rocimc )

plot.roc( rocimc, legacy.axes = TRUE, print.thres = "best", print.auc = TRUE,
          auc.polygon = FALSE, max.auc.polygon = FALSE, auc.polygon.col = "gainsboro", col = 2, grid = TRUE )


peso <- read.table( "http://gauss.inf.um.es/datos/peso.dat", sep = " ", header = TRUE )
rocpeso <- roc(peso$attr, peso$vals, auc = TRUE )
print(rocpeso)

plot.roc( rocpeso, add = TRUE, print.thres = "best", print.auc = FALSE, col = 3 )