9 Spatial and temporal aspects

9.1 Mantel tests

Tests correlations between distance matrices (e.g. distance in species composition vs. environmental or geographical distance). Significance can be calculated by randomizations since the distance matrix is not a set of independent values.

library(vegan)
load("community.rda")

De=dist(soil.data$pH.KCl)  ## calculates Euclidean distance beteen soil pH values

Dv=vegdist(log1p(vas.plants),method="euclidean") # vegetation distance

plot(De,Dv) # all pairwise combinations, for n sites n*(n-1)/2

 
mantel(De,Dv) # p value from randomizations
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = De, ydis = Dv) 
## 
## Mantel statistic r: 0.2568 
##       Significance: 0.008 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.136 0.181 0.220 0.251 
## Permutation: free
## Number of permutations: 999


library(sp)  # for spatial analyses, real metric distance between lon-lat
Dxy=as.dist(spDists(as.matrix(xy[,1:2]),longlat=T))


# correllograms with space
m.c1=mantel.correlog(De,Dxy)
m.c2=mantel.correlog(Dv,Dxy)

plot(m.c1)

plot(m.c2)


# Correlation vs. distance. Empty symbols are not significant.

plot(mantel.correlog(Dv,De)) # Distance in soil pH vs vegetation


# We can use partial mantel tests to take account the third distance matrix
 
mantel.partial(De,Dv,Dxy)
## 
## Partial Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel.partial(xdis = De, ydis = Dv, zdis = Dxy) 
## 
## Mantel statistic r: 0.2578 
##       Significance: 0.005 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.127 0.172 0.200 0.228 
## Permutation: free
## Number of permutations: 999

Make a graph to explore how real distance between points differ from Euclidean distance calculated from raw latitude and longitude.

Answer 
plot(Dxy, dist(xy[,-3]))

9.2 Variance partitioning

A useful method to see how much variation in species composition is described by some other variable groups. Can also calculate how much variation different variables describe jointly.


# Making a table space with squares and products. This can describe spatial pattern
space=xy[,-3]
space$x2=space$lon**2
space$y2=space$lat**2
space$xy=space$lon*space$lat

vp.o=varpart(vegdist(vas.plants,"bray"),soil.data$pH.KCl,soil.data[,2:4],space)
showvarparts(3)

plot(vp.o,bg=2:4,Xnames=c("pH","Nutrients","Space"))

9.3 Temporal data and ordination

Using data from DroughtNet, an international collaboration. Three years 6 permanent plots are recorded. Half of them are control, another half drought treament (rain shadow excludes ca 50% or rainfall)

temp=read.table("droughtnet.txt",stringsAsFactors = F) # long data format
# Plots 1,3,5 are drought treatment
temp[1:10,]
##    plot                  taxa cover     date
## 1     1  Achillea_millefolium     3 01-07-15
## 2     1 Anthoxanthum_odoratum     1 01-07-15
## 3     1     Avenula_pubescens     5 01-07-15
## 4     1           Briza_media    20 01-07-15
## 5     1      Campanula_patula     2 01-07-15
## 6     1     Cerastium_arvense     1 01-07-15
## 7     1    Dactylis_glomerata     5 01-07-15
## 8     1     Equisetum_arvense     1 01-07-15
## 9     1     Festuca_pratensis     1 01-07-15
## 10    1         Festuca_rubra     2 01-07-15
temp$year_plot=paste0(substring(temp$date,7,8),"-",temp$plot)
temp=as.data.frame.matrix(xtabs(cover~year_plot+taxa,data=temp))

ord=metaMDS(log1p(temp),"manhattan",2)
## Run 0 stress 0.210736 
## Run 1 stress 0.2161235 
## Run 2 stress 0.2458476 
## Run 3 stress 0.2265871 
## Run 4 stress 0.2189891 
## Run 5 stress 0.2320847 
## Run 6 stress 0.210736 
## ... Procrustes: rmse 6.886996e-05  max resid 0.0002606918 
## ... Similar to previous best
## Run 7 stress 0.2746193 
## Run 8 stress 0.2106909 
## ... New best solution
## ... Procrustes: rmse 0.02954776  max resid 0.1347751 
## Run 9 stress 0.2400829 
## Run 10 stress 0.2170666 
## Run 11 stress 0.2106907 
## ... New best solution
## ... Procrustes: rmse 0.0001966733  max resid 0.0008297866 
## ... Similar to previous best
## Run 12 stress 0.2434451 
## Run 13 stress 0.2170667 
## Run 14 stress 0.2351685 
## Run 15 stress 0.2562188 
## Run 16 stress 0.210736 
## ... Procrustes: rmse 0.02953941  max resid 0.1339546 
## Run 17 stress 0.2545797 
## Run 18 stress 0.2106907 
## ... New best solution
## ... Procrustes: rmse 3.933951e-05  max resid 0.0001613128 
## ... Similar to previous best
## Run 19 stress 0.210736 
## ... Procrustes: rmse 0.02955329  max resid 0.1339921 
## Run 20 stress 0.2332111 
## *** Best solution repeated 1 times
plot(scores(ord)$sites,type="n")
points(ord,"sites",pch=16,col=c("red","blue"))

row.names(temp)
##  [1] "15-1" "15-2" "15-3" "15-4" "15-5" "15-6" "16-1" "16-2" "16-3" "16-4"
## [11] "16-5" "16-6" "17-1" "17-2" "17-3" "17-4" "17-5" "17-6" "18-1" "18-2"
## [21] "18-3" "18-4" "18-5" "18-6" "19-1" "19-2" "19-3" "19-4" "19-5" "19-6"
for(i in 1:6) {
  n=which(grepl(paste0("-",i),row.names(temp)))
  for(j in 1:(length(n)-1)) {
  arrows(scores(ord)$sites[n[j],1],scores(ord)$sites[n[j],2],
         scores(ord)$sites[n[j+1],1],scores(ord)$sites[n[j+1],2],length=0.1,angle=20,
         col=ifelse(i%%2==1,"red","blue"))}
# i%%2 is looking if i is odd or even  
  
}

9.4 Temporal and spatial variation

Dv=vegdist(log1p(temp),method="manhattan")

treatment=as.numeric(substring(row.names(temp),4,4))%%2==1
year=as.numeric(substring(row.names(temp),1,2))
space=as.numeric(substring(row.names(temp),4,4))


ev=envfit(ord,data.frame(drought=as.numeric(treatment),year,space))
ev
## 
## ***VECTORS
## 
##            NMDS1    NMDS2     r2 Pr(>r)    
## drought -0.43570  0.90009 0.0809  0.314    
## year     0.99723 -0.07438 0.4217  0.002 ** 
## space   -0.23126 -0.97289 0.5315  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
plot(ord,type="n")
points(ord,"sites",pch=16,col="red")
plot(ev,col="green")

## Let's study more with distances!

De=dist(treatment)  ## calculates treatment distance
Dt=dist(year) ## temporal distance
Dp=log1p(dist(space)) ## spatial distance (plots 1...6 in row)
plot(De,Dv) # all pairwise combinations, for n sites n*(n-1)/2

mantel(Dv,De) # p value from randomizations
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = Dv, ydis = De) 
## 
## Mantel statistic r: 0.1696 
##       Significance: 0.002 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0448 0.0581 0.0788 0.0975 
## Permutation: free
## Number of permutations: 999
mantel(Dv,Dt) # p value from randomizations
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = Dv, ydis = Dt) 
## 
## Mantel statistic r: 0.2629 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0864 0.1165 0.1451 0.1653 
## Permutation: free
## Number of permutations: 999
mantel(Dv,Dp) # p value from randomizations
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = Dv, ydis = Dp) 
## 
## Mantel statistic r: 0.3485 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0682 0.0898 0.1057 0.1227 
## Permutation: free
## Number of permutations: 999
mantel.partial(Dv,De,Dp)
## 
## Partial Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel.partial(xdis = Dv, ydis = De, zdis = Dp) 
## 
## Mantel statistic r: 0.1322 
##       Significance: 0.005 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0449 0.0611 0.0849 0.1049 
## Permutation: free
## Number of permutations: 999
mantel.partial(Dv,De,Dt)
## 
## Partial Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel.partial(xdis = Dv, ydis = De, zdis = Dt) 
## 
## Mantel statistic r: 0.1892 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0489 0.0670 0.0852 0.0945 
## Permutation: free
## Number of permutations: 999
# Variance partitioning using distances
plot(varpart(Dv,De,Dp,Dt),bg=2:4,Xnames=c("Treatment","Space","Time"))
## Warning: collinearity detected in X1: mm = 30, m = 1
## Warning: collinearity detected in X2: mm = 30, m = 5
## Warning: collinearity detected in X3: mm = 30, m = 4
## Warning: collinearity detected in cbind(X1,X2): mm = 60, m = 5
## Warning: collinearity detected in cbind(X1,X3): mm = 60, m = 5
## Warning: collinearity detected in cbind(X2,X3): mm = 60, m = 9
## Warning: collinearity detected in cbind(X1,X2,X3): mm = 90, m = 9

## Constrained ordination cca
o.temp=cca(temp~treatment+year)
anova(o.temp,by="mar")
## Permutation test for cca under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = temp ~ treatment + year)
##           Df ChiSquare      F Pr(>F)    
## treatment  1   0.05855 2.1117  0.007 ** 
## year       1   0.08305 2.9951  0.001 ***
## Residual  27   0.74866                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Conditioning out (removing) the effect of space
o.temp=cca(temp~treatment+year+Condition(space))

anova(o.temp,by="mar")
## Permutation test for cca under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
## 
## Model: cca(formula = temp ~ treatment + year + Condition(space))
##           Df ChiSquare      F Pr(>F)    
## treatment  1   0.07133 2.8329  0.001 ***
## year       1   0.08141 3.2333  0.001 ***
## Residual  26   0.65465                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(o.temp)

# Analysis of variance based on dissimilarity, treatment, year and interactions, 
adonis2(Dv~treatment*year*space)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = Dv ~ treatment * year * space)
##                      Df SumOfSqs      R2      F Pr(>F)    
## treatment             1    484.7 0.09801 4.7200  0.001 ***
## year                  1    598.5 0.12100 5.8274  0.001 ***
## space                 1    665.2 0.13449 6.4768  0.001 ***
## treatment:year        1     49.1 0.00992 0.4779  0.900    
## treatment:space       1    683.4 0.13818 6.6548  0.001 ***
## year:space            1     55.3 0.01117 0.5382  0.862    
## treatment:year:space  1    150.3 0.03040 1.4638  0.150    
## Residual             22   2259.4 0.45682                  
## Total                29   4945.8 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# same but randomizing within plot id ("space")
adonis2(Dv~treatment*year, strata=space)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks:  strata 
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = Dv ~ treatment * year, strata = space)
##                Df SumOfSqs      R2      F Pr(>F)    
## treatment       1    484.7 0.09801 3.3048  0.001 ***
## year            1    598.5 0.12100 4.0802  0.001 ***
## treatment:year  1     49.1 0.00992 0.3346  0.800    
## Residual       26   3813.6 0.77106                  
## Total          29   4945.8 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

9.5 Some GIS with R

R is a very powerful gis program.


library(maps) # global maps

map("world")


map("world",fill=T,col="grey",xlim=c(-25,40),ylim=c(35,70))


map("world",fill=T,col="grey",xlim=c(-25,40),ylim=c(35,70),border=F)
box()



## R reading mapinfo and shape (ArcInfo) files

library(rgdal)
maa=readOGR("eesti reljeef.tab")
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## OGR data source with driver: MapInfo File 
## Source: "C:\Users\Carlos Perez Carmona\Dropbox\Courses\DataAnalysisCommunityEcology\2023-2024\Data\eesti reljeef.tab", layer: "eesti reljeef"
## with 17 features
## It has 2 fields
plot(maa)
## Warning in is.projected(p4str): Package sf not available

plot(maa,col="grey",border=F)
## Warning in is.projected(p4str): Package sf not available

plot(maa,col=terrain.colors(17)[rank(maa$ID)],border=F)
## Warning in is.projected(p4str): Package sf not available
box()

# latitude and longitude need to be transformed to x and y
coord=xy[,-3]
coordinates(coord)= c("lon","lat")
proj4string(coord) = CRS("+proj=longlat")
## Warning in CRS("+proj=longlat"): sf required for evolution_status==2L
coord=spTransform(coord,CRS(proj4string(maa)))
## Warning in CRS(proj4string(maa)): sf required for evolution_status==2L
## Warning in spTransform(coord, CRS(proj4string(maa))): NULL source CRS comment,
## falling back to PROJ string
## Warning in spTransform(coord, CRS(proj4string(maa))): NULL target CRS comment,
## falling back to PROJ string
## Warning in is.projected(p4str): Package sf not available

points(coord,cex=0.5)


# Another gis package for raster data

library(raster)

est=raster::getData("alt", country="EST", mask=T)
## Warning in raster::getData("alt", country = "EST", mask = T): getData will be removed in a future version of raster
## . Please use the geodata package instead
## Warning in sp::CRS(...): sf required for evolution_status==2L

## Warning in sp::CRS(...): sf required for evolution_status==2L

plot(est,col="green")

hist(est) # histogram of raster values
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 44%
## of the raster cells were used. 100000 values used.

plot(est>20) # Estonia land >20 m above the sea level


# Baltic countries

lat=raster::getData("alt", country="LVA", mask=T)
## Warning in raster::getData("alt", country = "LVA", mask = T): getData will be removed in a future version of raster
## . Please use the geodata package instead
lit=raster::getData("alt", country="LTU", mask=T)
## Warning in raster::getData("alt", country = "LTU", mask = T): getData will be removed in a future version of raster
## . Please use the geodata package instead

plot(merge(est,lat,lit),col=terrain.colors(286)) # merging all three countries

points(xy[,-3],cex=0.5) ## this raster map uses lat and lon



elev=raster::extract(est,xy[,-3]) # extracting values for our coordinates

elev
##  [1]  66  65  77  78  51  55  91  33  56  42  40  47  50 106 100  42  99 109  46
## [20]  42  44 114  96  53  79  87 110  80  54 131

Raster package has also current climate data and future climate predictions. It is a powerful GIS software.