Vignette ecospat package

In ecospat: Spatial Ecology Miscellaneous Methods

========================================================= Valeria Di Cola, Olivier Broennimann, Blaise Petitpierre, Manuela D'Amen, Frank Breiner & Antoine Guisan

08 November, 2017

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.

Load data

library(ecospat)
citation("ecospat")

Test data for the ecospat library

ecospat.testData()

data(ecospat.testData)
names(ecospat.testData)

Test data for the Niche Overlap Analysis

ecospat.testNiche.inv()

data(ecospat.testNiche.inv)
names(ecospat.testNiche.inv)

ecospat.testNiche.nat()

data(ecospat.testNiche.nat)
names(ecospat.testNiche.nat)

Test tree for Phylogenetic Diversity Analysis

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)
}

Pre-Modelling Analysis

Spatial Auto-correlation

Mantel Correlogram with ecospat.mantel.correlogram()

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).

Predictor Variable Selection

Number of Predictors with Pearson Correlation ecospat.npred()

colvar <- ecospat.testData[c(4:8)]
x <- cor(colvar, method="pearson")
ecospat.npred (x, th=0.75)

Number of Predictors with Spearman Correlation ecospat.npred()

x <- cor(colvar, method="spearman")
ecospat.npred (x, th=0.75)

Climate Analogy Tools

Climate Analogy with ecospat.climan()

x <- ecospat.testData[c(4:8)]
p<- x[1:90,] #A projection dataset.
ref<- x[91:300,] # A reference dataset

ecospat.climan(ref,p)

Extrapolation detection, creating a MESS object with ecospat.mess()

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.

Phylogenetic Diversity Measures

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)
}

Niche Quantification and Comparison with Ordination techniques

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

Calculate the Occurrence Densities Grid with ecospat.grid.clim.dyn()

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) 

Calculate Niche Overlap with ecospat.niche.overlap()

# 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%.

Perform the Niche Equivalency Test with ecospat.niche.equivalency.test() according to Warren et al. (2008)

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")

Perform the Niche Similarity Test with ecospat.niche.similarity.test()

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")

Delimiting niche categories and quantifying niche dynamics in analogue climates with ecospat.niche.dyn.index()

niche.dyn <- ecospat.niche.dyn.index (grid.clim.nat, grid.clim.inv, intersection = 0.1)

Visualizing niche categories, niche dynamics and climate analogy between ranges with ecospat.plot.niche.dyn()

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")

Plot the niche dynamics along one gradient (here temperature) with ecospat.plot.niche.dyn()

# 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")

Biotic Interactions

Species Co-occurrences Analysis with a Presence-absence matrix using the function ecospat.co_occurrences()

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)

Pairwise co-occurrence Analysis with calculation of the C-score index using the function ecospat.Cscore()

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 Preparation

Correlation Plot of Variables with ecospat.cor.plot()

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).

Calibration And Evaluation Dataset

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.

Core Niche Modelling

Model Evaluation

Presence-only Evaluation Indices- Boyce Index

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.

Accuracy of Community Prediction

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)

Spatial Predictions and Projections

ESM Ensemble of Small Models

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)

Spatial prediction of communities

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))

SESAM framework with ecospat.SESAM.prr()

prr<-ecospat.SESAM.prr(proba, sr)
head(prr)[,1:4]

Post-Modelling

Spatial Predictions of species assamblages

Co-occurrence analysis & Environmentally Constrained Null Models

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).

ecospat documentation built on Oct. 18, 2023, 1:19 a.m.