# Script for Qiongdao et al: Local-scale determinants of arboreal spider beta diversity in a temperate forest
# prepared by Qiongdao Zhang and Dong He
# 24 February 2018

## DATA MANAGEMENT

spider.abd<-read.table(file="spider abd.csv",header=T, row.names=1)
spider.trans<-decostand(t(spider.abd),"hellinger")

spider.env<-read.table("spider env.csv", header=T, row.names=1)

spider.xy<-read.csv("spider.xy.csv", header=T, row.names=1)

spider.trait<-read.csv("spider trait.csv", header=T, row.names=1)

### Species compositional dissimilarity among samples

spider.abd.rand<-array(0,dim=c(289,289,1000))
for(i in 1:1000){
      print(1:i)
      rand<-randomizeMatrix(spider.abd,null.model="richness",iterations=1) 
      beta.value<-as.matrix(bray.part(t(rand))$bray)
      diag(beta.value)<-0
      beta.value[is.nan(beta.value)]=0
      spider.abd.rand[,,i]<-beta.value
}

spider.beta.rand<-array(0,dim=c(1000,41616))
for(i in 1:1000){
      beta<-spider.abd.rand[,,i][lower.tri(spider.abd.rand[,,i])]
      spider.beta.rand[i,]<-beta
}

df<-apply(spider.beta.rand,2,function(x) quantile(x,probs = c(2.5,5,95,97.5)/100))



spider.beta.data<-cbind(spider.abd.beta.obs,obs.ses)

dat<-array(0,dim=c(41616,1))
for(i in 1:41616){
      if(spider.beta.data[i,2]<spider.beta.data[i,1] & 
           spider.beta.data[i,1]<spider.beta.data[i,5])
      {dat[i,]<-"Null"
}
      else {
            if(spider.beta.data[i,1]>spider.beta.data[i,5])
            dat[i,]<-"Positive"
            if(spider.beta.data[i,1]<spider.beta.data[i,2])
            dat[i,]<-"Negative"
}
               }

subset(spider.beta.data,obs.ses<-2)
dat<-subset(spider.beta.data,obs.ses<-2,type="neg")
factor(dat[,7])

## Transformation of Standardized Effect Size 
spider.abd.beta.obs<-as.matrix(bray.part(t(spider.abd))$bray)
diag(spider.abd.beta.obs)<-0
spider.abd.beta.obs[is.nan(spider.abd.beta.obs)]=0
spider.abd.beta.obs<-spider.abd.beta.obs[lower.tri(spider.abd.beta.obs)]

rand.mean<-apply(spider.abd.rand,c(1,2),mean,na.rm=T)
rand.sd<-apply(spider.abd.rand,c(1,2),sd,na.rm=T)
obs.ses<-(spider.abd.beta.obs-rand.mean)/rand.sd

diag(obs.ses)<-0
obs.ses[is.nan(obs.ses)]=0
obs.ses<-obs.ses[lower.tri(obs.ses)]

# plotting
hist(spider.beta.data$obs.ses,xlab=" Bray-Curtis dissimilarity S.E.S. ",ylab=" Number of pairwise"
     , breaks=seq(-6,4,by=0.5))
abline(v=0,col=2)

png(file="beta diversity.png",height=7,width=7,res=1000, units = "in")

ggplot(spider.beta.data, aes(obs,fill=Types)) +
  geom_histogram() +
    xlab("Bray-Curtis Dissimilarity") +
    ylab("Number of pairwise")  + 
    labs(fill="Significant levels") +     
    theme(axis.title.x=element_text(size = 15),
          axis.title.y=element_text(size = 15),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          axis.line = element_line(colour = "black"),
          legend.title=element_text( size =15),
          legend.text = element_text( size =15))
   
dev.off()
          
p2<-ggplot(spider.env, aes(Height,fill=Tree.Age)) +
  geom_histogram() +
    xlab("Total Height") +
    ylab("Number of pairwise")  + 
    labs(fill="Levels") +     
    theme(axis.title.x=element_text(size = 15),
          axis.title.y=element_text(size = 15),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          axis.line = element_line(colour = "black"),
          legend.title=element_text( size =15),
          legend.text = element_text( size =15))


###Inventory completeness in vegan
specpool(t(spider.abd))
pool<-apply(t(spider.abd.C),2,sum)
estimateR(pool)

### RDA triplot of the spider species
test<-spider.env[,c(5,12,13)]
spider.rda<-rda(spider.trans,test)

png(file="Triplot spider RDA.png",height=7,width=7,res=600, units = "in")

plot(spider.rda,xlab="RDA1=0.5964",ylab="RDA2=0.2581",scaling=2)
plot(spider.rda,display=c("wa","cn"),xlab="RDA1(56.4%)",ylab="RDA2(29.1%)",scaling=2)
legend("topright",c("R=0.1385"))

dev.off()

### PCNM Analysis and Variation Partitioning

#1.Test trend. If significant, forword selection of coordinates

spider.trans<-decostand(t(spider.abd),"hellinger")
spider.XY.rda<-rda(spider.trans~.,spider.xy)
anova.cca(spider.XY.rda)
spider.XY.R2a<-RsquareAdj(spider.XY.rda)$adj.r.squared
spider.XY.fwd<-forward.sel(as.data.frame.matrix(spider.trans),as.matrix(spider.xy),
adjR2thresh=spider.XY.R2a)
{
ordistep(rda(spider.trans~1,data=spider.xy),
adjR2thresh=spider.XY.R2a,scope=formula(spider.XY.rda),
direction="forward",pstep=10000)
XY.sign<-sort(spider.XY.fwd$terminfo$ordered)
}

XY.sign<-sort(spider.XY.fwd$order)
XY.red<-spider.xy[,c(XY.sign)]

#2.Test and forword selection of environmental variables

pc.vegetation<-rda(spider.env[,c(5:11,13,14)],scale=T)
ev<-pc.vegetation$CA$eig
ev[ev>mean(ev)]
rda(spider.trans~.,as.data.frame(pc.vegetation$CA$u[,1:4]))

pc.env<-as.data.frame(cbind(spider.env[,c(15,18:19)],pc.vegetation$CA$u[,1:4]))

spider.env.rda<-rda(spider.trans~.,pc.env)
spider.env.R2a<-RsquareAdj(spider.env.rda)$adj.r.squared
spider.env.fwd<-

            ordiR2step(rda(spider.trans~1,data=pc.env),
            adjR2thresh=spider.env.R2a,scope=formula(spider.env.rda),
            direction="forward",Pin = 0.1, Pout = 0.1,,pstep=10000)

forward.sel(spider.trans,pc.env,alpha = 0.1,adjR2thresh=spider.env.R2a)

#3.Test and forword selection of PCNM variables

library(spdep)
library(boot)

spider.xy.dist<-dist(spider.xy)
spider.PCNM.auto<-PCNM(spider.xy.dist)
select<-which(spider.PCNM.auto$Moran_I$Positive==TRUE)
spider.PCNM.pos<-as.data.frame(spider.PCNM.auto$vectors)[,select]

spider.PCNM.rda<-rda(spider.det,spider.PCNM.pos)
anova.cca(spider.PCNM.rda)

spider.undet.PCNM.rda<-rda(spider.trans~.,spider.PCNM.pos)
anova.cca(spider.undet.PCNM.rda,by="axis")

spider.undet.PCNM.R2a<-RsquareAdj(spider.undet.PCNM.rda)$adj.r.squared
spider.undet.PCNM.fwd<-forward.sel(as.data.frame.matrix(spider.trans),as.matrix(spider.PCNM.pos),
,alpha=0.1,adjR2thresh=spider.undet.PCNM.R2a)
nrow(spider.undet.PCNM.fwd)

PCNM.sign<-sort(spider.undet.PCNM.fwd$order)
PCNM.red<-spider.PCNM.pos[,c(PCNM.sign)]

## Variation Partitioning

varpart(spider.trans,spider.xy,PCNM.red)

#Broad scale
PCNM.rda<-rda(spider.trans,spider.xy)
RsquareAdj(PCNM.rda)
PCNM.abd<-fitted.values(PCNM.rda)

test<-varpart(PCNM.abd,spider.env[,c(5,12,13,14)],PCNM.abd)

#Fine scale
PCNM.rda<-rda(spider.det,PCNM.red)
RsquareAdj(PCNM.rda)
PCNM.abd<-fitted.values(PCNM.rda)

test<-varpart(PCNM.abd,spider.env[,c(5,12,13,14)],PCNM.abd)

plot(test,bg=c(3,158),digits=3)


## Permutation test 

png(file="Venn.png",height=8.0,width=8,res=600, units = "in")
plot()
dev.off()

spider.bray<-vegdist(t(spider.abd),"bray")
spider.pcoa<-cmdscale(spider.bray,k=nrow(t(spider.abd))-2,eig=T,add=T)
spider.scores<-spider.pcoa$points

test<-varpart(spider.trans,spider.env[,c(5,12,13,14)],spider.xy,PCNM.red)

showvarparts(3,bg=c(4,7,10),Xnames=c("","",""))
plot(test,bg=c(4,7,10),cutoff=0,digits=2,cex=1,Xnames=c("","",""))
text(locator(1),"Vegetation structure",cex=1)

anova.cca(rda(spider.trans,spider.env[,c(5,12,13,14)],
cbind(PCNM.red,spider.xy)))

anova.cca(rda(spider.bray,spider.xy,
cbind(PCNM.red,spider.env2[,c(5,12,13,14)])))

anova.cca(rda(spider.bray,PCNM.red,
cbind(spider.xy,spider.env2[,c(5,12,13,14)])))

## Mantel and Partial Mantel Tests for Dissimilarity Matrices

mantel(vegdist(spider.trans,"euc"),vegdist(spider.env[,c(5,12:14)],"euc"))

### Relative influence of spatially structured environmental variation (Environment กษ Space) 
and pure spatial variation on the spider composition 
at the broad versus fine scales
#Percentage bar graph 

x<-data.frame(rbind(c(7553,1405),c(1054,7790)))
names(x)<-c("Broad scale","Finer scales")

test<-melt(cbind(x, class = c("Enviroment+Space","Space")), id.vars = c("class"))
colnames(test)[3]<-"variation"


ggplot(test,aes(x = variable, y = variation,width=0.3,fill = class)) + 
    geom_bar(position = "fill",stat = "identity") + 
    xlab(" ") +
    ylab("Spatial variation (%)") +
    scale_y_continuous(labels = percent_format()) +
    scale_y_continuous(labels = percent_format()) +
     theme_bw() + 
 theme(axis.title.x=element_text(size = 15),
          axis.title.y=element_text(size = 15),
          legend.justification=c(1,0),
legend.text = element_text( size =15),
          legend.position="bottom")




#######################################
Variation partitioning analyses of beta diversity along environmental and 
spatial gradients were conducted for each spider dispersal capacity class
###The spider group with high vagility

test<-spider.trait[rownames(spider.abd),]

sp<-subset(spider.trait,ballooning.ability=="Common"|ballooning.ability=="Less common"|ballooning.ability=="Uncommon")
spider.abd.uncom<-spider.abd[rownames(sp),]
site<-names(which(apply(spider.abd.uncom,2,sum)>0))

spider.abd.uncom<-spider.abd.uncom[,site]

spider.uncom.xy<-spider.xy[colnames(spider.abd.uncom),]

## PCNM Analysis and Variation Partitioning

#1.Test trend. If significant, forword selection of coordinates

spider.com.trans<-decostand(t(spider.abd.com),"hellinger")

spider.com.XY.rda<-rda(spider.com.trans~.,spider.com.xy)
anova.cca(spider.com.XY.rda)
spider.com.XY.R2a<-RsquareAdj(spider.com.XY.rda)$adj.r.squared
spider.com.XY.fwd<-forward.sel(as.data.frame.matrix(spider.com.trans),as.matrix(spider.com.xy),
adjR2thresh=spider.com.XY.R2a)
{
ordistep(rda(spider.trans~1,data=test.xy),
adjR2thresh=spider.XY.R2a,scope=formula(spider.XY.rda),
direction="forward",pstep=10000)
XY.sign<-sort(spider.XY.fwd$terminfo$ordered)
}

com.XY.sign<-sort(spider.com.XY.fwd$order)
com.XY.red<-spider.com.xy[,c(com.XY.sign)]

#2.Test and forword selection of environmental variables

com.pc.vegetation<-rda(decostand(spider.com.env[,c(5:11,13,14)],"hel"))
ev<-com.pc.vegetation$CA$eig
ev[ev>mean(ev)]
rda(spider.com.trans~.,as.data.frame(com.pc.vegetation$CA$u[,1:2]))


#3.Test and forword selection of PCNM variables

spider.com.undet.PCNM.rda<-rda(spider.com.trans~.,spider.com.PCNM.pos)
anova.cca(spider.com.undet.PCNM.rda,by="axis")

spider.com.undet.PCNM.R2a<-RsquareAdj(spider.com.undet.PCNM.rda)$adj.r.squared
spider.com.undet.PCNM.fwd<-forward.sel(as.data.frame.matrix(spider.com.trans),as.matrix(spider.com.PCNM.pos),
adjR2thresh=spider.com.undet.PCNM.R2a)
nrow(spider.com.undet.PCNM.fwd)

com.PCNM.sign<-sort(spider.com.undet.PCNM.fwd$order)
com.PCNM.red<-spider.com.PCNM.pos[,c(com.PCNM.sign)]

#4.Variation Partitioning

spider.com.env<-spider.env[rownames(spider.com.trans),]


test<-varpart(spider.com.trans,spider.com.env[,c(5,12:14)],com.PCNM.red)

plot(test,bg=c(4,7,10),cutoff=0,digits=2,cex=1,Xnames=c("","",""))
text(locator(1),"Vegetation structure")

anova.cca(rda(spider.com.trans,spider.com.env[,c(5:11,15,18,19)],
com.PCNM.red)))

anova.cca(rda(spider.trans,PCNM.red[,2:6],
cbind(PCNM.red[,1],pc.env)))

anova.cca(rda(spider.trans,pc.env[,3:4],
cbind(pc.env[,1:2],PCNM.red)))

### The spider group with low vagility

#1.Test trend. If significant, forword selection of coordinates

spider.uncom.trans<-decostand(t(spider.abd.uncom),"hellinger")

spider.uncom.XY.rda<-rda(spider.uncom.trans~.,spider.uncom.xy)
anova.cca(spider.uncom.XY.rda)
spider.uncom.XY.R2a<-RsquareAdj(spider.uncom.XY.rda)$adj.r.squared
spider.uncom.XY.fwd<-forward.sel(as.data.frame.matrix(spider.uncom.trans),
as.matrix(spider.uncom.xy),adjR2thresh=spider.uncom.XY.R2a)

{
ordistep(rda(spider.trans~1,data=test.xy),
adjR2thresh=spider.XY.R2a,scope=formula(spider.XY.rda),
direction="forward",pstep=10000)
XY.sign<-sort(spider.XY.fwd$terminfo$ordered)
}

uncom.XY.sign<-sort(spider.uncom.XY.fwd$order)
uncom.XY.red<-spider.uncom.xy[,c(uncom.XY.sign)]

#2.Test and forword selection of environmental variables

uncom.pc.vegetation<-rda(decostand(spider.uncom.env2[,c(5:11,13,14)],"hel"))
ev<-uncom.pc.vegetation$CA$eig
ev[ev>mean(ev)]
rda(spider.uncom.trans~.,as.data.frame(uncom.pc.vegetation$CA$u[,1:2]))

uncom.pc.env<-as.data.frame(cbind(uncom.pc.climate$CA$u[,1:2],
uncom.pc.vegetation$CA$u[,1:2]))

spider.uncom.env.rda<-rda(spider.uncom.trans~.,uncom.pc.env)
spider.uncom.env.R2a<-RsquareAdj(spider.uncom.env.rda)$adj.r.squared
spider.uncom.env.fwd<-forward.sel(as.data.frame.matrix(spider.uncom.trans),
as.matrix(uncom.pc.env),adjR2thresh=spider.uncom.env.R2a)


ordiR2step(rda(spider.com.trans~1,data=com.pc.env),
adjR2thresh=spider.com.env.R2a,scope=formula(spider.com.env.rda),
direction="forward",pstep=10000)

#3.Test and forword selection of PCNM variables

spider.uncom.undet.PCNM.rda<-rda(spider.uncom.trans~.,spider.uncom.PCNM.pos)
anova.cca(spider.uncom.undet.PCNM.rda,by="axis")

spider.uncom.undet.PCNM.R2a<-RsquareAdj(spider.uncom.undet.PCNM.rda)$adj.r.squared
spider.uncom.undet.PCNM.fwd<-forward.sel(as.data.frame.matrix(spider.uncom.trans),
as.matrix(spider.uncom.PCNM.pos),adjR2thresh=spider.uncom.undet.PCNM.R2a)

nrow(spider.uncom.undet.PCNM.fwd)

uncom.PCNM.sign<-sort(spider.uncom.undet.PCNM.fwd$order)
uncom.PCNM.red<-spider.uncom.PCNM.pos[,c(uncom.PCNM.sign)]


#4.Variation Partitioning

spider.uncom.env<-spider.env[rownames(spider.uncom.trans),]

test<-varpart(spider.uncom.trans,spider.uncom.env[,c(5,12:14)],uncom.PCNM.red)


plot(test,bg=c(4,7,10),cutoff=-Inf,digits=2)
text(locator(1),"PC1-climate")

anova.cca(rda(spider.trans,pc.env,
cbind(PCNM.red,PCNM.red)))

anova.cca(rda(spider.trans,PCNM.red[,2:6],
cbind(PCNM.red[,1],pc.env)))

anova.cca(rda(spider.trans,pc.env[,3:4],
cbind(pc.env[,1:2],PCNM.red)))