###------------------------------------------------------------------------### ### R code to construct figures 1a and 1b, including underlying modelling ### Knut Rydgren, April 2018 ###------------------------------------------------------------------------### #--------------------------------------------------------------------------------------------------# #### Load data & R-packages #### #--------------------------------------------------------------------------------------------------# rm(list=ls(all=TRUE)) setwd("") # add the directory where you have saved the data files, that below is where I have placed the files setwd("c:\\Dta\\prj\\Release\\Art1_VegDyn\\Fig1_script") # use the path where you have saved the data files to be used library(lmerTest) library(lme4) dta.ordscore<-read.table("Rydgren&al_dtafig1a.txt",header=T) #imports the data with the ordination scores (DCA axes 1 and 2) for the restored plots (restoration plots from three different time periods) and the reference plots. In this case we had a static reference from the last time period. attach(dta.ordscore) names(dta.ordscore) head(dta.ordscore) dta.sucdist<-read.table("Rydgren&al_dtafig1b.txt",header=T) # imports the data with the successional distances attach(dta.sucdist) names(dta.sucdist) head(dta.sucdist) # Extract subsets for plotting purposes restored.dta <- subset(dta.ordscore,type=="restored") reference.dta<- subset(dta.ordscore,type=="reference") #--------------------------------------------------------------------------------------------------# #### Figure 1a #### #### USING DTA.ORDSCORE #--------------------------------------------------------------------------------------------------# par(mar=(c(4,4,3,2)+0.1)) # Changes the plotting margins par(mfrow=c(2,1)) # makes space for figures 1a and 1b on the same page plot(DCA1,DCA2,type="n",xlab="Ordination axis 1",ylab="Ordination axis 2", main="",xlim=c(0,5.0),ylim=c(-1.0,3.5)) points(DCA1[type=="reference" & yr=="2015"],DCA2[type=="reference" & yr=="2015"],col="green",pch=22, cex=1.0) points(DCA1[type=="restored" & yr=="2015"],DCA2[type=="restored" & yr=="2015"],col="darkorchid2",pch=15, cex=1.0) points(DCA1[type=="restored" & yr=="2008"],DCA2[type=="restored" & yr=="2008"],col="darkgoldenrod4",pch=15, cex=1.0) points(DCA1[type=="restored" & yr=="1991"],DCA2[type=="restored" & yr=="1991"],col="darkorange3",pch=15, cex=1.0) CentroidRef1<-3.890 # the centroid position of the reference along the first ordination axis CentroidRef2<-1.388 # the centroid position of the reference along the second ordination axis points(CentroidRef1, CentroidRef2, col="green",pch=16,cex=1.5) # These arrows indicate the distances from three chosen restoration plots from three different time points, along the first ordination axis to the centroid of the static reference arrows(0.18, 1.43, 0.18, 0.35,length=0.1,angle=20,col="darkorange3",lwd=0.3) arrows(0.18, 0.35, 3.89, 0.35,lty=2,length=0.2,angle=20,col="darkorange3",lwd=1.0) arrows(1.75, 1.26, 1.75,-0.05,length=0.1,angle=20,col="darkgoldenrod4",lwd=0.3) arrows(1.75,-0.05, 3.89,-0.05,lty=2,col="darkgoldenrod4",length=0.2,angle=20,lwd=1,code=2) arrows(2.74, 0.85, 2.74,-0.40,length=0.1,angle=20,col="darkorchid2",lwd=0.3) arrows(2.74,-0.40, 3.89,-0.40,lty=2,col="darkorchid2",length=0.2,angle=20,lwd=1,code=2) text(-0.0, -0.5,expression("(1)"), cex=0.9) text(0.42, -0.5, expression(d[paste("jt,0")] == x[0] - r[t]),cex=0.9) # Equation 1 in our ms. Adjust your plotting window so that you avoid overlapping text text(0.68, 0.5, expression(d[paste("121,1,0")] == "3.89" - "0.18"),cex=0.7) # Equation 1 has been applied to the chosen restoration plot from time 1; in this way we obtain the successional distance, i.e. the distance between a restored plot and the reference along the successional gradient (to be done for all restoration plots) text(2.22, 0.07, expression(d[paste("121,2,0")] == "3.89" - "1.75"),cex=0.7) # Similar to the line above except that it is from the chosen restored plot from time 2 text(3.22, -0.25, expression(d[paste("121,3,0")] == "3.89" - "2.74"),cex=0.7) # Similar to the line above except that it is from the chosen restored plot from time 3 legend("topright",legend=c("Restored (time 1)","Restored (time 2)","Restored (time 3)", "Static reference"),pch=c(15,15,15,22),col=c("darkorange3", "gold3", "darkorchid2", "green"), bty="n", cex=0.9) legend("topleft",c("(a)"), bty="n", cex=1.0) #---------------------------------------------------------------------------------------------------------------------# #### Figure 1b #### #### Make the models for the estimation of time-to-recovery and draw the model lines, USING DTA.SUCDIST #---------------------------------------------------------------------------------------------------------------------# # Plot the successional distances for all restored plots from the three different time points, i.e. 7, 24 and 31 years after disturbance plot(SucAge,DCA1s,xlim=c(1,80),ylim=c(-0.2,5.0), col=1,pch=1, main="",xlab="Time since disturbance (years)", ylab=expression("Successional distance" ~ (italic(d[jt][","][0])))) abline(h=0.0,lty=2, col="black") # Mark the three plots in colour for which the calculations in Fig 1A are shown points(7,3.710, pch=16, col="darkorange3", cex=1.4) points(24,2.140, pch=16, col="darkgoldenrod4", cex=1.4) points(31,1.150, pch=16, col="darkorchid2", cex=1.4) legend("topleft",c("(b)"), bty="n", cex=1.0) # Accounting for temporal (ftemporal_pseudo) and spatial pseudoreplication (fblock); first make factor variables ftemporal_pseudo<-as.factor(temporal_pseudo) fblock<-as.factor(block) # Modelling the successional distance as a function of years after disturbance; applying a linear model and logarithmic (asymptotic) model mod.lmer1<-lmer(DCA1s~SucAge+(1|ftemporal_pseudo)+(1|fblock), REML="TRUE", data=dta.sucdist) summary (mod.lmer1) mod.log.lmer1<-lmer(log(DCA1s)~SucAge+(1|ftemporal_pseudo)+(1|fblock), REML="TRUE", data=dta.sucdist) summary (mod.log.lmer1) # Plotting the linear model (mod.lmer1 obtained above), with confidence intervals for predicting time to recovery xv <- 1:80 yv <- expand.grid(SucAge = xv) X <- model.matrix(~ SucAge,data = yv) yv$fit <- X %*% fixef(mod.lmer1) yv$SE <- sqrt(diag(X %*% vcov(mod.lmer1) %*% t(X))) yv$lo <- yv$fit - (1.96 * yv$SE) yv$up <- yv$fit + (1.96 * yv$SE) confint(mod.lmer1) # The bottom two lines give the usual 95% confidence intervals for the fixed intercept and random effects. I The top two lines are approximate confidence intervals for the variance components parameters ?? and ?? polygon(x=c(xv,rev(xv)),y=c(yv$lo,rev(yv$up)),col=rgb(0.2, 0.3, 0.4, alpha=0.3),border=rgb(0.2, 0.3, 0.4, alpha=0.3)) lines(xv,yv$fit,col="black",lty=1) # Plotting the the asymptotic (logarithmic) model (mod.log.lmer1 obtained above), with confidence intervals for predicting time to recovery yv <- expand.grid(SucAge = xv) X <- model.matrix(~ SucAge,data = yv) yv$fit <- X %*% fixef(mod.log.lmer1) yv$SE <- sqrt(diag(X %*% vcov(mod.log.lmer1) %*% t(X))) yv$lo <- yv$fit - (1.96 * yv$SE ) yv$up <- yv$fit + (1.96 * yv$SE ) polygon(x=c(xv,rev(xv)),y=exp(c(yv$lo,rev(yv$up))),col=rgb(1, 0, 0, alpha=0.3),border=rgb(1, 0, 0, alpha=0.3)) lines(xv,exp(yv$fit),col="red") legend("topright",bty="n",legend=c("Linear model","Asymptotic model", "Static reference"),lty=c(1,1,1),col=c("black","red","green"), cex=0.9) # Plotting the threshold value for the successional distance djt,0 at which restoration is regarded as successful # Here we used +1 S.D. calcualted from the ordination scores for the first ordination axes of the static reference plots. In this case it was= 0.192 polygon(x=c(xv,rev(xv)),y=c(rep(-0.192,length(xv)),rep(0.192,length(xv))),col=rgb(0, 1, 0, alpha=0.3),border=rgb(0, 1, 0, alpha=0.3)) #---------------------------------------------------------------# #### The end #### #---------------------------------------------------------------#