####################################################################################################### # # # FUNCTIONAL DIVERSITY MESOCOSM MANIPULATION # # # ####################################################################################################### #Author: Jon Lefcheck #Last updated: 2015-05-18 ####################################################################################################### # TABLE OF CONTENTS # # # # Line 35: Importing and formatting the data # # Line 240: Exploratory plots # # Line 288: Linear mixed effects models # # Line 508: Importance of individual traits # # Line 624: Modeling individual species contributions # # Line 704: Structural equation models # # # ####################################################################################################### library(FD) #Calls: dbFD, gowdis library(fpc) #Calls: pamk library(ggplot2) #Calls: ggplot library(gridExtra) #Calls: grid.arrange library(nlme) #Calls: lme, lmeControl library(piecewiseSEM) #Calls: sem.fit, sem.coefs, sem.model.fits library(plotrix) #Calls: std.error library(plyr) #Calls: ddply library(reshape2) #Calls: melt, dcast #Read in modified dbFD() function from FD package source("dbFD_new.R") ####################################################################################################### # IMPORTING AND FORMATTING THE DATA # ####################################################################################################### #Import initial abundance data fd=read.csv("initial_species_data.csv") #Remove some replicates because of contamination fd=fd[!fd$id %in% c("mono7-1","mono9-5","poly-5"),] #Read in final bug sieve data sieve=read.csv("final_species_data.csv") #Read in final macrophyte and invert AFDM data afdm=read.csv("final_afdm_data.csv") #Import traits traits=read.csv("trait_data.csv") traits=traits[order(traits$species),] rownames(traits)=traits$species; traits=traits[,-1] #Change month of max abundance to ordered factor traits$month.max.abund=ordered(traits$month.max.abund,levels=c("May","July","August","October")) ####################################################################################################### #Estimate initial grazer abundance and initial biomass (g) for predators based on length-weight equations grazers=c("ampithoe","bittium","erichsonella","gammarus","hippolyte") predators=c("callinectes","fundulus","palaemonetes","syngnathus") fd=cbind(fd, initial.grazer.abund=rowSums(fd[,grepl(paste("initial.",grazers,".abund",sep="",collapse="|"),colnames(fd))]), initial.predator.abund=rowSums(fd[,grepl(paste("initial.",predators,".abund",sep="",collapse="|"),colnames(fd))]), #From Cadman & Weinstein 1985 J Crust Biol initial.callinectes.biomass=exp(-11.247 + 2.918*log(fd$initial.callinectes.length)), #From http://www.fishbase.org/summary/Syngnathus-fuscus.html initial.fundulus.biomass=0.00049*(fd$initial.fundulus.length/10)^3.12, #From http://www.fishbase.org/summary/Fundulus-heteroclitus.html initial.syngnathus.biomass=0.01122*(fd$initial.syngnathus.length/10)^3.04) #Estimate initial total predator biomass fd$initial.predator.biomass=rowSums(fd[,grepl("callinectes.biomass|fundulus.biomass|palaemonetes.biomass|syngnathus.biomass",colnames(fd))]) ####################################################################################################### #Calculate final abundance of stocked species abund.final=join_all(list( #Grazers dcast(cbind(sieve,abund=rowSums(sieve[,grepl("mm",colnames(sieve))],na.rm=T)),id~species,fun.aggregate=sum,value.var="abund"), #Predators dcast(ddply(subset(afdm,species %in% c("Blue crab","Mummichog","Pipefish")),c("id","species"),nrow), id~species,value.var="V1"), #Shrimp dcast(ddply(subset(afdm,species %in% c("Hippolyte","Grass shrimp")),c("id","species"),function(x) sum(as.numeric(as.character(x$number.length)))),id~species,value.var="V1") ),by="id") colnames(abund.final)=c("id","ampithoe","barnacle","bittium","bubble.snail","caprella","clam","corophium", "doridella","edotea","erichsonella","gammarus","larval.fish","melita","rhithropanopeus", "geukensia","nereid","spionid","tunicates","unk.amphipod","callinectes", "fundulus","syngnathus","palaemonetes","hippolyte") colnames(abund.final)[-1]=paste("final.",colnames(abund.final)[-1],".abund",sep="") #Summarize final lengths of predators lengths.final=dcast(subset(afdm,species %in% c("Blue crab","Mummichog","Pipefish")),id~species,value.var="number.length") colnames(lengths.final)=c("id","final.callinectes.length","final.fundulus.length","final.syngnathus.length") #Convert all lengths to numeric lengths.final[,grepl(".length",colnames(lengths.final))]=apply(lengths.final[,grepl(".length",colnames(lengths.final))],2,function(x) as.numeric(as.character(x))) ####################################################################################################### #Prepare final biomass of ashed species afdm[afdm$species=="(2) Gracilaria","species"]="Gracilaria" afdm.cast=dcast(afdm,id~species,fun.aggregate=sum,value.var="AFDM") colnames(afdm.cast)=c("id","barnacle","callinectes","bubble.snail","gracilaria","palaemonetes","green.algae", "hippolyte","fundulus","syngnathus","red.algae","tunicate") #Set estimated mg AFDM individual-1 from Edgar 1990 JEMBE multiples= outer(sieve$group=="Crustacean",c(37.007,14.659,5.807,2.300,0.911,0.361,0.143,0.058,0.023))+ outer(sieve$group=="Mollusc",c(47.015,18.076,6.950,2.672,1.027,0.395,0.152,0.060,0.023))+ outer(sieve$group=="Polychaete",c(27.959,11.780,4.963,2.091,0.881,0.371,0.156,0.067,0.028))+ outer(sieve$group=="Caprellid",c(4.382,2.267,1.173,0.607,0.314,0.163,0.084,0.044,0.023)) #Multiply by abundances in each sieve class, sum by row, and divide by 1000 to get g AFDM sieve=sieve[!sieve$group=="",] edgar.biomass=cbind.data.frame(id=sieve$id,species=sieve$species,biomass=rowSums(sieve[,grepl("mm",colnames(sieve))]*multiples,na.rm=T)/1e3) #Cast longways edgar.biomass=dcast(edgar.biomass,id~species,fun.aggregate=sum,value.var="biomass") #Rename columns colnames(edgar.biomass)=c("id","ampithoe","bittium","bubble.snail","caprella","clam","corophium","doridella","edotea", "erichsonella","gammarus","melita","geukensia","nereid","spionid","unk.amphipod") #Merge biomass datasets and sum identical colnames biomass.final=rbind(melt(afdm.cast,id.vars=c(1),measure.vars=c(2:ncol(afdm.cast))),melt(edgar.biomass,id.vars=c(1),measure.vars=c(2:ncol(edgar.biomass)))) biomass.final=dcast(biomass.final,id~variable,fun.aggregate=sum,value.var="value") #Append colnames colnames(biomass.final)[-1]=paste("final.",colnames(biomass.final)[-1],".biomass",sep="") #Get final biomass of grazers biomass.final$final.grazer.biomass=rowSums(biomass.final[,paste("final.",grazers,".biomass",sep="")]) #Get final biomass of predators biomass.final$final.predator.biomass=rowSums(biomass.final[,paste("final.",predators,".biomass",sep="")]) #Get final biomass of recruiting invertebrates recruit.inverts=c("barnacle","bubble.snail","nereid","tunicate") biomass.final$final.recruit.invert.biomass=rowSums(biomass.final[,paste("final.",recruit.inverts,".biomass",sep="")]) #Get final total biomass of stocked species biomass.final$final.biomass=rowSums(biomass.final[,-1]) #Merge into original dataset fd=join_all(list(fd,abund.final,lengths.final,biomass.final),by="id") #Replace NAs with zeros fd[,grepl("final.",colnames(fd))][is.na(fd[,grepl("final.",colnames(fd))])]=0 #Replace any negative numbers with zero (amount too small to account for measurement error) fd[,grepl("final.",colnames(fd))][fd[,grepl("final.",colnames(fd))]<0]=0 ####################################################################################################### #Calculate species & functional richness #Initial richness stocked=c("ampithoe","bittium","callinectes","erichsonella","fundulus","gammarus","hippolyte","palaemonetes","syngnathus") fd$initial.richness=rowSums(fd[,paste("initial.",stocked,".abund",sep="")]>0) #Initial grazer richness fd$initial.grazer.richness=rowSums(fd[,paste("initial.",grazers,".abund",sep="")]>0) #Initial predator richness fd$initial.predator.richness=rowSums(fd[,paste("initial.",predators,".abund",sep="")]>0) #Calculate initial FD for all species all.initial.abund=fd[,paste("initial.",stocked,".abund",sep="")] colnames(all.initial.abund)=rownames(traits) all.initial.FD=do.call(data.frame,dbFD(traits,all.initial.abund,stand.FRic=T,calc.CWM=F,messages=F))[,"FRic",drop=F] colnames(all.initial.FD)=paste("initial.",colnames(all.initial.FD),sep="") all.initial.FD=cbind.data.frame(id=fd[match(rownames(all.initial.FD),rownames(fd)),"id"],all.initial.FD) #Initial grazer FD grazers.initial.abund=fd[,paste("initial.",grazers,".abund",sep="")] colnames(grazers.initial.abund)=rownames(traits)[pmatch(gsub("\\b([a-z])([a-z]+)","\\U\\1\\L\\2",substr(grazers,1,nchar(grazers)-1),perl=TRUE),rownames(traits))] grazers.initial.abund=grazers.initial.abund[rowSums(grazers.initial.abund)!=0,] grazers.initial.FD=do.call(data.frame,dbFD(subset(traits,rownames(traits) %in% colnames(grazers.initial.abund)),grazers.initial.abund,stand.FRic=T,calc.CWM=F,messages=F))[,"FRic",drop=F] colnames(grazers.initial.FD)=paste("initial.grazer.",colnames(grazers.initial.FD),sep="") grazers.initial.FD=cbind.data.frame(id=fd[match(rownames(grazers.initial.FD),rownames(fd)),"id"],grazers.initial.FD) #Initial predator FD predators.initial.abund=fd[,paste("initial.",predators,".abund",sep="")] colnames(predators.initial.abund)=rownames(traits)[pmatch(gsub("\\b([a-z])([a-z]+)", "\\U\\1\\L\\2",substr(predators,1,nchar(predators)-1),perl=TRUE),rownames(traits))] predators.initial.abund=predators.initial.abund[rowSums(predators.initial.abund)!=0,] predators.initial.FD=do.call(data.frame,dbFD(subset(traits,rownames(traits) %in% colnames(predators.initial.abund)),predators.initial.abund,stand.FRic=T,calc.CWM=F,messages=F))[,"FRic",drop=F] colnames(predators.initial.FD)=paste("initial.predator.",colnames(predators.initial.FD),sep="") predators.initial.FD=cbind.data.frame(id=fd[match(rownames(predators.initial.FD),rownames(fd)),"id"],predators.initial.FD) ####################################################################################################### #Final richness of stocked species fd$final.richness=rowSums(fd[,paste("final.",stocked,".abund",sep="")]>0) #Final grazer richness fd$final.grazer.richness=rowSums(fd[,paste("final.",grazers,".abund",sep="")]>0) #Final predator richness fd$final.predator.richness=rowSums(fd[,paste("final.",predators,".abund",sep="")]>0) #Final recruiting invert richness fd$final.recruit.invert.richness=rowSums(ddply(fd,"id",function(x) { #Subset on recruiting inverts x=x[,grepl(paste("final.",recruit.inverts,collapse="|",sep=""),colnames(x))] #Remove everything but species name colnames(x)=gsub("final.|.abund|.biomass","",colnames(x)) sapply(unique(colnames(x)),function(y) sum(x[,y]) ) } )[,-1]>0) #Final total richness fd$final.total.richness=rowSums(ddply(fd,"id",function(x) { x=x[,grepl("final.*.abund|final.*.biomass",colnames(x))] colnames(x)=gsub("final.|.abund|.biomass","",colnames(x)) x2=sapply(unique(colnames(x)),function(y) sum(x[,y]) ) x2[!grepl("grazer|predator|biomass|recruit.invert",names(x2))] } )[,-1]>0) #Final FD for all species all.final.abund=fd[,paste("final.",stocked,".abund",sep="")] colnames(all.final.abund)=rownames(traits) all.final.abund=all.final.abund[rowSums(all.final.abund)!=0 & !is.na(rowSums(all.final.abund)),] all.final.FD=do.call(data.frame,dbFD(traits,all.final.abund,stand.FRic=T,calc.CWM=F,messages=F))[,"FRic",drop=F] colnames(all.final.FD)=paste("final.",colnames(all.final.FD),sep="") all.final.FD=cbind.data.frame(id=fd[match(rownames(all.final.FD),rownames(fd)),"id"],all.final.FD) #Final grazer FD grazers.final.abund=fd[,paste("final.",grazers,".abund",sep="")] colnames(grazers.final.abund)=rownames(traits)[pmatch(gsub("\\b([a-z])([a-z]+)", "\\U\\1\\L\\2",substr(grazers,1,nchar(grazers)-1),perl=TRUE),rownames(traits))] grazers.final.abund=grazers.final.abund[rowSums(grazers.final.abund)!=0 & !is.na(rowSums(grazers.final.abund)),] grazers.final.FD=do.call(data.frame,dbFD(subset(traits,rownames(traits) %in% colnames(grazers.final.abund)),grazers.final.abund,stand.FRic=T,calc.CWM=F,messages=F))[,"FRic",drop=F] colnames(grazers.final.FD)=paste("final.grazer.",colnames(grazers.final.FD),sep="") grazers.final.FD=cbind.data.frame(id=fd[match(rownames(grazers.final.FD),rownames(fd)),"id"],grazers.final.FD) #Final predator FD predators.final.abund=fd[,paste("final.",predators,".abund",sep="")] colnames(predators.final.abund)=rownames(traits)[pmatch(gsub("\\b([a-z])([a-z]+)", "\\U\\1\\L\\2",substr(predators,1,nchar(predators)-1),perl=TRUE),rownames(traits))] predators.final.abund=predators.final.abund[rowSums(predators.final.abund)!=0 & !is.na(rowSums(predators.final.abund)),] predators.final.FD=do.call(data.frame,dbFD(subset(traits,rownames(traits) %in% colnames(predators.final.abund)),predators.final.abund,stand.FRic=T,calc.CWM=F,messages=F))[,"FRic",drop=F] colnames(predators.final.FD)=paste("final.predator.",colnames(predators.final.FD),sep="") predators.final.FD=cbind.data.frame(id=fd[match(rownames(predators.final.FD),rownames(fd)),"id"],predators.final.FD) ####################################################################################################### #Merge data into master datafile fd.final=join_all(list(fd,all.initial.FD,grazers.initial.FD,predators.initial.FD,all.final.FD, grazers.final.FD,predators.final.FD),by="id") #Replace NAs with zeros fd.final[,grepl("initial.|final.",colnames(fd.final))][is.na(fd.final[,grepl("initial.|final.",colnames(fd.final))])]=0 #Calculate summary responses fd.final$final.algal.biomass=fd.final$final.green.algae.biomass+fd.final$final.red.algae.biomass fd.final$final.plant.biomass=fd.final$final.gracilaria.biomass+fd.final$final.green.algae.biomass+fd.final$final.red.algae.biomass fd.final$scaled.initial.richness=scale(fd.final$initial.richness) fd.final$scaled.initial.FRic=scale(fd.final$initial.FRic) ####################################################################################################### # EXPLORATORY PLOTS # ####################################################################################################### #Calculate Gower distance matrix traits.gowdis=gowdis(traits,ord="podani") #Perform multidimensional scaling on distance matrix traits.pcoa=pcoa(traits.gowdis) #Bind into dataframe traits.plot.df=data.frame( traits.pcoa$vectors[,1:2], species=rownames(traits.pcoa$vectors) ) #Get the convex hulls of each unique point set #hulls.df=ddply(traits.plot.df,"groups",function(x) x[chull(x$MDS1,x$MDS2),]) hulls.df=data.frame(traits.pcoa$vectors[chull(traits.pcoa$vectors[,1],traits.pcoa$vectors[,2]),]) #Rearrange position of text labels manually text.position.df=data.frame( Axis.1=c(-0.2,-0.15,0.11,-0.21,0.62,-0.24,-0.25,-0.08,0.4), Axis.2=c(0.17,-0.35,0.-0.08,0.115,-0.1,0.26,-0.155,0.01,0.22), species=c("Ampithoid","Bittiolum","Callinectes","Erichsonella","Fundulus","Gammarus","Hippolyte","Palaemonetes","Syngnathus")) #Plot results ggplot()+ geom_polygon(data=hulls.df,aes(x=Axis.1,y=Axis.2),alpha=0.5,lwd=0.8,fill="grey90",col="grey70")+ geom_point(data=traits.plot.df,aes(x=Axis.1,y=Axis.2),size=4,col="black")+ geom_text(data=text.position.df,aes(x=Axis.1,y=Axis.2,label=species),size=6)+ coord_cartesian(xlim=c(-0.45,0.75),ylim=c(-0.4,0.3))+ labs(x="PCO1",y="PCO2")+ theme_bw(base_size=18)+ theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank()) ####################################################################################################### #Look at correlations between SR and FD ggplot(fd.final,aes(x=initial.richness,y=initial.FRic))+ geom_point(size=2.5,alpha=0.75,col="grey50")+ stat_smooth(method="lm",se=F,lwd=1,col="black")+ geom_text(data=data.frame( label=paste("r = ",round(cor(fd.final$initial.richness,fd.final$initial.FRic),2),sep="")), aes(x=8,y=0.1,label=label),size=8)+ labs(x="Initial species richness",y="Initial functional richness")+ scale_x_continuous(breaks=c(1,3,6,9))+ theme_bw(base_size=18)+ theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank()) ####################################################################################################### # LINEAR MIXED EFFECTS MODELS # ####################################################################################################### #Fit models against only richness richmodels.list=list( recruit.invert.biomass=lme(final.recruit.invert.biomass~scaled.initial.richness, random=~1|block, weights=varIdent(form=~1|scaled.initial.richness), method="ML", data=fd.final), grazer.biomass=lme(final.grazer.biomass~scaled.initial.richness+initial.grazer.abund, random=~1|block, weights=varIdent(form=~1|scaled.initial.richness), method="ML", data=fd.final), predator.biomass=lme(final.predator.biomass~scaled.initial.richness+initial.predator.biomass, random=~1|block, weights=varIdent(form=~1|scaled.initial.richness), method="ML", data=fd.final), algal.biomass=lme(final.algal.biomass~scaled.initial.richness, random=~1|block, method="ML", data=fd.final), gracilaria.biomass=lme(final.gracilaria.biomass~scaled.initial.richness, random=~1|block, method="ML", data=fd.final), final.richness=lme(final.richness~scaled.initial.richness, random=~1|block, method="ML", data=fd.final), final.FD=lme(final.FRic~scaled.initial.richness, random=~1|block, method="ML", data=fd.final) ) #Fit models against only FD fdmodels.list=list( recruit.invert.biomass=lme(final.recruit.invert.biomass~scaled.initial.FRic, random=~1|block, weights=varIdent(form=~1|scaled.initial.richness), method="ML", data=fd.final), grazer.biomass=lme(final.grazer.biomass~scaled.initial.FRic+initial.grazer.abund, random=~1|block, weights=varIdent(form=~1|scaled.initial.richness), method="ML", control=lmeControl(opt="optim"), data=fd.final), predator.biomass=lme(final.predator.biomass~scaled.initial.FRic+initial.predator.biomass, random=~1|block, weights=varIdent(form=~1|scaled.initial.richness), method="ML", data=fd.final), algal.biomass=lme(final.algal.biomass~scaled.initial.FRic, random=~1|block, method="ML", data=fd.final), gracilaria.biomass=lme(final.gracilaria.biomass~scaled.initial.FRic, random=~1|block, method="ML", data=fd.final), final.richness=lme(final.richness~scaled.initial.FRic, random=~1|block, method="ML", data=fd.final), final.FD=lme(final.FRic~scaled.initial.FRic, random=~1|block, method="ML", data=fd.final) ) #Calculate AIC and conditional and marginal R2 for each list sem.model.fits(richmodels.list) sem.model.fits(fdmodels.list) #Also run interaction models (for supplements) intmodels.list=list( recruit.invert.biomass=lme(final.recruit.invert.biomass~scaled.initial.richness*scaled.initial.FRic, random=~1|block, weights=varIdent(form=~1|scaled.initial.richness), data=fd.final), grazer.biomass=lme(final.grazer.biomass~scaled.initial.richness*scaled.initial.FRic+initial.predator.biomass, random=~1|block, weights=varIdent(form=~1|scaled.initial.richness), data=fd.final), predator.biomass=lme(final.predator.biomass~scaled.initial.richness*scaled.initial.FRic+initial.predator.biomass, random=~1|block, weights=varIdent(form=~1|scaled.initial.richness), data=fd.final), algal.biomass=lme(final.algal.biomass~scaled.initial.richness*scaled.initial.FRic, random=~1|block, weights=varIdent(form=~1|scaled.initial.richness), control=lmeControl(opt="optim"), data=fd.final), gracilaria.biomass=lme(final.gracilaria.biomass~scaled.initial.richness*scaled.initial.FRic, random=~1|block, weights=varIdent(form=~1|scaled.initial.richness), data=fd.final), final.richness=lme(final.richness~scaled.initial.richness*scaled.initial.FRic, random=~1|block, weights=varIdent(form=~1|scaled.initial.richness), data=fd.final), final.FD=lme(final.FRic~scaled.initial.richness*scaled.initial.FRic, random=~1|block, weights=varIdent(form=~1|scaled.initial.richness), data=fd.final) ) #Examine fit statistics sem.model.fits(intmodels.list) #Extract coefficients (if significant) ldply(intmodels.list,function(i) { coefs=ifelse(summary(i)$tTable[-1,5]<=0.05,round(summary(i)$tTable[-1,1],2),"NS") t(data.frame(coefs)) } ) #Get predicted fit across richness levels newdata=ldply(intmodels.list,function(i) { i=update(i,fixed=.~scaled.initial.FRic,control=lmeControl(opt="optim")) newdata=expand.grid( scaled.initial.FRic=seq(range(i$data$scaled.initial.FRic)[1], range(i$data$scaled.initial.FRic)[2],0.01), initial.predator.biomass=mean(i$data[i$data$initial.richness==3,"initial.predator.biomass"],na.rm=T) ) newdata$response.name=as.character(formula(i)[[2]]) newdata$response=predict(i,newdata,level=0) return(newdata) } ) #Get predicted fits by richness level (3 or 6) newdata.richness=ldply(intmodels.list,function(i) { newdata=rbind( expand.grid( scaled.initial.richness=unique(i$data[i$data$initial.richness==3,"scaled.initial.richness"]), scaled.initial.FRic=seq(range(i$data[i$data$initial.richness==3,"scaled.initial.FRic"])[1], range(i$data[i$data$initial.richness==3,"scaled.initial.FRic"])[2],0.01), initial.predator.biomass=mean(i$data[i$data$initial.richness==3,"initial.predator.biomass"],na.rm=T)), expand.grid( scaled.initial.richness=unique(i$data[i$data$initial.richness==6,"scaled.initial.richness"]), scaled.initial.FRic=seq(range(i$data[i$data$initial.richness==6,"scaled.initial.FRic"])[1], range(i$data[i$data$initial.richness==6,"scaled.initial.FRic"])[2],0.01), initial.predator.biomass=mean(i$data[i$data$initial.richness==6,"initial.predator.biomass"],na.rm=T)) ) newdata$response.name=as.character(formula(i)[[2]]) newdata$response=predict(i,newdata,level=0) return(newdata) } ) #Melt dataframe for plotting fd.melt=fd.final[,colnames(fd.final) %in% c("scaled.initial.FRic","scaled.initial.richness",unique(newdata$response.name))] fd.melt=melt(fd.melt,id.vars=c(ncol(fd.melt)-1,ncol(fd.melt)),measure.vars=c(1:(ncol(fd.melt)-2))) names(fd.melt)[3:4]=c("response.name","response") #Plot predicted fits #Reorder for better plotting fd.melt$response.name=factor(fd.melt$response.name,levels=c("final.grazer.biomass","final.recruit.invert.biomass","final.predator.biomass", "final.algal.biomass","final.gracilaria.biomass","final.richness","final.FRic")) newdata$response.name=factor(newdata$response.name,levels=c("final.grazer.biomass","final.recruit.invert.biomass","final.predator.biomass", "final.algal.biomass","final.gracilaria.biomass","final.richness","final.FRic")) newdata.richness$response.name=factor(newdata.richness$response.name,levels=c("final.grazer.biomass","final.recruit.invert.biomass","final.predator.biomass", "final.algal.biomass","final.gracilaria.biomass","final.richness","final.FRic")) #Rename levels to add units levels(fd.melt$response.name)=c("Grazer biomass\n(mg AFDM)","Recruit invert biomass\n(g AFDM)","Predator biomass\n(g AFDM)", "Algal biomass\n(g AFDM)","Gracilaria biomass\n(g AFDM)","Final richness", "Final functional\ndiversity") levels(newdata$response.name)=c("Grazer biomass\n(mg AFDM)","Recruit invert biomass\n(g AFDM)","Predator biomass\n(g AFDM)", "Algal biomass\n(g AFDM)","Gracilaria biomass\n(g AFDM)","Final richness", "Final functional\ndiversity") levels(newdata.richness$response.name)=c("Grazer biomass\n(mg AFDM)","Recruit invert biomass\n(g AFDM)","Predator biomass\n(g AFDM)", "Algal biomass\n(g AFDM)","Gracilaria biomass\n(g AFDM)","Final richness", "Final functional\ndiversity") #And plot without DO & final diversity: 10" x 6.5" ggplot()+ geom_point(data=subset(fd.melt,!response.name %in% c("Final richness","Final functional\ndiversity")), aes(x=scaled.initial.FRic,y=response,shape=factor(scaled.initial.richness,labels=c("1","3","6","9")), col=factor(scaled.initial.richness,labels=c("1","3","6","9"))),size=3,alpha=0.85)+ geom_line(data=subset(newdata.richness,!response.name %in% c("Final richness","Final functional\ndiversity")), aes(x=scaled.initial.FRic,y=response,group=factor(scaled.initial.richness,labels=c("3","6")), shape=factor(scaled.initial.richness,labels=c("3","6")), col=factor(scaled.initial.richness,labels=c("3","6"))),lwd=1.2,lty=5)+ geom_line(data=subset(newdata,!response.name %in% c("Final richness","Final functional\ndiversity")), aes(x=scaled.initial.FRic,y=response),col="black",lwd=1.2)+ geom_text(data=data.frame( response.name=c("Grazer biomass\n(mg AFDM)","Recruit invert biomass\n(g AFDM)","Predator biomass\n(g AFDM)", "Algal biomass\n(g AFDM)","Gracilaria biomass\n(g AFDM)"), labels=letters[1:5]), aes(x=-0.3,y=Inf,label=labels),vjust=1.5,col="black",fontface="bold",size=6)+ scale_shape_manual(values=c(15:17,3),name="Richness")+ scale_color_manual(values=c("grey80","grey70","grey50","grey35"),name="Richness")+ labs(x="Functional richness",y="Ecosystem functioning")+ facet_wrap(~response.name,scales="free_y",nrow=2)+ theme_bw(base_size=18)+ theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(), legend.position=c(0.8,0.25))+ guides(color=guide_legend(override.aes=list(alpha=1,size=4,linetype=0))) #Repeat for final diversity: 10" x 4.5" ggplot()+ geom_point(data=subset(fd.melt,response.name %in% c("Final richness","Final functional\ndiversity")), aes(x=scaled.initial.FRic,y=response,shape=factor(scaled.initial.richness,labels=c("1","3","6","9")), col=factor(scaled.initial.richness,labels=c("1","3","6","9"))),size=3,alpha=0.85)+ geom_line(data=subset(newdata,response.name %in% c("Final richness","Final functional\ndiversity")), aes(x=scaled.initial.FRic,y=response),col="black",lwd=1.2)+ geom_line(data=subset(newdata.richness,response.name %in% c("Final richness","Final functional\ndiversity")), aes(x=scaled.initial.FRic,y=response,group=factor(scaled.initial.richness,labels=c("3","6")), col=factor(scaled.initial.richness,labels=c("3","6"))),lwd=1.2,lty=5,alpha=0.8)+ geom_text(data=data.frame( response.name=c("Final richness","Final functional\ndiversity"), labels=letters[1:2]), aes(x=-0.3,y=Inf,label=labels),vjust=1.5,col="black",fontface="bold",size=6)+ scale_shape_manual(values=c(15:17,3),name="Richness")+ scale_color_manual(values=c("grey80","grey70","grey50","grey35"),name="Richness")+ labs(x="Functional richness",y="Final diversity")+ facet_wrap(~response.name,scales="free_y",nrow=1)+ theme_bw(base_size=18)+ theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank())+ guides(color=guide_legend(override.aes=list(alpha=1,size=4,linetype=0))) ####################################################################################################### # IMPORTANCE OF INDIVIDUAL TRAITS # ####################################################################################################### #Calculate initial FD for each trait for all species all.jackknifed.initial.FD=do.call(data.frame,lapply(colnames(traits),function(i) { x=do.call(data.frame,dbFD(traits[,colnames(traits)!=i,drop=F],all.initial.abund,stand.FRic=T,calc.CWM=F,messages=F))[,"FRic",drop=F] colnames(x)=paste(i,"jackknife.all",colnames(x),sep=".") x[is.na(x)]=0 return(x) } ) ) all.jackknifed.initial.FD=cbind.data.frame(id=fd[match(rownames(all.jackknifed.initial.FD),rownames(fd)),"id"],all.jackknifed.initial.FD) #Bind into final data fd.final=join_all(list(fd.final,all.jackknifed.initial.FD),by="id") #Regress each response against multivariate FD and jack-knifed FD jackknifed.AIC=do.call(data.frame,lapply(c("initial.FRic",names(all.jackknifed.initial.FD)[-1]),function(i) { #Replace predictors with jackknifed FD indices models.list=lapply(fdmodels.list,function(j) update(j,fixed=formula(paste(formula(j)[[2]],"~",i))) ) #Extract AIC scores for each model AIC.scores=ldply(models.list,AIC) rownames(AIC.scores)=AIC.scores[,1]; AIC.scores=AIC.scores[,-1,drop=F] colnames(AIC.scores)=gsub(".jackknife.all.FRic","",i) return(AIC.scores) } ) ) #Calculate delta AIC scores for each response for(i in 1:nrow(jackknifed.AIC)) jackknifed.AIC[i,]=round(jackknifed.AIC[i,]-jackknifed.AIC[i,1],2) #Look at results jackknifed.AIC[,-1] ####################################################################################################### #Calculate initial FD for each trait for all species all.individual.initial.FD=do.call(data.frame,lapply(colnames(traits),function(i) { x=do.call(data.frame,dbFD(traits[,colnames(traits)==i,drop=F],all.initial.abund,stand.FRic=T,calc.CWM=F,messages=F))[,"FRic",drop=F] colnames(x)=paste(i,"all",colnames(x),sep=".") return(x) } ) ) all.individual.initial.FD=cbind.data.frame(id=fd[match(rownames(all.individual.initial.FD),rownames(fd)),"id"],all.individual.initial.FD) #Bind into final data fd.final=join_all(list(fd.final,all.individual.initial.FD),by="id") #Regress each response against multivariate and univariate FD individual.coefs=ldply(c("initial.FRic",names(all.individual.initial.FD)[-1]),function(i) { #Replace predictors with individual traits indices models.list=lapply(fdmodels.list,function(j) update(j,fixed=formula(paste(formula(j)[[2]],"~",i))) ) #Extract regression coefficients for each model coefs.df=do.call(rbind,lapply(models.list,function(j) data.frame(estimate=round(summary(j)$tTable[2,1],2), se=round(summary(j)$tTable[2,2],2), sig=ifelse(summary(j)$tTable[2,5] <= 0.05,"yes","no") ) ) ) return(data.frame(trait=i,response=rownames(coefs.df),coefs.df,row.names=NULL)) } ) #Relevel for better plotting levels(individual.coefs$trait)=c("All","Armor","Body plan","Trophic level","Max. biomass","Mean biomass","Mobility","Reprod. mode","Month max. abund") individual.coefs$trait=factor(individual.coefs$trait,levels=c("All","Armor","Body plan","Trophic level","Max. biomass","Mean biomass","Mobility","Reprod. mode","Month max. abund")) levels(individual.coefs$response)=c( "Algal biomass","Final functional\ndiversity","Final richness","Gracilaria biomass","Grazer biomass", "Predator biomass","Recruit invert biomass") individual.coefs$response=factor(individual.coefs$response,levels=c( "Grazer biomass","Recruit invert biomass","Predator biomass","Algal biomass", "Gracilaria biomass","Final richness","Final functional\ndiversity")) #Plot results: 12.5" x 5" sub.individual.coefs=subset(individual.coefs,!response %in% c("Final richness","Final functional\ndiversity")) sub.individual.coefs$trait=factor(sub.individual.coefs$trait,levels=rev(levels(sub.individual.coefs$trait))) ggplot()+ geom_point(data=subset(sub.individual.coefs,trait!="All"),aes(y=trait,x=estimate,fill=sig),size=0)+ geom_rect(data=subset(sub.individual.coefs,trait=="All"), aes(ymax=Inf,ymin=-Inf,xmax=estimate+2*se,xmin=estimate-2*se), fill="grey90")+ geom_vline(data=subset(sub.individual.coefs,trait=="All"),aes(xintercept=estimate),lwd=0.8)+ geom_errorbarh(data=subset(sub.individual.coefs,trait!="All"),aes(y=trait,x=estimate,xmax=estimate+2*se,xmin=estimate-2*se),height=0.2)+ geom_point(data=subset(sub.individual.coefs,trait!="All"),aes(y=trait,x=estimate,fill=sig),size=5,shape=21)+ geom_text(data=data.frame( response=c("Grazer biomass","Recruit invert biomass","Predator biomass", "Algal biomass","Gracilaria biomass"), labels=letters[1:5]), aes(x=-Inf,y=Inf,label=labels,group=1),vjust=1.5,hjust=-1,col="black",fontface="bold",size=6)+ scale_fill_manual(values=c("black","white"))+ facet_grid(~response,scales="free_x")+ labs(x="Standardized regression coefficient",y="")+ theme_bw(base_size=18)+ theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(), legend.position="none") #Look at how variance changes with for univariate traits: 10" x 8" ggplot()+ geom_hline(data=subset(sub.individual.coefs,trait=="All"),aes(yintercept=se),lwd=0.8)+ geom_point(data=subset(sub.individual.coefs,trait!="All"),aes(x=trait,y=se),size=4)+ geom_text(data=data.frame( response=c("Grazer biomass","Recruit invert biomass","Predator biomass", "Algal biomass","Gracilaria biomass"), labels=letters[1:5], y=ddply(subset(sub.individual.coefs,trait!="All"),c("response"),summarize,max(abs(se)))[,2]+ ddply(subset(sub.individual.coefs,trait!="All"),c("response"),summarize,diff(range(se)))[,2]*0.2), aes(x=-Inf,y=y,label=labels,group=1),hjust=-1,vjust=1,col="black",fontface="bold",size=6)+ facet_wrap(~response,scales="free_y")+ labs(x="",y="Standard error of linear coefficient")+ theme_bw(base_size=18)+ theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(), axis.text.x=element_text(angle=270,hjust=0,vjust=0.5), legend.position="none") #Quantify % change in standard errors ddply(sub.individual.coefs,c("response"),function(x) data.frame(response=unique(x$response), trait=x[x$trait!="All","trait"], se.reduction=1-x[x$trait=="All","se"]/x[x$trait!="All","se"]) ) ####################################################################################################### # MODELING SPECIES-SPECIFIC CONTRIBUTIONS # ####################################################################################################### #Model individual contributions (weighted by presence-absence) #Bind presence-absence data all.initial.presabs=all.initial.abund all.initial.presabs[all.initial.presabs>0]=1 colnames(all.initial.presabs)=paste("initial.",stocked,".presabs",sep="") fd.final=cbind(fd.final,all.initial.presabs) #Define fixed formula as each individual stocked species fixed.formula=paste("initial.",stocked,".presabs",sep="",collapse="+") speciesmods.list=list( recruit.invert.biomass=lme(formula(paste("scale(final.recruit.invert.biomass)~",fixed.formula)),random=~1|block,data=fd.final), grazer.biomass=lme(formula(paste("scale(final.grazer.biomass)~",fixed.formula)),random=~1|block,data=fd.final), predator.biomass=lme(formula(paste("scale(final.predator.biomass)~",fixed.formula)),random=~1|block,data=fd.final), algal.biomass=lme(formula(paste("scale(final.algal.biomass)~",fixed.formula)),random=~1|block,data=fd.final), gracilaria.biomass=lme(formula(paste("scale(final.gracilaria.biomass)~",fixed.formula)),random=~1|block,data=fd.final), final.richness=lme(formula(paste("scale(final.richness)~",fixed.formula)),random=~1|block, data=fd.final), final.FD=lme(formula(paste("scale(final.FRic)~",fixed.formula)),random=~1|block,data=fd.final) ) #Extract coefficients speciescoefs.df=ldply(speciesmods.list,function(x) data.frame( species=gsub("initial.(.*).presabs","\\1",rownames(summary(x)$tTable)[-1]), round(summary(x)$tTable[-1,],2)) ) #Cast species as columns species.coefs.long=dcast(speciescoefs.df,.id~species,value.var="Value") species.coefs.long[c(5:7,1,4),] #Same order as tables in paper #Cast p-values in similar fashion species.coefs.pvalues=dcast(speciescoefs.df,.id~species,value.var="p.value") species.coefs.long.sig=species.coefs.long[c(5:7,1,4),-1]*(species.coefs.pvalues[c(5:7,1,4),-1]<=0.05)*1 species.coefs.long.sig[species.coefs.long.sig==0]="" species.coefs.long.sig ####################################################################################################### #Look at functional distinctness against effect size #Calculate average distance to all other species dist2species=colMeans(as.matrix(traits.gowdis)) names(dist2species)=unique(speciescoefs.df$species) #Average predator distinctness vs average grazer distinctness mean(dist2species[c(3,5,8:9)]); std.error(dist2species[c(3,5,8:9)]) mean(dist2species[c(1:2,4,6:7)]); std.error(dist2species[c(1:2,4,6:7)]) #Merge to coefs dataset speciescoefs.df$pairwise=dist2species[match(speciescoefs.df$species,names(dist2species))] #Reorder for better plotting speciescoefs.df$.id=as.factor(as.character(speciescoefs.df$.id)) levels(speciescoefs.df$.id)=c( "Algal biomass\n(g AFDM)","Final functional\ndiversity","Final richness","Gracilaria biomass\n(g AFDM)", "Grazer biomass\n(mg AFDM)","Predator biomass\n(g AFDM)","Recruit invert biomass\n(g AFDM)") speciescoefs.df$.id=factor(speciescoefs.df$.id,levels=c( "Grazer biomass\n(mg AFDM)","Recruit invert biomass\n(g AFDM)","Predator biomass\n(g AFDM)","Algal biomass\n(g AFDM)", "Gracilaria biomass\n(g AFDM)","Final richness","Final functional\ndiversity")) levels(speciescoefs.df$species)=c("Amp","Bitt","Call","Erich","Fund","Gamm","Hippo","Pal","Syn") #And plot without DO ggplot(subset(speciescoefs.df,!.id %in% c("Final richness","Final functional\ndiversity")), aes(y=Value,x=pairwise,group=.id,shape=species))+ geom_point(size=3,alpha=0.75)+ stat_smooth(method="lm",se=F,lwd=1,col="black")+ facet_wrap(~.id,scales="free",nrow=2)+ scale_shape_manual(values=c(15:18,0:2,5,8),name="")+ #scale_color_manual(values=paste("grey",seq(50,90,5),sep=""),name="")+ labs(x="Average pairwise functional distance",y="Linear coefficient",nrow=2)+ geom_text(data=data.frame( .id=c("Grazer biomass\n(mg AFDM)","Recruit invert biomass\n(g AFDM)","Predator biomass\n(g AFDM)","Algal biomass\n(g AFDM)", "Gracilaria biomass\n(g AFDM)"), species=NA,#unique(speciescoefs.df$species), labels=letters[1:5]), aes(x=-Inf,y=Inf,label=labels),hjust=-1,vjust=1.5,col="black",fontface="bold",size=6)+ theme_bw(base_size=18)+ guides(shape=guide_legend(ncol=2))+ theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank(), legend.position=c(0.8,0.25)) ####################################################################################################### # STRUCTURAL EQUATION MODELS # ####################################################################################################### #Generate piecewise models and store in a list SEM.FDmodels.list=list( final.grazer.fd=lme(final.grazer.FRic~initial.grazer.FRic+final.predator.biomass, #weights=varIdent(form=~1|initial.richness), control=lmeControl(opt="optim"), random=~1|block,data=fd.final), final.predator.fd=lme(final.predator.FRic~initial.predator.FRic, #weights=varIdent(form=~1|initial.richness), control=lmeControl(opt="optim"), random=~1|block,data=fd.final), final.predator.biomass=lme(final.predator.biomass~final.predator.FRic+initial.predator.FRic+ initial.predator.biomass+initial.grazer.FRic, random=~1|block,data=fd.final), final.grazer.biomass=lme(final.grazer.biomass~initial.grazer.abund+initial.predator.FRic+ initial.grazer.FRic+final.grazer.FRic+final.predator.FRic+ initial.predator.biomass+final.predator.biomass, weights=varIdent(form=~1|initial.richness), control=lmeControl(opt="optim",msTol=1e-4), random=~1|block,data=fd.final), final.plant.biomass=lme(final.plant.biomass~final.grazer.biomass+final.predator.biomass, random=~1|block,data=fd.final), final.recruit.invert.biomass=lme(final.recruit.invert.biomass~final.predator.biomass+final.grazer.biomass+ final.grazer.FRic+final.predator.FRic, random=~1|block,data=fd.final) ) #Run Shipley's test of d-separation sem.fit(SEM.FDmodels.list,fd.final,corr.errors=c("initial.grazer.FRic~~initial.predator.FRic")) #Extract coefficients and visually assess model fit lapply(SEM.FDmodels.list,plot) sem.coefs(SEM.FDmodels.list,fd.final,standardize="scale",corr.errors=c("initial.grazer.FRic~~initial.predator.FRic")) #Double check inferences using vcov approach sem.lavaan(SEM.FDmodels.list, fd.final) #Generate piecewise models and store in a list SEM.SRmodels.list=list( final.grazer.richness=lme(final.grazer.richness~initial.grazer.richness+final.predator.biomass, #weights=varIdent(form=~1|initial.richness), control=lmeControl(opt="optim"), random=~1|block,data=fd.final), final.predator.richness=lme(final.predator.richness~initial.predator.richness, #weights=varIdent(form=~1|initial.richness), control=lmeControl(opt="optim"), random=~1|block,data=fd.final), final.predator.biomass=lme(final.predator.biomass~final.predator.richness+initial.predator.richness+ initial.predator.biomass+initial.grazer.richness, random=~1|block,data=fd.final), final.grazer.biomass=lme(final.grazer.biomass~initial.grazer.abund+initial.predator.richness+ initial.grazer.richness+final.grazer.richness+final.predator.richness+ initial.predator.biomass+final.predator.biomass, weights=varIdent(form=~1|initial.richness), control=lmeControl(opt="optim",msTol=1e-4), random=~1|block,data=fd.final), final.plant.biomass=lme(final.plant.biomass~final.grazer.biomass+final.predator.biomass, random=~1|block,data=fd.final), final.recruit.invert.biomass=lme(final.recruit.invert.biomass~final.predator.biomass+final.grazer.biomass+ final.grazer.richness+final.predator.richness, random=~1|block,data=fd.final) ) #Get model fit and AIC sem.fit(SEM.SRmodels.list,fd.final,corr.errors=c("initial.grazer.richness~~initial.predator.richness")) #Examine standardized regression coefficients sem.coefs(SEM.SRmodels.list,fd.final,standardize="scale",corr.errors=c("initial.grazer.richness~~initial.predator.richness")) ####################################################################################################### #Look at partial correlations between grazer FD and grazer biomass gbiomass.resids.df=partial.resid(final.grazer.biomass~final.grazer.FRic,SEM.FDmodels.list, plotit = F) #And plot ggplot(gbiomass.resids.df,aes(x=final.grazer.FRic,y=final.grazer.biomass))+ geom_point(size=2.5,alpha=0.75,col="grey50")+ stat_smooth(method="lm",se=F,lwd=1,col="black")+ scale_color_manual(values=paste("grey",seq(50,90,5),sep=""),name="")+ labs(x="Final grazer functional richness | others",y="Final grazer biomass | others",nrow=2)+ theme_bw(base_size=18)+ guides(col=guide_legend(ncol=2))+ theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank()) #Look at grazer biomass for only replicates that had grazers sem.coefs(list(lme(final.grazer.biomass~initial.grazer.abund+initial.predator.FRic+ initial.grazer.FRic+final.grazer.FRic+final.predator.FRic+ initial.predator.biomass+final.predator.biomass, control=lmeControl(opt="optim",msTol=1e-4), random=~1|block,data=subset(fd.final,initial.predator.richness>0))), data=subset(fd.final,initial.predator.richness>0), standardize="scale") #Look at predator biomass for only replicates that had predators sem.coefs(list(lme(final.predator.biomass~final.predator.FRic+initial.predator.FRic+ initial.predator.biomass+initial.grazer.FRic, random=~1|block,data=subset(fd.final,initial.predator.richness>0))), data=subset(fd.final,initial.predator.richness>0), standardize="scale") #Look at bottom-up paths of diversity on predator biomass predbiomass.FD.resids.df=partial.resid(final.predator.biomass~initial.grazer.FRic,SEM.FDmodels.list,plotit=F) colnames(predbiomass.FD.resids.df)=c("response","predictor") predbiomass.SR.resids.df=partial.resid(final.predator.biomass~initial.grazer.richness,SEM.SRmodels.list,plotit=F) colnames(predbiomass.SR.resids.df)=c("response","predictor") ggplot(rbind(data.frame(group="FRic",predbiomass.FD.resids.df),data.frame(group="Richness",predbiomass.SR.resids.df)), aes(x=predictor,y=response))+ geom_point()+ geom_point(size=2.5,alpha=0.75,col="grey50")+ stat_smooth(method="lm",se=F,lwd=1,col="black")+ scale_color_manual(values=paste("grey",seq(50,90,5),sep=""),name="")+ labs(y="Final predator biomass | others",x="Grazer diversity | others",nrow=2)+ facet_grid(~group,scales="free_x")+ theme_bw(base_size=18)+ guides(col=guide_legend(ncol=2))+ theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank()) #Look at partial correlations between grazer FD and grazer biomass gbiomass.resids.df=partial.resid(final.grazer.biomass~final.grazer.FRic,SEM.FDmodels.list, plotit = F) #And plot ggplot(gbiomass.resids.df,aes(x=final.grazer.FRic,y=final.grazer.biomass))+ geom_point(size=2.5,alpha=0.75,col="grey50")+ stat_smooth(method="lm",se=F,lwd=1,col="black")+ scale_color_manual(values=paste("grey",seq(50,90,5),sep=""),name="")+ labs(x="Final grazer functional richness | others",y="Final grazer biomass | others",nrow=2)+ theme_bw(base_size=18)+ guides(col=guide_legend(ncol=2))+ theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank()) #Look at bottom-up paths of diversity on predator biomass predbiomass.FD.resids.df=partial.resid(final.predator.biomass~initial.grazer.FRic,SEM.FDmodels.list,plotit=F) colnames(predbiomass.FD.resids.df)=c("response","predictor") predbiomass.SR.resids.df=partial.resid(final.predator.biomass~initial.grazer.richness,SEM.SRmodels.list,plotit=F) colnames(predbiomass.SR.resids.df)=c("response","predictor") #And plot: 8" x 4.5" ggplot(rbind(data.frame(group="Functional\nrichness",predbiomass.FD.resids.df),data.frame(group="Species\nrichness",predbiomass.SR.resids.df)), aes(x=predictor,y=response))+ geom_point()+ geom_point(size=2.5,alpha=0.75,col="grey50")+ stat_smooth(method="lm",se=F,lwd=1,col="black")+ scale_color_manual(values=paste("grey",seq(50,90,5),sep=""),name="")+ geom_text(data=data.frame( group=c("Functional\nrichness","Species\nrichness"), label=c("a","b")), aes(x=-Inf,y=Inf,label=label),hjust=-1,vjust=1.5,col="black",fontface="bold",size=6)+ labs(y="Final predator biomass | others",x="Grazer diversity | others",nrow=2)+ facet_grid(~group,scales="free_x")+ theme_bw(base_size=18)+ guides(col=guide_legend(ncol=2))+ theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank())