## ---- echo=FALSE, include=FALSE----------------------------------------------- library(maptools) ## ---- sdm46a------------------------------------------------------------------ library(dismo) library(maptools) data(wrld_simpl) predictors <- stack(list.files(file.path(system.file(package="dismo"), 'ex'), pattern='grd$', full.names=TRUE )) 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) absvals <- extract(predictors, backgr) pb <- c(rep(1, nrow(presvals)), rep(0, nrow(absvals))) sdmdata <- data.frame(cbind(pb, rbind(presvals, absvals))) sdmdata[,'biome'] <- as.factor(sdmdata[,'biome']) ## ---- sdm46b------------------------------------------------------------------ pred_nf <- dropLayer(predictors, 'biome') ## ---- sdm47------------------------------------------------------------------- set.seed(0) group <- kfold(bradypus, 5) pres_train <- bradypus[group != 1, ] pres_test <- bradypus[group == 1, ] ## ---- sdm48------------------------------------------------------------------- ext <- extent(-90, -32, -33, 23) ## ---- sdm49------------------------------------------------------------------- set.seed(10) backg <- randomPoints(pred_nf, n=1000, ext=ext, extf = 1.25) colnames(backg) = c('lon', 'lat') group <- kfold(backg, 5) backg_train <- backg[group != 1, ] backg_test <- backg[group == 1, ] ## ---- sdm50------------------------------------------------------------------- r <- raster(pred_nf, 1) plot(!is.na(r), col=c('white', 'light grey'), legend=FALSE) plot(ext, add=TRUE, col='red', lwd=2) points(backg_train, pch='-', cex=0.5, col='yellow') points(backg_test, pch='-', cex=0.5, col='black') points(pres_train, pch= '+', col='green') points(pres_test, pch='+', col='blue') ## ---- sdm60------------------------------------------------------------------- bc <- bioclim(pred_nf, pres_train) plot(bc, a=1, b=2, p=0.85) ## ---- sdm61a------------------------------------------------------------------ e <- evaluate(pres_test, backg_test, bc, pred_nf) e ## ---- sdm61b------------------------------------------------------------------ tr <- threshold(e, 'spec_sens') tr ## ---- sdm62, fig.width=9, fig.height=6---------------------------------------- pb <- predict(pred_nf, bc, ext=ext, progress='') pb par(mfrow=c(1,2)) plot(pb, main='Bioclim, raw values') plot(wrld_simpl, add=TRUE, border='dark grey') plot(pb > tr, main='presence/absence') plot(wrld_simpl, add=TRUE, border='dark grey') points(pres_train, pch='+') ## ---- sdm63, fig.width=9, fig.height=6---------------------------------------- dm <- domain(pred_nf, pres_train) e <- evaluate(pres_test, backg_test, dm, pred_nf) e pd = predict(pred_nf, dm, ext=ext, progress='') par(mfrow=c(1,2)) plot(pd, main='Domain, raw values') plot(wrld_simpl, add=TRUE, border='dark grey') tr <- threshold(e, 'spec_sens') plot(pd > tr, main='presence/absence') plot(wrld_simpl, add=TRUE, border='dark grey') points(pres_train, pch='+') ## ---- sdm64, fig.width=9, fig.height=6---------------------------------------- mm <- mahal(pred_nf, pres_train) e <- evaluate(pres_test, backg_test, mm, pred_nf) e pm = predict(pred_nf, mm, ext=ext, progress='') par(mfrow=c(1,2)) pm[pm < -10] <- -10 plot(pm, main='Mahalanobis distance') plot(wrld_simpl, add=TRUE, border='dark grey') tr <- threshold(e, 'spec_sens') plot(pm > tr, main='presence/absence') plot(wrld_simpl, add=TRUE, border='dark grey') points(pres_train, pch='+') ## ---- sdm65------------------------------------------------------------------- train <- rbind(pres_train, backg_train) pb_train <- c(rep(1, nrow(pres_train)), rep(0, nrow(backg_train))) envtrain <- extract(predictors, train) envtrain <- data.frame( cbind(pa=pb_train, envtrain) ) envtrain[,'biome'] = factor(envtrain[,'biome'], levels=1:14) head(envtrain) testpres <- data.frame( extract(predictors, pres_test) ) testbackg <- data.frame( extract(predictors, backg_test) ) testpres[ ,'biome'] = factor(testpres[ ,'biome'], levels=1:14) testbackg[ ,'biome'] = factor(testbackg[ ,'biome'], levels=1:14) ## ---- sdm66------------------------------------------------------------------- # logistic regression: gm1 <- glm(pa ~ bio1 + bio5 + bio6 + bio7 + bio8 + bio12 + bio16 + bio17, family = binomial(link = "logit"), data=envtrain) summary(gm1) coef(gm1) gm2 <- glm(pa ~ bio1+bio5 + bio6 + bio7 + bio8 + bio12 + bio16 + bio17, family = gaussian(link = "identity"), data=envtrain) evaluate(testpres, testbackg, gm1) ge2 <- evaluate(testpres, testbackg, gm2) ge2 ## ---- sdm67, fig.width=9, fig.height=6---------------------------------------- pg <- predict(predictors, gm2, ext=ext) par(mfrow=c(1,2)) plot(pg, main='GLM/gaussian, raw values') plot(wrld_simpl, add=TRUE, border='dark grey') tr <- threshold(ge2, 'spec_sens') plot(pg > tr, main='presence/absence') plot(wrld_simpl, add=TRUE, border='dark grey') points(pres_train, pch='+') points(backg_train, pch='-', cex=0.25) ## ---- sdm68a, fig.width=9, fig.fig.height=6----------------------------------- maxent() xm <- maxent(predictors, pres_train, factors='biome') plot(xm) ## ---- sdm68b, fig.width=9, fig.height=6--------------------------------------- response(xm) ## ---- sdm69, fig.width=9, fig.height=6---------------------------------------- e <- evaluate(pres_test, backg_test, xm, predictors) e px <- predict(predictors, xm, ext=ext, progress='') par(mfrow=c(1,2)) plot(px, main='Maxent, raw values') plot(wrld_simpl, add=TRUE, border='dark grey') tr <- threshold(e, 'spec_sens') plot(px > tr, main='presence/absence') plot(wrld_simpl, add=TRUE, border='dark grey') points(pres_train, pch='+') ## ---- sdm80, fig.width=9, fig.height=6---------------------------------------- library(randomForest) model <- pa ~ bio1 + bio5 + bio6 + bio7 + bio8 + bio12 + bio16 + bio17 rf1 <- randomForest(model, data=envtrain) model <- factor(pa) ~ bio1 + bio5 + bio6 + bio7 + bio8 + bio12 + bio16 + bio17 rf2 <- randomForest(model, data=envtrain) rf3 <- randomForest(envtrain[,1:8], factor(pb_train)) erf <- evaluate(testpres, testbackg, rf1) erf pr <- predict(predictors, rf1, ext=ext) par(mfrow=c(1,2)) plot(pr, main='Random Forest, regression') plot(wrld_simpl, add=TRUE, border='dark grey') tr <- threshold(erf, 'spec_sens') plot(pr > tr, main='presence/absence') plot(wrld_simpl, add=TRUE, border='dark grey') points(pres_train, pch='+') points(backg_train, pch='-', cex=0.25) ## ---- sdm81, fig.width=9, fig.height=6---------------------------------------- library(kernlab) svm <- ksvm(pa ~ bio1+bio5+bio6+bio7+bio8+bio12+bio16+bio17, data=envtrain) esv <- evaluate(testpres, testbackg, svm) esv ps <- predict(predictors, svm, ext=ext) par(mfrow=c(1,2)) plot(ps, main='Support Vector Machine') plot(wrld_simpl, add=TRUE, border='dark grey') tr <- threshold(esv, 'spec_sens') plot(ps > tr, main='presence/absence') plot(wrld_simpl, add=TRUE, border='dark grey') points(pres_train, pch='+') points(backg_train, pch='-', cex=0.25) ## ---- sdm82, fig.width=9, fig.height=6---------------------------------------- models <- stack(pb, pd, pm, pg, pr, ps) names(models) <- c("bioclim", "domain", "mahal", "glm", "rf", "svm") plot(models) ## ---- sdm83, fig.width=9, fig.height=6---------------------------------------- m <- mean(models) plot(m, main='average score') ## ---- sdm84, fig.width=9, fig.height=6---------------------------------------- auc <- sapply(list(ge2, erf, esv), function(x) x@auc) w <- (auc-0.5)^2 m2 <- weighted.mean( models[[c("glm", "rf", "svm")]], w) plot(m2, main='weighted mean of three models')