Script for ROC
--------------

> Xroc = rnorm(300)
> Yroc = rbinom(300,1,plogis(-1+0.4*Xroc))
> fitLg = glm(Yroc ~ Xroc, family=binomial)
> SensFP = function(pi0, Y, fit) c(mean(fit[Y==0]>pi0),
     mean(fit[Y==1]>pi0))
plot(1,1, xlab="1-Spec", ylab="Sens", type="n", xlim=c(0,1),
       ylim=c(0,1), main="ROC Plot")
 abline(0,1, col="red")

> Yroc2 = rbinom(300,1, plogis(-2+0.6*Xroc))
  fitLg2 = glm(Yroc2 ~ Xroc, family=binomial)

> ROCplot = function(Y,fit,pch) {
     for(i in 1:99) {
       aux = SensFP(i/100, Y, fit)
       points(aux[1],aux[2], pch=pch) } }
 
> ROCplot(Yroc,fitLg$fit,1)
  ROCplot(Yroc2,fitLg2$fit,5)

> Yroc3 = rbinom(300,1, plogis(-2+2*Xroc))
  fitLg3 = glm(Yroc3 ~ Xroc, family=binomial)
  ROCplot(Yroc3,fitLg3$fit,18)

CALCULATE AUC USING TRAPEZOID RULE:

SenFParr = array(1, c(101,2,3), dimnames=list(NULL, 
     c("FP","Sens"), paste0("fit",1:3)))
SenFParr[101,,]=0
for(i in 1:99) {
    SenFParr[i+1,,1] = SensFP(i/100,Yroc,fitLg$fit)
    SenFParr[i+1,,2] = SensFP(i/100,Yroc2,fitLg2$fit)
    SenFParr[i+1,,3] = SensFP(i/100,Yroc3,fitLg3$fit) }


> 0.5*sum((SenFParr[-101,1,1]-SenFParr[-1,1,1])*(SenFParr[-101,2,1]+SenFParr[-1,2,1]))
[1] 0.6684605                #### COMPARE with 0.5

> 0.5*sum((SenFParr[-101,1,3]-SenFParr[-1,1,3])*(SenFParr[-101,2,3]+SenFParr[-1,2,3]))
[1] 0.8608454