## ----sdm41a-------------------------------------------------------------------
sdmdata <- readRDS("sdm.Rds")
presvals <- readRDS("pvals.Rds")


## ----sdm25--------------------------------------------------------------------
m1 <- glm(pb ~ bio1 + bio5 + bio12, data=sdmdata)
class(m1)
summary(m1)

m2 = glm(pb ~ ., data=sdmdata)
m2


## ----sdm26--------------------------------------------------------------------
library(predicts)
bc <- envelope(presvals[,c("bio1", "bio5", "bio12")])
bc


## ----sdm27a-------------------------------------------------------------------
bio1 <- c(40, 150, 200)
bio5 <- c(60, 115, 290)
bio12 <- c(600, 1600, 1700)
pd <- data.frame(cbind(bio1, bio5, bio12))
pd
predict(m1, pd)


## ----sdm27b, fig.width=9, fig.height=6----------------------------------------
pr <- partialResponse(bc, presvals, "bio1")
plot(pr, type="l")


## ----sdm27c, fig.width=9, fig.height=6----------------------------------------
library(predicts)
f <- system.file("ex/bio.tif", package="predicts")
predictors <- rast(f)
names(predictors)
p <- predict(predictors, m1)
plot(p)


## ----sdm28, width=7, height=5-------------------------------------------------
p <- rnorm(50, mean=0.7, sd=0.3)
a <- rnorm(50, mean=0.4, sd=0.4)
par(mfrow=c(1, 2))
plot(sort(p), col='red', pch=21)
points(sort(a), col='blue', pch=24)
legend(1, 0.95 * max(a,p), c('presence', 'absence'),
          pch=c(21,24), col=c('red', 'blue'))
comb <- c(p,a)
group <- c(rep('presence', length(p)), rep('absence', length(a)))
boxplot(comb~group, col=c('blue', 'red'))


## ----sdm29--------------------------------------------------------------------
group = c(rep(1, length(p)), rep(0, length(a))) 
cor.test(comb, group)$estimate
mv <- wilcox.test(p,a)
auc <- as.numeric(mv$statistic) / (length(p) * length(a))
auc


## ----sdm40, width=7,height=5--------------------------------------------------
library(predicts)
e <- pa_evaluate(p=p, a=a)
class(e)
e
par(mfrow=c(1, 2))
plot(e, "density")
plot(e, "boxplot", col=c('blue', 'red'))


## ----sdm41--------------------------------------------------------------------
samp <- sample(nrow(sdmdata), round(0.75 * nrow(sdmdata)))
traindata <- sdmdata[samp,]
traindata <- traindata[traindata[,1] == 1, 2:9]
testdata <- sdmdata[-samp,]
bc <- envelope(traindata)
e <- pa_evaluate(testdata[testdata==1,], testdata[testdata==0,], bc)
e
plot(e, "ROC")


## ----sdm42--------------------------------------------------------------------
pres <- sdmdata[sdmdata[,1] == 1, 2:9]
back <- sdmdata[sdmdata[,1] == 0, 2:9]


## ----sdm43--------------------------------------------------------------------
k <- 5
group <- folds(pres, k)
group[1:10]
unique(group)
table(group)


## ----sdm44--------------------------------------------------------------------
e <- list()
for (i in 1:k) {
	train <- pres[group != i,]
	test <- pres[group == i,]
	bc <- envelope(train)
	p <- predict(bc, test)
	a <- predict(bc, back)
	e[[i]] <- pa_evaluate(p, a)
}	


## ----sdm45a-------------------------------------------------------------------
stats <- lapply(e, function(x){x@stats})
stats <- do.call(rbind, stats)
colMeans(stats)


## ----sdm45b-------------------------------------------------------------------
fsp <- system.file("/ex/bradypus.csv", package="predicts")
bradypus <- read.csv(fsp)
bradypus <- bradypus[,-1]
presvals <- extract(predictors, bradypus)
set.seed(0)
backgr <- spatSample(predictors, 500, na.rm=TRUE, xy=TRUE, values=FALSE)

nr <- nrow(bradypus)
s <- sample(nr, 0.25 * nr)
pres_train <- bradypus[-s, ]
pres_test <- bradypus[s, ]

nr <- nrow(backgr)
set.seed(9)
s <- sample(nr, 0.25 * nr)
back_train <- backgr[-s, ]
back_test <- backgr[s, ]


## ----sdm45c-------------------------------------------------------------------
library(dismo)
sb <- dismo::ssb(pres_test, back_test, pres_train)
sb[,1] / sb[,2]


## ----sdm45d-------------------------------------------------------------------
i <- pwd_sample(pres_test, back_test, pres_train, n=1, tr=0.1)
pres_test_pwd <- pres_test[!is.na(i[,1]), ]
back_test_pwd <- back_test[na.omit(as.vector(i)), ]
sb2 <- dismo::ssb(pres_test_pwd, back_test_pwd, pres_train)
sb2[1]/ sb2[2]


## ----sdm45e-------------------------------------------------------------------
bc <- envelope(predictors, pres_train)
ptst <- predict(bc, extract(predictors, pres_test))
btst <- predict(bc, extract(predictors, back_test))
pa_evaluate(ptst, btst)@stats

pwdptst <- predict(bc, extract(predictors, pres_test_pwd))
pwdbtst <- predict(bc, extract(predictors, back_test_pwd))
pa_evaluate(pwdptst, pwdbtst)@stats