#=======================================================
# Illustrate use of ROC curves with ungrouped beetle data
#===============================================================

# Read ungrouped beetle data from the web:
data="http://www.stat.ufl.edu/~aa/glm/data/Beetles.dat"
beetle.u = read.table(data, header = T)
head(beetle.u)

# Will use the ROCR library
library(ROCR)

#============================================
# ROC curve from a logistic regression model
#============================================ 
fit.logistic=glm(y~x,family=binomial,data=beetle.u)
pred.logistic=prediction(fitted(fit.logistic),beetle.u$y)
perf.logistic=performance(pred.logistic,"tpr","fpr")
plot(perf.logistic,main="Logistic",col="blue")
#Area under curve (AUC):
auc.logistic <- performance(pred.logistic,"auc")
auc.logistic@y.values[[1]]

#==================================
# ROC curve from using cloglog-link
#==================================
fit.cloglog=glm(y~x,family=binomial(link=cloglog),data=beetle.u)
pred.cloglog=prediction(fitted(fit.cloglog),beetle.u$y)
perf.cloglog=performance(pred.cloglog,"tpr","fpr")
plot(perf.cloglog,main="cloglog",add=T,col="red")
#Area under curve (AUC):
auc.cloglog <- performance(pred.cloglog,"auc")
auc.cloglog@y.values[[1]]
# For these two models the ROC curves are the same

#==============================================
# Correlation between observed and fitted model
# see section 5.2.5 in the text book
#==============================================

cor(fitted(fit.logistic),beetle.u$y)
cor(fitted(fit.cloglog),beetle.u$y)

# Here the cloglog-link gives a slightly better result