ENMeval
Rob Harbert and Peter Galante
5/15/2018
ENMeval
Load libraries
require(dismo)## Loading required package: dismo## Loading required package: raster## Loading required package: sprequire(ENMeval)## Loading required package: ENMeval## Loading required package: rJava## Loading required package: parallelrequire(phyloclim)## Loading required package: phyloclim## Loading required package: ape##
## Attaching package: 'ape'## The following objects are masked from 'package:raster':
##
## rotate, zoomrequire(sp)
require(rgdal)## Loading required package: rgdal## Warning: package 'rgdal' was built under R version 3.4.3## rgdal: version: 1.2-20, (SVN revision 725)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
## Path to GDAL shared files: /usr/local/anaconda3/share/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: (autodetected)
## Linking to sp version: 1.2-7require(rgeos)## Loading required package: rgeos## Warning: package 'rgeos' was built under R version 3.4.3## rgeos version: 0.3-26, (SVN revision 560)
## GEOS runtime version: 3.6.2-CAPI-1.10.2 0
## Linking to sp version: 1.2-7
## Polygon checking: TRUErequire(ENMTools)## Loading required package: ENMToolsrequire(spocc)## Loading required package: spoccGet data from GBIF
#occ = occ('Ambystoma opacum', from = 'gbif', limit=10000)
occ = occ('Quercus virginiana', from = 'gbif')## Warning in C_valid_tz(tzone): System timezone name is unknown. Please set
## environment variable TZ.occdf = occ2df(occ)
summary(occ)## <source> gbif
## <time> 2018-05-16 15:56:45
## <found> 1353
## <returned> 500
## <query options>
## scientificName: Quercus virginiana
## limit: 500
##
## <source> bison
## <time>
## <found>
## <returned>
## <query options>
##
## <source> inat
## <time>
## <found>
## <returned>
## <query options>
##
## <source> ebird
## <time>
## <found>
## <returned>
## <query options>
##
## <source> ecoengine
## <time>
## <found>
## <returned>
## <query options>
##
## <source> antweb
## <time>
## <found>
## <returned>
## <query options>
##
## <source> vertnet
## <time>
## <found>
## <returned>
## <query options>
##
## <source> idigbio
## <time>
## <found>
## <returned>
## <query options>
##
## <source> obis
## <time>
## <found>
## <returned>
## <query options>
##
## <source> ala
## <time>
## <found>
## <returned>
## <query options>
## loc = cbind(occdf$longitude, occdf$latitude)
loc = loc[loc[,1]<= -60,]
loc = na.omit(loc)Load WorldClim
# Env = stack("/data/spbio/climgrids/bio.gri")
Env = stack("~/bio.gri")
ext = extent(c(-100, - 60, 25, 55))
Env = crop(Env, ext)Maxent
max = maxent(Env, loc)
pred = predict(Env, max)ENMeval
require(ENMeval)
extr = extract(Env, loc)
loc = loc[!is.na(extr[,1]),]
res = ENMevaluate(occ=loc, env = Env, method='block', parallel=T, numCores=4, fc=c("L", "LQ", "H"), RMvalues=seq(0.5,4,0.5), rasterPreds=F)## Doing evaluations using spatial blocks...## Of 4 total cores using 4## Running in parallel...## Warning in calc.aicc(nparm, occ, predictive.maps): Can't calculate AICc
## when rasterPreds=F - returning NA's.## ENMeval completed in 3 minutes 36.8 seconds.print(res@results)## settings features rm full.AUC Mean.AUC Var.AUC Mean.AUC.DIFF
## 1 L_0.5 L 0.5 0.9697 0.9569684 0.003429567 0.02230974
## 2 LQ_0.5 LQ 0.5 0.9738 0.9574020 0.002265401 0.02123469
## 3 H_0.5 H 0.5 0.5266 0.8696230 0.003337644 0.05431994
## 4 L_1 L 1.0 0.9696 0.9544510 0.003443974 0.02226292
## 5 LQ_1 LQ 1.0 0.9726 0.9570222 0.002139333 0.02159246
## 6 H_1 H 1.0 0.9801 0.9160309 0.017145005 0.04864855
## 7 L_1.5 L 1.5 0.9694 0.9559225 0.003515027 0.02161554
## 8 LQ_1.5 LQ 1.5 0.9721 0.9538468 0.001923180 0.02499574
## 9 H_1.5 H 1.5 0.9800 0.9219303 0.012742220 0.04610627
## 10 L_2 L 2.0 0.9686 0.9570515 0.003648001 0.02137069
## 11 LQ_2 LQ 2.0 0.9718 0.9560926 0.001845243 0.02267767
## 12 H_2 H 2.0 0.9786 0.9517008 0.006418386 0.03275595
## 13 L_2.5 L 2.5 0.9684 0.9572677 0.003656116 0.02095777
## 14 LQ_2.5 LQ 2.5 0.9715 0.9541847 0.001920868 0.02468450
## 15 H_2.5 H 2.5 0.9776 0.9588992 0.003928384 0.02599059
## 16 L_3 L 3.0 0.9681 0.9576261 0.003651724 0.02068371
## 17 LQ_3 LQ 3.0 0.9712 0.9562190 0.001798412 0.02249986
## 18 H_3 H 3.0 0.9762 0.9583676 0.004006370 0.02607509
## 19 L_3.5 L 3.5 0.9679 0.9579271 0.003632170 0.02043099
## 20 LQ_3.5 LQ 3.5 0.9710 0.9573940 0.001786568 0.02130192
## 21 H_3.5 H 3.5 0.9754 0.9576215 0.004131148 0.02638563
## 22 L_4 L 4.0 0.9675 0.9584640 0.003651334 0.02010122
## 23 LQ_4 LQ 4.0 0.9709 0.9587027 0.001725562 0.01977046
## 24 H_4 H 4.0 0.9747 0.9549281 0.003912834 0.02708876
## Var.AUC.DIFF Mean.OR10 Var.OR10 Mean.ORmin Var.ORmin AICc
## 1 0.003117867 0.2752809 0.12342718 0.002808989 3.156167e-05 NA
## 2 0.002317543 0.2584270 0.09030846 0.002808989 3.156167e-05 NA
## 3 0.007989614 0.3174157 0.09328578 0.000000000 0.000000e+00 NA
## 4 0.003333256 0.2921348 0.11825106 0.002808989 3.156167e-05 NA
## 5 0.002094054 0.2612360 0.09631570 0.002808989 3.156167e-05 NA
## 6 0.019312808 0.2696629 0.11118125 0.000000000 0.000000e+00 NA
## 7 0.003290342 0.2808989 0.12851913 0.002808989 3.156167e-05 NA
## 8 0.001677747 0.3061798 0.09724151 0.002808989 3.156167e-05 NA
## 9 0.015494708 0.3033708 0.13138072 0.022471910 2.019947e-03 NA
## 10 0.003267382 0.2696629 0.13407398 0.002808989 3.156167e-05 NA
## 11 0.001584966 0.2921348 0.10411143 0.002808989 3.156167e-05 NA
## 12 0.006402846 0.3005618 0.12484745 0.028089888 3.156167e-03 NA
## 13 0.003193726 0.2668539 0.13031814 0.002808989 3.156167e-05 NA
## 14 0.001584951 0.3089888 0.10053444 0.002808989 3.156167e-05 NA
## 15 0.003567986 0.3005618 0.13158061 0.011235955 5.049867e-04 NA
## 16 0.003098966 0.2612360 0.12989732 0.002808989 3.156167e-05 NA
## 17 0.001489960 0.2949438 0.10330135 0.002808989 3.156167e-05 NA
## 18 0.003631498 0.2949438 0.13132812 0.011235955 5.049867e-04 NA
## 19 0.003070202 0.2471910 0.12994992 0.002808989 3.156167e-05 NA
## 20 0.001470447 0.2893258 0.10877204 0.002808989 3.156167e-05 NA
## 21 0.003751117 0.3005618 0.13360056 0.011235955 5.049867e-04 NA
## 22 0.003024709 0.2331461 0.13158061 0.002808989 3.156167e-05 NA
## 23 0.001456968 0.2724719 0.11382191 0.002808989 3.156167e-05 NA
## 24 0.003892553 0.3146067 0.14619366 0.019662921 1.041535e-03 NA
## delta.AICc w.AIC nparam
## 1 NA NA 18
## 2 NA NA 29
## 3 NA NA 106
## 4 NA NA 13
## 5 NA NA 24
## 6 NA NA 76
## 7 NA NA 13
## 8 NA NA 18
## 9 NA NA 77
## 10 NA NA 11
## 11 NA NA 17
## 12 NA NA 62
## 13 NA NA 10
## 14 NA NA 15
## 15 NA NA 47
## 16 NA NA 11
## 17 NA NA 13
## 18 NA NA 42
## 19 NA NA 11
## 20 NA NA 13
## 21 NA NA 32
## 22 NA NA 11
## 23 NA NA 12
## 24 NA NA 28require(spThin)## Loading required package: spThin## Loading required package: spam## Loading required package: dotCall64## Loading required package: grid## Spam version 2.1-4 (2018-04-12) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.##
## Attaching package: 'spam'## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve## Loading required package: fields## Loading required package: maps## See www.image.ucar.edu/~nychka/Fields for
## a vignette and other supplements.## Loading required package: knitrdf = data.frame(occdf)
extr = extract(Env, cbind(df$longitude, df$latitude))
df = df[!is.na(extr[,1]),]
thin<-thin(loc.data = df,
lat.col = "latitude",
long.col = "longitude",
spec.col = "name",
thin.par = 10,
reps = 10,
locs.thinned.list.return = T,
write.files = T,
max.files = 2,
out.dir = "Thinlocs/",
write.log.file = T)## **********************************************
## Beginning Spatial Thinning.
## Script Started at: Wed May 16 16:00:30 2018
## lat.long.thin.count
## 203
## 10
## [1] "Maximum number of records after thinning: 203"
## [1] "Number of data.frames with max records: 10"
## [1] "Writing new *.csv files"## Warning in thin(loc.data = df, lat.col = "latitude", long.col =
## "longitude", : Writing new *.csv files to output directory: Thinlocs/## Warning in thin(loc.data = df, lat.col = "latitude", long.col = "longitude", : Thinlocs/thinned_data_thin1_new.csv ' already exists. Renaming file
## to avoid overwriting.## Warning in thin(loc.data = df, lat.col = "latitude", long.col = "longitude", : Thinlocs/thinned_data_thin1_new_new.csv ' already exists. Renaming file
## to avoid overwriting.## [1] "Writing file: Thinlocs/thinned_data_thin1_new_new.csv"## Warning in thin(loc.data = df, lat.col = "latitude", long.col = "longitude", : Thinlocs/thinned_data_thin2_new.csv ' already exists. Renaming file
## to avoid overwriting.## Warning in thin(loc.data = df, lat.col = "latitude", long.col = "longitude", : Thinlocs/thinned_data_thin2_new_new.csv ' already exists. Renaming file
## to avoid overwriting.## [1] "Writing file: Thinlocs/thinned_data_thin2_new_new.csv"newthin = read.csv('Thinlocs/thinned_data_thin1.csv')Re-run ENMeval
thinres = ENMevaluate(occ=newthin[,2:3], env = Env, method='block', parallel=T, numCores=4, fc=c("L", "LQ", "H"), RMvalues=seq(1,4,1), rasterPreds=F)## Doing evaluations using spatial blocks...## Of 4 total cores using 4## Running in parallel...## Warning in calc.aicc(nparm, occ, predictive.maps): Can't calculate AICc
## when rasterPreds=F - returning NA's.## ENMeval completed in 2 minutes 8.7 seconds.print(thinres@results)## settings features rm full.AUC Mean.AUC Var.AUC Mean.AUC.DIFF
## 1 L_1 L 1 0.9669 0.9561117 0.002658811 0.02010466
## 2 LQ_1 LQ 1 0.9705 0.9613613 0.001222391 0.01485223
## 3 H_1 H 1 0.9748 0.7980157 0.300175882 0.12336765
## 4 L_2 L 2 0.9659 0.9564613 0.003093781 0.02004019
## 5 LQ_2 LQ 2 0.9699 0.9619874 0.001190614 0.01407087
## 6 H_2 H 2 0.9722 0.9633446 0.001978232 0.01761102
## 7 L_3 L 3 0.9651 0.9568635 0.003175750 0.01946517
## 8 LQ_3 LQ 3 0.9695 0.9617705 0.001184982 0.01378699
## 9 H_3 H 3 0.9710 0.9623671 0.002170074 0.01777951
## 10 L_4 L 4 0.9642 0.9583281 0.002747695 0.01756281
## 11 LQ_4 LQ 4 0.9688 0.9618434 0.001262957 0.01360482
## 12 H_4 H 4 0.9693 0.9610099 0.001978685 0.01762392
## Var.AUC.DIFF Mean.OR10 Var.OR10 Mean.ORmin Var.ORmin AICc
## 1 0.001905733 0.2270588 0.08073264 0.005000000 0.0001000000 NA
## 2 0.001234402 0.1473529 0.05285950 0.005000000 0.0001000000 NA
## 3 0.119443025 0.2278431 0.06961077 0.009901961 0.0001307574 NA
## 4 0.002027161 0.1871569 0.07187547 0.005000000 0.0001000000 NA
## 5 0.001181087 0.1473529 0.06567509 0.005000000 0.0001000000 NA
## 6 0.001594909 0.2024510 0.06198481 0.009901961 0.0001307574 NA
## 7 0.001995562 0.1722549 0.06744517 0.005000000 0.0001000000 NA
## 8 0.001150405 0.1473529 0.06567509 0.005000000 0.0001000000 NA
## 9 0.001713479 0.1624510 0.05797958 0.014705882 0.0008650519 NA
## 10 0.001594124 0.1673529 0.06261626 0.005000000 0.0001000000 NA
## 11 0.001173703 0.1571569 0.07604802 0.005000000 0.0001000000 NA
## 12 0.001676887 0.1624510 0.05797958 0.014705882 0.0008650519 NA
## delta.AICc w.AIC nparam
## 1 NA NA 12
## 2 NA NA 18
## 3 NA NA 54
## 4 NA NA 12
## 5 NA NA 16
## 6 NA NA 45
## 7 NA NA 10
## 8 NA NA 11
## 9 NA NA 31
## 10 NA NA 8
## 11 NA NA 13
## 12 NA NA 26Check out the pattern of model fitting:
eval.plot(thinres@results, "Mean.AUC", )
“Best” models
#minimize ommission rate AND optimize AUC
setsort = res@results[order(res@results[,'Mean.ORmin']),]
setsort2 = setsort[order(setsort[,'Mean.AUC'], decreasing=TRUE),]
top = setsort2[1,]
best = which(as.character(res@results[,1]) == as.character(setsort2[1,1]))
pred.raw = predict(Env, res@models[[best]])
plot(pred.raw, col=viridis::viridis(99))
#and for the thinned model...
setsort = thinres@results[order(thinres@results[,'Mean.ORmin']),]
setsort2 = setsort[order(setsort[,'Mean.AUC'], decreasing=TRUE),]
top = setsort2[1,]
best.thin = which(as.character(thinres@results[,1]) == as.character(setsort2[1,1]))
pred.thin = predict(Env, thinres@models[[best.thin]])
plot(pred.raw, col=viridis::viridis(99))
plot(pred.thin, col=viridis::viridis(99))
Did you get different models parameters? Output? Which was more complex? Which performed better?
print(res@results[best,])## settings features rm full.AUC Mean.AUC Var.AUC Mean.AUC.DIFF
## 15 H_2.5 H 2.5 0.9776 0.9588992 0.003928384 0.02599059
## Var.AUC.DIFF Mean.OR10 Var.OR10 Mean.ORmin Var.ORmin AICc
## 15 0.003567986 0.3005618 0.1315806 0.01123596 0.0005049867 NA
## delta.AICc w.AIC nparam
## 15 NA NA 47print(thinres@results[best.thin,])## settings features rm full.AUC Mean.AUC Var.AUC Mean.AUC.DIFF
## 6 H_2 H 2 0.9722 0.9633446 0.001978232 0.01761102
## Var.AUC.DIFF Mean.OR10 Var.OR10 Mean.ORmin Var.ORmin AICc
## 6 0.001594909 0.202451 0.06198481 0.009901961 0.0001307574 NA
## delta.AICc w.AIC nparam
## 6 NA NA 45Binary Range – Thresholding
#via dismo
#Get straight from ENMeval output:
ev.set <- evaluate(newthin[,2:3], thinres@bg.pts, thinres@models[[best.thin]], Env)
th1 = threshold(ev.set)
p1.nomit = pred.thin >= th1$no_omission
p1.equal = pred.thin >= th1$equal_sens_spec
#OR: calculate evaluation from predicted model output:
occ <- na.omit(extract(pred.thin, newthin[,2:3]));
bg = randomPoints(Env[[1]], n=20000)
bg <- na.omit(extract(pred.thin, bg));
ev2 <- evaluate(occ, bg)
th2 = threshold(ev2)
p2.nomit = pred.thin >= th2$no_omission
p2.equal = pred.thin >= th2$equal_sens_spec
#p1 and p2 should be virtually identicalLook at object ‘th’ for other threshold options.
Niche overlap –
Warren, D.L., R.E. Glor, M. Turelli, and D. Funk. 2009. Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evolution 62:2868-2883; Erratum: Evolution 65: 1215
#Values closer to 1 indicate stronger niche overlap
nicheOverlap(pred.raw, pred.thin, "D", mask=T)## [1] 0.8786465Get background data from minimum convex polygon of occurrence data.
#Minimum convex polygon function
mcp <- function (xy) {
xy <- as.data.frame(coordinates(xy))
coords.t <- chull(xy[, 1], xy[, 2])
xy.bord <- xy[coords.t, ]
xy.bord <- rbind(xy.bord[nrow(xy.bord), ], xy.bord)
return(SpatialPolygons(list(Polygons(list(Polygon(as.matrix(xy.bord))), 1))))
}
MCP.thinlocs = mcp(df[,2:3])
plot(Env[[1]], col=viridis::viridis(99))
plot(MCP.thinlocs, add=T)
envPoly <- rasterToPolygons(Env[[1]], fun=NULL, na.rm=T)
#Get bg
bg.area.thinlocs <- gIntersection(envPoly, MCP.thinlocs)## Warning in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td,
## unaryUnion_if_byid_false, : spgeom1 and spgeom2 have different proj4
## stringsMCP.raster.thinlocs <- rasterize(bg.area.thinlocs, Env[[1]])
bg.points.thinlocs <- randomPoints(mask = MCP.raster.thinlocs, n = 5000)
plot(Env[[1]], col=viridis::viridis(99))
plot(bg.area.thinlocs, add=T)
Use ‘bg.points.thinlocs’ to build a new best model
ENVIREM
Try the Envirem dataset instead of Worldclim
enirem = stack('/data/spbio/climgrids/envirem.gri')Paleoclimate models
http://worldclim.org/paleo-climate1
The climgrids folder has several of the GCM paleoclimate outputs for the Mid Holocene and LGM scenarios provided by Worldclim and Envirem.
Use at least one of these to project a model into the past.
lgm = stack('/data/spbio/climgrids/miroc/lgm_wc.gri')
names(lgm) = names(Env) #Make names match
lgm = crop(lgm, extent(Env))
lgm.pred.thin = predict(lgm, thinres@models[[best.thin]]) Homework 1
- Pick two (or more) new species from a single genus and use ENMeval to find the best model for each.
- Make maps of best model output for each as continuous and binary ranges.
- Project into LGM for all 3 GCMs. (2 sp. x 3 GCMs)
- Calculate range overlap between each LGM projection and the modern as a proxy for geographic range shift.