Insight into Niche overlap using ECOSPAT

[Abid Ali]

Dear Maxent Team,

My question is not about Maxent however, I would post an insight into the niche overlap under the Ecospat package. Before that, I would like to discuss following points. 

1. I have a species with three major lineages, I'd go to see niche overlap of these lineages with each other however, before that, I ran ENM and we chose variables like BIO2, BIO4, BI013, BI014 and BIO19 for modelling of each single lineage from the same extent study area called we called it M-ZONE. I would use the same set of variables for niche overlap located in a folder named "bios_sp1", but see the codes below!

According to DI COLA et al 2017, we want to test for niche conservatism, so we use the alternative 'greater', i.e. the niche overlap is more similar than random expectations. Note that in the case of invasive species, the native niche is a reference and shifted only to the invasive niche (rand.type 2). However, in our case where there are no assumptions about a reference niche (e.g.sister species), both niches can be simultaneously shifted, so we used (rand.type 1). Thus, when alternative="greater" my null hypothesis is that the niches are more similar than expected by chance.

However, as seen in step 15 of my analysis, the code encounters an error. I'm unsure how to resolve this issue.

Moving forward, I have a few questions:

1. Is it appropriate to select the same set of environmental variables for different lineages of the same species from slightly different geographical locations within the same M-ZONE used for the previous ENM?
2. How can I address the coding error in step 15, where the 'greater' option is not functioning properly? It works fine when I use the default option.
3. The given codes lack sp3, I wonder how can I include sp3 to run the same analysis for three lineages instead of only two.
4. Additionally, I'd like to plot the niche dynamics against a single gradient using ecospat.plot.niche.dyn(). Could you please assist me in rewriting this query so I can create the plot with a clear theme?
5. I would like to see if the script (see below) is fine to go with.

I appreciate your assistance and guidance in resolving these issues.

 

rm(list = ls())

library(ecospat)

library(raster)

library(dismo)

library(ecospat)

library(ade4)

library(raster)

 

 

# 1. Load the records of the species.

sp1<-read.csv("spg_a.csv", header=T, sep=",")

sp2<-read.csv("spg_d.csv", header=T, sep=",")

sp3<-read.csv("spg_c.csv", header=T, sep=",") I have the sp3 data loaded however, will only run for sp1 and sp2 since I not aware of including more than two species.lineages.

# 2. Load in a stack for each species the variables (We are using the same set of env's for both species sp1 and sp2, since they are sibling from the same M zone and these variables were used for ENM previously). Stored in directory "bios_sp1"

env_sp1<- list.files("bios_sp1/", pattern = "*.tif$", full.names = TRUE)

env_sp1<- stack(env_sp1)

# 3. From raster to points

clim_M_sp1 <- as.data.frame(rasterToPoints(env_sp1))

clim_M_sp2 <- as.data.frame(rasterToPoints(env_sp1))

clim_M_all <- as.data.frame(rasterToPoints(env_sp1))

clim_occ_sp1<- extract(env_sp1, sp1[,c(2,3)])

clim_occ_sp2<- extract(env_sp1, sp2[,c(2,3)])

xy_sp1<-sp1[,2:3]

colnames(xy_sp1)[1]="x"

colnames(xy_sp1)[2]="y"

 

xy_sp2<-sp2[,2:3]

colnames(xy_sp2)[1]="x"

colnames(xy_sp2)[2]="y"

 

clim_p_sp1<-cbind(xy_sp1, clim_occ_sp1)

clim_p_sp2<-cbind(xy_sp2, clim_occ_sp2)

 

 

# 5. dataset for the PCA that includes all sites for the study area and all species occurrences.

#In total I have seven columns, the first two coordinated and the rest five for corresponding env (BIO2, BIO4, BI013, BI014 and BIO19), Thus, [3:7] takes env data only. 

data<- rbind(clim_M_all[,3:7], clim_M_sp1[,3:7],clim_M_sp2[,3:7],clim_p_sp1[,3:7], clim_p_sp2[,3:7])

 

# 6. Weight vectors 0 and 1 for no repeated climates. 0 for occurrences and 1 for study area sites.

w<-c(rep(1,nrow(clim_M_all)) ,rep(0,nrow(clim_M_sp1)),rep(0,nrow(clim_M_sp2)),rep(0,nrow(clim_p_sp1)),rep(0,nrow(clim_p_sp2)))

 

# 7. PCA generation.... The PCA is performed with all the data from the study area. Occurrences are not used for PCA calibration (weight of 0) but their coordinates in the PCA are calculated.

 pca.cal <- dudi.pca(data, row.w = w, center = T, scale = T, scannf = F, nf = 2)

head(pca.cal)

 

# 8. rows containing clim and individual species data

row.clim.all<-1:nrow(clim_M_all)

row.clim.sp1<-(1+nrow(clim_M_all)):(nrow(clim_M_all) + nrow(clim_M_sp1))

row.clim.sp2<-(1+nrow(clim_M_all)+nrow(clim_M_sp1)):(nrow(clim_M_all)+nrow(clim_M_sp1)+nrow(clim_M_sp2))

row.sp1<-(1+nrow(clim_M_all)+nrow(clim_M_sp1)+nrow(clim_M_sp2)):(nrow(clim_M_all)+nrow(clim_M_sp1)+nrow(clim_M_sp2)+nrow(clim_p_sp1))

 

row.sp2 <-(1+nrow(clim_M_all)+nrow(clim_M_sp1)+nrow(clim_M_sp2)+nrow(clim_p_sp1)) : (nrow(clim_M_all)+nrow(clim_M_sp1)+nrow(clim_M_sp2)+nrow(clim_p_sp1)+ nrow(clim_p_sp2))

 

# 9. coordinates on each of the PCA axes for the data for the entire study area and for each of the species.

scores.clim.all    <- pca.cal$li[row.clim.all,] #coordinates for all cells

scores.clim.sp1 <- pca.cal$li[row.clim.sp1,] #coordinates for all cells in m of D.barbatus

scores.clim.sp2 <- pca.cal$li[row.clim.sp2,] #coordinates for all cells in m of D.macroura

scores.sp1 <- pca.cal$li[row.sp1,]   #coordinates for sp1

scores.sp2 <- pca.cal$li[row.sp2,]   #coordinates for sp2

 

# 10. Contribution of each variable to each component of the GWP.

ecospat.plot.contrib(contrib=pca.cal$co, eigen=pca.cal$eig)

 

# 11. Generate occurrence density surfaces in environmental space (maximum two axes). For testing, only 100 replications are considered.

z1<-ecospat.grid.clim.dyn (scores.clim.all, scores.clim.sp1, scores.sp1, R=100)

z2<-ecospat.grid.clim.dyn (scores.clim.all, scores.clim.sp2, scores.sp2, R=100)

 

# 12. Observed niche overlap metrics (D - Schoener's metric and I - Warren's metric)

D_overlap<-ecospat.niche.overlap (z1=z1, z2=z2, cor=TRUE)

D_overlap$D  # The overlap is of the proportion

D_overlap$I

# It is recommended to use at least 1000 replications for the equivalency test. Here are 10 to finish it fast

a.dyn<-ecospat.niche.equivalency.test(z1=z1 , z2=z2, rep=10)

 

# 14. CHART FOR NICHE EQUIVALENCE TEST

ecospat.plot.overlap.test(a.dyn,"D","Equivalency")

 The codes working till this phase

# We use Two sibling species and we consider that their niches are more similar, therefore, we set it "greater".

# 15.  b.dyn<-ecospat.niche.similarity.test(z1=z1 , z2=z2, alternative = "greater", rep=100, rand.type=1)

I am getting the following error.

Error in ecospat.niche.similarity.test(z1 = z1, z2 = z2, alternative = "greater",  :

  unused argument (alternative = "greater")

Any help on what this may mean is much appreciated!

Thanks