rocr.r

gistfile1.txt

args <- commandArgs(TRUE)
dat =  read.table(file=args[1])

library(ROCR)

predictions = dat$V2
labels = dat$V3

# perform prediction
pred <- prediction(predictions, labels)

# performance
roc.perf <- performance(pred, measure = "tpr", x.measure = "fpr")


# draw roc curve
plot(roc.perf, col=rainbow(10))
abline(a=0, b= 1)

# precise and recall
rec.perf <- performance(pred, "prec", "rec")
plot(rec.perf, col=rainbow(10))

# accuracy
acc.perf = performance(pred, "acc")
plot(acc.perf, avg= "vertical", spread.estimate="boxplot", show.spread.at= seq(0.1, 0.9, by=0.1))

ind = which.max( slot(acc.perf, "y.values")[[1]] )
acc = slot(acc.perf, "y.values")[[1]][ind]
cutoff = slot(acc.perf, "x.values")[[1]][ind]
print(c(accuracy= acc, cutoff = cutoff))

# sensitivity and specificity
sens.perf<- performance(pred, "sens", "spec")
plot(sens.perf)

# check calibration
cal.perf <- performance(pred, "cal", window.size=50)
plot(cal.perf)

# cutoff density plot
plot(0,0,type="n", xlim= c(0,1), ylim=c(0,7), xlab="Cutoff", ylab="Density",
main="How well do the predictions separate the classes?")

for (runi in 1:length(pred@predictions)) {
  lines(density(pred@predictions[[runi]][pred@labels[[runi]]=="1"]), col= "red")
  lines(density(pred@predictions[[runi]][pred@labels[[runi]]=="0"]), col="green")
}

# find optimal cutoff
opt.cut = function(perf, pred){
    cut.ind = mapply(FUN=function(x, y, p){
        d = (x - 0)^2 + (y-1)^2
        ind = which(d == min(d))
        c(sensitivity = y[[ind]], specificity = 1-x[[ind]], cutoff = p[[ind]])
    }, [email protected], [email protected], pred@cutoffs)
}

print(opt.cut(roc.perf, pred))

# calculate AUC value
auc.perf = performance(pred, measure = "auc")
[email protected]