################################################################################## ## R code for analyses reported in Buckley, H. L., T. E. Miller, A. M. Ellison, ## ## and N. J. Gotelli. 2010. Local to continental-scale variation in the richness## ## and composition of an aquatic food web. Global Ecology and Biogeography ## ################################################################################## #setwd("") to your directory and uncomment library(MASS) library(labdsv) library(vegan) library(BRugs) library(arm) # read in raw pa data for each taxon # row numbers are pitchseq inv.raw <- read.table("invpa.txt",header=T,sep="\t") pro.raw <- read.table("propa.txt",header=T,sep="\t") bac.raw <- read.table("bacpa.txt",header=T,sep="\t") # read in raw ab data for each taxon inv.ab <- read.table("invab.txt",header=T,sep="\t") pro.ab <- read.table("proab.txt",header=T,sep="\t") bac.ab <- read.table("bacab.txt",header=T,sep="\t") # read in environmental data site <- read.table("siteenv.txt",header=T,sep="\t") pitcher <- read.table("pitchenv.txt",header=T,sep="\t") # read in species pitcher richness data richness <- read.table("pitchrich.txt",header=T,sep="\t") # read in PCA results pca.pitcher <- read.table("PCA pitcher morphology.txt",header=T) pca.plant <- read.table("PCA plant morphology.txt",header=T) pca.pp <- read.table("PCA pitcher and plant morphology.txt",header=T) pca.vege <- read.table("PCA vegetation.txt",header=T) pca.site <- read.table("PCA site.txt",header=T) # merge site vars with pitcher to get repeated at pitcher level allvars <- merge(pitcher,site,by="site",all.x=TRUE) allvars <- merge(allvars,pca.site,by="site",all.x=TRUE) allvars <- allvars[order(allvars$siteno.x),] #write.table(allvars,"F://Hannah//sarracenia inquilines//August 2008//allvars.csv",sep=",",row.names=FALSE) # calculate distance matrices inv.dis <- dsvdis(inv.raw,index="steinhaus") pro.dis <- dsvdis(pro.raw,index="steinhaus") bac.dis <- dsvdis(bac.raw,index="steinhaus") # perform a principal coordinates analysis inv.pco <- pco(inv.dis,k=30) pro.pco <- pco(pro.dis,k=30) bac.pco <- pco(bac.dis,k=30) # extract scores for the first two axes i1<-sqrt(inv.pco$points[,1]+1) i2<-sqrt(inv.pco$points[,2]+1) p1<-sqrt(pro.pco$points[,1]+1) p2<-sqrt(pro.pco$points[,2]+1) b1<-sqrt(bac.pco$points[,1]+1) b2<-sqrt(bac.pco$points[,2]+1) # plot pco for species patterns plot(inv.pco) surf(inv.pco,inv.raw$MOSQUI>0,family=binomial) # calculate proportion of variance explained (inv.pco$eig[1]+inv.pco$eig[2])/sum(inv.pco$eig) (pro.pco$eig[1]+pro.pco$eig[2])/sum(pro.pco$eig) (bac.pco$eig[1]+bac.pco$eig[2])/sum(bac.pco$eig) # perform a principal components analysis summary(inv.pca <- pca(inv.raw)) summary(pro.pca <- pca(pro.raw)) summary(bac.pca <- pca(bac.raw)) # define variables veg1 <- as.vector(pca.vege$X1) veg2 <- as.vector(pca.vege$X2) pp1 <- as.vector(pca.pp$X1) pp2 <- as.vector(pca.pp$X2) pp3 <- as.vector(pca.pp$X3) # additional environmental variables and site codes # pitcher level age <- as.vector(pitcher$LFYR) standing <- ifelse(pitcher$stwater>0,1,0) mosquito <- log(pitcher$MOSQUI+1) mosquito.pa <- as.vector(ifelse(pitcher$MOSQUI>0,1,0)) midge <- log(pitcher$CHIRON+1) midge.pa <- as.vector(ifelse(pitcher$CHIRON>0,1,0)) mite <- log(pitcher$MITESS+1) mite.pa <- as.vector(ifelse(pitcher$MITESS>0,1,0)) rothab <- log(pitcher$ROTHAB+1) rothab.pa <- as.vector(ifelse(pitcher$ROTHAB>0,1,0)) fletch.pa <- as.vector(ifelse(inv.ab$BLAESO>0,1,0)) bodo.pa <- as.vector(ifelse(pro.ab$BODSP1>0,1,0)) total.s <- richness$ts inv.s <- log(richness$is+1) pro.s <- log(richness$ps+1) bac.s <- log(richness$bs+1) # site level bog.area <- allvars$log_area site1 <- allvars$ax1 site2 <- allvars$ax2 s1 <- pca.site$ax1 s2 <- pca.site$ax2 # other parameters N1 <- 592 N2 <- 30 sno <- as.vector(richness$siteseq) ## Set dependent variable dep <- fletch.pa ## HBM ############################# ## models # unconditional gaussian model sarracenia0.bug <- " model{ for(i in 1:N1) { dep[i] ~ dnorm(pitcher.mu[i],pitcher.tau) pitcher.mu[i] <- a0[sno[i]] } for(j in 1:N2) { a0[j] ~ dnorm(site.mu[j],site.tau) site.mu[j] <- a1 } # Priors pitcher.sigma ~ dunif(0,100) pitcher.tau <- pow(pitcher.sigma,-2) pitcher.sigma2 <- pow(pitcher.sigma,2) site.sigma ~ dunif(0,100) site.tau <- pow(site.sigma,-2) a1 ~ dnorm(0,0.001) } " write(sarracenia0.bug, "sarracenia0.bug") # unconditional binomial model sarracenia0.bug <- " model{ for(i in 1:N1) { dep[i] ~ dbern(pitcher.mu[i]) logit(pitcher.mu[i]) <- a0[sno[i]] } for(j in 1:N2) { a0[j] ~ dnorm(site.mu[j],site.tau) site.mu[j] <- a1 #logit(site.pred.pa[j]) <- a0[j] } # Priors site.sigma ~ dunif(0,100) site.tau <- pow(site.sigma,-2) a1 ~ dnorm(0,0.001) } " write(sarracenia0.bug, "sarracenia0.bug") # conditional models # invertebrates sarracenia1.bug <- " model{ for(i in 1:N1) { dep[i] ~ dnorm(pitcher.mu[i],pitcher.tau) pitcher.mu[i] <- a0[sno[i]] + beta.p[1]*age[i] + beta.p[2]*pp1[i] + beta.p[3]*pp2[i] + beta.p[4]*pp3[i] + beta.p[5]*veg1[i] + beta.p[6]*veg2[i] e.dep[i] <- dep[i] - pitcher.mu[i] # R^2 variance explained } for(j in 1:N2) { a0[j] ~ dnorm(site.mu[j],site.tau) site.mu[j] <- a1 + beta.s[1]*s1[j] + beta.s[2]*s2[j] } # Priors for(i in 1:6){ beta.p[i] ~ dnorm(0,.001) } for(i in 1:2){ beta.s[i] ~ dnorm(0,.001) } pitcher.sigma ~ dunif(0,100) pitcher.tau <- pow(pitcher.sigma,-2) pitcher.sigma2 <- pow(pitcher.sigma,2) site.sigma ~ dunif(0,100) site.tau <- pow(site.sigma,-2) a1 ~ dnorm(0,0.001) } " write(sarracenia1.bug, "sarracenia1.bug") # protozoa and bacteria sarracenia2.bug <- " model{ for(i in 1:N1) { dep[i] ~ dnorm(pitcher.mu[i],pitcher.tau) pitcher.mu[i] <- a0[sno[i]] + beta.p[1]*age[i] + beta.p[2]*pp1[i] + beta.p[3]*pp2[i] + beta.p[4]*pp3[i] + beta.p[5]*veg1[i] + beta.p[6]*veg2[i] + beta.p[7]*mosquito.pa[i] + beta.p[8]*midge.pa[i] e.dep[i] <- dep[i] - pitcher.mu[i] # R^2 variance explained } for(j in 1:N2) { a0[j] ~ dnorm(site.mu[j],site.tau) site.mu[j] <- a1 + beta.s[1]*s1[j] + beta.s[2]*s2[j] } # Priors for(i in 1:8){ beta.p[i] ~ dnorm(0,.001) } for(i in 1:2){ beta.s[i] ~ dnorm(0,.001) } pitcher.sigma ~ dunif(0,100) pitcher.tau <- pow(pitcher.sigma,-2) pitcher.sigma2 <- pow(pitcher.sigma,2) site.sigma ~ dunif(0,100) site.tau <- pow(site.sigma,-2) a1 ~ dnorm(0,0.001) } " write(sarracenia2.bug, "sarracenia2.bug") # Presence of obligates sarracenia3.bug <- " model{ for(i in 1:N1) { dep[i] ~ dbern(pitcher.mu[i]) logit(pitcher.mu[i]) <- a0[sno[i]] + beta.p[1]*age[i] + beta.p[2]*pp1[i] + beta.p[3]*pp2[i] + beta.p[4]*pp3[i] + beta.p[5]*veg1[i] + beta.p[6]*veg2[i] e.dep[i] <- dep[i] - pitcher.mu[i] # R^2 variance explained } for(j in 1:N2) { a0[j] ~ dnorm(site.mu[j],site.tau) site.mu[j] <- a1 + beta.s[1]*s1[j] + beta.s[2]*s2[j] } # Priors for(i in 1:6){ beta.p[i] ~ dnorm(0,.001) } for(i in 1:2){ beta.s[i] ~ dnorm(0,.001) } site.sigma ~ dunif(0,100) site.tau <- pow(site.sigma,-2) a1 ~ dnorm(0,0.001) } " write(sarracenia3.bug, "sarracenia3.bug") ############################################################# ## data files ## unconditional model bugsData(list(dep=dep,sno=sno,N1=N1,N2=N2),"data0.txt") ## conditional models ## pitcher-level predictors: ## pp1,pp2,pp3,age,veg1,veg2,standing,mosquito.pa,midge.pa ## site-level predictors: ## s1,s2,bog.area bugsData(list(dep=dep,sno=sno,N1=N1,N2=N2,pp1=pp1,pp2=pp2,pp3=pp3, age=age,veg1=veg1,veg2=veg2,s1=s1,s2=s2),"data1.txt") # invertebrates bugsData(list(dep=dep,sno=sno,N1=N1,N2=N2,pp1=pp1,pp2=pp2,pp3=pp3, age=age,veg1=veg1,veg2=veg2,s1=s1,s2=s2, mosquito.pa=mosquito.pa, midge.pa=midge.pa),"data2.txt") # protozoa and bacteria ############################################################# ## initial values # unconditional gaussian model bugsData(list(a0=rnorm(N2,0,1), a1=rnorm(1,0,1), pitcher.sigma=rgamma(1,1,1), site.sigma=rgamma(1,1,1) ),"sarracenia0.inits") # unconditional binomial model bugsData(list(a0=rnorm(N2,0,1), a1=rnorm(1,0,1), site.sigma=rgamma(1,2,0.5) ),"sarracenia0.inits") # conditional models bugsData(list(a0=rnorm(N2,0,1), a1=rnorm(1,0,1), pitcher.sigma=rgamma(1,1,1), site.sigma=rgamma(1,1,1), beta.p=rnorm(6,0,1), beta.s=rnorm(2,0,1) ),"sarracenia1.inits") # invertebrate bugsData(list(a0=rnorm(N2,0,1), a1=rnorm(1,0,1), pitcher.sigma=rgamma(1,1,1), site.sigma=rgamma(1,1,1), beta.p=rnorm(8,0,1), beta.s=rnorm(2,0,1) ),"sarracenia2.inits") # protozoa and bacteria bugsData(list(a0=rnorm(N2,0,1), a1=rnorm(1,0,1), site.sigma=rgamma(1,1,1), beta.p=rnorm(6,0,1), beta.s=rnorm(2,0,1) ),"sarracenia3.inits") # obligates ########################################################### ########################################################### ## run the obligate models modelCheck("sarracenia3.bug") modelData("data1.txt") modelCompile(numChains=3) modelInits("sarracenia3.inits",1) modelInits("sarracenia3.inits",2) modelInits("sarracenia3.inits",3) modelGenInits() modelUpdate(10) samplesSet(c("beta.p","beta.s","site.sigma2","a1")) dicSet() modelUpdate(5000) modelUpdate(10000) modelUpdate(30000) modelUpdate(50000) samplesSetBeg(5000) samplesSetBeg(10000) samplesHistory("beta.p") samplesDensity("beta.p") samplesHistory("beta.s") samplesDensity("beta.s") samplesHistory("site.sigma2") samplesHistory("a1") #samplesDensity("*") #samplesAutoC("*",chain=1) #samplesBgr("*") #ss <- samplesStats("*") dic <- dicStats() dic param.p <- samplesStats("beta.p") write.table(param.p,"bodo_full_beta_p.csv",sep=",") param.s <- samplesStats("beta.s") write.table(param.s,"bodo_full_beta_s.csv",sep=",") # param.e.dep <- samplesStats("e.dep") # write.table(param.e.dep,"mosqui_full_e_dep.csv",sep=",") # ## run the invertebrate model modelCheck("sarracenia1.bug") modelData("data1.txt") modelCompile(numChains=3) modelInits("sarracenia1.inits",1) modelInits("sarracenia1.inits",2) modelInits("sarracenia1.inits",3) modelGenInits() modelUpdate(10) samplesSet(c("e.dep","beta.p","beta.s","pitcher.sigma2","site.sigma2","a1")) dicSet() modelUpdate(5000) modelUpdate(10000) modelUpdate(30000) modelUpdate(50000) samplesSetBeg(5000) samplesSetBeg(10000) samplesHistory("beta.p") samplesDensity("beta.p") samplesHistory("beta.s") samplesDensity("beta.s") samplesHistory("pitcher.sigma2") samplesHistory("site.sigma2") samplesHistory("a1") #samplesDensity("*") #samplesAutoC("*",chain=1) #samplesBgr("*") #ss <- samplesStats("*") dic <- dicStats() dic param.p <- samplesStats("beta.p") write.table(param.p,"is_full_beta_p_test.csv",sep=",") param.s <- samplesStats("beta.s") write.table(param.s,"is_full_beta_s_test.csv",sep=",") param.e.dep <- samplesStats("e.dep") write.table(param.e.dep,"is_full_e_dep_test.csv",sep=",") ## run the protozoan or bacterial model modelCheck("sarracenia2.bug") modelData("data2.txt") modelCompile(numChains=3) modelInits("sarracenia2.inits",1) modelInits("sarracenia2.inits",2) modelInits("sarracenia2.inits",3) modelGenInits() modelUpdate(10) samplesSet(c("e.dep","beta.p","beta.s","pitcher.sigma2","site.sigma2","a1")) dicSet() modelUpdate(5000) modelUpdate(10000) modelUpdate(30000) modelUpdate(50000) modelUpdate(60000) samplesSetBeg(5000) samplesHistory("beta.p") samplesDensity("beta.p") samplesHistory("beta.s") samplesDensity("beta.s") samplesHistory("pitcher.sigma2") samplesHistory("site.sigma2") samplesHistory("a1") samplesDensity("*") samplesAutoC("*",chain=1,thin=10) samplesBgr("*") #ss <- samplesStats("*") dic <- dicStats() dic param.p <- samplesStats("beta.p",thin=10) write.table(param.p,"F:\\Hannah\\Sarracenia inquilines\\August 2008\\rothab_full_beta_p.csv",sep=",") param.s <- samplesStats("beta.s",thin=10) write.table(param.s,"F:\\Hannah\\Sarracenia inquilines\\August 2008\\rothab_full_beta_s.csv",sep=",") param.e.dep <- samplesStats("e.dep",thin=10) write.table(param.e.dep,"F:\\Hannah\\Sarracenia inquilines\\August 2008\\rothab_full_e_dep.csv",sep=",") ## Calculate R2 ep2 <- as.data.frame(samplesHistory("e.dep",plot=FALSE,thin=100,firstChain=1,lastChain=1)) ep<-ep2[,1:N1] rsquared.abund <- 1 - mean (apply (ep, 1, var)) / var(dep,na.rm=TRUE) rsquared.abund # = lambda.abund <- 1 - var (apply (ep, 2, mean)) / mean (apply (ep, 1, var)) lambda.abund # = ########################################################### ########################################################### ## run the unconditional gaussian model modelCheck("sarracenia0.bug") modelData("data0.txt") modelCompile(numChains=3) modelInits("sarracenia0.inits",1) modelInits("sarracenia0.inits",2) modelInits("sarracenia0.inits",3) modelGenInits() modelUpdate(10) samplesSet(c("pitcher.sigma2","site.sigma2","a1")) dicSet() modelUpdate(30000) modelUpdate(5000) modelUpdate(10000) samplesSetBeg(5000) samplesHistory("*",thin=10) samplesDensity("*",thin=10) samplesAutoC("*",chain=1,thin=10) samplesBgr("*",thin=10) ss <- samplesStats("*",thin=10) write.table(ss,"b2_uncon_ss.csv",sep=",") ss dic <- dicStats() dic ## run the unconditional binomial model modelCheck("sarracenia0.bug") modelData("data0.txt") modelCompile(numChains=3) modelInits("sarracenia0.inits",1) modelInits("sarracenia0.inits",2) modelInits("sarracenia0.inits",3) modelGenInits() modelUpdate(10) samplesSet(c("site.sigma2","a1")) dicSet() modelUpdate(30000) modelUpdate(5000) modelUpdate(10000) samplesSetBeg(5000) samplesHistory("*",thin=10) samplesDensity("*",thin=10) samplesAutoC("*",chain=1,thin=10) samplesBgr("*",thin=10) ss <- samplesStats("*",thin=10) write.table(ss,"habrot.pa_uncon_ss.csv",sep=",") ss dic <- dicStats() dic # Variance partitioning pitcher.sigma2 <- samplesStats("pitcher.sigma2")$mean site.sigma2 <- samplesStats("site.sigma2")$mean tot <- pitcher.sigma2+site.sigma2 v<-array() v[1] <- pitcher.sigma2/tot v[2] <- site.sigma2/tot signif(v,3) # i1: 0.66 0.34 # i2: 0.80 0.20 # is: 0.70 0.30 # p1: 0.84 0.16 # plot results param <- rbind(param.p,param.s) p <- length(param$mean) # number of parameters xl <- -0.15 xr <- 0.3 pp <- p+1 linew <- 1.4 pl = c("age","pp1","pp2","pp3","veg1","veg2","s1","s2") #pl = c("age","pp1","pp2","pp3","veg1","veg2","mosq","midg","s1","s2") x11(width=9, height=4) par(mar=c(6,13,3,4)) # set plot margins for yaxis labels #mar=c(bottom,left,top,right) default is (5,4,4,2) plot(range(c(xl,xr)),range(c(.8,p+.5)),axes=F,xlab='',ylab='', type='n', main='Arthropod richness') axis(1,tick=T,labels=T,at=seq(xl,xr,by=0.1)) axis(2,at=seq(1,p,by=1),tick=T,labels=pl,las=1) segments(x0=0,y0=0,x1=0,y1=p+.2,lty=2) # center dashed vertical line for (i in 1:p){ segments(x0=param$val2.5pc[i],y0=i,x1=param$val97.5pc[i],y1=i,lwd=linew) points(param$mean[i],i,pch=16,cex=1) } segments(x0=-2,y0=6.5,x1=2,y1=6.5,lty=3,lwd=linew) mtext('parameter values',side=1,line=2.5,cex=1,font=1,at=0) text(1.1,3.75,"pitcher level",adj=c(0,0),cex=1,font=2) text(1.1,7.5,"site level",adj=c(0,0),cex=1,font=2)