setwd("/media/DATA/Dropbox/MyR/EpiTools/APS")

# Cargar los datos. Los dias deben estar ordenados en forma ascendente
data <- read.table(file="abcpe.txt", sep="",header=T, dec=".")

# Crear el gr?fico de la severidad de la enfermedad a trav?s del tiempo
plot(sev~dia,data=data,
   xlab="Time",
   ylim=c(0,max(data$sev)+1),
   xlim=c(0,max(data$dia)+0.5),
   ylab="Disease Severity (%)",
   type="o",
   pch=19,
   col="mediumblue")

# Agregar un t?tulo y subt?tulo a nuestro gr?fico 
title(main="Illustration of AUDPC Calculation")

# Ilustrar el AUDPC con resct?ngulos 
sev.promedio <- ((data$sev[1:(nrow(data)-1)])+(data$sev[2:nrow(data)]))/2
rect(data$dia[1:(nrow(data)-1)],
     0,
     data$dia[2:nrow(data)],
     sev.promedio,border="orange")

# Sumatoria de todos los trapezios y escritura en grafico
ABCPE <- sum(diff(data$dia)*sev.promedio)
ABCPEst<-ABCPE/max(data$dia)-min(data$dia)
text(data$dia[2],sev.promedio[length(sev.promedio)], paste("ABCPE =",as.character(round(ABCPE,2)),sep=" "))
text(data$dia[2],(0.85*sev.promedio[length(sev.promedio)]), paste("ABCPEst =",as.character(round(ABCPEst,2)),sep=" "))

require(pracma)
## Loading required package: pracma
(AUC = trapz(data$dia,data$sev))
## [1] 206
text(data$dia[6],sev.promedio[length(sev.promedio)], paste("AUC(){pracma}=",as.character(round(AUC,2)),sep=" "))