###################################################
### chunk number 1: load data
###################################################
library("PresenceAbsence")
data("SPDATA")


###################################################
### chunk number 2: examine data 1
###################################################
head(SPDATA)


###################################################
### chunk number 3: define variables
###################################################
species<-as.character(unique(SPDATA$SPECIES))
model.names<-as.character(names(SPDATA)[-c(1,2)])
N.models<-ncol(SPDATA)-2
N.sp<-length(species)
N.obs<-length(SPDATA$SPECIES[SPDATA$SPECIES==species[1]])
Obs.prev<-table(SPDATA$SPECIES,SPDATA$OBSERVED)[,2]/ N.obs
Obs.prev<-round(Obs.prev,digits=2)


###################################################
### chunk number 4: fig1plot
###################################################
par(oma=c(0,0,3,0),mar=c(4,4,4,1),mfrow=c(1,3),cex=1)
sp<-1
DATA<-SPDATA[SPDATA$SPECIES==species[sp],]
for(mod in 1:N.models){
	presence.absence.hist(	DATA,
	which.model=mod,
	legend.cex=1,
	N.bars=20)}
mtext(species[sp],side=3,line=.5,cex=1.5,outer=TRUE)


###################################################
### chunk number 5: fig1
###################################################
par(oma=c(0,0,3,0),mar=c(4,4,4,1),mfrow=c(1,3),cex=1)
sp<-1
DATA<-SPDATA[SPDATA$SPECIES==species[sp],]
for(mod in 1:N.models){
	presence.absence.hist(	DATA,
	which.model=mod,
	legend.cex=1,
	N.bars=20)}
mtext(species[sp],side=3,line=.5,cex=1.5,outer=TRUE)


###################################################
### chunk number 6: fig2plot
###################################################
par(oma=c(0,0,3,0),mar=c(4,4,4,1),mfrow=c(1,3),cex=1)
sp<-1
DATA<-SPDATA[SPDATA$SPECIES==species[sp],]
for(mod in 1:N.models){
	presence.absence.hist(	DATA,
	which.model=mod,
	legend.cex=1,
	N.bars=15,
	truncate.tallest=TRUE)}
mtext(species[sp],side=3,line=.5,cex=1.5,outer=TRUE)


###################################################
### chunk number 7: fig2
###################################################
par(oma=c(0,0,3,0),mar=c(4,4,4,1),mfrow=c(1,3),cex=1)
sp<-1
DATA<-SPDATA[SPDATA$SPECIES==species[sp],]
for(mod in 1:N.models){
	presence.absence.hist(	DATA,
	which.model=mod,
	legend.cex=1,
	N.bars=15,
	truncate.tallest=TRUE)}
mtext(species[sp],side=3,line=.5,cex=1.5,outer=TRUE)


###################################################
### chunk number 8: fig3plot
###################################################
par(oma=c(0,0,3,0),mar=c(4,4,4,1),mfrow=c(1,3),cex=1)
mod<-2
for(sp in c("ACGR3","PIEN","POTR5")){
	presence.absence.hist(	DATA=SPDATA[SPDATA$SPECIES==sp,],
	which.model=mod,
	legend.cex=1,
	N.bars=15,
	main=sp)}
mtext(model.names[mod],side=3,line=0,cex=1.5,outer=TRUE)


###################################################
### chunk number 9: fig3
###################################################
par(oma=c(0,0,3,0),mar=c(4,4,4,1),mfrow=c(1,3),cex=1)
mod<-2
for(sp in c("ACGR3","PIEN","POTR5")){
	presence.absence.hist(	DATA=SPDATA[SPDATA$SPECIES==sp,],
	which.model=mod,
	legend.cex=1,
	N.bars=15,
	main=sp)}
mtext(model.names[mod],side=3,line=0,cex=1.5,outer=TRUE)


###################################################
### chunk number 10: fig4plot
###################################################
par(oma=c(0,0,3,0),mar=c(4,4,4,1),mfrow=c(1,3),cex=1)
mod<-2
for(sp in c("JUSC2","ABCO","PICO")){
	presence.absence.hist(	DATA=SPDATA[SPDATA$SPECIES==sp,],
	which.model=mod,
	legend.cex=1,
	N.bars=15,
	truncate.tallest=TRUE,
	main=sp)}
mtext(model.names[mod],side=3,line=0,cex=1.5,outer=T)


###################################################
### chunk number 11: fig4
###################################################
par(oma=c(0,0,3,0),mar=c(4,4,4,1),mfrow=c(1,3),cex=1)
mod<-2
for(sp in c("JUSC2","ABCO","PICO")){
	presence.absence.hist(	DATA=SPDATA[SPDATA$SPECIES==sp,],
	which.model=mod,
	legend.cex=1,
	N.bars=15,
	truncate.tallest=TRUE,
	main=sp)}
mtext(model.names[mod],side=3,line=0,cex=1.5,outer=T)


###################################################
### chunk number 12: predprev 1
###################################################
DATA<-SPDATA[SPDATA$SPECIES==species[1],]
pred.prev<-predicted.prevalence(DATA=DATA)
for(sp in 2:N.sp){
	DATA<-SPDATA[SPDATA$SPECIES==species[sp],]
	pred.prev<-rbind(pred.prev,predicted.prevalence(DATA=DATA))}
pred.prev<-cbind(species=species,pred.prev)
pred.prev[,3:6]<-round( pred.prev[,3:6], digits = 2)
pred.prev


###################################################
### chunk number 13: predprev 2
###################################################
sp<-1
DATA<-SPDATA[SPDATA$SPECIES==species[1],]
pred.prev<-predicted.prevalence(DATA=DATA,threshold=11)
pred.prev[,2:5]<-round( pred.prev[,2:5], digits = 2)
print(paste("Species:",species[sp]))
pred.prev


###################################################
### chunk number 14: accu 1
###################################################
sp<-1
DATA<-SPDATA[SPDATA$SPECIES==species[sp],]
accu<-presence.absence.accuracy(DATA)
accu[,-c(1,2)]<-signif(accu[,-c(1,2)], digits =2)
print(paste("Species:",species[sp]))
accu


###################################################
### chunk number 15: accu 2
###################################################
sp<-1
DATA<-SPDATA[SPDATA$SPECIES==species[sp],]
accu<-presence.absence.accuracy(	DATA, 
  which.model=3,
  threshold=11,
  st.dev=FALSE)
accu[,-c(1,2)]<-signif(accu[,-c(1,2)], digits =2)
print(paste("Species:",species[sp]))
accu


###################################################
### chunk number 16: accu 3
###################################################
accu<-presence.absence.accuracy(	DATA,
  which.model=3,
  threshold=c(.2,.4,.5,.89),
  st.dev=FALSE)
accu[,-c(1,2)]<-signif(accu[,-c(1,2)], digits =2)
print(paste("Species:",species[sp]))
accu


###################################################
### chunk number 17: fig5plot
###################################################
par(oma=c(0,5,0,0),mar=c(3,3,3,1),mfrow=c(2,3),cex=.7,cex.lab=1.4,mgp=c(2,.5,0))
sp<-c("JUSC2","PICO")
for(i in 1:2){
	DATA=SPDATA[SPDATA$SPECIES==sp[i],]
	for(mod in 1:3){
		if(i==2 & mod==3){legend.flag<-TRUE}else{legend.flag<-FALSE}
		error.threshold.plot(	DATA, 
		which.model=mod,
		color=TRUE,
		add.legend=legend.flag,
		legend.cex=1.3)
		if(mod==1){
			mtext(sp[i],side=2,line=6.3,cex=1.6)
			mtext(paste(	"(",Obs.prev[names(Obs.prev)==sp[i]],"Prevalence )"),
				side=2,line=4.3,cex=1.3)}}}


###################################################
### chunk number 18: fig5
###################################################
par(oma=c(0,5,0,0),mar=c(3,3,3,1),mfrow=c(2,3),cex=.7,cex.lab=1.4,mgp=c(2,.5,0))
sp<-c("JUSC2","PICO")
for(i in 1:2){
	DATA=SPDATA[SPDATA$SPECIES==sp[i],]
	for(mod in 1:3){
		if(i==2 & mod==3){legend.flag<-TRUE}else{legend.flag<-FALSE}
		error.threshold.plot(	DATA, 
		which.model=mod,
		color=TRUE,
		add.legend=legend.flag,
		legend.cex=1.3)
		if(mod==1){
			mtext(sp[i],side=2,line=6.3,cex=1.6)
			mtext(paste(	"(",Obs.prev[names(Obs.prev)==sp[i]],"Prevalence )"),
				side=2,line=4.3,cex=1.3)}}}


###################################################
### chunk number 19: fig6plot
###################################################
par(oma=c(0,0,0,0),mfrow=c(1,2),cex=.7,cex.lab=1.5)
sp<-c("ACGR3","POTR5")
for(i in 1:2){
	DATA<-SPDATA[SPDATA$SPECIES==sp[i],]
	auc.roc.plot(	DATA, 
				color=TRUE,
				legend.cex=1.4,
				main="")
	mtext(sp[i],side=3,line=2.5,cex=1.6)
	mtext(paste(	Obs.prev[names(Obs.prev)==sp[i]],"Prevalence"),side=3,line=.5,cex=1.3)}


###################################################
### chunk number 20: fig6
###################################################
par(oma=c(0,0,0,0),mfrow=c(1,2),cex=.7,cex.lab=1.5)
sp<-c("ACGR3","POTR5")
for(i in 1:2){
	DATA<-SPDATA[SPDATA$SPECIES==sp[i],]
	auc.roc.plot(	DATA, 
				color=TRUE,
				legend.cex=1.4,
				main="")
	mtext(sp[i],side=3,line=2.5,cex=1.6)
	mtext(paste(	Obs.prev[names(Obs.prev)==sp[i]],"Prevalence"),side=3,line=.5,cex=1.3)}


###################################################
### chunk number 21: fig7plot
###################################################
par(oma=c(0,0,0,0),mfrow=c(1,2),cex=.7,cex.lab=1.5)
sp<-c("ACGR3","POTR5")
for(i in 1:2){
	DATA<-SPDATA[SPDATA$SPECIES==sp[i],]
	auc.roc.plot(	DATA, 
				color=TRUE,
				legend.cex=1.4,
				mark=5,
				which.model=3,
				main="")
	mtext(sp[i],side=3,line=2.5,cex=1.6)
	mtext(paste(	Obs.prev[names(Obs.prev)==sp[i]],"Prevalence"),side=3,line=.5,cex=1.3)}


###################################################
### chunk number 22: fig7
###################################################
par(oma=c(0,0,0,0),mfrow=c(1,2),cex=.7,cex.lab=1.5)
sp<-c("ACGR3","POTR5")
for(i in 1:2){
	DATA<-SPDATA[SPDATA$SPECIES==sp[i],]
	auc.roc.plot(	DATA, 
				color=TRUE,
				legend.cex=1.4,
				mark=5,
				which.model=3,
				main="")
	mtext(sp[i],side=3,line=2.5,cex=1.6)
	mtext(paste(	Obs.prev[names(Obs.prev)==sp[i]],"Prevalence"),side=3,line=.5,cex=1.3)}


###################################################
### chunk number 23: fig8plot
###################################################
par(oma=c(0,0,2,0),mar=c(4,4,4,1),mfrow=c(1,3),cex=.7,cex.lab=1.4,mgp=c(2,.5,0))
sp<- c("JUSC2","ABCO","PICO") 
for(i in 1:3){
	DATA<-SPDATA[SPDATA$SPECIES==sp[i],]
	auc.roc.plot(	DATA, 
				color=TRUE,
				legend.cex=1.2,
				main="")
	mtext(sp[i],side=3,line=2,cex=1.6)
	mtext(paste(	Obs.prev[names(Obs.prev)==sp[i]],"Prevalence"),side=3,line=.5,cex=1.3)}


###################################################
### chunk number 24: fig8
###################################################
par(oma=c(0,0,2,0),mar=c(4,4,4,1),mfrow=c(1,3),cex=.7,cex.lab=1.4,mgp=c(2,.5,0))
sp<- c("JUSC2","ABCO","PICO") 
for(i in 1:3){
	DATA<-SPDATA[SPDATA$SPECIES==sp[i],]
	auc.roc.plot(	DATA, 
				color=TRUE,
				legend.cex=1.2,
				main="")
	mtext(sp[i],side=3,line=2,cex=1.6)
	mtext(paste(	Obs.prev[names(Obs.prev)==sp[i]],"Prevalence"),side=3,line=.5,cex=1.3)}


###################################################
### chunk number 25: opt 1
###################################################
sp<-"QUGA"
DATA<-SPDATA[SPDATA$SPECIES==sp,]
print(paste("Species:",sp))
optimal.thresholds(DATA,opt.methods=1:12)


###################################################
### chunk number 26: opt 2
###################################################
sp<-"QUGA"
DATA<-SPDATA[SPDATA$SPECIES==sp,]
print(paste("Species:",sp))
optimal.thresholds(DATA,opt.methods=1:12, req.sens=0.90, req.spec=0.90, FPC=2, FNC=1)


###################################################
### chunk number 27: fig9plot
###################################################
par(oma=c(0,0,3,0),mar=c(4,4,4,1),mfrow=c(1,3),cex=1)
sp<-c("JUSC2","ABCO","PICO")
mod<-2
for(i in 1:3){
	presence.absence.hist(	DATA=SPDATA[SPDATA$SPECIES==sp[i],],
					which.model=mod,
					N.bars=15,
					truncate.tallest=TRUE,
					legend.cex=.60,
					opt.legend.cex=.60,
					opt.thresholds=TRUE,
					opt.methods=c(1,2,10,4),
					req.sens=0.85,
					main=sp[i])
	mtext(paste(Obs.prev[names(Obs.prev)==sp[i]],"Prevalence"),side=3,line=.5,cex=1.1)}
mtext(model.names[mod],side=3,line=0,cex=1.5,outer=T)


###################################################
### chunk number 28: fig9
###################################################
par(oma=c(0,0,3,0),mar=c(4,4,4,1),mfrow=c(1,3),cex=1)
sp<-c("JUSC2","ABCO","PICO")
mod<-2
for(i in 1:3){
	presence.absence.hist(	DATA=SPDATA[SPDATA$SPECIES==sp[i],],
					which.model=mod,
					N.bars=15,
					truncate.tallest=TRUE,
					legend.cex=.60,
					opt.legend.cex=.60,
					opt.thresholds=TRUE,
					opt.methods=c(1,2,10,4),
					req.sens=0.85,
					main=sp[i])
	mtext(paste(Obs.prev[names(Obs.prev)==sp[i]],"Prevalence"),side=3,line=.5,cex=1.1)}
mtext(model.names[mod],side=3,line=0,cex=1.5,outer=T)


###################################################
### chunk number 29: fig10plot
###################################################
par(oma=c(0,5,0,0),mar=c(3,3,3,1),mfrow=c(2,3),cex=.7,cex.lab=1.4,mgp=c(2,.5,0))
sp<-c("JUSC2","PICO")
for(i in 1:2){
	DATA=SPDATA[SPDATA$SPECIES==sp[i],]
	for(mod in 1:3){
		if(i==2 & mod==2){
			legend.flag<-TRUE}else{legend.flag<-FALSE}
		if(i==2 & mod==3){
			opt.legend.flag<-TRUE}else{opt.legend.flag<-FALSE}
		error.threshold.plot(	DATA, 
					which.model=mod, 
					color=TRUE,
					opt.thresholds=TRUE,
					opt.methods=c(1,2,10,4),
					req.sens=0.9,
					add.legend=legend.flag,
					add.opt.legend= opt.legend.flag,
					legend.cex=1.1,
					opt.legend.cex=1.1,
					vert.lines=FALSE)
		if(mod==1){
			mtext(sp[i],side=2,line=6.3,cex=1.6)
			mtext(paste("(",Obs.prev[names(Obs.prev)==sp[i]],"Prevalence )"),
				side=2,line=4.3,cex=1.3)}}}


###################################################
### chunk number 30: fig10
###################################################
par(oma=c(0,5,0,0),mar=c(3,3,3,1),mfrow=c(2,3),cex=.7,cex.lab=1.4,mgp=c(2,.5,0))
sp<-c("JUSC2","PICO")
for(i in 1:2){
	DATA=SPDATA[SPDATA$SPECIES==sp[i],]
	for(mod in 1:3){
		if(i==2 & mod==2){
			legend.flag<-TRUE}else{legend.flag<-FALSE}
		if(i==2 & mod==3){
			opt.legend.flag<-TRUE}else{opt.legend.flag<-FALSE}
		error.threshold.plot(	DATA, 
					which.model=mod, 
					color=TRUE,
					opt.thresholds=TRUE,
					opt.methods=c(1,2,10,4),
					req.sens=0.9,
					add.legend=legend.flag,
					add.opt.legend= opt.legend.flag,
					legend.cex=1.1,
					opt.legend.cex=1.1,
					vert.lines=FALSE)
		if(mod==1){
			mtext(sp[i],side=2,line=6.3,cex=1.6)
			mtext(paste("(",Obs.prev[names(Obs.prev)==sp[i]],"Prevalence )"),
				side=2,line=4.3,cex=1.3)}}}


###################################################
### chunk number 31: fig11plot
###################################################
par(oma=c(0,0,2,0),mar=c(4,4,4,1),mfrow=c(1,3),cex=.7,cex.lab=1.4,mgp=c(2,.5,0))
sp<-c("ACGR3","QUGA","JUOS")
mod=1
for(i in 1:3){
	DATA=SPDATA[SPDATA$SPECIES==sp[i],]
	if(i==2){legend.flag<-TRUE}else{legend.flag<-FALSE}
	if(i==3){opt.legend.flag<-TRUE}else{opt.legend.flag<-FALSE}
	error.threshold.plot(	DATA, 
	which.model=mod, 
	color=TRUE,
	opt.thresholds=TRUE,
	opt.methods=c(	"Default", "MaxPCC", "Sens=Spec"),
	add.legend=legend.flag,
	add.opt.legend=opt.legend.flag,
	legend.cex=1.1,
	opt.legend.cex=1.1,
	vert.lines=TRUE,
	main="")
	mtext(sp[i],side=3,line=2,cex=1.6)
	mtext(paste(	Obs.prev[names(Obs.prev)==sp[i]],"Prevalence"),side=3,line=.5,cex=1.3)}
#mtext(model.names[mod],side=3,line=-1,cex=1.6,outer=TRUE)


###################################################
### chunk number 32: fig11
###################################################
par(oma=c(0,0,2,0),mar=c(4,4,4,1),mfrow=c(1,3),cex=.7,cex.lab=1.4,mgp=c(2,.5,0))
sp<-c("ACGR3","QUGA","JUOS")
mod=1
for(i in 1:3){
	DATA=SPDATA[SPDATA$SPECIES==sp[i],]
	if(i==2){legend.flag<-TRUE}else{legend.flag<-FALSE}
	if(i==3){opt.legend.flag<-TRUE}else{opt.legend.flag<-FALSE}
	error.threshold.plot(	DATA, 
	which.model=mod, 
	color=TRUE,
	opt.thresholds=TRUE,
	opt.methods=c(	"Default", "MaxPCC", "Sens=Spec"),
	add.legend=legend.flag,
	add.opt.legend=opt.legend.flag,
	legend.cex=1.1,
	opt.legend.cex=1.1,
	vert.lines=TRUE,
	main="")
	mtext(sp[i],side=3,line=2,cex=1.6)
	mtext(paste(	Obs.prev[names(Obs.prev)==sp[i]],"Prevalence"),side=3,line=.5,cex=1.3)}
#mtext(model.names[mod],side=3,line=-1,cex=1.6,outer=TRUE)


###################################################
### chunk number 33: fig12plot
###################################################
par(oma=c(0,5,0,0),mar=c(3,3,3,1),mfrow=c(2,3),cex=.7,cex.lab=1.4,mgp=c(2,.5,0))
sp<-c("JUSC2","PICO")
for(i in 1:2){
	DATA<-SPDATA[SPDATA$SPECIES==sp[i],]
	for(mod in 1:3){
		auc.roc.plot(	DATA, 
		color=mod+1,
		which.model=mod,
		main=model.names[mod],
		opt.thresholds=TRUE,
		opt.methods=c(1,4,6),
		mark.numbers=TRUE,
		mark.color=FALSE,
		legend.cex=1.1,
		opt.legend.cex=1.1)
	if(mod==1){
			mtext(sp[i],side=2,line=6.3,cex=1.6)
			mtext(paste("(",Obs.prev[names(Obs.prev)==sp[i]],"Prevalence )"),
				side=2,line=4.3,cex=1.3)}}}


###################################################
### chunk number 34: fig12
###################################################
par(oma=c(0,5,0,0),mar=c(3,3,3,1),mfrow=c(2,3),cex=.7,cex.lab=1.4,mgp=c(2,.5,0))
sp<-c("JUSC2","PICO")
for(i in 1:2){
	DATA<-SPDATA[SPDATA$SPECIES==sp[i],]
	for(mod in 1:3){
		auc.roc.plot(	DATA, 
		color=mod+1,
		which.model=mod,
		main=model.names[mod],
		opt.thresholds=TRUE,
		opt.methods=c(1,4,6),
		mark.numbers=TRUE,
		mark.color=FALSE,
		legend.cex=1.1,
		opt.legend.cex=1.1)
	if(mod==1){
			mtext(sp[i],side=2,line=6.3,cex=1.6)
			mtext(paste("(",Obs.prev[names(Obs.prev)==sp[i]],"Prevalence )"),
				side=2,line=4.3,cex=1.3)}}}


###################################################
### chunk number 35: fig13plot
###################################################
par(oma=c(0,5,0,0),mar=c(3,3,3,1),mfrow=c(2,3),cex=.7,cex.lab=1.4,mgp=c(2,.5,0))
sp<-c("JUSC2","PICO")
for(i in 1:2){
	DATA=SPDATA[SPDATA$SPECIES==sp[i],]
	for(mod in 1:3){
		calibration.plot(	DATA, 
				which.model=mod, 
				color=mod+1,
				main=model.names[mod],
				xlab="",ylab="")
	if(mod==1){
			mtext(sp[i],side=2,line=6.3,cex=1.6)
			mtext(paste("(",Obs.prev[names(Obs.prev)==sp[i]],"Prevalence )"),
				side=2,line=4.3,cex=1.3)}}}
mtext("Predicted Probability of Occurrence",side=1,line=-1,cex=1.4,outer=TRUE)
mtext("Observed Occurrence as Proportion of Sites Surveyed",side=2,line=-1,cex=1.4,outer=TRUE)


###################################################
### chunk number 36: fig13
###################################################
par(oma=c(0,5,0,0),mar=c(3,3,3,1),mfrow=c(2,3),cex=.7,cex.lab=1.4,mgp=c(2,.5,0))
sp<-c("JUSC2","PICO")
for(i in 1:2){
	DATA=SPDATA[SPDATA$SPECIES==sp[i],]
	for(mod in 1:3){
		calibration.plot(	DATA, 
				which.model=mod, 
				color=mod+1,
				main=model.names[mod],
				xlab="",ylab="")
	if(mod==1){
			mtext(sp[i],side=2,line=6.3,cex=1.6)
			mtext(paste("(",Obs.prev[names(Obs.prev)==sp[i]],"Prevalence )"),
				side=2,line=4.3,cex=1.3)}}}
mtext("Predicted Probability of Occurrence",side=1,line=-1,cex=1.4,outer=TRUE)
mtext("Observed Occurrence as Proportion of Sites Surveyed",side=2,line=-1,cex=1.4,outer=TRUE)


###################################################
### chunk number 37: fig14plot
###################################################
par(mar=c(4,4,4,1))
sp<- "JUSC2"
DATA=SPDATA[SPDATA$SPECIES==sp,]
MAIN<- paste(	sp,"(",Obs.prev[names(Obs.prev)==sp],"Prevalence ) - ",model.names[mod])
presence.absence.summary(	DATA,
					opt.methods=c("MaxKappa","Default"), 
					which.model=mod,
					N.bins=10,N.bars=15,
					main=MAIN,
					legend.cex=.6,
					opt.legend.cex=.6,
					vert.lines=TRUE)


###################################################
### chunk number 38: fig14
###################################################
par(mar=c(4,4,4,1))
sp<- "JUSC2"
DATA=SPDATA[SPDATA$SPECIES==sp,]
MAIN<- paste(	sp,"(",Obs.prev[names(Obs.prev)==sp],"Prevalence ) - ",model.names[mod])
presence.absence.summary(	DATA,
					opt.methods=c("MaxKappa","Default"), 
					which.model=mod,
					N.bins=10,N.bars=15,
					main=MAIN,
					legend.cex=.6,
					opt.legend.cex=.6,
					vert.lines=TRUE)