Load the data
### Community dataset
comm <- read.table("Community.txt")
# Create a presence-absence matrix
idx <- rep(1:39,each=3)
comm <- sapply(unique(idx), function(x) rowSums(comm[,idx==x]))
colnames(comm) <- seq(1,ncol(comm),1)
comm <- t(comm)
commPA <- ifelse(comm > 0, yes = 1, no = 0)
### Environmental dataset
env <- read.table("Env.txt",h=T)
# Create a coordinates matrix and a sedimentary variables matrix
coord <- env[,2:3]
sed <- env[,-c(1:4)]
# Note how sedimentary variables are highly correlated
cor(sed)
## Mz (mm) Sorting Skewness Sand (%) Fine (%) Ht index WC (%) TOC (%)
## Mz (mm) 1.0000000 -0.8556151 -0.0382096 0.923477 -0.930812 -0.8808441 -0.877664 -0.806925
## Sorting -0.8556151 1.0000000 0.0167622 -0.861804 0.858400 0.9674139 0.898430 0.822186
## Skewness -0.0382096 0.0167622 1.0000000 0.241447 -0.230248 0.0922166 -0.190792 -0.307952
## Sand (%) 0.9234769 -0.8618041 0.2414471 1.000000 -0.999599 -0.8752916 -0.915433 -0.860178
## Fine (%) -0.9308116 0.8584001 -0.2302482 -0.999599 1.000000 0.8740202 0.915649 0.859712
## Ht index -0.8808441 0.9674139 0.0922166 -0.875292 0.874020 1.0000000 0.886271 0.792710
## WC (%) -0.8776639 0.8984303 -0.1907923 -0.915433 0.915649 0.8862712 1.000000 0.892737
## TOC (%) -0.8069246 0.8221856 -0.3079524 -0.860178 0.859712 0.7927100 0.892737 1.000000
# Calculate a PCA to summarize the sedimentary variation
require(stats)
pca <- princomp (sed, cor = T)
summary(pca)
## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## Standard deviation 2.514458 1.073999 0.4850163 0.3920144 0.27516449 0.19517920 0.14558255 1.02937e-02
## Proportion of Variance 0.790312 0.144184 0.0294051 0.0192094 0.00946444 0.00476186 0.00264928 1.32452e-05
## Cumulative Proportion 0.790312 0.934497 0.9639018 0.9831112 0.99257561 0.99733747 0.99998675 1.00000e+00
peso.pc1 <- pca$loadings[,1] * pca$sdev[1]
peso.pc2 <- pca$loadings[,2] * pca$sdev[2]
print(rbind(peso.pc1, peso.pc2))
## Mz (mm) Sorting Skewness Sand (%) Fine (%) Ht index WC (%) TOC (%)
## peso.pc1 0.941223 -0.939290 0.147454 0.972228 -0.9723293 -0.939661 -0.9629611 -0.911966
## peso.pc2 0.170676 -0.180328 -0.982596 -0.102431 0.0925434 -0.253637 0.0532297 0.200272
# Take salinity and sedimentary variation as predictor variables
PC1 <- as.vector(pca$scores[,1])
env <- cbind(Sal=env[,4], PC1)
Redundancy Analysis (RDA)
### RDA
require(vegan)
require(ade4)
# Analysis of the Y, F and R tables (PCA, RDA, PRA)
dist <- vegdist(commPA, method = "bray")
pca.Y <- dudi.pco(dist, scannf = F)
rda.F <- pcaiv(pca.Y, env, scannf = F)
testR2.rda.F <- randtest(rda.F, nrepet = 999)
testR2.rda.F # test of the R2 (portion of community data explained by environmental variables)
## Monte-Carlo test
## Call: randtest.pcaiv(xtest = rda.F, nrepet = 999)
##
## Observation: 0.381911
##
## Based on 999 replicates
## Simulated p-value: 0.001
## Alternative hypothesis: greater
##
## Std.Obs Expectation Variance
## 1.59155e+01 5.35672e-02 4.25616e-04
rda.F$eig / sum(rda.F$eig) * 100 # percentage of variations explained by the first two axes
## [1] 76.8689 23.1311
rda.F$cor # correlation with environmental variables
## RS1 RS2
## Sal 0.977205 0.212298
## PC1 0.247526 -0.968881
pra.R <- pcaivortho(pca.Y, env, scannf = F) # variation not explained by environmental variables
### Spatial predictors (Moran Eigenvector Maps - MEMs)
# Gabriel graph
require(spdep)
nb1 <- graph2nb(gabrielneigh(as.matrix(coord)), sym = T)
lw1 <- nb2listw(nb1)
# MEMs
require(spacemakeR)
U <- scores.listw(lw1)
## listw not symmetric, (w+t(w)) used in the place of w
# Plot scalograms for the axes of the different analyses
nG <- 19 # number of groups
nC <- 2 # number of spatial components
orthosmooth.randtest <- function(x, MEM, nrepet = 999)
{
R2 <- as.vector(cor(x, MEM)^2)
debut <- seq(1,by = nC, length = nG)
fin <- seq(nC,by = nC, length = nG)
R2.smooth <- colSums(sapply(1:nG, function(x) R2[debut[x]:fin[x]]))
sim <- matrix(0, nrepet, nG)
for(i in 1:nrepet)
{
R2.sim <- as.vector(cor(sample(x), MEM)^2)
sim[i, ] <- colSums(sapply(1:nG, function(x) R2.sim[debut[x]:fin[x]]))
}
res <- as.krandtest(sim, R2.smooth)
return(list(res=res,sim=sim))
}
testrda.F1 <- orthosmooth.randtest(rda.F$li[,1], U$vectors)
testrda.F2 <- orthosmooth.randtest(rda.F$li[,2], U$vectors)
testpra.R1 <- orthosmooth.randtest(pra.R$li[,1], U$vectors)
testpra.R2 <- orthosmooth.randtest(pra.R$li[,2], U$vectors)
plotsmooth <- function(obj)
{
par(mar = c(3,2,1.5,0.5))
themax <- which.max(obj$res$obs)
thecol <- rep("grey80", 19)
thecol[themax] <- 'grey40'
debut <- seq(1,by = nC, length = nG)
fin <- seq(nC,by = nC, length = nG)
mp <- barplot(obj$res$obs, col = thecol, ylim = c(0,1), axes = FALSE)
axis(2,at=c(0,1),las=1,tcl=-0.3)
axis(1, at = mp[,1], label = 1:nG, tcl = 0, padj = -1)
mtext(1, text = "Smoothed MEMs", padj = 3)
mtext(2, text = expression(R^2), padj = -1)
qu95 <- apply(obj$sim, 2, quantile, 0.95)
points(mp[,1], qu95, ty = 'b', pch = 3, lty = 3)
text(mp[themax], obj$res$obs[themax], obj$res$pvalue[themax], pos = 3)
}
plotsmooth(testrda.F1)
plotsmooth(testrda.F2)
plotsmooth(testpra.R1)
plotsmooth(testpra.R2)
# Test the significance of Moran's I for MEMs
mctests <- test.scores(U, lw1, 1000)
colnames(U$vectors) <- paste("MEM", 1:ncol(U$vectors), sep="")
U.sel <- U
U.sel$vectors <- U.sel$vectors[ ,mctests[,2]<0.05] # keep only MEMS with signif values
# Forward selection
Usel.frame <- as.data.frame(U.sel$vectors)
rda.0 <- rda(pca.Y$li[,1:2] ~ 1, Usel.frame)
rda.full <- rda(pca.Y$li[,1:2] ~ ., Usel.frame)
rda1.fw.vegan <- ordiR2step(rda.0, formula(rda.full), perm.max = 10000, trace = F)
RsquareAdj(rda1.fw.vegan) # community variation explained by MEMs
## $r.squared
## [1] 0.642113
##
## $adj.r.squared
## [1] 0.63244
### Variation partitioning
index <- match(attr(rda1.fw.vegan$terms,"term.labels"),colnames(U$vectors))
val.sel <- U$values[index]
vec.sel <- U$vectors[,index, drop=F]
broad <- vec.sel[,val.sel>0, drop=F]
fine <- vec.sel[,val.sel<0, drop=F]
colnames(broad) # MEMs selected as spatial predictors
## [1] "MEM1"
colnames(fine) # MEMs selected as spatial predictors
## NULL
pc.varpart <- varpart(pca.Y$li[,1:2],env[,1], env[,2], broad, fine)
plot(pc.varpart, bg="grey80",Xnames=labv)
Components of β-diversity
require(betapart)
require(ecodist)
# Calculate the components of β-diversity
beta <- beta.pair(commPA)
# Calculate the geographic and environmental distance matrices
geodist <- dist(coord)
saldist <- dist(env[,1])
seddist <- dist(env[,2])
# Multiple regression on distance matrices
MRM(1-beta$beta.sor~.,data=decostand(as.data.frame(cbind(geodist,saldist,seddist)),"standardize"))
## $coef
## 1 - beta$beta.sor pval
## Int 0.53868346 0.001
## geodist -0.09828244 0.001
## saldist 0.00411635 0.667
## seddist -0.02380433 0.001
##
## $r.squared
## R2 pval
## 0.485229 0.001000
##
## $F.test
## F F.pval
## 231.568 0.001
MRM(1-beta$beta.sim~.,data=decostand(as.data.frame(cbind(geodist,saldist,seddist)),"standardize"))
## $coef
## 1 - beta$beta.sim pval
## Int 0.70249431 0.001
## geodist -0.08072323 0.001
## saldist -0.00362489 0.838
## seddist 0.00750761 0.255
##
## $r.squared
## R2 pval
## 0.251721 0.001000
##
## $F.test
## F F.pval
## 82.6421 0.0010
MRM(1-beta$beta.sne~.,data=decostand(as.data.frame(cbind(geodist,saldist,seddist)),"standardize"))
## $coef
## 1 - beta$beta.sne pval
## Int 0.83618914 0.002
## geodist -0.01755922 0.267
## saldist 0.00774124 0.627
## seddist -0.03131194 0.002
##
## $r.squared
## R2 pval
## 0.073122 0.002000
##
## $F.test
## F F.pval
## 19.3808 0.0020
Components of β-diversity by estuarine sectors
### Split the sectors
# Community
InnerComm <- commPA[1:15,]
MiddleComm <- commPA[16:29,]
OuterComm <- commPA[30:39,]
# Environment
InnerEnv <- env[1:15,]
MiddleEnv <- env[16:29,]
OuterEnv <- env[30:39,]
# Coordinates
InnerCoord <- coord[1:15,]
MiddleCoord <- coord[16:29,]
OuterCoord <- coord[30:39,]
### Calculate the components of β-diversity for each sector
# Here an example using the inner sector
beta <- beta.pair(InnerComm)
geodist <- dist(InnerCoord)
saldist <- dist(InnerEnv[,1])
seddist <- dist(InnerEnv[,2])
# Multiple regression on distance matrices
MRM(1-beta$beta.sor~.,data=decostand(as.data.frame(cbind(geodist,saldist,seddist)),"standardize"))
## $coef
## 1 - beta$beta.sor pval
## Int 0.68337981 0.001
## geodist -0.00809516 0.511
## saldist -0.01824340 0.171
## seddist -0.05331592 0.001
##
## $r.squared
## R2 pval
## 0.296683 0.001000
##
## $F.test
## F F.pval
## 14.2017 0.0010
MRM(1-beta$beta.sim~.,data=decostand(as.data.frame(cbind(geodist,saldist,seddist)),"standardize"))
## $coef
## 1 - beta$beta.sim pval
## Int 0.86591390 0.001
## geodist 0.00404359 0.823
## saldist -0.03256782 0.087
## seddist 0.04081437 0.002
##
## $r.squared
## R2 pval
## 0.171943 0.001000
##
## $F.test
## F F.pval
## 6.99078 0.00100
MRM(1-beta$beta.sne~.,data=decostand(as.data.frame(cbind(geodist,saldist,seddist)),"standardize"))
## $coef
## 1 - beta$beta.sne pval
## Int 0.8174659 0.001
## geodist -0.0121387 0.400
## saldist 0.0143244 0.324
## seddist -0.0941303 0.001
##
## $r.squared
## R2 pval
## 0.498212 0.001000
##
## $F.test
## F F.pval
## 33.4267 0.0010