#Load required libraries
library(ROCR)
library(PresenceAbsence)

#Set working directory
setwd("D:/Area2/MGET/V2")

#Read table (ID, observed class, predicted probability)
table <- read.table("rf2000_2_export_for_PresenceAbsence.txt", header = TRUE, sep = ",")
pred <- prediction( table$Prob, table$Obs)
ACC <- performance(pred, "acc","cutoff")

#Define required sensitivity and specificity
ReqSens <- 0.95
ReqSpec <- 0.95

#Calculate optimal thresholds
optimal.thresholds(table, opt.methods = 10:11, req.sens = ReqSens, req.spec = ReqSpec)

#Define thresholds
PA <- 0.619 #Maximizing the Youden index (see rf2000_2_model_stats.txt)
SensTh <- 0.20 #Sensitivity threshold
SpecTh <- 0.80 #Specificity threshold

#Presence-Absence summary plot
jpeg("Summary plot.jpeg", width = 16, height = 10, res = 500, units = "cm", pointsize = 8, quality = 100)
par(mar=c(4,4,1,1))
par(mfrow=c(2,2))

auc.roc.plot(table, xlab = "1 - Specificity", ylab = "Sensitivity", main = "", line.type = FALSE, lwd=2, legend.text = "")

plot(ACC, lwd=2, xlab = "Threshold", ylab = "PCC")
abline(v=PA)

error.threshold.plot(table, lwd=2, main = "", opt.thresholds = TRUE, opt.methods = 10:11, req.sens = ReqSens, req.spec = ReqSpec, add.legend = FALSE, add.opt.legend = FALSE)
abline(v=SensTh)
abline(v=SpecTh)
abline(v=PA)

presence.absence.hist(table, main = "", xlab = "Predicted probability", ylab = "Frequency", N.bars = 10, opt.thresholds = TRUE, opt.methods = 10:11, req.sens = ReqSens, req.spec = ReqSpec, add.legend = FALSE, add.opt.legend = FALSE)
abline(v=PA)
dev.off()
