#These code for the statistical freeware package R were used in the data analyses and graph drawing in ‘Bradford, M.A., Berg, B., Maynard, D.S., Wieder, W.R. and Wood, S.A. Understanding the dominant controls on litter decomposition. Journal of Ecology, in press.’ #Data and code, along with a ReadMe file (BergDecompReadMe.txt), were submitted to the Dryad data repository October 28th 2015 (the day of the manuscript’s acceptance). #The persons associated with this code have dedicated the work to the public domain by waiving all of their rights to the code worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law. You can copy, modify, distribute and use the code, even for commercial purposes, all without asking permission. dat<-read.csv() #### this is the main data library(lme4) library(sampling) library(LMERConvenienceFunctions) dat$MAT2<-dat$MAT*dat$MAT dat$Siteno<-as.character(dat$Siteno) dat.o<-dat ### some misc functions used throughout strat<-function(tdat,nsamp){ n<-nrow(tdat) tdat[sample(1:n,nsamp),] } use.norm<-function(tnorm.vec,nsamp){ tv<-as.matrix(rnorm(nsamp, mean=tnorm.vec$mean,sd=tnorm.vec$std)) apply(tv,1,min.max) } min.max<-function(t){ min(max(t,0),100) } #stratify the data based on temperature dat.high<-dat.o[dat.o$MAT>15.2-(15.2+1.7)/2,] dat.low<-dat.o[dat.o$MAT<=15.2-(15.2+1.7)/2,] ### the main bootstrapping function boot.fun<-function(tdat,norm.vec,nsamp=4,niter=1000){ out.dat<-data.frame( lm.inter=rep(-9,niter),lm.slope=rep(-9,niter),lm.X2=rep(-9,niter), lm.r2.fix=rep(-9,niter),lm.r2.rand=rep(-9,niter),lm.m.inter=rep(-9,niter), lm.m.slope=rep(-9,niter),lm.m.X2=rep(-9,niter),lm.m.r2=rep(-9,niter)) for(i in 1:niter){ tsdat<-by(cbind(tdat$Percloss,dat$MAT),tdat$Siteno,strat,nsamp=nsamp) keep.mass<-rep(c(rep(TRUE,nsamp),rep(FALSE,nsamp)),length(names(tsdat))) adat<-unlist(as.list(tsdat)) ttdat<-data.frame(Percloss=as.numeric(adat[keep.mass]),MAT=as.numeric(adat[!keep.mass]),Siteno=rep(names(tsdat),each=nsamp)) ttdat.m<-data.frame(Percloss=as.numeric(unlist(as.list(by(ttdat$Percloss,ttdat$Siteno,mean)))),MAT=as.numeric(unlist(as.list(by(ttdat$MAT,ttdat$Siteno,mean)))),Siteno=names(by(ttdat$MAT,ttdat$Siteno,mean))) ttdat$MAT2<-ttdat$MAT*ttdat$MAT ttdat.m$MAT2<-ttdat.m$MAT*ttdat.m$MAT tf.lm.x2<-lmer(Percloss~MAT+MAT2+(1|Siteno),data=ttdat) tf.m.x2<-lm(Percloss~MAT+MAT2,data=ttdat.m) tf.lm<-lmer(Percloss~MAT+(1|Siteno),data=ttdat) tf.m<-lm(Percloss~MAT,data=ttdat.m) if(anova(tf.lm,tf.lm.x2)$"Pr(>Chisq)"[2]<.05){ out.dat$lm.inter[i]<-fixef(tf.lm.x2)[1] out.dat$lm.slope[i]<-fixef(tf.lm.x2)[2] out.dat$lm.X2[i]<-fixef(tf.lm.x2)[3] out.dat$lm.r2.fix[i]<-r2.mixed(tf.lm.x2)$fR2 out.dat$lm.r2.rand[i]<-r2.mixed(tf.lm.x2)$rfR2 } else{ out.dat$lm.inter[i]<-fixef(tf.lm)[1] out.dat$lm.slope[i]<-fixef(tf.lm)[2] out.dat$lm.r2.fix[i]<-r2.mixed(tf.lm)$fR2 out.dat$lm.r2.rand[i]<-r2.mixed(tf.lm)$rfR2 } if(anova(tf.m,tf.m.x2)$"Pr(>F)"[2]<.05){ out.dat$lm.m.inter[i]<-tf.m.x2$coef[1] out.dat$lm.m.slope[i]<-tf.m.x2$coef[2] out.dat$lm.m.X2[i]<-tf.lm.x2$coef[3] out.dat$lm.m.r2[i]<-summary(tf.m.x2)$adj.r.squared } else{ out.dat$lm.m.inter[i]<-tf.m$coef[1] out.dat$lm.m.slope[i]<-tf.m$coef[2] out.dat$lm.m.r2[i]<-summary(tf.m)$adj.r.squared } } out.dat } ############ This function calculates a pseudo-R2 for mixed effects models, and is modified from: # Nakagawa and Schielzeth, 2013. "A general and simple method for obtaining R 2 from generalized linear mixed-effects models". Methods in Ecology and Evolution, 4(2):133-142 r2.mixed<-function(mF){ #extract the design matrix for the fixed effects mFX<-model.matrix(mF) #calculate the variance of the fixed effects VarF <- var(as.vector(fixef(mF) %*% t(mFX))) #calculate the variance of the random effects VarR<-sum(as.numeric(VarCorr(mF))) #calculate the residual variance VarResid<-attr(VarCorr(mF), "sc")^2 #calculate the fraction of total variance explained by the fixed effect variance fR2<-VarF/(VarF + VarR + VarResid) #calculate the fraction of total variance explained by the fixed effect AND random effect variance rfR2<-(VarF + VarR)/(VarF + VarR + VarResid) list(fR2=fR2,rfR2=rfR2) } ###### plot the mean and all data using temp^2 dat.o.m<-data.frame(Percloss=as.numeric(unlist(as.list(by(dat.o$Percloss,dat.o$Siteno,mean)))),MAT=as.numeric(unlist(as.list(by(dat.o$MAT,dat.o$Siteno,mean)))),Siteno=names(by(dat.o$Percloss,dat.o$Siteno,mean))) dat.o.m$MAT2<-dat.o.m$MAT*dat.o.m$MAT summary(f1<-lm(Percloss~MAT+MAT2,data=dat.o)) summary(f2<-lm(Percloss~MAT+MAT2,data=dat.o.m)) summary(f3<-lmer(Percloss~MAT+MAT2+(1|Siteno),data=dat.o)) pamer.fnc(f3) r2.mixed(f3) pred.lmer<-function(tfit,newvals){ tcoefs<-fixef(tfit) tcoefs[1]*newvals[,1]+tcoefs[2]*newvals[,2]+tcoefs[3]*newvals[,3] } xvals<-seq(-10,30,.01) xvals2<-xvals^2 newMM<-as.matrix(data.frame(int=rep(1,length(xvals)),MAT=xvals,MAT2=xvals2)) yvals<-pred.lmer(f3,newMM) newdat<-data.frame(MAT=xvals) newdat$y<-newMM%*%fixef(f3) fit.band2<-data.frame(predict(f2,interval="confidence",newdata=data.frame(MAT=seq(-10,20,.1),MAT2=(seq(-10,20,.1))^2))) #bootstrap CI for lmer output predFun<-function(f3) newMM%*%fixef(f3) bb<-bootMer(f3,FUN=predFun,nsim=1000) ### occasionally a random selection will not converge. just rerun if needed. bb_se<-apply(bb$t,2,function(x) tout<-quantile(x,probs=c(0.025,0.975))) newdat$blo<-bb_se[1,] newdat$bhi<-bb_se[2,] ### plot the data plot(dat.o$Percloss~dat.o$MAT, main="",xlab="",ylab="",xlim=c(-2,15.5),ylim=c(0,60),pch=16,cex=1,family="serif",xaxs="i",yaxs="i",cex.axis=1.5) lines(newdat$y~newdat$MAT,lwd=1.2) points(dat.o.m$Percloss~dat.o.m$MAT,pch=4,cex=2.5) polygon(c(newdat$MAT,rev(newdat$MAT)),c(newdat$blo,rev(newdat$bhi)),border="red",lty=1,lwd=1.8) lines(fit.band2$lwr~seq(-10,20,.1),col="red",lwd=1.8,lty=2) lines(fit.band2$upr~seq(-10,20,.1),col="red",lwd=1.8,lty=2) ########### look at n=4 vs 20 in the high, low regions ############################# t.niter<-10000 dat<-dat.low dat$Siteno<-as.character(dat$Siteno) norm.vec<-data.frame(Siteno=names(by(dat$Percloss,dat$Siteno,mean)),mean=as.numeric(unlist(as.list(by(dat$Percloss,dat$Siteno,mean)))),std=as.numeric(unlist(as.list(by(dat$Percloss,dat$Siteno,sd))))) norm.vec<-norm.vec[!is.na(norm.vec$mean),] norm.vec$Siteno<-as.character(norm.vec$Siteno) boot.vals.low.4<-boot.fun(dat,norm.vec,nsamp=4,niter=t.niter) boot.vals.low.20<-boot.fun(dat,norm.vec,nsamp=20,niter=t.niter) dat<-dat.high dat$Siteno<-as.character(dat$Siteno) norm.vec<-data.frame(Siteno=names(by(dat$Percloss,dat$Siteno,mean)),mean=as.numeric(unlist(as.list(by(dat$Percloss,dat$Siteno,mean)))),std=as.numeric(unlist(as.list(by(dat$Percloss,dat$Siteno,sd))))) norm.vec<-norm.vec[!is.na(norm.vec$mean),] norm.vec$Siteno<-as.character(norm.vec$Siteno) boot.vals.high.4<-boot.fun(dat,norm.vec,nsamp=4,niter=t.niter) boot.vals.high.20<-boot.fun(dat,norm.vec,nsamp=20,niter=t.niter) #################### get bootstrap CI, low values ## low, n=4 lval<-seq(-2,7,.01) lval2<-lval^2 mboundl4<-matrix(-9,nrow=length(lval),ncol=2) tpred<-rep(-9,t.niter) for(i in 1:length(lval)){ for(j in 1:t.niter){ tpred[j]<-boot.vals.low.4$lm.inter[j] + boot.vals.low.4$lm.slope[j]*lval[i] } mboundl4[i,]<-as.numeric(quantile(tpred,probs=c(0.025,0.975))) } ## low,n=20 mboundl20<-matrix(-9,nrow=length(lval),ncol=2) tpred<-rep(-9,t.niter) for(i in 1:length(lval)){ for(j in 1:t.niter){ tpred[j]<-boot.vals.low.20$lm.inter[j] + boot.vals.low.20$lm.slope[j]*lval[i] } mboundl20[i,]<-as.numeric(quantile(tpred,probs=c(0.025,0.975))) } ################ get bootstrap CI, high values #high, n=4 lval<-seq(6,16,.01) lval2<-lval^2 mboundh4<-matrix(-9,nrow=length(lval),ncol=2) tpred<-rep(-9,t.niter) for(i in 1:length(lval)){ for(j in 1:t.niter){ tpred[j]<-boot.vals.high.4$lm.inter[j] + boot.vals.high.4$lm.slope[j]*lval[i] } mboundh4[i,]<-as.numeric(quantile(tpred,probs=c(0.025,0.975))) } #high,n=20 mboundh20<-matrix(-9,nrow=length(lval),ncol=2) tpred<-rep(-9,t.niter) for(i in 1:length(lval)){ for(j in 1:t.niter){ tpred[j]<-boot.vals.high.20$lm.inter[j] + boot.vals.high.20$lm.slope[j]*lval[i] } mboundh20[i,]<-as.numeric(quantile(tpred,probs=c(0.025,0.975))) } ############## plot the data par(mfrow=c(1,2)) lval<-seq(-2,7,.01) lval2<-lval^2 plot(Percloss~MAT,data=dat.low,main="",xlab="",ylab="",xlim=c(-2,15.5),ylim=c(0,60),pch=1,cex=1,family="serif",xaxs="i",yaxs="i",cex.axis=1.5) flow<-lm(Percloss~MAT,data=dat.low) lines(x=lval,y=mboundl4[,1],col="red",lwd=1.5,lty=2) lines(x=lval,y=mboundl4[,2],col="red",lwd=1.5,lty=2) lines(x=lval,y=flow$coef[1]+flow$coef[2]*lval,col="black",lwd=1.5) lines(x=lval,y=mboundl20[,1],col="red",lwd=1.5,lty=1) lines(x=lval,y=mboundl20[,2],col="red",lwd=1.5,lty=1) lval<-seq(6,16,.01) lval2<-lval^2 plot(Percloss~MAT,data=dat.high,main="",xlab="",ylab="",xlim=c(-2,15.5),ylim=c(0,60),pch=1,cex=1,family="serif",xaxs="i",yaxs="i",cex.axis=1.5) flow<-lm(Percloss~MAT,data=dat.high) lines(x=lval,y=mboundh4[,1],col="red",lwd=1.5,lty=2) lines(x=lval,y=mboundh4[,2],col="red",lwd=1.5,lty=2) lines(x=lval,y=flow$coef[1]+flow$coef[2]*lval,col="black",lwd=1.5) lines(x=lval,y=mboundh20[,1],col="red",lwd=1.5,lty=1) lines(x=lval,y=mboundh20[,2],col="red",lwd=1.5,lty=1) ### get summary tables and stats boot.vals.low.4$lm.rand<-boot.vals.low.4$lm.r2.rand-boot.vals.low.4$lm.r2.fix boot.vals.low.20$lm.rand<-boot.vals.low.20$lm.r2.rand-boot.vals.low.20$lm.r2.fix boot.vals.high.4$lm.rand<-boot.vals.high.4$lm.r2.rand-boot.vals.high.4$lm.r2.fix boot.vals.high.20$lm.rand<-boot.vals.high.20$lm.r2.rand-boot.vals.high.20$lm.r2.fix calc.quantiles<-function(ttvec,qvec=c(0.025,0.975),tkeep=c("lm.slope","lm.m.slope","lm.r2.fix","lm.rand","lm.r2.rand","lm.m.r2")){ toutd<-data.frame(matrix(nrow=length(tkeep),ncol=length(qvec)+2)) row.names(toutd)<-tkeep names(toutd)<-c(qvec,"mean","median") ln<-length(qvec) for(i in 1:length(tkeep)){ toutd[i,1:ln]<-quantile(ttvec[,names(ttvec)==tkeep[i]],probs=qvec,na.rm=T) toutd$mean[i]<-mean(ttvec[,names(ttvec)==tkeep[i]],na.rm=T) toutd$median[i]<-median(ttvec[,names(ttvec)==tkeep[i]],na.rm=T) } toutd } qlow4<-calc.quantiles(boot.vals.low.4) qlow20<-calc.quantiles(boot.vals.low.20) qhigh4<-calc.quantiles(boot.vals.high.4) qhigh20<-calc.quantiles(boot.vals.high.20) toutd<-data.frame(model=rep(c("qlow4","qlow20","qhigh4","qhigh20"),each=nrow(qlow4)),rbind(qlow4,qlow20,qhigh4,qhigh20)) toutd