#RECREATING FIGURE 8 FROM:
#Tomašových A., and Kidwell S.M. 2009. Preservation of spatial and environmental gradients by death assemblages. Paleobiology 35:122-148
#2) BETWEEN-SAMPLE FIDELITY - HOW STRONGLY ARE COMPOSITIONAL RELATIONSHIPS BETWEEN MULTIPLE SAMPLES IN ONE REGION AFFECTED BY RANDOM SPECIES DEGRADATION?
#IN OTHER WORDS, THIS ALLOWS TESTING HOW MANY SPECIES NEED TO BE REMOVED IN ORDER TO LOSE ORIGINAL SIMILARITIES AMONG SAMPLES,
#SLOW DECAY IN SIMILARITY CAN IMPLY STRUCTURAL REDUNDANCY IN LOCAL SAMPLES, SO THAT EVEN WHEN SOME SPECIES ARE LOST,
#SIMILARITIES AMONG SAMPLES CAN BE STILL PRESERVED IF OTHER SPECIES RESPOND TO ENVIRONMENTAL GRADIENTS IN A SIMILAR WAY
#COMPUTING SPEARMAN RANK CORRELATION IN DISSIMILARITIES BETWEEN COMPLETE AND DEGRADED SET OF MULTIPLE LOCAL COMMUNITIES
#example file - dataset with 10 living mollusk-brachiopod assemblages sampled in subtidal habitats of San Juan Islands
#Kowalewski et al. 2003, Journal of Taphonomy
#load(file="Living assemblages-San Juan Islands.Rdata") will load "abslivematrix" (includes also undetermined species)
library(vegan) #vegan needs to be installed from R packages
#INPUT
abund<-abslivematrix #input table with species abundance or presence-absence in rows and samples in columns
#DOUBLECHECKING THAT ALL SPECIES HAVE AT LEAST ONE INDIVIDUAL
abund<-t(abund[apply(abund, MARGIN=1, sum)>0,])
#INITIAL ALLOCATION
pairwise<-(nrow(abund)*(nrow(abund)-1))/2 #number of pairwise dissimilarities per dataset
simulations<-100 #number of simulations
species<-ncol(abund) #total number of species
ninetypercentofspecies<-ceiling(species-(species/100)*10) #WARNING: 90% of species are removed, to avoid the possibility that only one sample with non-zero abundances remains
bcrankcor<-array(NA, dim=c(ninetypercentofspecies,simulations)) #output of correlation values per removal of 1, 2...N species per simulation for Bray-Curtis
sorrankcor<-array(NA, dim=c(ninetypercentofspecies,simulations)) #output of correlation values per removal of 1, 2...N species per simulation for Sorenson
bcsignif<-array(NA, dim=c(ninetypercentofspecies,simulations)) #output of p values per removal of 1, 2...N species per simulation for Bray-Curtis
sorsignif<-array(NA, dim=c(ninetypercentofspecies,simulations)) #output of p values per removal of 1, 2...N species per simulation for Sorenson
averagebcsignif<-numeric();averagesorsignif<-numeric()
averagebcrankcor<-numeric();minusbcrankcor<-numeric()
plusbcrankcor<-numeric();averagesorrankcor<-numeric()
minussorrankcor<-numeric(); plussorrankcor<-numeric()
permutations<-100 #number of permutations in Mantel test
#DISSIMILARITIES BASED ON SUBSETS OF SPECIES, RANDOMLY REMOVING (DEGRADING) SPECIES FROM THE FULL DATASET
for (l in 1:ninetypercentofspecies) {
for (k in 1:simulations) {
random<-sample(c(1:species), l, replace=F)
abund2<-t(abund[,-random])
#detection and removal of samples without species
zerosamples<-apply(abund2,MARGIN=2,sum)
newabund2<-abund2[,zerosamples>0]
Qmodeabund<-decostand(newabund2, method="total", MARGIN=2)
degradedbc<-vegdist(t(Qmodeabund), method = "bray")
degradedsor<-vegdist(t(Qmodeabund), method = "bray", binary=TRUE)
Qmodeabund<-decostand(abund[zerosamples>0,], method="total", MARGIN=1)
completebc<-vegdist(Qmodeabund, method = "bray")
completesor<-vegdist(Qmodeabund, method = "bray", binary=TRUE)
bcsignif[l,k]<- mantel(completebc, degradedbc, permutations=permutations, method = "spearman")$signif
sorsignif[l,k]<- mantel(completesor, degradedsor, permutations=permutations, method = "spearman")$signif
bcrankcor[l,k]<-cor(completebc, degradedbc, method="spearman", use="complete.obs")
sorrankcor[l,k]<-cor(completesor, degradedsor, method="spearman", use="complete.obs")
}
pvaluebcvector<-bcsignif[l,]
pvaluebcvector<-pvaluebcvector[!is.na(pvaluebcvector)]
pvaluesorvector<-sorsignif[l,]
pvaluesorvector<-pvaluesorvector[!is.na(pvaluesorvector)]
averagebcsignif[l]<-length(pvaluebcvector[pvaluebcvector<0.05])/length(pvaluebcvector)
averagesorsignif[l]<-length(pvaluesorvector[pvaluesorvector<0.05])/length(pvaluesorvector)
averagebcrankcor[l]<-mean(bcrankcor[l,], na.rm=TRUE)
minusbcrankcor[l]<-quantile(bcrankcor[l,], c(0.975), na.rm=TRUE)
plusbcrankcor[l]<-quantile(bcrankcor[l,], c(0.025), na.rm=TRUE)
averagesorrankcor[l]<-mean(sorrankcor[l,], na.rm=TRUE)
minussorrankcor[l]<-quantile(sorrankcor[l,], c(0.975), na.rm=TRUE)
plussorrankcor[l]<-quantile(sorrankcor[l,], c(0.025), na.rm=TRUE)
print(l)
}
#PLOT
par(pty="s")
par(cex=1.4)
plot(1:ninetypercentofspecies,averagebcrankcor, type="l",ylim=c(-0.4,1), ylab="Sperman rank correlation (r)",
main="Correlation between dissimilarities in complete and degraded dataset",
xlab="Number of degraded species", lwd=3, cex.main=0.8)
lines(1:ninetypercentofspecies,minusbcrankcor, col="black", lwd=1)
lines(1:ninetypercentofspecies,plusbcrankcor, col="black",lwd=1)
lines(1:ninetypercentofspecies,averagesorrankcor, col="grey41", lwd=3)
lines(1:ninetypercentofspecies,minussorrankcor, col="grey41", lwd=1)
lines(1:ninetypercentofspecies,plussorrankcor, col="grey41",lwd=1)
lines(1:ninetypercentofspecies,averagebcsignif, col="black", lwd=1, lty=2)
lines(1:ninetypercentofspecies,averagesorsignif, col="grey41",lwd=1, lty=2)
legend(1,0, c("Bray-Curtis", "Sorenson", "proportion of p<0.05 [BC]","proportion of p<0.05 [Sor]"),
col=c("black", "grey41","black", "grey41"), lwd=c(2,2,1,1), lty=c(1,1,2,2), cex=0.8, bty="n")
########################################################################################################################################
#SIMULATION OF RANDOM SPECIES DEGRADATION AND ITS EFFECT ON WITHIN-SAMPLE FIDELITY MEASURED BY PAIRWISE SPEARMAN RANK CORRELATION IN
#SPECIES PROPORTIONAL ABUNDANCES BETWEEN TWO SAMPLES (COMPLETE AND DEGRADED)
#1) WITHIN-SAMPLE FIDELITY - HOW STRONGLY ARE SPECIES ABUNDANCES within ONE SAMPLE AFFECTED BY RANDOM SPECIES DEGRADATION?
#COMPUTING SPEARMAN RANK CORRELATION IN ABUNDANCES BETWEEN COMPLETE AND DEGRADED LOCAL SAMPLE (LOCAL COMMUNITY),
#WITH THE RESULTING CORRELATION BEING AVERAGED ACROSS THE NUMBER OF LOCAL SAMPLES IN THE FULL REGIONAL DATASET
#INPUT
abund<-abslivematrix #input table with species abundance or presence-absence in rows and samples in columns
#DOUBLECHECKING THAT ALL SPECIES HAVE AT LEAST ONE INDIVIDUAL
abund<-t(abund[apply(abund, MARGIN=1, sum)>0,])
#SPEARMAN WITHIN-HABITAT FIDELITY IN ABUNDANCES - SIMULATION OF REDUCED RICHNESS
simulations=100
species<-ncol(abund) #total number of species
fidelity<-array(NA, dim=c(nrow(abund), simulations))
signif<-array(NA, dim=c(nrow(abund), simulations))
meanfidelity<-array(NA, dim=c(species-1, simulations))
meansignif<-array(NA, dim=c(species-1, simulations))
matrix1<-array(NA, dim=c(nrow(abund), ncol(abund)))
matrix2<-array(NA, dim=c(nrow(abund), ncol(abund)))
grandmeanfidelity<-rep(NA, species-1)
grandplusfidelity<-rep(NA, species-1)
grandminusfidelity<-rep(NA, species-1)
grandsignif<-rep(NA, species-1)
ninetypercentofspecies<-ceiling(species-(species/100)*10)
for (l in 1:ninetypercentofspecies) {
for (k in 1:simulations) {
random<-sample(c(1:species), l, replace=F)
abund2<-t(abund[,-random])
liveabund<-decostand(abund,method="total")
deadabund<-decostand(abund,method="total")
deadabund[,random]<-0
#loop for Spearman rank correlation so that species with zeros in two samples are avoided
for (m in 1:nrow(abund)) { #nrow - SAMPLES
for (n in 1:ncol(abund)) { #ncol - SPECIES
if (liveabund[m,n] == 0 & deadabund[m,n] == 0)
{matrix1[m,n]<-NA; matrix2[m,n]<-NA}
else
{matrix1[m,n]<-liveabund[m,n]; matrix2[m,n]<-deadabund[m,n]
}
}
fidelity[m,k]<-cor(matrix1[m,], matrix2[m,], use="complete.obs",method="spearman")
if (!is.na(fidelity[m,k])) {signif[m,k]<-cor.test(matrix1[m,], matrix2[m,], use="complete.obs",method="spearman")$p.value}}
meanfidelity[l,k]<-mean(fidelity[,k], na.rm=TRUE) #fidelity - 100 simulations for m samples
meansignif[l,k]<-mean(signif[,k], na.rm=TRUE)
}
grandmeanfidelity[l]<-mean(meanfidelity[l,], na.rm=TRUE)
grandplusfidelity[l]<-quantile(meanfidelity[l,],c(0.975), na.rm=TRUE)
grandminusfidelity[l]<-quantile(meanfidelity[l,],c(0.025), na.rm=TRUE)
grandsignif[l]<-mean(meansignif[l,], na.rm=TRUE)
}
#PLOT
par(pty="s")
par(cex=1.4)
plot(1:(species-1),grandmeanfidelity, type="l",ylim=c(-0.5,1),
ylab="Average Spearman rank correlation", xlab="Number of degraded species", lwd=3)
lines(1:(species-1),grandminusfidelity, col="black", lwd=1)
lines(1:(species-1),grandplusfidelity, col="black",lwd=1)