## ---- echo=FALSE, include=FALSE----------------------------------------------- library(knitr) opts_chunk$set(fig.width = 6, fig.height = 6, fig.cap='', collapse = TRUE) library(dismo) ## ---- sdm41a------------------------------------------------------------------ library(dismo) 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------------------------------------------------------------------- bc <- bioclim(presvals[,c('bio1', 'bio5', 'bio12')]) class(bc) bc pairs(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) predict(bc, pd) ## ---- sdm27b, fig.width=9, fig.height=6--------------------------------------- response(bc) ## ---- sdm27c, fig.width=9, fig.height=6--------------------------------------- predictors <- stack(list.files(file.path(system.file(package="dismo"), 'ex'), pattern='grd$', full.names=TRUE )) 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------------------------------------------------- e <- evaluate(p=p, a=a) class(e) e par(mfrow=c(1, 2)) density(e) boxplot(e, 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 <- bioclim(traindata) e <- 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 <- kfold(pres, k) group[1:10] unique(group) ## ---- sdm44------------------------------------------------------------------- e <- list() for (i in 1:k) { train <- pres[group != i,] test <- pres[group == i,] bc <- bioclim(train) e[[i]] <- evaluate(p=test, a=back, bc) } ## ---- sdm45a------------------------------------------------------------------ auc <- sapply(e, function(x){x@auc}) auc mean(auc) sapply( e, function(x){ threshold(x)['spec_sens'] } ) ## ---- sdm45b------------------------------------------------------------------ file <- file.path(system.file(package="dismo"), "ex/bradypus.csv") bradypus <- read.table(file, header=TRUE, sep=',') bradypus <- bradypus[,-1] presvals <- extract(predictors, bradypus) set.seed(0) backgr <- randomPoints(predictors, 500) 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------------------------------------------------------------------ sb <- ssb(pres_test, back_test, pres_train) sb[,1] / sb[,2] ## ---- sdm45d------------------------------------------------------------------ i <- pwdSample(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 <- ssb(pres_test_pwd, back_test_pwd, pres_train) sb2[1]/ sb2[2] ## ---- sdm45e------------------------------------------------------------------ bc <- bioclim(predictors, pres_train) evaluate(bc, p=pres_test, a=back_test, x=predictors) evaluate(bc, p=pres_test_pwd, a=back_test_pwd, x=predictors)