rm(list=ls()) ###### required packages library(MASS) library(nlme)## to fit linear mixed-effects models that are the components of the multi-level confirmatory path analyses library(lme4)## to fit negative binomials library(reshape) ## data formatting (though more recent packages could streamline this) library(vegan)####to calculate diversity metrics from abudance data.frames library(splitstackshape) ## data formatting library(piecewiseSEM)###to conduct confirmatory path analysis (piecewiseSEM) library(car) library(shape) library(pscl) library(reshape2) ## data formatting (though more recent packages could streamline this) library(DHARMa)###for negative binomial diagnostics library(vcd)###for negative binomial diagnostics ## change to match your computer setwd("C:/Users/Isaac Towers/Desktop/Proc B/Proc B Resubmission 4_08_2018") #######make function for to put correlation value in the top right panel ## define a funciton to standardise variables - to use later std.func<-function(x) as.numeric((x-mean(x, na.rm=T))/sd(x, na.rm=T)) ##########define a function to replace missing paths that use the wrong degrees of freedom (AICc uses 12 because lowest number of observations in dataset, same process used in piecewiseSEM) replace.path.func<-function(old.paths, row.to.replace, new.mod){ if(class(new.mod)=="lm"){ old.paths$missing.paths[,6]<-ifelse(old.paths$missing.paths[,6]==0, 1e-9, old.paths$missing.paths[,6]) ## this is necessary if some of the p vals for missing paths are really small and get rounded to zero in the missing path data.frame old.paths$missing.paths[row.to.replace,c(2,3,4,5,6)]<-c(summary(new.mod)$coef[2,1], summary(new.mod)$coef[2,2], summary(new.mod)$df[2], summary(new.mod)$coef[2,3], summary(new.mod)$coef[2,4]) old.paths$missing.paths[row.to.replace,7]<-ifelse((summary(new.mod)$coef[2,4])<0.05, "*","") missing.paths<-old.paths$missing.paths fisher.c<--2*sum(log(old.paths$missing.paths$p.value)) dfs<-2*(dim(old.paths$missing.paths)[1]) p.val<-1-pchisq(q=fisher.c, df=dfs) K<-old.paths$AIC$K AIC<-fisher.c+2*old.paths$AIC$K AICC<-fisher.c+2*old.paths$AIC$K*(12/(12-K-1)) return(list(missing.paths=missing.paths, AIC=list(K=K, AIC=AIC, AICC=AICC), stats=data.frame(fisher.c, dfs, p.val))) } } ##standard error function std.error<-function (x) sd(x)/sqrt(length(!is.na(x))) ###function to predict linear curves only lm.predict<-function(mod, newdat){ pred = predict(mod, newdata = newdat, interval = "confidence", level = 0.95) upper<-pred[,3] lower<-pred[,2] return(data.frame(pred=pred[,1], upper, lower)) } ###function to predict linear mixed effects models lme.predict<-function(mod, newdat, int.type=NA, se.mult){ pred <- as.vector(predict(mod, newdat, level = 0)) Designmat <- model.matrix(formula(mod)[-2], newdat) predvar <- diag(Designmat %*% vcov(mod) %*% t(Designmat)) SE <- sqrt(predvar) SE2 <- sqrt(predvar+mod$sigma^2) if(int.type=="CI"){ upper<-pred + (se.mult*SE) lower<-pred - (se.mult*SE) } if(int.type=="PI"){ upper<-pred + (se.mult*SE2) lower<-pred - (se.mult*SE2) } return(data.frame(pred, upper, lower)) } ####function to plot the curve and associated condfidence interval plot.CI.func<- function(x.for.plot, pred, upper, lower, env.colour, env.trans=NA, line.colour, line.weight){ colour.rgb<-col2rgb(col=env.colour) polygon.coords<-data.frame(rbind(cbind(x.for.plot[1], lower[1]), cbind(x.for.plot, upper), cbind(x.for.plot, lower)[rev(order(x.for.plot)),])) names(polygon.coords)<-c("x", "y") polygon(polygon.coords$x, polygon.coords$y, col=rgb(red=colour.rgb["red",],blue=colour.rgb["blue",], green=colour.rgb["green",] , alpha=env.trans, maxColorValue = 255), border=NA) lines(x.for.plot, pred, col=line.colour, lwd=line.weight) } ###function to inverse logits inv.logit<-function(x)(exp(x)/(1+exp(x))) ##creates a sequence of values frmo min to max by 100 seq.func<-function(x)(seq(min(x, na.rm=T), max(x, na.rm=T), length.out=100)) ### plot with 95% CI using lmer or lme objects lmer.predict<-function(mod, newdat, se.mult, binom=NULL, poisson=NULL){ pvar1 <- diag(as.matrix(newdat) %*% tcrossprod(vcov(mod),as.matrix(newdat))) newdat$y<- as.matrix(newdat) %*% fixef(mod) newdat <- data.frame(newdat, plo = newdat$y-(se.mult*sqrt(pvar1)), phi = newdat$y+(se.mult*sqrt(pvar1))) ## if you have used binomial errors then this will back transform logits to the probability scale if(binom==T) { newdat$y<-plogis(newdat$y); newdat$plo<-plogis(newdat$plo); newdat$phi<-plogis(newdat$phi) } else ## if you have used poisson errors or have log-transformed your response, then this will back transform to the original scale (e.g. abundance) ## not used in this code if(poisson==T) { newdat$y<-exp(newdat$y); newdat$plo<-exp(newdat$plo); newdat$phi<-exp(newdat$phi) } return(with(newdat, data.frame(y, phi, plo))) } plot.CI.func<- function(x.for.plot, pred, upper, lower, env.colour, env.trans=NA, line.colour, line.weight, line.type){ colour.rgb<-col2rgb(col=env.colour) polygon.coords<-data.frame(rbind(cbind(x.for.plot[1], lower[1]), cbind(x.for.plot, upper), cbind(x.for.plot, lower)[rev(order(x.for.plot)),])) names(polygon.coords)<-c("x", "y") polygon(polygon.coords$x, polygon.coords$y, col=rgb(red=colour.rgb["red",],blue=colour.rgb["blue",], green=colour.rgb["green",] , alpha=env.trans, maxColorValue = 255), border=NA) lines(x.for.plot, pred, col=line.colour, lwd=line.weight, lty=line.type) } #####function to standardise variables standardised.variable.func<-function(original.var){ if(((max(original.var)-min(original.var))/5>1)){ round(seq(from=(min(original.var)), to=(max(original.var)), by=((max(original.var)-min(original.var))/5))) } if(((max(original.var)-min(original.var))/5>0.5 && max(original.var)-min(original.var))/5<1){ round(seq(from=(min(original.var)), to=(max(original.var)), by=((max(original.var)-min(original.var))/5))/0.5)*0.5 } } axis.text.func <- function(side.n, variable, increment, sq.root) { ifelse(sq.root, axis(side=side.n, at=(sqrt(seq(floor(min(variable^2, na.rm=T)), ceiling(max(variable^2, na.rm=T)), increment)) - mean(variable, na.rm=T)) / sd(variable, na.rm=T), labels=seq(floor(min(variable^2, na.rm=T)), ceiling(max(variable^2, na.rm=T)), increment)), axis(side=side.n, at=(seq(floor(min(variable, na.rm=T)), ceiling(max(variable, na.rm=T)), increment) - mean(variable, na.rm=T)) / sd(variable, na.rm=T), labels=seq(floor(min(variable, na.rm=T)), ceiling(max(variable, na.rm=T)), increment))) } ## species in quadrat scale data ## one row for every species recorded in every quadrat ## this file includes a bunch of trait values for each species in each quadrat plants.core.2011<-read.csv("plants.core.2011.csv", header=T) # now read in the quadrat scale data (qs stands for quadrat scale) qs.core.2011<-read.csv("quadrat_scale_data.csv", header=T) ## first just look at 2011 plants.core.2011<-droplevels(subset(plants.core.2011, year==2011)) qs.core.2011<-droplevels(subset(qs.core.2011, year==2011)) ##just for native plants - so that we can work out native species richness native.plants.core.2011<-droplevels(subset(plants.core.2011, exotic==0)) native.species<-native.plants.core.2011$species exotic.plants.core.2011<-droplevels(subset(plants.core.2011, exotic==1)) exotic.species<-exotic.plants.core.2011$species ## create "community data" objects that VEGAN likes ## quadrat level data (quadrats as rows and species as columns) ("interaction neighbourhood scale") comm.data.2011<-with(plants.core.2011, tapply(count, list(yrsq, species), sum))[qs.core.2011$yrsq,] ## this bit at the end maintains the correct order comm.data.2011<-ifelse(is.na(comm.data.2011), 0,comm.data.2011) ##and the same for the native plants native.comm.data.2011<-with(native.plants.core.2011, tapply(count, list(yrsq, species), sum))[qs.core.2011$yrsq,] native.comm.data.2011<-ifelse(is.na(native.comm.data.2011), 0,native.comm.data.2011) ## from this object we can calculate diversity measures for each quadrat! qs.core.2011$richness<-specnumber(comm.data.2011) qs.core.2011$native.richness<-specnumber(native.comm.data.2011) qs.core.2011$exotic.richness<-qs.core.2011$richness-qs.core.2011$native.richness ## make some transformed environemntal variables qs.core.2011$sqrt.p<-sqrt(qs.core.2011$colwell.p) qs.core.2011$sqrt.k<-sqrt(qs.core.2011$colwell.k) qs.core.2011$ln.tot.n<-log(qs.core.2011$tot.n) qs.core.2011$ln.plants.per.quadrat<-log(qs.core.2011$plants.per.quadrat) ##standardise each of these variables and paste std. on the front std.data.frame<-with(qs.core.2011, data.frame(crust,sqrt.p, sqrt.k, ln.tot.n, mean.cover.n.w.e, litter.cover, log, relative.compact, dist.to.edge, native.richness, exotic.richness, ln.plants.per.quadrat)) std.data<-data.frame(apply(std.data.frame, 2, std.func)) colnames(std.data)<-paste("std.", colnames(std.data), sep="") ## add these to the qs.core.2011 dataset qs.core.2011<-data.frame(qs.core.2011, std.data) ###correlations between abiotic variables at quadrat scale?? abiotic.correlations<-cor(std.data, method="pearson") ## now we do the same for each site ("community" in our case) comm.site.data.2011<-with(plants.core.2011, tapply(count, list(yrs, species), sum)) comm.site.data.2011<-ifelse(is.na(comm.site.data.2011), 0,comm.site.data.2011) ##as well as the native comm.site.data.2011.native<-with(native.plants.core.2011, tapply(count, list(yrs, species), sum)) comm.site.data.2011.native<-ifelse(is.na(comm.site.data.2011.native), 0,comm.site.data.2011.native) ##create site-scale dataset site.data<-data.frame(site.richness=specnumber(comm.site.data.2011), site.native.richness=specnumber(comm.site.data.2011.native)) site.data$site.exotic.richness<-site.data$site.richness-site.data$site.native.richness ####find the mean value of neighbourhood variables and assign them to the site scale means.of.nei.var<-with(qs.core.2011, data.frame(exotic.richness, std.crust, std.dist.to.edge, std.log,std.sqrt.p,std.sqrt.k,std.ln.tot.n, std.relative.compact, std.mean.cover.n.w.e, std.litter.cover, std.ln.plants.per.quadrat)) means.of.nei.var.melt<-melt(data.frame(yrs=qs.core.2011$yrs, means.of.nei.var)) means.of.nei.var<-cast(means.of.nei.var.melt, yrs~variable, mean, na.rm=T) site.data<-data.frame(site.data, means.of.nei.var) ####standardise some variables that originated at the site scale - i.e. diversity measures std.data.frame.site<-with(site.data, data.frame(site.native.richness,site.exotic.richness)) std.data.site<-data.frame(apply(std.data.frame.site, 2, std.func)) colnames(std.data.site)<-paste("std.", colnames(std.data.site), sep="") site.data<-data.frame(site.data, std.data.site) #####all the heterogeneity metrics soil.nut.dist<-dist(with(qs.core.2011, data.frame(std.sqrt.p, std.sqrt.k, std.ln.tot.n))) soil.nut.disper1<-betadisper(soil.nut.dist, group=qs.core.2011$yrs, type="median") site.data$std.site.soil.nut.disper<-std.func(tapply(soil.nut.disper1$distances, qs.core.2011$yrs, mean)) phys.dist<-dist(with(qs.core.2011, data.frame(std.crust, std.relative.compact, std.log, std.litter.cover))) phys.disper1<-betadisper(phys.dist, group=qs.core.2011$yrs, type="median") site.data$std.site.phys.disper<-std.func(tapply(phys.disper1$distances, qs.core.2011$yrs, mean)) ### just the cover heterogeneity variable - using the distance from average dist.from.avg.func<-function(x) mean(abs(x-mean(x, na.rm=T)), na.rm=T) site.data$std.canopy.cover.disper<-std.func(tapply(qs.core.2011$mean.cover.n.w.e, qs.core.2011$yrs, dist.from.avg.func)) ###.... and for each remnant comm.remnant.data.2011<-with(plants.core.2011, tapply(count, list(yr, species), sum)) comm.remnant.data.2011<-ifelse(is.na(comm.remnant.data.2011), 0,comm.remnant.data.2011) comm.native.remnant.data.2011<-with(native.plants.core.2011, tapply(count, list(yr,species), sum)) comm.native.remnant.data.2011<-ifelse(is.na(comm.native.remnant.data.2011), 0,comm.native.remnant.data.2011) ###add yr to the site scale data frame site.data$yr<-qs.core.2011$yr[match(site.data$yrs, qs.core.2011$yrs)] ###contstruct a dataset for remnants remnant.data<-data.frame(remnant.richness=specnumber(comm.remnant.data.2011), remnant.native.richness=specnumber(comm.native.remnant.data.2011)) remnant.data$remnant.exotic.richness<-remnant.data$remnant.richness-remnant.data$remnant.native.richness remnant.data$yr=rownames(remnant.data) #####for each remnant find the remnant-scale mean of each neighbourhood variable means.of.nei.var.rem<-with(qs.core.2011, data.frame(std.crust, exotic.richness, std.dist.to.edge, std.log,std.sqrt.p,std.sqrt.k,std.ln.tot.n, std.relative.compact, std.mean.cover.n.w.e, std.litter.cover, std.ln.plants.per.quadrat,mean.gstmin,gsai)) means.of.nei.var.melt.rem<-melt(data.frame(yr=qs.core.2011$yr, means.of.nei.var.rem)) means.of.nei.var.rem<-cast(means.of.nei.var.melt.rem, yr~variable, mean, na.rm=T) remnant.data<-data.frame(remnant.data, means.of.nei.var.rem) #######Jackknife richness estimates at the community scale #######subset the species data to only exotic species###### ###extract the names of exotic species exotic.plants.core.2011<-plants.core.2011[plants.core.2011$exotic==1,] exotic.names<-(droplevels(unique(exotic.plants.core.2011$species))) exotic.names<-as.character(exotic.names) ####convert community matrix into a dataframe comm.data.2011<-as.data.frame(comm.data.2011) ###subset community matrix based on exotic species names exotic.comm.data.2011<-comm.data.2011[(exotic.names)] exotic.species.richness.estimates<-specpool(exotic.comm.data.2011, pool=qs.core.2011$yrs) ####set communities with NA for richness to 1, maximum species richess estimate in that remnant (see Towers and Dwyer methods) exotic.species.richness.estimates[7:8,]<-1 site.data$jack1exotic<-exotic.species.richness.estimates$jack1 site.data$std.jack1exotic<-std.func(exotic.species.richness.estimates$jack1) qs.core.2011$std.jack1exotic<-rep(site.data$std.jack1exotic, each=15) qs.core.2011$jack1exotic<-rep(site.data$jack1exotic, each=15) native.species.richness.estimates<-specpool(native.comm.data.2011, pool=qs.core.2011$yrs) site.data$jack1native<-native.species.richness.estimates$jack1 site.data$std.jack1native<-std.func(native.species.richness.estimates$jack1) qs.core.2011$std.jack1native<-rep(site.data$std.jack1native, each=15) qs.core.2011$jack1native<-rep(site.data$jack1native, each=15) ###find the mean and standard errors for the native and exotic esimtates at the the community scale mean(exotic.species.richness.estimates$jack1) mean(native.species.richness.estimates$jack1) std.error(exotic.species.richness.estimates$jack1) std.error(native.species.richness.estimates$jack1) ###find the mean and standard errors for the native and exotic esimtates at the the neighbourhood scale mean(qs.core.2011$native.richness) mean(qs.core.2011$exotic.richness) std.error(qs.core.2011$native.richness) std.error(qs.core.2011$exotic.richness) ##########to find the richness estimates at the metacommuntiy scale richness.estimators.metacommunity.exotic<-specpool(exotic.comm.data.2011, pool=qs.core.2011$yr) remnant.data$jack1remnant.exotic<-richness.estimators.metacommunity.exotic$jack1 richness.estimators.metacommunity.native<-specpool(native.comm.data.2011, pool=qs.core.2011$yr) remnant.data$jack1remnant.native<-richness.estimators.metacommunity.native$jack1 remnant.data$std.jack1remnant.exotic<-std.func(remnant.data$jack1remnant.exotic) remnant.data$std.jack1remnant.native<-std.func(remnant.data$jack1remnant.native) mean(richness.estimators.metacommunity.native$jack1) mean(richness.estimators.metacommunity.exotic$jack1) std.error(richness.estimators.metacommunity.native$jack1) std.error(richness.estimators.metacommunity.exotic$jack1) ####################### to find richness accumulation curves for each metacommunity rems<-levels(plants.core.2011$remnant) rem.temp<-with(qs.core.2011, tapply(mean.gstmin, remnant, mean)) rems.sorted.by.temp<-rems[order(rem.temp)] jpeg("Figure_S1_asymptotic_estimators.jpeg", width=7000, height=5000, quality=100, res=300) par(mfrow=c(3,4)) par(mar=c(5.1,6,4.1,2.1)) for (i in 1:length(rems.sorted.by.temp)){ rem.data<-droplevels(subset(plants.core.2011, remnant==rems.sorted.by.temp[i])) rem.natives<-unique(rem.data$species[rem.data$exotic==0]) rem.exotics<-unique(rem.data$species[rem.data$exotic==1]) rem.qs.data<-droplevels(subset(qs.core.2011, remnant==rems.sorted.by.temp[i])) rem.native.comm.data<-as.data.frame.matrix(with(rem.data, table(yrsq, species)))[,rem.natives] rem.native.accum<-specaccum(rem.native.comm.data, permutations=999) rem.native.upper<-rem.native.accum$richness + rem.native.accum$sd rem.native.lower<-rem.native.accum$richness - rem.native.accum$sd plot(c(0:45), seq(0,70, length.out=46), type="n", xlab="Quadrats", ylab="Species richness", xlim=c(0,45), cex.lab=2.5, cex.axis=2) plot.CI.func(x.for.plot=c(0:45), pred=c(0,rem.native.accum$richness), upper=c(0,rem.native.upper), lower=c(0,rem.native.lower), env.colour="black", env.trans=40, line.colour="black", line.weight=2, line.type=1) rem.exotic.comm.data<-as.data.frame.matrix(with(rem.data, table(yrsq, species)))[,rem.exotics] if(length(rem.exotics)>3){ rem.exotic.accum<-specaccum(rem.exotic.comm.data, permutations=999) rem.exotic.upper<-rem.exotic.accum$richness + rem.exotic.accum$sd rem.exotic.lower<-rem.exotic.accum$richness - rem.exotic.accum$sd plot.CI.func(x.for.plot=c(0:45), pred=c(0,rem.exotic.accum$richness), upper=c(0,rem.exotic.upper), lower=c(0,rem.exotic.lower), env.colour="goldenrod", env.trans=40, line.colour="goldenrod", line.weight=2, line.type=1) } } dev.off() #### to create the variable describing the distance between communities - remnant scale variable #find the mean x and y location of each quadrat within a site site.data$x.dist<-tapply(qs.core.2011$x_m, FUN=mean, qs.core.2011$yrs) site.data$y.dist<-tapply(qs.core.2011$y_m, FUN=mean, qs.core.2011$yrs) ###find the mean distance between communities for each remnant mean.distance<-vector("numeric", 12) for(i in 1:length(unique(site.data$yr))){ subset.site<-site.data[site.data$yr==unique(site.data$yr)[i],] mean.distance[i]<-mean(with(subset.site,(dist(cbind(x.dist,y.dist))))) } ###allocate this variable to the metacommunity scale remnant.data$mean.distance<-mean.distance #####for each remnant find the remnant-scale mean of each site variable means.of.nei.var.site<-with(site.data, data.frame(std.jack1native, std.jack1exotic, std.canopy.cover.disper, std.site.soil.nut.disper,std.site.phys.disper)) means.of.nei.var.melt.site<-melt(data.frame(yr=site.data$yr, means.of.nei.var.site)) means.of.nei.var.site<-cast(means.of.nei.var.melt.site, yr~variable, mean, na.rm=T) remnant.data<-data.frame(remnant.data, means.of.nei.var.site) ###standardise some variables that originated in the remnant dataset std.data.frame.remnant<-with(remnant.data, data.frame(mean.distance, mean.gstmin, remnant.exotic.richness, remnant.native.richness, gsai)) std.data.remnant<-data.frame(apply(std.data.frame.remnant, 2, std.func)) colnames(std.data.remnant)<-paste("std.", colnames(std.data.remnant), sep="") remnant.data<-data.frame(remnant.data, std.data.remnant) #####match all the required variables down from the remnant to the site and quadrat level ####required to do this for piecewiseSEM remnant.2.match<-with(remnant.data, data.frame(std.mean.distance, std.jack1remnant.native, std.jack1remnant.exotic, std.remnant.native.richness, std.remnant.exotic.richness, yr, std.mean.gstmin, std.gsai) ) qs.core.2011<-merge(qs.core.2011,remnant.2.match, "yr") site.data<-merge(site.data,remnant.2.match, "yr") #####match all the required variables down from the site to the quadrat level site.2.quadrat<-with(site.data,data.frame(std.canopy.cover.disper,std.site.exotic.richness,std.site.native.richness, yrs, std.site.soil.nut.disper,std.site.phys.disper)) qs.core.2011<-merge(qs.core.2011, site.2.quadrat, "yrs") ########override the neighbourhood scale value of distance to the remnant edge ########so that the d-sep analysis will use the site scale value qs.core.2011$std.dist.to.edge<-rep(site.data$std.dist.to.edge, each=15) #####diagnostics - should we use a negative binomial distribution? Ord_plot(qs.core.2011$exotic.richness)####looks good P>0.05 fit <- goodfit(qs.core.2011$exotic.richness, type="nbinom") summary(fit) #######goodness of fit to negative binomial indicates good fit rootogram(fit) #####good fit = when bars touch the line, doesn't look too bad ############## Now for the Analysis ################################ ########## define the component models for the native-driven model i.e. native-> exotic modlist.template.nativedriven=list( lm(std.jack1remnant.native~std.gsai+std.mean.gstmin+std.mean.distance, na.action=na.omit, remnant.data), lm(std.jack1remnant.exotic~std.jack1remnant.native+std.gsai+std.mean.gstmin+std.mean.distance, na.action=na.omit, remnant.data), lme(std.jack1native~std.jack1remnant.native+std.site.soil.nut.disper+std.site.phys.disper+ std.canopy.cover.disper, random=~1|yr, na.action=na.omit, site.data), lme(std.jack1exotic~std.dist.to.edge+std.jack1native+std.jack1remnant.exotic+std.site.soil.nut.disper+std.site.phys.disper+ std.canopy.cover.disper, random=~1|yr, na.action=na.omit, site.data), lme(std.native.richness~std.crust+std.ln.plants.per.quadrat+std.jack1native+std.log+std.sqrt.p+std.sqrt.k+std.mean.cover.n.w.e+std.litter.cover+std.relative.compact+std.ln.tot.n, random=~1|yr/yrs, na.action=na.omit, qs.core.2011), glmer.nb(exotic.richness~std.crust+std.native.richness+std.ln.plants.per.quadrat+std.jack1exotic+std.ln.tot.n+std.log+std.sqrt.p+std.sqrt.k+std.mean.cover.n.w.e+std.litter.cover+std.relative.compact+(1|yr/yrs), na.action=na.omit, qs.core.2011), lme(std.ln.plants.per.quadrat~std.crust+std.mean.gstmin+std.gsai+std.log+std.sqrt.p+std.sqrt.k+std.mean.cover.n.w.e+std.litter.cover+std.relative.compact+std.ln.tot.n, random=~1|yr/yrs, na.action=na.omit, qs.core.2011) ) ###retrieve the model coefficients from the component models sem.coefs(modlist.template.nativedriven) ###extract the r2 values for each model rsquared(modlist.template.nativedriven) ######diagnostics for the glmer model fittedModel<-glmer.nb(exotic.richness~std.crust+std.native.richness+std.ln.plants.per.quadrat+std.jack1exotic+std.ln.tot.n+std.log+std.sqrt.p+std.sqrt.k+std.mean.cover.n.w.e+std.litter.cover+std.relative.compact+(1|yr/yrs), na.action=na.omit, qs.core.2011) simulationOutput <- simulateResiduals(fittedModel = fittedModel, n = 250, use.u=T) ############we expect to see no deviation from qqnorm line and that residuals are flat if model fit perfect plotSimulatedResiduals(simulationOutput = simulationOutput) ##not bad, some deviation from horizantal is ok ###perform a d-sep test modtest<-sem.fit(modlist.template.nativedriven, qs.core.2011, grouping.vars=c("yrs","yr"), condition=T) ####at the time of writing piecewiseSEM is unable to handle 3-level hiearchies automatically ####need to manually override the d-sep tests for testing the effect of neighbourhood ####and site level exogenous variables on remnant native and exotic richness ####in other words, need to make sure that sem.fit is using 12 observations for ####neigbourhood scale exogenous variables rather than 36 modtest<-replace.path.func(modtest, 14, lm( std.jack1remnant.native ~ std.site.soil.nut.disper + std.gsai + std.mean.gstmin+std.mean.distance,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 16, lm( std.jack1remnant.exotic ~ std.site.soil.nut.disper + std.gsai + std.mean.gstmin+std.mean.distance + std.jack1remnant.native,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 19, lm( std.jack1remnant.native ~ std.site.phys.disper + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 21, lm( std.jack1remnant.exotic ~ std.site.phys.disper + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.native,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 24, lm( std.jack1remnant.native ~ std.canopy.cover.disper + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 26, lm( std.jack1remnant.exotic ~ std.canopy.cover.disper + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.native,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 29, lm( std.jack1remnant.native ~ std.dist.to.edge + std.gsai+std.mean.distance + std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 31, lm( std.jack1remnant.exotic ~ std.dist.to.edge + std.gsai+std.mean.distance + std.mean.gstmin + std.jack1remnant.native,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 35, lm( std.jack1remnant.native ~ std.crust + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 36, lm( std.jack1remnant.exotic ~ std.crust + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.native,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 39, lm( std.jack1remnant.native ~ std.log + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 40, lm( std.jack1remnant.exotic ~ std.log + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.native,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 43, lm( std.jack1remnant.native ~ std.sqrt.p + std.gsai+std.mean.distance + std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 44, lm( std.jack1remnant.exotic ~ std.sqrt.p + std.gsai+std.mean.distance + std.mean.gstmin + std.jack1remnant.native,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 47, lm( std.jack1remnant.native ~ std.sqrt.k + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 48, lm( std.jack1remnant.exotic ~ std.sqrt.k + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.native,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 51, lm( std.jack1remnant.native ~ std.mean.cover.n.w.e + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 52, lm( std.jack1remnant.exotic ~ std.mean.cover.n.w.e + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.native,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 55, lm( std.jack1remnant.native ~ std.litter.cover + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 56, lm( std.jack1remnant.exotic ~ std.litter.cover + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.native,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 59, lm( std.jack1remnant.native ~ std.relative.compact + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 60, lm( std.jack1remnant.exotic ~ std.relative.compact + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.native,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 63, lm( std.jack1remnant.native ~ std.ln.tot.n + std.gsai+std.mean.distance + std.mean.gstmin,na.action=na.omit, remnant.data)) modtest.nativedriven<-replace.path.func(modtest, 64, lm( std.jack1remnant.exotic ~ std.ln.tot.n + std.gsai+std.mean.distance + std.mean.gstmin + std.jack1remnant.native,na.action=na.omit, remnant.data)) ####view the results modtest.nativedriven ########## define the component models for the exotic-driven model i.e. exotic -> native modlist.template.exoticdriven=list( lm(std.jack1remnant.native~std.gsai+std.mean.gstmin+std.jack1remnant.exotic+std.mean.distance, na.action=na.omit, remnant.data), lm(std.jack1remnant.exotic~std.gsai+std.mean.gstmin+std.mean.distance, na.action=na.omit, remnant.data), lme(std.jack1native~std.jack1exotic+std.jack1remnant.native+std.site.soil.nut.disper+std.site.phys.disper+std.canopy.cover.disper, random=~1|yr, na.action=na.omit, site.data), lme(std.jack1exotic~std.dist.to.edge+std.jack1remnant.exotic+std.site.soil.nut.disper+std.site.phys.disper+ std.canopy.cover.disper, random=~1|yr, na.action=na.omit, site.data), lme(std.native.richness~exotic.richness+std.crust+std.ln.plants.per.quadrat+std.jack1native+std.ln.tot.n+std.log+std.sqrt.p+std.sqrt.k+std.mean.cover.n.w.e+std.litter.cover+std.relative.compact, random=~1|yr/yrs, na.action=na.omit, qs.core.2011), glmer.nb(exotic.richness~std.crust+std.ln.plants.per.quadrat+std.jack1exotic+std.ln.tot.n+std.log+std.sqrt.p+std.sqrt.k+std.mean.cover.n.w.e+std.relative.compact+std.litter.cover+(1|yr/yrs), na.action=na.omit, qs.core.2011), lme(std.ln.plants.per.quadrat~std.crust+std.mean.gstmin+std.gsai+std.log+std.sqrt.p+std.sqrt.k+std.mean.cover.n.w.e+std.litter.cover+std.relative.compact+std.ln.tot.n, random=~1|yr/yrs, na.action=na.omit, qs.core.2011) ) sem.coefs(modlist.template.exoticdriven) ######diagnostics for the glmer model fittedModel<-glmer.nb(exotic.richness~std.crust+std.ln.plants.per.quadrat+std.jack1exotic+std.ln.tot.n+std.log+std.sqrt.p+std.sqrt.k+std.mean.cover.n.w.e+std.relative.compact+std.litter.cover+(1|yr/yrs), na.action=na.omit, qs.core.2011) simulationOutput <- simulateResiduals(fittedModel = fittedModel, n = 250, use.u=T) ############we expect to see no deviation from qqnorm line and that residuals are flat if model fit perfect plotSimulatedResiduals(simulationOutput = simulationOutput) ##not bad, some deviation from horizantal is ok ###conduct d-sep tests modtest<-sem.fit(modlist.template.exoticdriven, qs.core.2011, grouping.vars=c("yrs","yr"), condition=T) ####at the time of writing piecewiseSEM is unable to handle 3-level hiearchies automatically ####need to manually override the d-sep tests for testing the effect of neighbourhood ####and site level exogenous variables on remnant native and exotic richness####in other words, need to make sure that sem.fit is using 12 observations for ####neigbourhood scale exogenous variables rather than 36 modtest<-replace.path.func(modtest, 14, lm( std.jack1remnant.exotic ~ std.site.soil.nut.disper + std.gsai + std.mean.gstmin+std.mean.distance,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest,16, lm(std.jack1remnant.native ~ std.site.soil.nut.disper+std.gsai + std.mean.gstmin+std.mean.distance + std.jack1remnant.exotic,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 19, lm( std.jack1remnant.exotic ~ std.site.phys.disper + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 21, lm( std.jack1remnant.native ~ std.site.phys.disper + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.exotic,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 24, lm( std.jack1remnant.exotic ~ std.canopy.cover.disper + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 26, lm( std.jack1remnant.native ~ std.canopy.cover.disper + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.exotic,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 29, lm( std.jack1remnant.exotic ~ std.dist.to.edge + std.gsai+std.mean.distance + std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 31, lm( std.jack1remnant.native ~ std.dist.to.edge + std.gsai+std.mean.distance + std.mean.gstmin + std.jack1remnant.exotic,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 35, lm( std.jack1remnant.exotic ~ std.crust + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 36, lm( std.jack1remnant.native ~ std.crust + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.exotic,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 39, lm( std.jack1remnant.exotic ~ std.ln.tot.n + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 40, lm( std.jack1remnant.native ~ std.ln.tot.n + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.exotic,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 43, lm( std.jack1remnant.exotic ~ std.log + std.gsai+std.mean.distance + std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 44, lm( std.jack1remnant.native ~ std.log + std.gsai+std.mean.distance + std.mean.gstmin + std.jack1remnant.exotic,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 47, lm( std.jack1remnant.exotic ~ std.sqrt.p + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 48, lm( std.jack1remnant.native ~ std.sqrt.p + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.exotic,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 51, lm( std.jack1remnant.exotic ~ std.sqrt.k + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 52, lm( std.jack1remnant.native ~ std.sqrt.k + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.exotic,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 55, lm( std.jack1remnant.exotic ~ std.mean.cover.n.w.e + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 56, lm( std.jack1remnant.native ~ std.mean.cover.n.w.e + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.exotic,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 59, lm( std.jack1remnant.exotic ~ std.litter.cover + std.gsai +std.mean.distance+ std.mean.gstmin,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 60, lm( std.jack1remnant.native ~ std.litter.cover + std.gsai +std.mean.distance+ std.mean.gstmin + std.jack1remnant.exotic,na.action=na.omit, remnant.data)) modtest<-replace.path.func(modtest, 63, lm( std.jack1remnant.exotic ~ std.relative.compact + std.gsai+std.mean.distance + std.mean.gstmin,na.action=na.omit, remnant.data)) modtest.exotic.driven<-replace.path.func(modtest, 64, lm( std.jack1remnant.native ~ std.relative.compact + std.gsai+std.mean.distance + std.mean.gstmin + std.jack1remnant.exotic,na.action=na.omit, remnant.data)) ####view the results modtest.exotic.driven ##### this is the best model, lets extract the coefficients and variation explained for each component model (Table S1) ###have to re-run the model with std.exotic.richness to extract standardised coefficients modlist.template.exoticdriven=list( lm(std.jack1remnant.native~std.gsai+std.mean.gstmin+std.jack1remnant.exotic+std.mean.distance, na.action=na.omit, remnant.data), lm(std.jack1remnant.exotic~std.gsai+std.mean.gstmin+std.mean.distance, na.action=na.omit, remnant.data), lme(std.jack1native~std.jack1exotic+std.jack1remnant.native+std.site.soil.nut.disper+std.site.phys.disper+std.canopy.cover.disper, random=~1|yr, na.action=na.omit, site.data), lme(std.jack1exotic~std.dist.to.edge+std.jack1remnant.exotic+std.site.soil.nut.disper+std.site.phys.disper+ std.canopy.cover.disper, random=~1|yr, na.action=na.omit, site.data), lme(std.native.richness~std.exotic.richness+std.crust+std.ln.plants.per.quadrat+std.jack1native+std.ln.tot.n+std.log+std.sqrt.p+std.sqrt.k+std.mean.cover.n.w.e+std.litter.cover+std.relative.compact, random=~1|yr/yrs, na.action=na.omit, qs.core.2011), glmer.nb(exotic.richness~std.crust+std.ln.plants.per.quadrat+std.jack1exotic+std.ln.tot.n+std.log+std.sqrt.p+std.sqrt.k+std.mean.cover.n.w.e+std.relative.compact+std.litter.cover+(1|yr/yrs), na.action=na.omit, qs.core.2011), lme(std.ln.plants.per.quadrat~std.crust+std.mean.gstmin+std.gsai+std.log+std.sqrt.p+std.sqrt.k+std.mean.cover.n.w.e+std.litter.cover+std.relative.compact+std.ln.tot.n, random=~1|yr/yrs, na.action=na.omit, qs.core.2011) ) sem.coefs(modlist.template.exoticdriven) ##### this is the best model, lets extract the variation explained for each component model rsquared(modlist.template.exoticdriven) ####remove all figures to this point dev.off() ###pair-wise correlations - to find estimates and p values ##metacommunity summary(lm(remnant.data$std.jack1remnant.native~remnant.data$std.mean.gstmin)) summary(lm(remnant.data$std.jack1remnant.native~remnant.data$std.gsai)) summary(lm(remnant.data$std.jack1remnant.native~remnant.data$std.jack1remnant.exotic)) summary(lm(remnant.data$std.jack1remnant.native~remnant.data$std.mean.distance)) summary(lm(remnant.data$std.jack1remnant.exotic~remnant.data$std.mean.gstmin)) summary(lm(remnant.data$std.jack1remnant.exotic~remnant.data$std.gsai)) summary(lm(remnant.data$std.jack1remnant.exotic~remnant.data$std.jack1remnant.native)) summary(lm(remnant.data$std.jack1remnant.exotic~remnant.data$std.mean.distance)) ##community summary((lme(std.jack1native~std.jack1remnant.native, random=~1|yr, site.data))) summary(lme(std.jack1native~std.site.phys.disper, random=~1|yr, site.data)) summary(lme(std.jack1native~std.site.soil.nut.disper, random=~1|yr, site.data)) summary(lme(std.jack1native~std.canopy.cover.disper, random=~1|yr, site.data)) summary(lme(std.jack1native~std.jack1exotic, random=~1|yr, site.data)) summary((lme(std.jack1exotic~std.jack1remnant.exotic, random=~1|yr, site.data))) summary((lme(std.jack1exotic~std.site.phys.disper, random=~1|yr, site.data))) summary(lme(std.jack1exotic~std.site.soil.nut.disper, random=~1|yr, site.data)) summary(lme(std.jack1exotic~std.canopy.cover.disper, random=~1|yr, site.data)) summary(lme(std.jack1exotic~std.dist.to.edge, random=~1|yr, site.data)) summary(lme(std.jack1exotic~std.jack1native, random=~1|yr, site.data)) ##neighbourhood summary(lme(std.native.richness~std.jack1native, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.native.richness~std.ln.plants.per.quadrat, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.native.richness~std.sqrt.k, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.native.richness~std.sqrt.p, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.native.richness~std.ln.tot.n, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.native.richness~std.log, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.native.richness~std.crust, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.native.richness~std.mean.cover.n.w.e, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.native.richness~std.litter.cover, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.native.richness~std.relative.compact, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.native.richness~std.exotic.richness, random=~1|yr/yrs, qs.core.2011)) summary(glmer.nb(exotic.richness~std.jack1exotic+(1|yr/yrs), qs.core.2011)) summary(glmer.nb(exotic.richness~std.ln.plants.per.quadrat+(1|yr/yrs), qs.core.2011)) summary(glmer.nb(exotic.richness~std.sqrt.k+(1|yr/yrs), qs.core.2011)) summary(glmer.nb(exotic.richness~std.sqrt.p+(1|yr/yrs), qs.core.2011)) summary(glmer.nb(exotic.richness~std.ln.tot.n+(1|yr/yrs), qs.core.2011)) summary(glmer.nb(exotic.richness~std.log+(1|yr/yrs), qs.core.2011)) summary(glmer.nb(exotic.richness~std.crust+(1|yr/yrs), qs.core.2011)) summary(glmer.nb(exotic.richness~std.mean.cover.n.w.e+(1|yr/yrs), qs.core.2011)) summary(glmer.nb(exotic.richness~std.litter.cover+(1|yr/yrs), qs.core.2011)) summary(glmer.nb(exotic.richness~std.relative.compact+(1|yr/yrs), qs.core.2011)) summary(glmer.nb(exotic.richness~std.native.richness+(1|yr/yrs), qs.core.2011)) ##plant density summary(lme(std.ln.plants.per.quadrat~std.sqrt.p, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.ln.plants.per.quadrat~std.sqrt.k, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.ln.plants.per.quadrat~std.ln.tot.n, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.ln.plants.per.quadrat~std.log, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.ln.plants.per.quadrat~std.crust, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.ln.plants.per.quadrat~std.mean.cover.n.w.e, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.ln.plants.per.quadrat~std.litter.cover, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.ln.plants.per.quadrat~std.relative.compact, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.ln.plants.per.quadrat~std.mean.gstmin, random=~1|yr/yrs, qs.core.2011)) summary(lme(std.ln.plants.per.quadrat~std.gsai, random=~1|yr/yrs, qs.core.2011)) #### ####to find pair-wise correlations, using marginal ##metacommunity rsquared(lm(remnant.data$std.jack1remnant.native~remnant.data$std.mean.gstmin)) rsquared(lm(remnant.data$std.jack1remnant.native~remnant.data$std.gsai)) rsquared(lm(remnant.data$std.jack1remnant.native~remnant.data$std.jack1remnant.exotic)) rsquared(lm(remnant.data$std.jack1remnant.native~remnant.data$std.mean.distance)) rsquared(lm(remnant.data$std.jack1remnant.exotic~remnant.data$std.mean.gstmin)) rsquared(lm(remnant.data$std.jack1remnant.exotic~remnant.data$std.gsai)) rsquared(lm(remnant.data$std.jack1remnant.exotic~remnant.data$std.jack1remnant.native)) rsquared(lm(remnant.data$std.jack1remnant.exotic~remnant.data$std.mean.distance)) ##community rsquared((lme(std.jack1native~std.jack1remnant.native, random=~1|yr, site.data))) rsquared(lme(std.jack1native~std.site.phys.disper, random=~1|yr, site.data)) rsquared(lme(std.jack1native~std.site.soil.nut.disper, random=~1|yr, site.data)) rsquared(lme(std.jack1native~std.canopy.cover.disper, random=~1|yr, site.data)) rsquared(lme(std.jack1native~std.jack1exotic, random=~1|yr, site.data)) rsquared((lme(std.jack1exotic~std.jack1remnant.exotic, random=~1|yr, site.data))) rsquared((lme(std.jack1exotic~std.site.phys.disper, random=~1|yr, site.data))) rsquared(lme(std.jack1exotic~std.site.soil.nut.disper, random=~1|yr, site.data)) rsquared(lme(std.jack1exotic~std.canopy.cover.disper, random=~1|yr, site.data)) rsquared(lme(std.jack1exotic~std.dist.to.edge, random=~1|yr, site.data)) rsquared(lme(std.jack1exotic~std.jack1native, random=~1|yr, site.data)) ##neighbourhood rsquared(lme(std.native.richness~std.jack1native, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.native.richness~std.ln.plants.per.quadrat, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.native.richness~std.sqrt.k, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.native.richness~std.sqrt.p, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.native.richness~std.ln.tot.n, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.native.richness~std.log, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.native.richness~std.crust, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.native.richness~std.mean.cover.n.w.e, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.native.richness~std.litter.cover, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.native.richness~std.relative.compact, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.native.richness~std.exotic.richness, random=~1|yr/yrs, qs.core.2011)) rsquared(glmer.nb(exotic.richness~std.jack1exotic+(1|yr/yrs), qs.core.2011)) rsquared(glmer.nb(exotic.richness~std.ln.plants.per.quadrat+(1|yr/yrs), qs.core.2011)) rsquared(glmer.nb(exotic.richness~std.sqrt.k+(1|yr/yrs), qs.core.2011)) rsquared(glmer.nb(exotic.richness~std.sqrt.p+(1|yr/yrs), qs.core.2011)) rsquared(glmer.nb(exotic.richness~std.ln.tot.n+(1|yr/yrs), qs.core.2011)) rsquared(glmer.nb(exotic.richness~std.log+(1|yr/yrs), qs.core.2011)) rsquared(glmer.nb(exotic.richness~std.crust+(1|yr/yrs), qs.core.2011)) rsquared(glmer.nb(exotic.richness~std.mean.cover.n.w.e+(1|yr/yrs), qs.core.2011)) rsquared(glmer.nb(exotic.richness~std.litter.cover+(1|yr/yrs), qs.core.2011)) rsquared(glmer.nb(exotic.richness~std.relative.compact+(1|yr/yrs), qs.core.2011)) rsquared(glmer.nb(exotic.richness~std.native.richness+(1|yr/yrs), qs.core.2011)) ##plant density rsquared(lme(std.ln.plants.per.quadrat~std.sqrt.p, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.ln.plants.per.quadrat~std.sqrt.k, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.ln.plants.per.quadrat~std.ln.tot.n, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.ln.plants.per.quadrat~std.log, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.ln.plants.per.quadrat~std.crust, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.ln.plants.per.quadrat~std.mean.cover.n.w.e, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.ln.plants.per.quadrat~std.litter.cover, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.ln.plants.per.quadrat~std.relative.compact, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.ln.plants.per.quadrat~std.mean.gstmin, random=~1|yr/yrs, qs.core.2011)) rsquared(lme(std.ln.plants.per.quadrat~std.gsai, random=~1|yr/yrs, qs.core.2011)) ############################### TO MAKE THE FIGURES ############################# ####################TO MAKE THE RAW NERR PLOTS################################ pdf(file="Figure_4_rawrelationships.pdf", width=11, height=3.9, useDingbats=F) par(mfrow=c(1,3), mar=c(4.5, 4, 4, 2), mgp=c(2.6,0.8,0)) ###neighbourhood level neighb.nerr<-(lme(native.richness~exotic.richness, random=~exotic.richness|yr/yrs, qs.core.2011, na.action=na.omit)) with(qs.core.2011, plot(native.richness ~ exotic.richness, type="n", ylab="Neighbourhood native richness", xlab="Neighbourhood exotic richness",cex.axis=1.1, cex.lab=1.4)) y.for.plot<-qs.core.2011$native.richness x.for.plot<-qs.core.2011$exotic.richness mod.for.plot<-lme(y.for.plot~x.for.plot, random=~exotic.richness|yr/yrs, qs.core.2011, na.action=na.omit) pred.data<-lme.predict(mod.for.plot, data.frame(y.for.plot=seq.func(y.for.plot), x.for.plot=seq.func(x.for.plot)), int.type = "CI", se.mult=1.96) ######### plot the fitted line and CI for the relationship between partials plot.CI.func(x.for.plot=seq.func(x.for.plot), pred=pred.data$pred, upper=pred.data$upper, lower=pred.data$lower, env.colour="black", env.trans=40, line.colour="black", line.type=1, line.weight=1) with(qs.core.2011, points(jitter(exotic.richness, 0.5), jitter(native.richness, 1), col="grey", pch=19, cex=1.25)) site.data$raw.n.v.e.int.site<-NA site.data$raw.n.v.e.slope.site<-NA for(i in 1:dim(site.data)[1]){ site.quadrats<-droplevels(subset(qs.core.2011, yrs==site.data$yrs[i])) curve(cbind(1,x)%*%as.numeric(coef(neighb.nerr)[i,]), add=T, col="grey40", from=min(site.quadrats$exotic.richness), to=max(site.quadrats$exotic.richness)) site.data$raw.n.v.e.int.site[i]<-coef(neighb.nerr)[i,1] site.data$raw.n.v.e.slope.site[i]<-coef(neighb.nerr)[i,2] } curve(cbind(1,x)%*%fixef(neighb.nerr), add=T, lwd=2.5) mtext("(a)",line=0.5, side=3, adj=0, cex=1.25) ####community level plot ##model not converging with 10 interactions so increase iteration ctrl <- lmeControl(opt='optim', maxIter= 300, msMaxIter = 400) community.nerr<-(lme(jack1native~jack1exotic, random=~jack1exotic|yr, site.data, control=ctrl,na.action=na.omit)) with(site.data, plot(jack1native~jack1exotic, type="n", ylab="Community native richness", xlab="Community exotic richness", cex.axis=1.1, cex.lab=1.4)) y.for.plot<-site.data$jack1native x.for.plot<-site.data$jack1exotic mod.for.plot<-(lme(y.for.plot~x.for.plot, random=~x.for.plot|yr, site.data, control=ctrl,na.action=na.omit)) pred.data<-lme.predict(mod.for.plot, data.frame(y.for.plot=seq.func(y.for.plot), x.for.plot=seq.func(x.for.plot)), int.type = "CI", se.mult=1.96) plot.CI.func(x.for.plot=seq.func(x.for.plot), pred=pred.data$pred, upper=pred.data$upper, lower=pred.data$lower, env.colour="black", env.trans=40, line.colour="black", line.type=1, line.weight=1) with(site.data, points(jack1native~jack1exotic, col="grey", pch=19, cex=1.25)) for(i in 1:dim(remnant.data)[1]){ remnant.sites<-droplevels(subset(site.data, yr==remnant.data$yr[i])) curve(cbind(1,x)%*%as.numeric(coef(community.nerr)[i,]), add=T, col="grey40", from=min(remnant.sites$jack1exotic)-0.5, to=max(remnant.sites$jack1exotic)+0.5) } curve(cbind(1,x)%*%fixef(community.nerr), add=T, lwd=2.5) mtext("(b)", line=0.5, side=3, adj=0, cex=1.25) ######remnant level metacommunity.nerr<-(lm(jack1remnant.native~jack1remnant.exotic, remnant.data,na.action=na.omit)) with(remnant.data, plot(jack1remnant.native~jack1remnant.exotic, type="n", ylab="Metacommunity native richness", xlab="Metacommunity exotic richness", cex.axis=1.1, cex.lab=1.4)) y.for.plot<-remnant.data$jack1remnant.native x.for.plot<-remnant.data$jack1remnant.exotic mod.for.plot<-(lm(y.for.plot~x.for.plot, site.data,na.action=na.omit)) pred.data<-as.data.frame(predict(object=mod.for.plot, newdat=data.frame(y.for.plot=seq.func(y.for.plot), x.for.plot=seq.func(x.for.plot)), interval = "confidence", level = 0.95)) plot.CI.func(x.for.plot=seq.func(x.for.plot), pred=pred.data$fit, upper=pred.data$upr, lower=pred.data$lwr, env.colour="black", env.trans=40, line.colour="black", line.type=1, line.weight=1) with(remnant.data, points(jack1remnant.native~jack1remnant.exotic, col="grey", pch=19, cex=1.25)) curve(cbind(1,x)%*%coef(metacommunity.nerr), add=T, lwd=2.5) mtext("(c)", line=0.5, side=3, adj=0, cex=1.25) dev.off() ########same, but with exotics as the response####### pdf(file="Figure_S2_rawrelationships_exotic_response.pdf", width=11, height=3.9, useDingbats=F) par(mfrow=c(1,3), mar=c(4.5, 4, 4, 2), mgp=c(2.6,0.8,0)) ###neighbourhood level neighb.nerr<-glmer.nb(exotic.richness ~ native.richness + (native.richness|yrs:yr) + (1|yr), qs.core.2011, na.action=na.omit) summary(neighb.nerr) with(qs.core.2011, plot(exotic.richness ~ native.richness, type="n", xlab="Neighbourhood native richness", ylab="Neighbourhood exotic richness",cex.axis=1.1, cex.lab=1.4)) x.for.plot<-qs.core.2011$native.richness y.for.plot<-qs.core.2011$exotic.richness mod.for.plot<-glmer.nb(y.for.plot~x.for.plot + (x.for.plot|yrs:yr) + (1|yr), qs.core.2011, na.action=na.omit) pred.data<-data.frame(intercept=rep(1, 100), x.for.plot=seq.func(x.for.plot)) pred<-lmer.predict(mod.for.plot, newdat=pred.data, se.mult=1.96, binom=F, poisson=T) plot.CI.func(x.for.plot=seq.func(x.for.plot), pred=pred$y, upper=pred$phi, lower=pred$plo, env.colour="black", env.trans=40, line.colour="black", line.type=1, line.weight=1) with(qs.core.2011, points(jitter(exotic.richness, 0.5) ~ jitter(native.richness, 1), col="grey", pch=19, cex=1.25)) for(i in 1:dim(site.data)[1]){ site.quadrats<-droplevels(subset(qs.core.2011, yrs==site.data$yrs[i])) curve(exp(cbind(1,x)%*%as.numeric(coef(neighb.nerr)$yrs[i,])), add=T, col="grey40", from=min(site.quadrats$native.richness), to=max(site.quadrats$native.richness)) } curve(exp(cbind(1,x)%*%fixef(neighb.nerr)), add=T, lwd=2.5) mtext("(a)", line=0.5, side=3, adj=0, cex=1.25) ####community level plot ##model not converging with 10 interations so increase iteration ctrl <- lmeControl(opt='optim', maxIter= 300, msMaxIter = 400) community.nerr<-(lme(jack1exotic~jack1native, random=~jack1native|yr, site.data, control=ctrl,na.action=na.omit)) with(site.data, plot(jack1exotic~jack1native, type="n", xlab="Community native richness", ylab="Community exotic richness", cex.axis=1.1, cex.lab=1.4)) x.for.plot<-site.data$jack1native y.for.plot<-site.data$jack1exotic mod.for.plot<-(lme(y.for.plot~x.for.plot, random=~x.for.plot|yr, site.data, control=ctrl,na.action=na.omit)) pred.data<-lme.predict(mod.for.plot, data.frame(y.for.plot=seq.func(y.for.plot), x.for.plot=seq.func(x.for.plot)), int.type = "CI", se.mult=1.96) plot.CI.func(x.for.plot=seq.func(x.for.plot), pred=pred.data$pred, upper=pred.data$upper, lower=pred.data$lower, env.colour="black", env.trans=40, line.colour="black", line.type=1, line.weight=1) with(site.data, points(jack1native, jack1exotic, col="grey", pch=19, cex=1.25)) for(i in 1:dim(remnant.data)[1]){ remnant.sites<-droplevels(subset(site.data, yr==remnant.data$yr[i])) curve(cbind(1,x)%*%as.numeric(coef(community.nerr)[i,]), add=T, col="grey40", from=min(remnant.sites$jack1native)-0.5, to=max(remnant.sites$jack1native)+0.5) } curve(cbind(1,x)%*%fixef(community.nerr), add=T, lwd=2.5) mtext("(b)", line=0.5, side=3, adj=0, cex=1.25) ######remnant level metacommunity.nerr<-(lm(jack1remnant.exotic~jack1remnant.native, remnant.data,na.action=na.omit)) with(remnant.data, plot(jack1remnant.exotic~jack1remnant.native, type="n", xlab="Metacommunity native richness", ylab="Metacommunity exotic richness", cex.axis=1.1, cex.lab=1.4)) x.for.plot<-remnant.data$jack1remnant.native y.for.plot<-remnant.data$jack1remnant.exotic mod.for.plot<-(lm(y.for.plot~x.for.plot, site.data,na.action=na.omit)) pred.data<-as.data.frame(predict(object=mod.for.plot, newdat=data.frame(y.for.plot=seq.func(y.for.plot), x.for.plot=seq.func(x.for.plot)), interval = "confidence", level = 0.95)) plot.CI.func(x.for.plot=seq.func(x.for.plot), pred=pred.data$fit, upper=pred.data$upr, lower=pred.data$lwr, env.colour="black", env.trans=40, line.colour="black", line.type=1, line.weight=1) with(remnant.data, points(jack1remnant.native, jack1remnant.exotic, col="grey", pch=19, cex=1.25)) curve(cbind(1,x)%*%coef(metacommunity.nerr), add=T, lwd=2.5) mtext("(c)", line=0.5, side=3, adj=0, cex=1.25) dev.off() #################FIGURE 6 #####prints the figure to a jpeg in your working directory pdf(file="Figure_6_Partialplotsmetacommunity.pdf",height=9,width=8.5,useDingbats=F) par(mfrow=c(2,2), mar=c(4.5, 4, 4, 2), mgp=c(2.6,0.8,0)) ########FIGURE 6A ######## #####extract partial residuals for the response and predictor variable r.n_Tmin<-residuals(lm(std.jack1remnant.native~std.jack1remnant.exotic+std.gsai+std.mean.distance, remnant.data)) Tmin_n<-residuals(lm(std.mean.gstmin~std.jack1remnant.exotic+std.gsai+std.mean.distance, remnant.data)) y.for.plot<-r.n_Tmin x.for.plot<-Tmin_n ########## plots residuals against each other expression('Mean'~'growing'~'season'~'Tmin'~degrees~C~''|''~'(others)') plot(y.for.plot~x.for.plot, ylab="Metacommunity native richness | (Others)", xlab=expression('Mean'~'growing'~'season'~'Tmin'~'('*degree*'C'*')'~'| (Others)'),pch=19, col="grey", cex=1.25,cex.axis=1, cex.lab=1.25, ylim=c(-1.7,1.6)) #########extract coefficients for the relationship between partial residuals mod.for.plot<-lm(y.for.plot~x.for.plot, na.action = na.exclude) pred.data<-lm.predict(mod.for.plot, data.frame(y.for.plot=seq.func(y.for.plot), x.for.plot=seq.func(x.for.plot))) ######### plot the fitted line and CI for the relationship between partials plot.CI.func(x.for.plot=seq.func(x.for.plot), pred=pred.data$pred, upper=pred.data$upper, lower=pred.data$lower, env.colour="black", env.trans=40, line.colour="black", line.type=1,line.weight=3) mtext("(a)", line=0.5, side=3, adj=0, cex=1.5) ########FIGURE 6B ######## #####extract partial residuals for the response and predictor variable r.n_GSAI<-residuals(lm(std.jack1remnant.native~std.jack1remnant.exotic+std.mean.gstmin+std.mean.distance, remnant.data)) GSAI_n<-residuals(lm(std.gsai~std.jack1remnant.exotic+std.mean.gstmin+std.mean.distance, remnant.data)) y.for.plot<-r.n_GSAI x.for.plot<-GSAI_n ########## plots residuals against each other plot(y.for.plot~x.for.plot, ylab="Metacommunity native richness | (Others)", xlab="Growing season moisture availability | (Others)",pch=19, col="grey", cex=1.25,cex.axis=1, cex.lab=1.25, ylim=c(-1.2,1.6)) #########extract coefficients for the relationship between partial residuals mod.for.plot<-lm(y.for.plot~x.for.plot, na.action = na.exclude) pred.data<-lm.predict(mod.for.plot, data.frame(y.for.plot=seq.func(y.for.plot), x.for.plot=seq.func(x.for.plot))) ######### plot the fitted line and CI for the relationship between partials plot.CI.func(x.for.plot=seq.func(x.for.plot), pred=pred.data$pred, upper=pred.data$upper, lower=pred.data$lower, env.colour="black", env.trans=40, line.colour="black", line.type=1,line.weight=3) mtext("(b)", line=0.5, side=3, adj=0, cex=1.5) ########FIGURE 6C ######## #####extract partial residuals for the response and predictor variable r.e_Tmin<-residuals(lm(std.jack1remnant.exotic~std.gsai+std.mean.distance, remnant.data)) Tmin_e<-residuals(lm(std.mean.gstmin~std.gsai+std.mean.distance, remnant.data)) y.for.plot<-r.e_Tmin x.for.plot<-Tmin_e ########## plots residuals against each other plot(y.for.plot~x.for.plot, ylab="Metacommunity exotic richness | (Others)", xlab=expression('Mean'~'growing'~'season'~'Tmin'~'('*degree*'C'*')'~'| (Others)'),pch=19, col="grey", cex=1.25,cex.axis=1, cex.lab=1.25, ylim=c(-1.2,2.2)) #########extract coefficients for the relationship between partial residuals mod.for.plot<-lm(y.for.plot~x.for.plot, na.action = na.exclude) pred.data<-lm.predict(mod.for.plot, data.frame(y.for.plot=seq.func(y.for.plot), x.for.plot=seq.func(x.for.plot))) ######### plot the fitted line and CI for the relationship between partials plot.CI.func(x.for.plot=seq.func(x.for.plot), pred=pred.data$pred, upper=pred.data$upper, lower=pred.data$lower, env.colour="black", env.trans=40, line.colour="black", line.type=1,line.weight=3) mtext("(c)", line=0.5, side=3, adj=0, cex=1.5) ########FIGURE 6D ######## #####extract partial residuals for the response and predictor variable r.e_GSAI<-residuals(lm(std.jack1remnant.exotic~std.mean.gstmin+std.mean.distance, remnant.data)) GSAI_e<-residuals(lm(std.gsai~std.mean.gstmin+std.mean.distance, remnant.data)) y.for.plot<-r.e_GSAI x.for.plot<-GSAI_e ########## plots residuals against each other plot(y.for.plot~x.for.plot, ylab="Metacommunity exotic richness | (Others)", xlab="Growing season moisture availability | (Others)",pch=19, col="grey", cex=1.25,cex.axis=1, cex.lab=1.25) #########extract coefficients for the relationship between partial residuals mod.for.plot<-lm(y.for.plot~x.for.plot, na.action = na.exclude) pred.data<-lm.predict(mod.for.plot, data.frame(y.for.plot=seq.func(y.for.plot), x.for.plot=seq.func(x.for.plot))) ######### plot the fitted line and CI for the relationship between partials plot.CI.func(x.for.plot=seq.func(x.for.plot), pred=pred.data$pred, upper=pred.data$upper, lower=pred.data$lower, env.colour="black", env.trans=40, line.colour="black", line.type=1,line.weight=3) mtext("(d)", line=0.5, side=3, adj=0, cex=1.5) dev.off() #################FIGURE 7 pdf(file="Figure_7_remaining_relationships.pdf", width=11, height=3.9, useDingbats=F) par(mfrow=c(1,3), mar=c(4.5, 4, 4, 2), mgp=c(2.6,0.8,0)) #### Neighborhood scale ####extract the residuals for the response variabe (Native richenss) #####cannot extract the residuals from a negative binomial for our purposes so we will ####plot exotic richenss as a raw variable n.r.residuals<-residuals(lme(std.native.richness~std.crust+std.ln.plants.per.quadrat+std.jack1native+std.ln.tot.n+std.log+std.sqrt.p+std.sqrt.k+std.mean.cover.n.w.e+std.litter.cover+std.relative.compact, random=~1|yr/yrs, qs.core.2011), level=0) y.for.plot<-n.r.residuals x.for.plot<-qs.core.2011$exotic.richness #######plot the effect of exotic richness on native richness residuals plot(y.for.plot~jitter(x.for.plot, amount=0.15), ylab="Neighbourhood native richness | (Others)", xlab="Neighbourhood exotic richness",pch=19, col="grey",cex=1.25, cex.axis=1.1, cex.lab=1.4) ####### model the effect of exotic richness on native richness residuals mod.for.plot<-lme(y.for.plot~x.for.plot, na.action = na.exclude, random=~1|yr/yrs, qs.core.2011) ######plot the fitted line and CI pred.data<-lme.predict(mod.for.plot, data.frame(y.for.plot=seq.func(y.for.plot), x.for.plot=seq.func(x.for.plot)), int.type = "CI", se.mult=1.96) plot.CI.func(x.for.plot=seq.func(x.for.plot), pred=pred.data$pred, upper=pred.data$upper, lower=pred.data$lower, env.colour="black", env.trans=40, line.colour="black", line.type=1, line.weight=2.5) mtext("(a)", line=0.5, side=3, adj=0, cex=1.25) ######extract the residuals at the community scale ##### community scale s.n.s.e<-residuals(lme(std.jack1native~std.jack1remnant.native+std.site.soil.nut.disper+ std.canopy.cover.disper+std.site.phys.disper, random=~1|yr, site.data), level=0) s.e.s.n<-residuals(lme(std.jack1exotic~std.jack1remnant.exotic+std.site.soil.nut.disper+ std.canopy.cover.disper+std.site.phys.disper+std.dist.to.edge, random=~1|yr, site.data), level=0) y.for.plot<-s.n.s.e x.for.plot<-s.e.s.n #####plot the residuals against each other plot(y.for.plot~x.for.plot, ylab="Community native richness | (Others)", xlab="Community exotic richness | (Others)",pch=19, col="grey",cex=1.25,cex.axis=1.1, cex.lab=1.4, ylim=c(-2.25,2)) ####model the relationship between residuals mod.for.plot<-lme(y.for.plot~x.for.plot, na.action = na.exclude, random=~1|yr, site.data) pred.data<-lme.predict(mod.for.plot, data.frame(y.for.plot=seq.func(y.for.plot), x.for.plot=seq.func(x.for.plot)), int.type = "CI", se.mult=1.96) ######### plot the fitted line and CI for the relationship between partials plot.CI.func(x.for.plot=seq.func(x.for.plot), pred=pred.data$pred, upper=pred.data$upper, lower=pred.data$lower, env.colour="black", env.trans=40, line.colour="black", line.type=1, line.weight=2.5) mtext("(b)", line=0.5, side=3, adj=0, cex=1.25) #####extract the residuals for the diversity variables at metacommunity scale ## metacommunity scale r.n.r.e<-residuals(lm(std.jack1remnant.native~std.gsai+std.mean.gstmin+std.mean.distance , remnant.data)) r.e.n.r<-residuals(lm(std.jack1remnant.exotic~std.gsai+std.mean.gstmin+std.mean.distance , remnant.data)) y.for.plot<-r.n.r.e x.for.plot<-r.e.n.r ########## plots residuals against each other plot(y.for.plot~x.for.plot, ylab="Metacommunity native richness | (Others)", xlab="Metacommunity exotic richness | (Others)",pch=19, col="grey", cex=1.25,cex.axis=1, cex.lab=1.25, ylim=c(-1.5,1.5)) #########extract coefficients for the relationship between partial residuals mod.for.plot<-lm(y.for.plot~x.for.plot, na.action = na.exclude) pred.data<-lm.predict(mod.for.plot, data.frame(y.for.plot=seq.func(y.for.plot), x.for.plot=seq.func(x.for.plot))) ######### plot the fitted line and CI for the relationship between partials plot.CI.func(x.for.plot=seq.func(x.for.plot), pred=pred.data$pred, upper=pred.data$upper, lower=pred.data$lower, env.colour="black", env.trans=40, line.colour="black", line.type=1,line.weight=2.5) mtext("(c)", line=0.5, side=3, adj=0, cex=1.25) dev.off()