========================================================= Valeria Di Cola, Olivier Broennimann, Blaise Petitpierre, Manuela D'Amen, Frank Breiner & Antoine Guisan
Miscellaneous methods and utilities for spatial ecology analysis, written by current and former members and collaborators of the ecospat group of Antoine Guisan, Department of Ecology and Evolution (DEE) & Institute of Earth Surface Dynamics (IDYST), University of Lausanne, Switzerland.
ecospat offers the possibility to perform Pre-modelling Analysis, such as Spatial autocorrelation analysis, MESS (Multivariate Environmental Similarity Surfaces) analyses, Phylogenetic diversity Measures, Biotic Interactions. It also provides functions to complement biomod2 in preparing the data, calibrating and evaluating (e.g. boyce index) and projecting the models. Complementary analysis based on model predictions (e.g. co-occurrences analyses) are also provided.
In addition, the ecospat package includes Niche Quantification and Overlap functions that were used in Broennimann et al. 2012 and Petitpierre et al. 2012 to quantify climatic niche shifts between the native and invaded ranges of invasive species.
library(ecospat) citation("ecospat")
ecospat.testData()
data(ecospat.testData) names(ecospat.testData)
ecospat.testNiche.inv()
data(ecospat.testNiche.inv) names(ecospat.testNiche.inv)
ecospat.testNiche.nat()
data(ecospat.testNiche.nat) names(ecospat.testNiche.nat)
ecospat.testTree()
if(requireNamespace("ape")){ fpath <- system.file("extdata", "ecospat.testTree.tre", package="ecospat") tree<-ape::read.tree(fpath) tree$tip.label plot(tree, cex=0.6) }
ecospat.mantel.correlogram(dfvar=ecospat.testData[c(2:16)],colxy=1:2, n=100, colvar=3:7, max=1000, nclass=10, nperm=100)
The graph indicates that spatial autocorrelation (SA) is minimal at a distance of 180 meters. Note however that SA is not significantly different than zero for several distances (open circles).
colvar <- ecospat.testData[c(4:8)] x <- cor(colvar, method="pearson") ecospat.npred (x, th=0.75)
x <- cor(colvar, method="spearman") ecospat.npred (x, th=0.75)
x <- ecospat.testData[c(4:8)] p<- x[1:90,] #A projection dataset. ref<- x[91:300,] # A reference dataset
ecospat.climan(ref,p)
x <- ecospat.testData[c(2,3,4:8)] proj<- x[1:90,] #A projection dataset. cal<- x[91:300,] #A calibration dataset mess.object<-ecospat.mess (proj, cal, w="default") ecospat.plot.mess (mess.object, cex=1, pch=15)
In the MESS plot pixels in red indicate sites where at least one environmental predictor has values outside of the range of that predictor in the calibration dataset. In the MESSw plot, same as previous plot but with weighted by the number of predictors.Finally, the MESSneg plot shows at each site how many predictors have values outside of their calibration range.
if(requireNamespace("ape")){ fpath <- system.file("extdata", "ecospat.testTree.tre", package="ecospat") tree <- ape::read.tree(fpath) data <- ecospat.testData[9:52] pd<- ecospat.calculate.pd(tree, data, method = "spanning", type = "species", root = TRUE, average = FALSE, verbose = FALSE ) plot(pd) }
First we load the test data for the niche dynamics analysis in invaded and native range. A PCA is calibrated on all the sites of the study area, including both native and invaded ranges (same as PCAenv in Broenniman et al. 2012). Finally, we plot the variables Contributions
library(ade4) inv <- ecospat.testNiche.inv nat <- ecospat.testNiche.nat pca.env <- ade4::dudi.pca(rbind(nat,inv)[,3:10],scannf=F,nf=2) ecospat.plot.contrib(contrib=pca.env$co, eigen=pca.env$eig)
The correlation circle indicate the contribution of original predictors to the PCA axes.
Now we can predict the scores on the axes
# PCA scores for the whole study area scores.globclim <- pca.env$li # PCA scores for the species native distribution scores.sp.nat <- ade4::suprow(pca.env,nat[which(nat[,11]==1),3:10])$li # PCA scores for the species invasive distribution scores.sp.inv <- ade4::suprow(pca.env,inv[which(inv[,11]==1),3:10])$li # PCA scores for the whole native study area scores.clim.nat <- ade4::suprow(pca.env,nat[,3:10])$li # PCA scores for the whole invaded study area scores.clim.inv <- ade4::suprow(pca.env,inv[,3:10])$li
For a species in the native range (North America)
# gridding the native niche grid.clim.nat <- ecospat.grid.clim.dyn(glob=scores.globclim, glob1=scores.clim.nat, sp=scores.sp.nat, R=100, th.sp=0)
For a species in the invaded range (Australia)
# gridding the invasive niche grid.clim.inv <- ecospat.grid.clim.dyn(glob=scores.globclim, glob1=scores.clim.inv, sp=scores.sp.inv, R=100, th.sp=0)
# Compute Schoener's D, index of niche overlap D.overlap <- ecospat.niche.overlap (grid.clim.nat, grid.clim.inv, cor = TRUE)$D D.overlap
The niche overlap between the native and the invaded range is 22%.
It is reccomended to use at least 1000 replications for the equivalency test. As an example we used rep = 10, to reduce the computational time.
eq.test <- ecospat.niche.equivalency.test(grid.clim.nat, grid.clim.inv,rep=10, intersection = 0.1, overlap.alternative = "higher", expansion.alternative = "lower", stability.alternative = "higher", unfilling.alternative = "lower")
Niche equivalency test H1: the observed overlap between the native and invaded niche is higher than if the two niches are randomized, the niche expansion of the invaded niche is lower, the niche stability is higher and the niche unfilling is lower than if the two niches are randomized.
Plot Equivalency test
ecospat.plot.overlap.test(eq.test, "D", "Equivalency")
Shifts randomly on niche (here the invasive niche) in the study area It is recomended to use at least 1000 replications for the similarity test. As an example we used rep = 10, to reduce the computational time.
sim.test <- ecospat.niche.similarity.test(grid.clim.nat, grid.clim.inv,rep=10, overlap.alternative = "higher", expansion.alternative = "lower", stability.alternative = "higher", unfilling.alternative = "lower", intersection = 0.1, rand.type=1)
Niche similarity test H1: the observed overlap between the native and invaded is higher than randomly shifted invasive niches in the invaded study area, the niche expansion of the invaded niche is lower, the niche stability is higher and the niche unfilling is lower than if the two niches are randomized.
Plot Similarity test
ecospat.plot.overlap.test(sim.test, "D", "Similarity")
niche.dyn <- ecospat.niche.dyn.index (grid.clim.nat, grid.clim.inv, intersection = 0.1)
Plot niche overlap
ecospat.plot.niche.dyn(grid.clim.nat, grid.clim.inv, quant=0.25, interest=2, title= "Niche Overlap", name.axis1="PC1", name.axis2="PC2") ecospat.shift.centroids(scores.sp.nat, scores.sp.inv, scores.clim.nat, scores.clim.inv)
Plot Similarity test for niche expansion, stability and unfilling
ecospat.plot.overlap.test(sim.test, "expansion", "Similarity") ecospat.plot.overlap.test(sim.test, "stability", "Similarity") ecospat.plot.overlap.test(sim.test, "unfilling", "Similarity")
# gridding the native niche grid.clim.t.nat <- ecospat.grid.clim.dyn(glob=as.data.frame(rbind(nat,inv)[,10]), glob1=as.data.frame(nat[,10]), sp=as.data.frame(nat[which(nat[,11]==1),10]), R=1000, th.sp=0) # gridding the invaded niche grid.clim.t.inv <- ecospat.grid.clim.dyn(glob=as.data.frame(rbind(nat,inv)[,10]), glob1=as.data.frame(inv[,10]), sp=as.data.frame(inv[which(inv[,11]==1),10]), R=1000, th.sp=0) t.dyn<-ecospat.niche.dyn.index (grid.clim.t.nat, grid.clim.t.inv, intersection=0.1) ecospat.plot.niche.dyn(grid.clim.t.nat, grid.clim.t.inv, quant=0, interest=2, title= "Niche Overlap", name.axis1="Average temperature")
data <- ecospat.testData[c(9:16,54:57)]
For each pair of species (sp1, sp2), the number (N) of plots where both species were present is divided by the number of plots where the rarest of the two species is present. This index ranges from 0 (no co-occurrence) to 1 (always in co-occurrence) as given in eq. 1.
where N(S1 intersects S2) is the number of times species S1 and S2 co-occur, while Min(NS1, NS2) is the number of times species S1 and S2 co-occur, while is the occurrence frequency of the rarest of the two species.
ecospat.co_occurrences (data)
This function allows to apply a pairwise null model analysis to a presence-absence community matrix to determine which species associations are significant across the study area. The strength of associations is quantified by the C-score index and a 'fixed-equiprobable' null model algorithm is applied.
It is recomended to use at least 10000 permutatiobns for the test. As an example we used nperm = 100, to reduce the computational time.
data<- ecospat.testData[c(53,62,58,70,61,66,65,71,69,43,63,56,68,57,55,60,54,67,59,64)] nperm <- 100 outpath <- getwd() ecospat.Cscore(data, nperm, outpath)
The function returns the C-score index for the observed community (ObsCscoreTot), p.value (PValTot) and standardized effect size (SES.Tot). It saves also a table in the working directory where the same metrics are calculated for each species pair (only the table with species pairs with significant p-values is saved in this version)
data <- ecospat.testData[,4:8] ecospat.cor.plot(data)
A scatter plot of matrices, with bivariate scatter plots below the diagonal, histograms on the diagonal, and the Pearson correlation above the diagonal. Useful for descriptive statistics of small data sets (better with less than 10 variables).
data <- ecospat.testData caleval <- ecospat.caleval (data = ecospat.testData[53], xy = data[2:3], row.num = 1:nrow(data), nrep = 2, ratio = 0.7, disaggregate = 0.2, pseudoabs = 100, npres = 10, replace = FALSE) head(caleval)
We obtained an evaluation and calibration dataset with a desired ratio of disaggregation.
The argument fit is a vector containing the predicted suitability values
fit <- ecospat.testData$glm_Saxifraga_oppositifolia
The argument obs is a vector containing the predicted suitability values of the validation points (presence records)
obs<-ecospat.testData$glm_Saxifraga_oppositifolia[which(ecospat.testData$Saxifraga_oppositifolia==1)]
Calculate and plot Boyce Index with ecospat.boyce
ecospat.boyce (fit, obs, nclass = 0, window.w = "default", res = 100, PEplot = TRUE)$cor
Here the boyce index is 0.91. If the rank of predicted expected ratio would be completely ordered along habitat suitability axis then boyce index would be 1.
Indices of accuracy of community predictions ecospat.CommunityEval()
eval<-ecospat.testData[c(53,62,58,70,61,66,65,71,69,43,63,56,68,57,55,60,54,67,59,64)] pred<-ecospat.testData[c(73:92)]
CommunityEval<-ecospat.CommunityEval (eval, pred, proba = TRUE, ntir=5,verbose = T)
library(biomod2)
# species # occurrences xy <- inv[,1:2] head(xy) sp_occ <- inv[11] # env current <- inv[3:7] head(current) ## BIOMOD t1 <- Sys.time() sp<-1
### Formating the data with the BIOMOD_FormatingData() function form the package biomod2 myBiomodData <- biomod2::BIOMOD_FormatingData( resp.var = as.numeric(sp_occ[,sp]), expl.var = current, resp.xy = xy, resp.name = colnames(sp_occ)[sp]) myBiomodOption <- biomod2::bm_DefaultModelingOptions() myBiomodOption@GLM$test = 'none' myBiomodOption@GBM$interaction.depth = 2
### Calibration of simple bivariate models # remove insivible(capture.output)) to print output in the console # this is just to keep the vignette short invisible(capture.output(my.ESM <- ecospat.ESM.Modeling( data=myBiomodData, models=c('GLM'), models.options=myBiomodOption, NbRunEval=2, DataSplit=70, weighting.score=c("AUC"), parallel=F) ) )
### Evaluation and average of simple bivariate models to ESMs my.ESM_EF <- ecospat.ESM.EnsembleModeling(my.ESM,weighting.score=c("SomersD"),threshold=0)
### Projection of simple bivariate models into new space my.ESM_proj_current <- ecospat.ESM.Projection(ESM.modeling.output=my.ESM, new.env=current)
### Projection of calibrated ESMs into new space my.ESM_EFproj_current <- ecospat.ESM.EnsembleProjection(ESM.prediction.output=my.ESM_proj_current, ESM.EnsembleModeling.output=my.ESM_EF)
Input data for the first argument (proba) as data frame of rough probabilities from SDMs for all species in columns in the considered sites in rows.
proba <- ecospat.testData[,73:92]
Input data for the second argument (sr) as data frame with richness value in the first column and sites.
sr <- as.data.frame(rowSums(proba))
prr<-ecospat.SESAM.prr(proba, sr) head(prr)[,1:4]
Input data as a matrix of plots (rows) x species (columns). Input matrices should have column names (species names) and row names (sampling plots).
presence<-ecospat.testData[c(53,62,58,70,61,66,65,71,69,43,63,56,68,57,55,60,54,67,59,64)] pred<-ecospat.testData[c(73:92)]
Define the number of permutations. It is recommended to use at least 10000 permutations for the test. As an example we used nperm = 100, to reduce the computational time. Then Define the outpath. Then we can run Run the function ecospat.cons_Cscore. The function tests for non-random patterns of species co-occurrence in a presence-absence matrix. It calculates the C-score index for the whole community and for each species pair. An environmental constraint is applied during the generation of the null communities.
nbpermut <- 100 outpath <- getwd() ecospat.cons_Cscore(presence, pred, nbpermut, outpath)
The function returns - the C-score index for the observed community (ObsCscoreTot), - the mean of C-score for the simulated communities (SimCscoreTot), - the p.values (PVal.less and PVal.greater) to evaluate the significance of the difference between the former two indices. - the standardized effect size for the whole community (SES.Tot). A SES that is greater than 2 or less than -2 is statistically significant with a tail probability of less than 0.05 (Gotelli & McCabe 2002 - Ecology). If a community is structured by competition, we would expect the C-score to be large relative to a randomly assembled community (positive SES). In this case the observed C-score is significantly lower than expected by chance, this meaning that the community is dominate by positive interactions (aggregated pattern).
A table is saved in the path specified where the same metrics are calculated for each species pair (only the table with species pairs with significant p.values is saved).
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.