1 Load data

1.1 Package

library(maptools)
## Warning: package 'maptools' was built under R version 3.3.3
## Loading required package: sp
## Checking rgeos availability: TRUE
library(maps)
library(mapdata)
## Warning: package 'mapdata' was built under R version 3.3.3
library(dismo)
## Loading required package: raster
library(raster)

1.2 Load csv file

pachy <- read.csv("pachystachya.csv")

1.3 Load map

data(wrld_simpl)

2 Plot distribution map

long_limit <- c(-76,-35)
lat_limit <- c(-40,15)
plot(wrld_simpl, xlim=long_limit, ylim=lat_limit, axes=TRUE, col="light yellow")

# plot the points
points(pachy$longitude, pachy$latitude, col="darkolivegreen4", pch=20, cex=0.8)

3 Make 1000 pseudo-absences points form the circles 50km

#' draw circle 50km diam.
cc = circles(pachy[,c("longitude","latitude")], d=50000, lonlat=T, dissolve=FALSE)

# draw random points that must fall within the circles in object cc
bg = spsample(cc@polygons, 1000, type='random', iter=1000)

Plot

plot(wrld_simpl, xlim=long_limit, ylim=lat_limit, axes=TRUE, col="light yellow")

# plot the presence points
points(pachy$longitude, pachy$latitude, col="darkolivegreen4", pch=20, cex=0.8)
# plot pseudo-absences points:
points(bg, col="orange", pch=17, cex=0.5)

4 Download data bioclimatic variables 2.5 minutes from worldclim

#library(raster) # package to handle raster-formatted spatial data
BClim = getData("worldclim", var="bio", res=2.5, path="data/")

Crop database for region Neotropic

region = extent(-90,-30,-40,15) # define the extent
BClim = crop(BClim, region)

Save raster

writeRaster(BClim, filename="data/neotropic.grd", overwrite=T)
BClim = brick("data/neotropic.grd")

Plot data climat

# this format plots the first (of 19) variables stored in BClim; change the 1 to 2-19 for the others
plot(BClim, 1, cex=0.5, legend=T, mar=par("mar"), xaxt="n", yaxt="n", main="Annual mean temperature (ºC x 10)", add=T)


# plot the presence points
points(pachy$longitude, pachy$latitude, col="darkgreen", pch=20, cex=0.5)
# and the pseudo-absence points
points(bg, cex=0.5, col="darkorange3")
# add axes
axis(1,las=1)
axis(2,las=1)

Or we can plot for all variable

plot(BClim)

plot(BClim, c(17,18,19))

Bioclimatic variables:

BIO1 = Annual Mean Temperature

BIO2 = Mean Diurnal Range (Mean of monthly (max temp - min temp))

BIO3 = Isothermality (BIO2/BIO7) (* 100)

BIO4 = Temperature Seasonality (standard deviation *100)

BIO5 = Max Temperature of Warmest Month

BIO6 = Min Temperature of Coldest Month

BIO7 = Temperature Annual Range (BIO5-BIO6)

BIO8 = Mean Temperature of Wettest Quarter

BIO9 = Mean Temperature of Driest Quarter

BIO10 = Mean Temperature of Warmest Quarter

BIO11 = Mean Temperature of Coldest Quarter

BIO12 = Annual Precipitation

BIO13 = Precipitation of Wettest Month

BIO14 = Precipitation of Driest Month

BIO15 = Precipitation Seasonality (Coefficient of Variation)

BIO16 = Precipitation of Wettest Quarter

BIO17 = Precipitation of Driest Quarter

BIO18 = Precipitation of Warmest Quarter

BIO19 = Precipitation of Coldest Quarter

5 Pulling bioclim values

# for the subsampled presence points
Ybrev_bc = extract(BClim, pachy[,c("longitude","latitude")]) 
# for the pseudo-absence points
bg_bc = extract(BClim, bg) 

#You’ll probably want to build some useful dataframes, with two columns of coordinates followed by the 19 bioclim variables. First, for the presence points:
Ybrev_bc = data.frame(lon=pachy$longitude, lat=pachy$latitude, Ybrev_bc)

#And then for the pseudo-absences:
bgpoints = bg@coords
colnames(bgpoints) = c("lon","lat")
bg_bc = data.frame(cbind(bgpoints,bg_bc))

# double-check for missing data
length(which(is.na(bg_bc$bio1)))
## [1] 73
# and pull out the missing lines
bg_bc = bg_bc[!is.na(bg_bc$bio1), ]

6 Prepare our database for test with MAXENT

group_p = kfold(Ybrev_bc, 5) # vector of group assignments splitting the Ybrev_bc into 5 groups
group_a = kfold(bg_bc, 5) # ditto for bg_bc

test = 3
train_p = Ybrev_bc[group_p!=test, c("lon","lat")]
train_a = bg_bc[group_a!=test, c("lon","lat")]
test_p = Ybrev_bc[group_p==test, c("lon","lat")]
test_a = bg_bc[group_a==test, c("lon","lat")]

7 MAXENT

We use function maxent in package dismo. Read help ?maxent

# run Maxent for train_p and train_a
me = maxent(BClim, p=train_p, a=train_a)
# test statistic
e = evaluate(test_p, test_a, me, BClim)
e
## class          : ModelEvaluation 
## n presences    : 89 
## n absences     : 185 
## AUC            : 0.6606134 
## cor            : 0.2286534 
## max TPR+TNR at : 0.6183459
plot(e, 'ROC')

plot(e, 'TPR')

boxplot(e)

density(e)

Generate the predictions

pred_me = predict(me, BClim)

8 make a nice plot

plot(pred_me, cex=0.5, legend=T, mar=par("mar"), xaxt="n", yaxt="n")
plot(wrld_simpl, xlim=long_limit, ylim=lat_limit, axes=TRUE, add=T)

# plot the presence points
points(pachy$longitude, pachy$latitude, col="red", pch=20, cex=0.5)
# and the pseudo-absence points
points(bg, cex=0.5, col="darkorange3")
# add axes
axis(1,las=1)
axis(2,las=1)