library("vegan") library("mvabund") # input and transform Thetford datasets arthropods=read.table("arthropods.txt", header=TRUE, sep="\t", row.names=1,na.strings="NA", dec=".", strip.white=TRUE) arthropods <- arthropods[,which(specnumber(arthropods,MARGIN=2)>0)] # remove spp without records arthropods.env=read.table("arthropods.env.txt", header=TRUE, sep="\t", row.names=1,na.strings="NA", dec=".", strip.white=TRUE) arthropodsB=read.table("MBCarthropods.txt", header=TRUE, sep="\t", row.names=1,na.strings="NA", dec=".", strip.white=TRUE) arthropodsB=data.frame(t(arthropodsB)) arthropodsB=decostand(arthropodsB, method = "log", logbase = Inf) # convert to presence/absence # arthropodsB=1*(arthropodsB>0) # alternative way to convert to presence/absence # remove Ride43, which was a Control site that got mistakenly mowed arthropods <- arthropods[-42,] arthropods.xy <- arthropods.xy[-42,] arthropods.env <- arthropods.env[-42,] arthropodsB <- arthropodsB[-42,] # create datasets without singleton species (OTUS present in only one site) arthropods.all <- arthropods # with singletons, for specpool analysis arthropods <- arthropods[,which(specnumber(arthropods,MARGIN=2)>1)] # remove singleton spp arthropods.dcs=decostand(arthropods, method = "log", logbase = Inf) # convert to presence/absence arthropodsB.all <- arthropodsB # with singletons, for specpool analysis arthropodsB <- arthropodsB[,which(specnumber(arthropodsB,MARGIN=2)>1)] # remove singleton spp # create datasets without singletons and without high-prevalence species arthropods_no45plusprevalents=arthropods[,which(specnumber(arthropods,MARGIN=2)<46)] # 45 is max number of MIDs: no OTUs present in 45+ MIDs arthropodsB_no45plusprevalents=arthropodsB[,which(specnumber(arthropodsB,MARGIN=2)<46)] # 45 is max number of MIDs: no OTUs present in 45+ MIDs # make a bunch of treatment variables colnames(arthropods.env) avpctbare <- arthropods.env[,1] pctcanopy <- arthropods.env[,2] avwidth <- arthropods.env[,3] # no values for Heath ph <- arthropods.env[,4] treatment <- arthropods.env[,5] # make datasets without Heath sites (to test the effects of the treatments only) arthropods.noheath <- arthropods[c(-61,-62,-63,-64,-65,-66,-67),] arthropods.noheath.all <- arthropods.all[c(-61,-62,-63,-64,-65,-66,-67),] arthropods.noheath.h <- arthropods.h[c(-61,-62,-63,-64,-65,-66,-67),] arthropodsB.noheath <- arthropodsB[c(-61,-62,-63,-64,-65,-66,-67),] arthropodsB.noheath.all <- arthropodsB.all[c(-61,-62,-63,-64,-65,-66,-67),] arthropods.noheath.env <- arthropods.env[c(-61,-62,-63,-64,-65,-66,-67),] arthropods.noheath.xy <- arthropods.xy[c(-61,-62,-63,-64,-65,-66,-67),] treatment.noheath <- treatment[c(-61,-62,-63,-64,-65,-66,-67)] avpctbare.noheath <- avpctbare[c(-61,-62,-63,-64,-65,-66,-67)] pctcanopy.noheath <- pctcanopy[c(-61,-62,-63,-64,-65,-66,-67)] avwidth.noheath <- avwidth[c(-61,-62,-63,-64,-65,-66,-67)] ph.noheath <- ph[c(-61,-62,-63,-64,-65,-66,-67)] ########## Mantel and Procrustes comparisons with Heath sites ########## arthropodsB.jdis <- vegdist(arthropodsB,method="jaccard",binary=TRUE) # binary for presence/absence arthropods.qjdis <- vegdist(arthropods,method="jaccard") # quantitative mantel(arthropods.qjdis,arthropodsB.jdis) # highly significantly correlated p = 0.001 (quantitative vs binary) # Using the dissimilarity matrices made previously for NMDS (binary jaccard) arthropods.qjmds <- metaMDS(arthropods.qjdis, wascores=FALSE, autotransform=FALSE,noshare=FALSE) arthropods.qjmds <- metaMDS(arthropods.qjdis, wascores=FALSE, autotransform=FALSE,noshare=FALSE, previous = arthropods.qjmds) plot(arthropods.qjmds,display="sites",type="n") with(arthropods.env,ordisurf(arthropods.qjmds,avpctbare,col="red",main="")) with(arthropods.env,ordiellipse(arthropods.qjmds, treatment,cex=.5, draw="polygon",col=c("red"),alpha=10,kind="se",conf=0.95,label=TRUE, show.groups=(c("Heath")))) with(arthropods.env,ordiellipse(arthropods.qjmds, treatment,cex=.5, draw="polygon",col=c("black"),alpha=100,kind="se",conf=0.95,label=TRUE, show.groups=(c("1Control")))) with(arthropods.env,ordiellipse(arthropods.qjmds, treatment,cex=.5, draw="polygon",col=c("blue"),alpha=10,kind="se",conf=0.95,label=FALSE, show.groups=(c("ForageHarvest","ForestPlough","Disc","Swipe")))) with(arthropods.env,ordiellipse(arthropods.qjmds, treatment,cex=.5, draw="polygon",col=c("green"),alpha=50,kind="se",conf=0.95,label=TRUE, show.groups=(c("TurfStrip","AgriPlough")))) setwd("~/Dropbox/Working Docs/Barcode soup censuses/FieldVal paper/figures") arthropodsB.jmds <- metaMDS(arthropodsB.jdis, wascores=FALSE, autotransform=FALSE,noshare=FALSE) arthropodsB.jmds <- metaMDS(arthropodsB.jdis, wascores=FALSE, autotransform=FALSE,noshare=FALSE, previous = arthropodsB.jmds) plot(arthropodsB.jmds,display="sites",type="n") with(arthropods.env,ordisurf(arthropodsB.jmds, avpctbare,col="red",main="")) with(arthropods.env,ordiellipse(arthropodsB.jmds, treatment,cex=.5, draw="polygon",col=c("red"),alpha=10,kind="se",conf=0.95,label=TRUE, show.groups=(c("Heath")))) with(arthropods.env,ordiellipse(arthropodsB.jmds, treatment,cex=.5, draw="polygon",col=c("black"),alpha=100,kind="se",conf=0.95,label=TRUE, show.groups=(c("1Control")))) with(arthropods.env,ordiellipse(arthropodsB.jmds, treatment,cex=.5, draw="polygon",col=c("blue"),alpha=10,kind="se",conf=0.95,label=FALSE, show.groups=(c("Disc","ForageHarvest","Swipe")))) with(arthropods.env,ordiellipse(arthropodsB.jmds, treatment,cex=.5, draw="polygon",col=c("green"),alpha=10,kind="se",conf=0.95,label=TRUE, show.groups=(c("AgriPlough")))) with(arthropods.env,ordiellipse(arthropodsB.jmds, treatment,cex=.5, draw="polygon",col=c("green"),alpha=70,kind="se",conf=0.95,label=TRUE, show.groups=(c("TurfStrip","ForestPlough")))) # Procrustes test proarthtest <- protest(arthropods.qjmds,arthropodsB.jmds,symmetric=TRUE) plot(proarthtest) proarthtest # highly significantly correlated, p=0.001 ####### Mantel and Procrustes comparisons without Heath sites ########## arthropodsB.noheath.jdis <- vegdist(arthropodsB.noheath,method="jaccard",binary=TRUE) # binary for presence/absence arthropods.noheath.qjdis <- vegdist(arthropods.noheath,method="jaccard") # quantitative mantel(arthropods.noheath.qjdis,arthropodsB.noheath.jdis) # highly significantly correlated p = 0.001 (quantitative vs binary) # Using the dissimilarity matrices made previously for NMDS (binary jaccard) arthropods.noheath.qjmds <- metaMDS(arthropods.noheath.qjdis, wascores=FALSE, autotransform=FALSE,noshare=FALSE) arthropods.noheath.qjmds <- metaMDS(arthropods.noheath.qjdis, wascores=FALSE, autotransform=FALSE,noshare=FALSE, previous = arthropods.noheath.qjmds) plot(arthropods.noheath.qjmds,display="sites",type="n") with(arthropods.noheath.env,ordisurf(arthropods.noheath.qjmds,avpctbare,col="red",main="Pedley")) with(arthropods.noheath.env,ordiellipse(arthropods.noheath.qjmds, treatment,cex=.5, draw="polygon",col=c("black"),alpha=100,kind="se",conf=0.95,label=TRUE, show.groups=(c("1Control")))) with(arthropods.noheath.env,ordiellipse(arthropods.noheath.qjmds, treatment,cex=.5, draw="polygon",col=c("blue"),alpha=10,kind="se",conf=0.95,label=FALSE, show.groups=(c("ForageHarvest","ForestPlough","Disc", "Swipe")))) with(arthropods.noheath.env,ordiellipse(arthropods.noheath.qjmds, treatment,cex=.5, draw="polygon",col=c("green"),alpha=50,kind="se",conf=0.95,label=TRUE, show.groups=(c("AgriPlough","TurfStrip")))) arthropodsB.noheath.jmds <- metaMDS(arthropodsB.noheath.jdis, wascores=FALSE, autotransform=FALSE,noshare=FALSE) arthropodsB.noheath.jmds <- metaMDS(arthropodsB.noheath.jdis, previous = arthropodsB.noheath.jmds) plot(arthropodsB.noheath.jmds,display="sites",type="n") with(arthropods.noheath.env,ordisurf(arthropodsB.noheath.jmds,avpctbare,col="red",main="MBC")) with(arthropods.noheath.env,ordiellipse(arthropodsB.noheath.jmds, treatment,cex=.5, draw="polygon",col=c("black"),alpha=100,kind="se",conf=0.95,label=TRUE, show.groups=(c("1Control")))) with(arthropods.noheath.env,ordiellipse(arthropodsB.noheath.jmds, treatment,cex=.5, draw="polygon",col=c("blue"),alpha=10,kind="se",conf=0.95,label=FALSE, show.groups=(c("Disc", "Swipe","ForageHarvest")))) with(arthropods.noheath.env,ordiellipse(arthropodsB.noheath.jmds, treatment,cex=.5, draw="polygon",col=c("green"),alpha=50,kind="se",conf=0.95,label=TRUE, show.groups=(c("AgriPlough","ForestPlough","TurfStrip")))) # Procrustes test proarthtest <- protest(arthropods.noheath.qjmds,arthropodsB.noheath.jmds,symmetric=TRUE) plot(proarthtest) proarthtest # highly significantly correlated, p=0.001 # If you have trouble installing mvabund from CRAN or want to install the development version, # scroll down to end of this file: ##### mvabund with Heath sites##### arthmvabund<-mvabund(arthropods) is.mvabund(arthmvabund) arthmvabundB<-mvabund(arthropodsB) is.mvabund(arthmvabundB) meanvar.plot(arthmvabund~arthropods.env$avpctbare,col=as.numeric(arthropods.env$treatment)) meanvar.plot(arthmvabundB~arthropods.env$avpctbare,col=as.numeric(arthropods.env$treatment)) plot(arthmvabund ~ arthropods.env$treatment) # predictor must be a factor plot(arthmvabundB ~ arthropods.env$treatment) # predictor must be a factor arth.nb <- manyglm(arthmvabund~arthropods.env$treatment, family="negative.binomial") plot(arth.nb, res.type="pit.norm") # residuals vs fitted, pretty good plot(arth.nb, res.type="pit.uniform") # residuals vs fitted, pretty good arthmvabund.summ.pit <- summary(arth.nb,test="wald",p.uni="adjusted",resamp="pit.trap", nBoot=999, show.time=TRUE) # nBoot=499 requires ~51 secs, time scales linearly with nBoot arthmvabund.summ.pit$coefficients # Some treatments significantly diff from 1Control pvalues=arthmvabund.summ.pit$coefficients[c(1:8),2] pvalues.corr.fdr<-p.adjust(pvalues, method = "fdr", n = length(pvalues)) pvalues.corr.fdr arthB.nb <- manyglm(arthmvabundB~arthropods.env$treatment, family="binomial") plot(arthB.nb, res.type="pit.norm") # residuals vs fitted, pretty good plot(arthB.nb, res.type="pit.uniform") # residuals vs fitted, pretty good arthmvabundB.summ5.pit <-summary(arthB.nb,test="wald",p.uni="adjusted",resamp="pit.trap", nBoot=999, show.time=TRUE) # nBoot=499 requires ~42 secs, time scales linearly with nBoot arthmvabundB.summ5.pit$coefficients # Some treatments are significantly diff from 1Control pvaluesB=arthmvabundB.summ5.pit$coefficients[c(1:8),2] pvaluesB.corr.fdr<-p.adjust(pvaluesB, method = "fdr", n = length(pvaluesB)) pvaluesB.corr.fdr ###### Alpha species diversity with Heath using incidence coverage estimates ########## (pool_arth<- specpool(arthropods.all)) # incidence based estimate (pool_MBC<- specpool(arthropodsB.all)) # incidence based estimate (pool_arth_by_treatment<- specpool(arthropods.all,treatment)) # incidence based estimate (pool_MBC_by_treatment<- specpool(arthropodsB.all,treatment)) # incidence based estimate # Installing mvabund from source # 1) and 2) are done only once per Mac OS installation # 1) Command Line Tools for Xcode (gcc and other simple stuff instead of all of Xcode) # https://developer.apple.com/downloads/index.action# # Install (automatically done with installer package) # 2) GNU Scientific Library # http://www.gnu.org/software/gsl # in Terminal: # > cd gsl-1.15 # > ./configure --disable-shared --disable-dependency-tracking # command for Mac OS X. See readme for other OSs # > make clean # > make # > sudo make install # 3) mvabund # http://r-forge.r-project.org/R/?group_id=1399 # download the Linux source code: (.tar.gz), not the Mac OS code # the downloaded code automatically decompresses. discard this # in Terminal: # cd to your downloads folder # > R --arch x86_64 CMD INSTALL mvabund_3.6.8.tar.gz # or whatever is the latest version # mvabund will show up in the packages list in R or RStudio