Gastropod diversification and community structuring processes in ancient Lake Ohrid: A metacommunity speciation perspective
This document showes how to reproduce the analyses of Hauffe et al. The file Hauffe.R does the same but is not resulting in a html.
Warning: Using one CPU core only, the whole distance-based STEPCAM will last months (almost a year)
Therefore, this script can perform the following tasks:
1.Run the whole analysis with the settings used in Hauffe et al. (months)
Run the whole analysis with inappropriate settings. Distance-based STEPCAM will not converge (about one day)
Just load the results of the spatial clustering and distance-based STEPCAM and will plot the figures and run the Bayesian GLM of Hauffe et al. (less than an hour)
Run <- 3 # Third option of the list above
# The geoclust function is only working on Linux and possibly Mac
# If this script is run on Windows, a already performed clusteringshould be loaded.
switch(Sys.info()[['sysname']],
Windows = {OS <- 0},
Linux = {OS <- 1},
Darwin = {OS <- 1})
if(Run == 1){
numParticles <- 500
stopRate <- 0.0001
SD <- 500
}
if(Run == 2){
numParticles <- 50
stopRate <- 0.1
SD <- 25
}
if(Run == 3){
OS <- 0
}library(STEPCAM)## Loading required package: vcd
## Loading required package: grid
## Loading required package: FD
## Loading required package: ade4
## Loading required package: ape
## Loading required package: geometry
## Loading required package: magic
## Loading required package: abind
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.2-1
##
## Attaching package: 'vegan'
##
## The following object is masked from 'package:ade4':
##
## ccalibrary(MCMCglmm)## Loading required package: Matrix
## Loading required package: codalibrary(vegan)
library(rgdal)## Loading required package: sp
## rgdal: version: 0.8-16, (SVN revision 498)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.10.1, released 2013/08/26
## Path to GDAL shared files: C:/Users/AGW-Torsten/Documents/R/win-library/3.0/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: C:/Users/AGW-Torsten/Documents/R/win-library/3.0/rgdal/projlibrary(raster)##
## Attaching package: 'raster'
##
## The following object is masked from 'package:magic':
##
## shift
##
## The following objects are masked from 'package:ape':
##
## rotate, zoom
##
## The following object is masked from 'package:vcd':
##
## mosaiclibrary(gdistance)## Loading required package: igraph
##
## Attaching package: 'igraph'
##
## The following object is masked from 'package:MCMCglmm':
##
## sir
##
## The following object is masked from 'package:ape':
##
## edgesif(OS != 0){
library(geoclust) # not on CRAN http://lbbe.univ-lyon1.fr/Download-5012.html?lang=fr
}
library(compositions)## Loading required package: tensorA
##
## Attaching package: 'tensorA'
##
## The following object is masked from 'package:Matrix':
##
## norm
##
## The following object is masked from 'package:base':
##
## norm
##
## Loading required package: robustbase
## Loading required package: energy
## Loading required package: bayesm
## Welcome to compositions, a package for compositional data analysis.
## Find an intro with "? compositions"
##
##
## Attaching package: 'compositions'
##
## The following object is masked from 'package:gdistance':
##
## normalize
##
## The following object is masked from 'package:raster':
##
## scale
##
## The following object is masked from 'package:magic':
##
## aplus
##
## The following object is masked from 'package:ape':
##
## balance
##
## The following object is masked from 'package:vcd':
##
## Kappa
##
## The following objects are masked from 'package:stats':
##
## cor, cov, dist, var
##
## The following objects are masked from 'package:base':
##
## %*%, scale, scale.defaultlibrary(gtools)##
## Attaching package: 'gtools'
##
## The following object is masked from 'package:bayesm':
##
## rdirichlet
##
## The following object is masked from 'package:permute':
##
## permutelibrary(knitr) # Not for the analysis but for nice html reportLoad the script for distance-based stepwise community assembly
This script uses a lot of code from the packages FD and STEPCAM. Please cite them correctly! citation(“FD”) citation(“STEPCAM”)
The code here (Hauffe et al., 2015) is consequently distributed under the same license: GPL-2 See https://cran.r-project.org/web/licenses/GPL-2
source("dbSTEPCAM.R")
source("auxiliaryFunctions.R")Spatially constrained clustering of ecological networks
OhridComm <- read.csv("OhridComm.csv", sep = ",", header = TRUE, row.names = 1)
TaxDist <- read.csv("TaxDistance.csv", sep = ",", header = TRUE, row.names = 1)
Localities <- read.csv("Localities.csv", sep = ",", header = TRUE, row.names = 1)
# This step is only working on Linux and possibly Mac (not tested)
if(Run == 1 & OS == 1){ # Slow version
# Ecological network
OhridCommDist <- as.matrix(vegdist(OhridComm))
# Spatial (structural) network
# 1. Depth differences:
OhridDepthDist <- as.matrix(dist(Localities$Depth))
MantelCorlg <- mantel.correlog(as.dist(OhridCommDist), as.dist(OhridDepthDist), nperm = 9999)
MantelCorlg
# 2. cost layer:
OhridShp <- readOGR(dsn = getwd(), layer="Ohrid")
OhridRas <- raster(xmn = bbox(OhridShp)[1,1], xmx = bbox(OhridShp)[1,2],
ymn = bbox(OhridShp)[2,1],ymx = bbox(OhridShp)[2,2],
crs = OhridShp@proj4string)
res(OhridRas) <- 10 # Grain size
# Lasts ~ 1h
OhridRas <- rasterize(OhridShp, OhridRas)
OhridTrans <- transition(OhridRas, mean, 16)
CostDist <- costDistance(OhridTrans, as.matrix(Localities[,c("Eastings","Northings")]))
CostDist <- as.matrix(CostDist)
# Normalize the [0,1] community dissimilarities
OhridCommDistAsin <- asin(as.matrix(OhridCommDist))
#boxplot(c(OhridCommDistAsin))
# Take the 8.2 m from the mantel correlogram
# Create structural network (as in Miele et al., 2014)
# by removing vertical (bathymetrical) links
StructNet <- CostDist * ifelse(OhridDepthDist > 8.2, 0, 1)
# Actual spatial clustering
# Explore different settings for the number of connected spatial neighbors
# in the structural network (i.e., modify the connectivity of the x and y dimension)
# Every locality is connected with at minimum 15 localities
K <- seq(15, 265, by = 10)
GC <- list() # Geographic clustering
for(k in 1:length(K)){ # ~1 day on my Atom CPU
print(K[k])
StructNetTmp <- apply(StructNet, 2, function(x) pruneK(x, K = K[k]))
ListWStructNet <- mat2listw(as.matrix(StructNetTmp), style = "M")
GC[[k]] <- geoclust(m = OhridCommDistAsin, lw = ListWStructNet,
k = 1:10, clump = 1.570796)
# Layers from Hauffe et al., 2011
# SL: Surface Layer
# IL: Intermediate Layer
# DL: Deep Layer
# FS: Feeder Springs
print(table(GC[[k]]$groups, Localities$Layer))
# Group number may differ between runs, though same settings are used!
}
}
# Fast version or if geoclust cannot be used (no Linux/Mac)
if(Run %in% c(2,3) | OS != 0 ){
GC <- readRDS("SpatiallyConstrainedClustering.rds")
}
# Check robustness of clustering
################################
# Spatially constrained clustering uses arbitrary numbers ('names') for the clusters
# and repeated runs with the same settings may also use different numbers.
# Here, cluster-numbers are replaced by the name of the eco-zones used in Hauffe et al.
EcoZones <- lapply(GC, function(x) nameGeoCluster(Loc = Localities, Groups = x$groups))
# Order all according to a referenz cluster (here 65 neighbors)
G <- 6
Localities <- data.frame(Localities, GeoClus = EcoZones[[G]])
RefClus <- order(EcoZones[[G]])
EcoZones <- lapply(EcoZones, function(x) x[RefClus])
# Color coded cluster assignment of all 264 localities
par(las = 2)
palette(c("#b2182b", "#ef8a62", "#fddbc7", "#f7f7f7", "#d1e5f0", "#67a9cf", "#2166ac"))
plot(1:length(EcoZones), 1:length(EcoZones), type = "n", ylim = c(1, length(RefClus)),
xlab = "Neighbors", ylab = "Locality", yaxt = "n", xaxt = "n")
axis(side = 1, at = 1:length(EcoZones), labels = seq(15, 265, by = 10))
for(i in 1:length(EcoZones)){
points(rep(i, length(RefClus)), 1:length(RefClus), pch = 15, col = as.factor(EcoZones[[i]]))
} 
# Figure 1 of Hauffe et al.
###########################
Ohrid0mShp <- readOGR(dsn = getwd(), layer="0m")## OGR data source with driver: ESRI Shapefile
## Source: "Z:/Staff/Torsten_Hauffe/Biogeosciences2015/HTML", layer: "0m"
## with 1 features and 2 fields
## Feature type: wkbPolygon with 2 dimensionslayout(matrix(1:3, nc = 3, nr = 1))
par(las = 2, mar = c(8.1, 4.5, 0.1, 0.1))
boxplot(Localities[Localities$GeoClus == "1 SE upper littoral zone 1", "Depth"]+0.1, #
Localities[Localities$GeoClus == "2 SE upper littoral zone 2", "Depth"]+0.1,
Localities[Localities$GeoClus == "3 SE upper littoral zone 3", "Depth"]+0.1,
Localities[Localities$GeoClus == "4 Non-SE upper littoral", "Depth"]+0.1,
Localities[Localities$GeoClus == "5 lower littoral", "Depth"]+0.1,
Localities[Localities$GeoClus == "6 upper sublittoral", "Depth"]+0.1,
Localities[Localities$GeoClus == "7 lower sublittoral", "Depth"]+0.1, # 4
col = c("#b2182b", "#ef8a62", "#fddbc7", "#f7f7f7", "#d1e5f0", "#67a9cf", "#2166ac"),
names = c("SE upper littoral\nzone 1","SE upper littoral\nzone 2","SE upper littoral\nzone 3",
"Non-SE upper\nlittoral","lower\nlittoral","upper\nsublittoral", "lower\nsublittoral"),
log = "y", ylab = "Depth (m)")
mtext("A", side = 2, line = 3.1, at = 60, cex = 1.4)
par(las = 2, mar = c(0.1, 0.1, 0.1, 0.1))
plot(Ohrid0mShp, col = "grey90")
# Different order to minimize issues with overplotting
points(Localities[Localities$GeoClus == "4 Non-SE upper littoral", 4:5],
pch = 21, cex = 1.5, bg = "#f7f7f7")
points(Localities[Localities$GeoClus == "3 SE upper littoral zone 3", 4:5],
pch = 21, cex = 1.5, bg = "#fddbc7")
points(Localities[Localities$GeoClus == "2 SE upper littoral zone 2", 4:5],
pch = 21, cex = 1.5, bg = "#ef8a62")
points(Localities[Localities$GeoClus == "1 SE upper littoral zone 1", 4:5],
pch = 21, cex = 1.5, bg = "#b2182b")
text(x = par("usr")[1] + par("usr")[1]/500, y = par("usr")[4] - par("usr")[4]/7000,
labels = "B", cex = 2)
plot(Ohrid0mShp, col = "grey90")
points(Localities[Localities$GeoClus == "5 lower littoral", 4:5],
pch = 21, cex = 1.5, bg = "#d1e5f0")
points(Localities[Localities$GeoClus == "6 upper sublittoral", 4:5],
pch = 21, cex = 1.5, bg = "#67a9cf")
points(Localities[Localities$GeoClus == "7 lower sublittoral", 4:5],
pch = 21, cex = 1.5, bg = "#2166ac")
text(x = par("usr")[1] + par("usr")[1]/500, y = par("usr")[4] - par("usr")[4]/7000,
labels = "C", cex = 2)
Distance-based STEPCAM
setwd("Z:/Staff/Torsten_Hauffe/Biogeosciences2015/Supplement")
if(Run %in% c(1,2)){ # Slow version
for(i in 1:nrow(OhridComm)){
O <- dbSTEPCAM_ABC(OhridComm, as.dist(TaxDist),
numParticles = numParticles,
plot_number = i, stopRate = stopRate,
SD = SD, ContinuePrevious = FALSE)
# Save output:
f <- list.files(pattern = "particles_t=")
f <- mixedsort(f)
d <- read.table(f[length(f)-1], header = FALSE) # Read last run
write.table(d, paste(getwd(),"/", i,"_t=",
length(f)-1,".txt", sep = ""),
row.names = FALSE, col.names = FALSE)
# Clean up:
unlink(paste(getwd(), f, sep = "/"))
}
}
# Posterior mean of STEPCAM
###########################
if(Run %in% c(1,2)){
OhridRuns <- list.files(getwd(), pattern = "_t=")
OhridRuns <- mixedsort(OhridRuns)
# Load the posterior distribution of process-importance for all 264 localities
StepcamRes <- sapply(OhridRuns, function(x)
read.table(paste(getwd(), "/", x, sep = ""), sep = "", header = FALSE),
simplify = FALSE)
# Process-importance in the interval [0,1]
ProImp <- lapply(StepcamRes, function(x) rowMeans(apply(x, 1, function(y) y[1:3]/sum(y[1:3]))) )
ProImp <- matrix(unlist(ProImp), nc = 3, byrow = TRUE)
colnames(ProImp) <- c("Dispersal", "Environment", "Interaction")
# Result:
Res <- data.frame(Localities, ProImp)
}
# STEPCAM results of Hauffe et al.
if(Run == 3){
Res <- read.csv("StepcamResultHauffe.csv", sep = ",", header = TRUE, row.names = 1)
}
# Figure 2 of Hauffe et al.
layout(matrix(1:3, nc = 3, nr = 1))
par(las = 2, mar = c(8.1, 4.1, 0.1, 0.1))
boxplot(Dispersal ~ GeoClus, data = Res,
#col = c("magenta", "red4", "green4", "cyan", "cornflowerblue", "blue", "blue4"),
col = c("#b2182b", "#ef8a62", "#fddbc7", "#f7f7f7", "#d1e5f0", "#67a9cf", "#2166ac"),
names = c("SE upper littoral\nzone 1","SE upper littoral\nzone 2","SE upper littoral\nzone 3",
"Non-SE upper\nlittoral","lower\nlittoral","upper\nsublittoral", "lower\nsublittoral"),
ylab = "Dispersal limitation importance", ylim = c(0,1))
mtext("A", side = 2, line = 2.5, at = 1, cex = 1.4)
boxplot(Environment ~ GeoClus, data = Res,
col = c("#b2182b", "#ef8a62", "#fddbc7", "#f7f7f7", "#d1e5f0", "#67a9cf", "#2166ac"),
names = c("SE upper littoral\nzone 1","SE upper littoral\nzone 2","SE upper littoral\nzone 3",
"Non-SE upper\nlittoral","lower\nlittoral","upper\nsublittoral", "lower\nsublittoral"),
ylab = "Environmental filtering importance", ylim = c(0,1))
mtext("B", side = 2, line = 2.5, at = 1, cex = 1.4)
boxplot(Interaction ~ GeoClus, data = Res,
col = c("#b2182b", "#ef8a62", "#fddbc7", "#f7f7f7", "#d1e5f0", "#67a9cf", "#2166ac"),
names = c("SE upper littoral\nzone 1","SE upper littoral\nzone 2","SE upper littoral\nzone 3",
"Non-SE upper\nlittoral","lower\nlittoral","upper\nsublittoral", "lower\nsublittoral"),
ylab = "Species interaction importance", ylim = c(0,1))
mtext("C", side = 2, line = 2.5, at = 1, cex = 1.4)
# Mean process-importance and bootstrapped confidence interval of the mean
colMeans(Res[, c("Dispersal", "Environment", "Interaction")])## Dispersal Environment Interaction
## 0.7865417 0.1220283 0.0914300MeansBoot <- matrix(NA, nc = 3, nr = 10000)
for(i in 1:nrow(MeansBoot)){
bsample <- sample(1:nrow(Res), nrow(Res), replace = TRUE )
ResBoot <- Res[bsample,]
MeansBoot[i,] <- colMeans(ResBoot[, c("Dispersal", "Environment", "Interaction")])
}
apply(MeansBoot, 2, function(x) quantile(x, c(0.025, 0.975)))## [,1] [,2] [,3]
## 2.5% 0.7552137 0.09983307 0.0785725
## 97.5% 0.8167957 0.14592331 0.1047689Bayesian generalized linear model
Test whether process importance changes with limnological characteristics of Lake Ohrid Two competing hypotheses: either lake depth or eco-zones are better correlated with processes. Operational criterion: deviance information criterion (Lower DIC is better).
# additive planar transformation
Res <- data.frame(Res, apt(Res[,9:7]))
# Correlation of process-importance with depth (Do processes change with lake depth?)
McmcDepth1 <- MCMCglmm(cbind(Interaction.1, Environment.1) ~ # Multivariate dependent variables
trait:Depth + trait:SpeciesRichness, # Predictors/independend
data = Res, rcov=~us(trait):units,
family = rep("gaussian", 2),
nitt = 30000, thin = 20, burnin = 5000,
verbose = FALSE)
summary(McmcDepth1)$DIC## [1] -641.1425kable(summary(McmcDepth1)$solutions, digits = 3) | post.mean | l-95% CI | u-95% CI | eff.samp | pMCMC | |
|---|---|---|---|---|---|
| (Intercept) | 0.132 | 0.099 | 0.163 | 1099.757 | 0.001 |
| traitInteraction.1:Depth | 0.002 | 0.001 | 0.003 | 1099.343 | 0.001 |
| traitEnvironment.1:Depth | 0.005 | 0.004 | 0.007 | 1140.510 | 0.001 |
| traitInteraction.1:SpeciesRichness | -0.007 | -0.011 | -0.004 | 950.650 | 0.001 |
| traitEnvironment.1:SpeciesRichness | -0.011 | -0.015 | -0.006 | 1287.682 | 0.001 |
# Check temporal autocorrelation of MCMC samples
autocorr.plot( McmcDepth1$Sol ) # QQ-plot (as in regular glm does not exist for BGLMs)
effectiveSize( McmcDepth1$Sol ) ## (Intercept) traitInteraction.1:Depth
## 1099.7567 1099.3428
## traitEnvironment.1:Depth traitInteraction.1:SpeciesRichness
## 1140.5098 950.6496
## traitEnvironment.1:SpeciesRichness
## 1287.6823# 2-way interaction between depth and species richness
# In the analyses of Hauffe et al., the interaction is significant
# and it does improve the model fit (DIC)
McmcDepth2 <- MCMCglmm(cbind(Interaction.1, Environment.1) ~
trait:Depth + trait:SpeciesRichness + trait:Depth:SpeciesRichness,
data = Res, rcov=~us(trait):units,
family = rep("gaussian", 2),
nitt = 60000, thin = 40, burnin = 10000,
verbose = FALSE)
summary(McmcDepth2)$DIC## [1] -662.7343kable(summary(McmcDepth2)$solutions, digits = 3)| post.mean | l-95% CI | u-95% CI | eff.samp | pMCMC | |
|---|---|---|---|---|---|
| (Intercept) | 0.135 | 0.098 | 0.170 | 1042.084 | 0.001 |
| traitInteraction.1:Depth | 0.000 | -0.002 | 0.002 | 1250.000 | 0.824 |
| traitEnvironment.1:Depth | 0.010 | 0.007 | 0.013 | 1250.000 | 0.001 |
| traitInteraction.1:SpeciesRichness | -0.008 | -0.013 | -0.004 | 1248.991 | 0.002 |
| traitEnvironment.1:SpeciesRichness | -0.009 | -0.014 | -0.004 | 1250.000 | 0.001 |
| traitInteraction.1:Depth:SpeciesRichness | 0.000 | 0.000 | 0.001 | 1250.000 | 0.093 |
| traitEnvironment.1:Depth:SpeciesRichness | -0.001 | -0.001 | 0.000 | 1250.000 | 0.002 |
autocorr.plot( McmcDepth2$Sol ) 
effectiveSize( McmcDepth2$Sol )## (Intercept)
## 1042.084
## traitInteraction.1:Depth
## 1250.000
## traitEnvironment.1:Depth
## 1250.000
## traitInteraction.1:SpeciesRichness
## 1248.991
## traitEnvironment.1:SpeciesRichness
## 1250.000
## traitInteraction.1:Depth:SpeciesRichness
## 1250.000
## traitEnvironment.1:Depth:SpeciesRichness
## 1250.000# Check spatial autocorrelation
# Prediction
# It is not possible to use predict with newdata (as in regular base:::glm)
# That is why I cannot plot a single regression line showing how process importance
# changes with depth or species richness
PredictMcmcDepth <- aptInv(matrix(predict(McmcDepth2), nc = 2, byrow = FALSE))[,3:1]
ResidualsMcmcDepth <- Res[,c("Dispersal", "Environment", "Interaction")] - PredictMcmcDepth
# No way to test three-dimensional spatial autocorrelation
# Spatial autocorrelation in x and y axes:
MantelCorlgDepth <- mantel.correlog(dist(ResidualsMcmcDepth),
dist(Localities[,c("Eastings","Northings")]), nperm = 999)
# Only first lag distance with significant
# but very little autocorrelation (rM = 0.02)
plot(MantelCorlgDepth)
# Spatial autocorrelation in z axis:
MantelCorlgDepth <- mantel.correlog(dist(ResidualsMcmcDepth),
dist(Localities[,"Depth"]), nperm = 999)
# Only second lag distance with significant autocorrelation (rM = 0.06)
plot(MantelCorlgDepth)
# Check modeled trend in process-importance
layout(matrix(1:6, nr = 2, nc = 3, byrow = TRUE))
plot(PredictMcmcDepth[,1] ~ Res$Depth, xlab = "Depth (m)",
ylab = "Dispersal limitation importance")
plot(PredictMcmcDepth[,2] ~ Res$Depth, xlab = "Depth (m)",
ylab = "Environmental filtering importance")
plot(PredictMcmcDepth[,3] ~ Res$Depth, xlab = "Depth (m)",
ylab = "Species interaction importance")
plot(PredictMcmcDepth[,1] ~ Res$SpeciesRichness, xlab = "Species Richness",
ylab = "Dispersal limitation importance")
plot(PredictMcmcDepth[,2] ~ Res$SpeciesRichness, xlab = "Species Richness",
ylab = "Environmental filtering importance")
plot(PredictMcmcDepth[,3] ~ Res$SpeciesRichness, xlab = "Species Richness",
ylab = "Species interaction importance")
# Correlation of process-importance with eco-zones (Do processes change with eco-zones?)
McmcEcoZone1 <- MCMCglmm(cbind(Interaction.1, Environment.1) ~
trait:GeoClus + trait:SpeciesRichness,
data = Res, rcov=~us(trait):units,
family = rep("gaussian", 2),
nitt = 30000, thin = 20, burnin = 5000,
verbose = FALSE)## Warning in MCMCglmm(cbind(Interaction.1, Environment.1) ~ trait:GeoClus +
## : some fixed effects are not estimable and have been removed. Use
## singular.ok=TRUE to sample these effects, but use an informative prior!summary(McmcEcoZone1)$DIC## [1] -723.8611kable(summary(McmcEcoZone1)$solutions, digits = 3)| post.mean | l-95% CI | u-95% CI | eff.samp | pMCMC | |
|---|---|---|---|---|---|
| (Intercept) | 0.449 | 0.385 | 0.515 | 1207.781 | 0.001 |
| traitInteraction.1:GeoClus1 SE upper littoral zone 1 | -0.284 | -0.369 | -0.203 | 1116.142 | 0.001 |
| traitEnvironment.1:GeoClus1 SE upper littoral zone 1 | 0.025 | -0.070 | 0.126 | 923.321 | 0.635 |
| traitInteraction.1:GeoClus2 SE upper littoral zone 2 | -0.376 | -0.450 | -0.291 | 1250.000 | 0.001 |
| traitEnvironment.1:GeoClus2 SE upper littoral zone 2 | -0.302 | -0.381 | -0.213 | 1250.000 | 0.001 |
| traitInteraction.1:GeoClus3 SE upper littoral zone 3 | -0.350 | -0.421 | -0.275 | 1250.000 | 0.001 |
| traitEnvironment.1:GeoClus3 SE upper littoral zone 3 | -0.224 | -0.297 | -0.151 | 1250.000 | 0.001 |
| traitInteraction.1:GeoClus4 Non-SE upper littoral | -0.388 | -0.452 | -0.314 | 1042.029 | 0.001 |
| traitEnvironment.1:GeoClus4 Non-SE upper littoral | -0.274 | -0.333 | -0.206 | 1196.814 | 0.001 |
| traitInteraction.1:GeoClus5 lower littoral | -0.309 | -0.375 | -0.233 | 1250.000 | 0.001 |
| traitEnvironment.1:GeoClus5 lower littoral | -0.206 | -0.274 | -0.135 | 1250.000 | 0.001 |
| traitInteraction.1:GeoClus6 upper sublittoral | -0.328 | -0.400 | -0.249 | 1250.000 | 0.001 |
| traitEnvironment.1:GeoClus6 upper sublittoral | -0.183 | -0.255 | -0.113 | 1250.000 | 0.001 |
| traitInteraction.1:GeoClus7 lower sublittoral | -0.270 | -0.341 | -0.212 | 955.481 | 0.001 |
| traitInteraction.1:SpeciesRichness | -0.002 | -0.007 | 0.002 | 1250.000 | 0.306 |
| traitEnvironment.1:SpeciesRichness | -0.019 | -0.025 | -0.012 | 1250.000 | 0.001 |
layout(matrix(1:16, 4, 4))
autocorr.plot( McmcEcoZone1$Sol, auto.layout = FALSE ) 
effectiveSize( McmcEcoZone1$Sol )## (Intercept)
## 1207.7815
## traitInteraction.1:GeoClus1 SE upper littoral zone 1
## 1116.1419
## traitEnvironment.1:GeoClus1 SE upper littoral zone 1
## 923.3209
## traitInteraction.1:GeoClus2 SE upper littoral zone 2
## 1250.0000
## traitEnvironment.1:GeoClus2 SE upper littoral zone 2
## 1250.0000
## traitInteraction.1:GeoClus3 SE upper littoral zone 3
## 1250.0000
## traitEnvironment.1:GeoClus3 SE upper littoral zone 3
## 1250.0000
## traitInteraction.1:GeoClus4 Non-SE upper littoral
## 1042.0290
## traitEnvironment.1:GeoClus4 Non-SE upper littoral
## 1196.8139
## traitInteraction.1:GeoClus5 lower littoral
## 1250.0000
## traitEnvironment.1:GeoClus5 lower littoral
## 1250.0000
## traitInteraction.1:GeoClus6 upper sublittoral
## 1250.0000
## traitEnvironment.1:GeoClus6 upper sublittoral
## 1250.0000
## traitInteraction.1:GeoClus7 lower sublittoral
## 955.4809
## traitInteraction.1:SpeciesRichness
## 1250.0000
## traitEnvironment.1:SpeciesRichness
## 1250.0000# 2-way interaction between eco-zones and species richness
# In the analyses of Hauffe et al., the interaction is significant
# and it does improve the model fit (DIC)
McmcEcoZone2 <- MCMCglmm(cbind(Interaction.1, Environment.1) ~
trait:GeoClus + trait:SpeciesRichness + trait:GeoClus:SpeciesRichness,
data = Res, rcov=~us(trait):units,
family = rep("gaussian", 2),
nitt = 60000, thin = 40, burnin = 10000,
verbose = FALSE)## Warning in MCMCglmm(cbind(Interaction.1, Environment.1) ~ trait:GeoClus +
## : some fixed effects are not estimable and have been removed. Use
## singular.ok=TRUE to sample these effects, but use an informative prior!summary(McmcEcoZone2)$DIC## [1] -740.7374kable(summary(McmcEcoZone2)$solutions, digits = 3)| post.mean | l-95% CI | u-95% CI | eff.samp | pMCMC | |
|---|---|---|---|---|---|
| (Intercept) | 0.626 | 0.477 | 0.758 | 1123.800 | 0.001 |
| traitInteraction.1:GeoClus1 SE upper littoral zone 1 | -0.657 | -0.863 | -0.460 | 1320.041 | 0.001 |
| traitEnvironment.1:GeoClus1 SE upper littoral zone 1 | -0.124 | -0.369 | 0.153 | 1250.000 | 0.342 |
| traitInteraction.1:GeoClus2 SE upper littoral zone 2 | -0.497 | -0.689 | -0.278 | 1250.000 | 0.001 |
| traitEnvironment.1:GeoClus2 SE upper littoral zone 2 | -0.442 | -0.699 | -0.175 | 1250.000 | 0.001 |
| traitInteraction.1:GeoClus3 SE upper littoral zone 3 | -0.506 | -0.660 | -0.333 | 1291.413 | 0.001 |
| traitEnvironment.1:GeoClus3 SE upper littoral zone 3 | -0.379 | -0.559 | -0.179 | 1231.418 | 0.001 |
| traitInteraction.1:GeoClus4 Non-SE upper littoral | -0.584 | -0.745 | -0.437 | 1250.000 | 0.001 |
| traitEnvironment.1:GeoClus4 Non-SE upper littoral | -0.520 | -0.679 | -0.341 | 1120.045 | 0.001 |
| traitInteraction.1:GeoClus5 lower littoral | -0.430 | -0.608 | -0.273 | 1250.000 | 0.001 |
| traitEnvironment.1:GeoClus5 lower littoral | -0.405 | -0.600 | -0.214 | 1250.000 | 0.001 |
| traitInteraction.1:GeoClus6 upper sublittoral | -0.446 | -0.600 | -0.275 | 1232.963 | 0.001 |
| traitEnvironment.1:GeoClus6 upper sublittoral | -0.360 | -0.544 | -0.168 | 1096.971 | 0.001 |
| traitInteraction.1:GeoClus7 lower sublittoral | -0.570 | -0.715 | -0.433 | 1129.481 | 0.001 |
| traitInteraction.1:SpeciesRichness | 0.043 | 0.010 | 0.077 | 1250.000 | 0.018 |
| traitEnvironment.1:SpeciesRichness | -0.025 | -0.074 | 0.021 | 1250.000 | 0.322 |
| traitInteraction.1:GeoClus2 SE upper littoral zone 2:SpeciesRichness | -0.058 | -0.104 | -0.012 | 1250.000 | 0.022 |
| traitEnvironment.1:GeoClus2 SE upper littoral zone 2:SpeciesRichness | -0.002 | -0.075 | 0.064 | 1250.000 | 0.960 |
| traitInteraction.1:GeoClus3 SE upper littoral zone 3:SpeciesRichness | -0.048 | -0.082 | -0.013 | 1250.000 | 0.014 |
| traitEnvironment.1:GeoClus3 SE upper littoral zone 3:SpeciesRichness | 0.004 | -0.045 | 0.057 | 1250.000 | 0.869 |
| traitInteraction.1:GeoClus4 Non-SE upper littoral:SpeciesRichness | -0.043 | -0.077 | -0.009 | 1250.000 | 0.026 |
| traitEnvironment.1:GeoClus4 Non-SE upper littoral:SpeciesRichness | 0.015 | -0.034 | 0.067 | 1250.000 | 0.546 |
| traitInteraction.1:GeoClus5 lower littoral:SpeciesRichness | -0.054 | -0.090 | -0.021 | 1250.000 | 0.005 |
| traitEnvironment.1:GeoClus5 lower littoral:SpeciesRichness | 0.010 | -0.043 | 0.062 | 1250.000 | 0.707 |
| traitInteraction.1:GeoClus6 upper sublittoral:SpeciesRichness | -0.051 | -0.082 | -0.014 | 1250.000 | 0.005 |
| traitEnvironment.1:GeoClus6 upper sublittoral:SpeciesRichness | 0.007 | -0.044 | 0.056 | 1250.000 | 0.794 |
| traitInteraction.1:GeoClus7 lower sublittoral:SpeciesRichness | -0.023 | -0.062 | 0.013 | 1149.098 | 0.222 |
| traitEnvironment.1:GeoClus7 lower sublittoral:SpeciesRichness | -0.026 | -0.080 | 0.026 | 1250.000 | 0.342 |
layout(matrix(1:30, 5, 6))
autocorr.plot( McmcEcoZone2$Sol, auto.layout = FALSE )
effectiveSize( McmcEcoZone2$Sol )## (Intercept)
## 1123.800
## traitInteraction.1:GeoClus1 SE upper littoral zone 1
## 1320.041
## traitEnvironment.1:GeoClus1 SE upper littoral zone 1
## 1250.000
## traitInteraction.1:GeoClus2 SE upper littoral zone 2
## 1250.000
## traitEnvironment.1:GeoClus2 SE upper littoral zone 2
## 1250.000
## traitInteraction.1:GeoClus3 SE upper littoral zone 3
## 1291.413
## traitEnvironment.1:GeoClus3 SE upper littoral zone 3
## 1231.418
## traitInteraction.1:GeoClus4 Non-SE upper littoral
## 1250.000
## traitEnvironment.1:GeoClus4 Non-SE upper littoral
## 1120.045
## traitInteraction.1:GeoClus5 lower littoral
## 1250.000
## traitEnvironment.1:GeoClus5 lower littoral
## 1250.000
## traitInteraction.1:GeoClus6 upper sublittoral
## 1232.963
## traitEnvironment.1:GeoClus6 upper sublittoral
## 1096.971
## traitInteraction.1:GeoClus7 lower sublittoral
## 1129.481
## traitInteraction.1:SpeciesRichness
## 1250.000
## traitEnvironment.1:SpeciesRichness
## 1250.000
## traitInteraction.1:GeoClus2 SE upper littoral zone 2:SpeciesRichness
## 1250.000
## traitEnvironment.1:GeoClus2 SE upper littoral zone 2:SpeciesRichness
## 1250.000
## traitInteraction.1:GeoClus3 SE upper littoral zone 3:SpeciesRichness
## 1250.000
## traitEnvironment.1:GeoClus3 SE upper littoral zone 3:SpeciesRichness
## 1250.000
## traitInteraction.1:GeoClus4 Non-SE upper littoral:SpeciesRichness
## 1250.000
## traitEnvironment.1:GeoClus4 Non-SE upper littoral:SpeciesRichness
## 1250.000
## traitInteraction.1:GeoClus5 lower littoral:SpeciesRichness
## 1250.000
## traitEnvironment.1:GeoClus5 lower littoral:SpeciesRichness
## 1250.000
## traitInteraction.1:GeoClus6 upper sublittoral:SpeciesRichness
## 1250.000
## traitEnvironment.1:GeoClus6 upper sublittoral:SpeciesRichness
## 1250.000
## traitInteraction.1:GeoClus7 lower sublittoral:SpeciesRichness
## 1149.098
## traitEnvironment.1:GeoClus7 lower sublittoral:SpeciesRichness
## 1250.000# Check spatial autocorrelation
PredictMcmcEcoZone <- aptInv(matrix(predict(McmcEcoZone2), nc = 2, byrow = FALSE))[,3:1]
ResidualsMcmcEcoZone <- Res[,c("Dispersal", "Environment", "Interaction")] - PredictMcmcEcoZone
# Spatial autocorrelation in x and y axes:
MantelCorlgEcoZone <- mantel.correlog(dist(ResidualsMcmcEcoZone),
dist(Localities[,c("Eastings","Northings")]), nperm = 999)
layout(matrix(1, 1, 1))
# No significant autocorrelation
plot(MantelCorlgEcoZone)
# Spatial autocorrelation in z axis:
MantelCorlgEcoZone <- mantel.correlog(dist(ResidualsMcmcEcoZone),
dist(Localities[,"Depth"]), nperm = 999)
# Only second and 6th lag distance with significant
# but very low autocorrelation (rM = 0.06)
plot(MantelCorlgEcoZone)
# Check modeled trend in process-importance
layout(matrix(1:6, nr = 2, nc = 3, byrow = TRUE))
par(las = 2, mar = c(8.1, 4.1, 0.1, 0.1))
boxplot(PredictMcmcEcoZone[,1] ~ Res$GeoClus,
ylab = "Dispersal limitation importance", ylim = c(0,1),
col = c("#b2182b", "#ef8a62", "#fddbc7", "#f7f7f7", "#d1e5f0", "#67a9cf", "#2166ac"))
boxplot(PredictMcmcEcoZone[,2] ~ Res$GeoClus,
ylab = "Environmental filtering importance", ylim = c(0,1),
col = c("#b2182b", "#ef8a62", "#fddbc7", "#f7f7f7", "#d1e5f0", "#67a9cf", "#2166ac"))
boxplot(PredictMcmcEcoZone[,3] ~ Res$GeoClus,
ylab = "Species interaction importance", ylim = c(0,1),
col = c("#b2182b", "#ef8a62", "#fddbc7", "#f7f7f7", "#d1e5f0", "#67a9cf", "#2166ac"))
par(las = 1, mar = c(4.1, 4.1, 3.1, 0.1))
plot(PredictMcmcEcoZone[,1] ~ Res$SpeciesRichness, xlab = "Species Richness",
ylab = "Dispersal limitation importance", pch = 16, col = as.factor(Res$GeoClus))
plot(PredictMcmcEcoZone[,2] ~ Res$SpeciesRichness, xlab = "Species Richness",
ylab = "Environmental filtering importance", pch = 16, col = as.factor(Res$GeoClus))
plot(PredictMcmcEcoZone[,3] ~ Res$SpeciesRichness, xlab = "Species Richness",
ylab = "Species interaction importance", pch = 16, col = as.factor(Res$GeoClus))