--- title: "Jasper Ridge Synchrony" author: "Jonathan Walter, Lauren Hallett, et al." date: "" output: pdf_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` ## Interspecific geographies of synchrony at Jasper Ridge This document organizes for completeness and reproducibility all analyses presented in the manuscript "Disturbance and competition structure spatial and interspecific synchrony in a serpentine plant community," by Jonathan Walter, Lauren Hallett, Lawrence Sheppard, Tom Anderson, Richard Hobbs, Katie Suding, and Dan Reuman. These analyses were performed by JW in R version 3.6.1, with 'wsyn' version 1.0.1, 'ecodist' 2.0.1, and 'tidyverse' 1.2.1, running on macOS Mojave version 10.14.6. ```{r data_import_cleaning, cache=TRUE, echo=FALSE, message=FALSE} rm(list=ls()) library(wsyn) library(ecodist) library(tidyverse) # ----------------------------------------------------------------------- # Data import, formatting, and cleaning setwd("~/Dropbox/JR_synchrony/EDIPackage") plantago<-read.csv("PLER_cov_1m.csv", header=TRUE) #normalize with square root bromus<-read.csv("BRHO_cov_1m.csv", header=TRUE) lasthenia<-read.csv("LACA_cov_1m.csv", header=TRUE) microseris<-read.csv("MIDO_cov_1m.csv", header=TRUE) vulpia<-read.csv("VUMI_cov_1m.csv", header=TRUE) calycadenia<-read.csv("CAMU_cov_1m.csv", header=TRUE) gophers<-read.csv("gopher_cov_1m.csv", header=TRUE) plot.coords<-read.csv("JR_plot_coordinates.csv", header=TRUE) soil.depth<-read.csv("JR_soil_reorder.csv", header=TRUE) PDSI.raw<-read.csv("NCO_PDSI_growingseason.csv") JR_aggcover <- read.csv("JR_cover_1mplot.csv") %>% tbl_df() agg.all<-JR_aggcover %>% group_by(year, uniqueID) %>% summarize(cover = sum(cover)) %>% spread(year, cover) JRprecip.raw<-read.csv("../JR_weather/JR_monthprecip_8315.csv") #precip2 has already been centered and scaled plot.ids<-plantago[,1] years<-1983:2015 pler.cln<-cleandat(as.matrix(plantago[,-1]),times=years,clev=5)$cdat brho.cln<-cleandat(as.matrix(bromus[,-1]),times=years,clev=5)$cdat laca.cln<-cleandat(as.matrix(lasthenia[,-1]),times=years,clev=5)$cdat mido.cln<-cleandat(as.matrix(microseris[,-1]),times=years,clev=5)$cdat vumi.cln<-cleandat(as.matrix(vulpia[,-1]),times=years,clev=5)$cdat camu.cln<-cleandat(as.matrix(calycadenia[,-1]),times=years,clev=5)$cdat goph.cln<-cleandat(as.matrix(gophers[,-1]),times=years,clev=3)$cdat pagg.cln<-cleandat(as.matrix(agg.all[-1]),times=years,clev=5)$cdat colnames(plot.coords)<-c("plotID","X","Y") PDSI<-matrix(PDSI.raw$PDSI[PDSI.raw$samplingYear %in% years], nrow=length(plot.ids), ncol=length(years), byrow=T) dat0 <- JRprecip.raw %>% tbl_df() falldat <- dat0 %>% filter(month == 10 | month == 11) %>% group_by(samplingYear) %>% summarize(precip = sum(allprecip)) winterdat <- dat0 %>% filter(month == 12 | month == 1 | month == 2) %>% group_by(samplingYear) %>% summarize(precip = sum(allprecip)) fallprcp<-matrix(falldat$precip, nrow=length(plot.ids), ncol=length(years), byrow = T) wintprcp<-matrix(winterdat$precip, nrow=length(plot.ids), ncol=length(years), byrow = T) PDSI.cln<-cleandat(PDSI, times=years, clev=5)$cdat fallprcp.cln<-cleandat(fallprcp, times=years, clev=5)$cdat wintprcp.cln<-cleandat(wintprcp, times=years, clev=5)$cdat ``` ## Plot the raw plant time series ```{R plot_raw, echo=FALSE, cache=TRUE} time<-1983:2015 line.plot.matrix<-function(inmat,drop,yrange,title,yaxislab){ plot(time, rep(NA, length(time)), ylim=yrange,main=title,xlab="",ylab=yaxislab,cex.lab=1.2) colors<-rainbow(nrow(inmat)) for(mm in 1:nrow(inmat)){ lines(time, inmat[mm,-drop], col=colors[mm]) } } tiff("./Fig2_raw_timeseries.tif", units="in", height=4.25, width=6.5, res=300) par(mfrow=c(2,3), mar=c(2.5,2.5,2,1), mgp=c(2.5,0.8,0),oma=c(1.5,1.5,0,0)) line.plot.matrix(plantago,drop=1,yrange=c(0,100),title="Plantago",yaxislab="Cover %") text(1984,97,"a)") line.plot.matrix(bromus,drop=1,yrange=c(0,100),title="Bromus",yaxislab="Cover %") text(1984,97,"b)") line.plot.matrix(lasthenia,drop=1,yrange=c(0,100),title="Lasthenia",yaxislab="Cover %") text(1984,97,"c)") line.plot.matrix(microseris,drop=1,yrange=c(0,100),title="Microseris",yaxislab="Cover %") text(1984,97,"d)") line.plot.matrix(vulpia, drop=1, yrange=c(0,100),title="Vulpia",yaxislab="Cover %") text(1984,97,"e)") line.plot.matrix(calycadenia, drop=1, yrange=c(0,100),title="Calycadenia",yaxislab="Cover %") text(1984,97,"f)") #line.plot.matrix(gophers, drop=1, yrange=c(0,5),title="Gopher Disturbance",yaxislab="Subquadrats disturbed") mtext("Year",1,outer=T) mtext("Cover %",2,outer=T) dev.off() ``` ## Geography of Synchrony Analyses Create coherence matrices and clustering; plot distance decay of synchrony and networks; matrix regression analyses ```{r geog_synch, echo=FALSE, cache=TRUE} # ----------------------------------------------------------------------- # 1. Geography of synchrony # ----------------------------------------------------------------------- clust.pler.st<-clust(pler.cln,years,plot.coords,method="ReXWT",tsrange=c(2,6),weighted=T) #short timescales clust.brho.st<-clust(brho.cln,years,plot.coords,method="ReXWT",tsrange=c(2,6),weighted=T) clust.laca.st<-clust(laca.cln,years,plot.coords,method="ReXWT",tsrange=c(2,6),weighted=T) clust.mido.st<-clust(mido.cln,years,plot.coords,method="ReXWT",tsrange=c(2,6),weighted=T) clust.vumi.st<-clust(vumi.cln,years,plot.coords,method="ReXWT",tsrange=c(2,6),weighted=T) clust.camu.st<-clust(camu.cln,years,plot.coords,method="ReXWT",tsrange=c(2,6),weighted=T) clust.goph.st<-clust(goph.cln,years,plot.coords,method="ReXWT",tsrange=c(2,6),weighted=T) clust.pler.lt<-clust(pler.cln,years,plot.coords,method="ReXWT",tsrange=c(6,Inf),weighted=T) #long timescales clust.brho.lt<-clust(brho.cln,years,plot.coords,method="ReXWT",tsrange=c(6,Inf),weighted=T) clust.laca.lt<-clust(laca.cln,years,plot.coords,method="ReXWT",tsrange=c(6,Inf),weighted=T) clust.mido.lt<-clust(mido.cln,years,plot.coords,method="ReXWT",tsrange=c(6,Inf),weighted=T) clust.vumi.lt<-clust(vumi.cln,years,plot.coords,method="ReXWT",tsrange=c(6,Inf),weighted=T) clust.camu.lt<-clust(camu.cln,years,plot.coords,method="ReXWT",tsrange=c(6,Inf),weighted=T) clust.goph.lt<-clust(goph.cln,years,plot.coords,method="ReXWT",tsrange=c(6,Inf),weighted=T) # ---------------------------------------------------------------------- # 1a. Distance-decay of synchrony geog.dist<-dist(plot.coords[,2:3]) #define function for doing spling fitting spline.fit<-function(dists,zmat,nresamp=1000,quantiles=c(0.025,0.975)){ xemp<-c(dists) triang<-lower.tri(zmat) yemp<-zmat[triang] dfs=sqrt(nrow(zmat)) distmat<-full(dists) out<-list() emp.spline<-smooth.spline(xemp,yemp,df=dfs) out$emp.spline<-emp.spline resamp.splines<-matrix(NA, nrow=nresamp, ncol=length(emp.spline$y)) for(ii in 1:nresamp){ shuffle<-sample(1:nrow(zmat), size=nrow(zmat), replace=TRUE) xres<-distmat[shuffle,shuffle][triang] yres<-zmat[shuffle,shuffle][triang] drop.NaNs<-!is.na(yres) xres<-xres[drop.NaNs]; yres<-yres[drop.NaNs] res.spline<-smooth.spline(xres,yres,df=dfs) resamp.splines[ii,]<-predict(res.spline, x=emp.spline$x)$y } out$resamp.splines<-resamp.splines out$spline.quantiles<-apply(resamp.splines,2,quantile,probs=quantiles) return(out) } #do spline fits spline.pler.st<-spline.fit(geog.dist,clust.pler.st$adj) #short timescales spline.brho.st<-spline.fit(geog.dist,clust.brho.st$adj) spline.laca.st<-spline.fit(geog.dist,clust.laca.st$adj) spline.mido.st<-spline.fit(geog.dist,clust.mido.st$adj) spline.vumi.st<-spline.fit(geog.dist,clust.vumi.st$adj) spline.camu.st<-spline.fit(geog.dist,clust.camu.st$adj) spline.goph.st<-spline.fit(geog.dist,clust.goph.st$adj) spline.pler.lt<-spline.fit(geog.dist,clust.pler.lt$adj) #long timescales spline.brho.lt<-spline.fit(geog.dist,clust.brho.lt$adj) spline.laca.lt<-spline.fit(geog.dist,clust.laca.lt$adj) spline.mido.lt<-spline.fit(geog.dist,clust.mido.lt$adj) spline.vumi.lt<-spline.fit(geog.dist,clust.vumi.lt$adj) spline.camu.lt<-spline.fit(geog.dist,clust.camu.lt$adj) spline.goph.lt<-spline.fit(geog.dist,clust.goph.lt$adj) colors=c("red","blue") #blue for in-module, red for out-of-module mod.in.out<-function(mod.ids){ xx<-matrix(1, nrow=length(mod.ids), ncol=length(mod.ids)) for(i in 1:length(mod.ids)){ for(j in 1:length(mod.ids)){ if(mod.ids[i]==mod.ids[j]){xx[i,j]<-2} } } return(xx) } inout.pler.st<-mod.in.out(clust.pler.st$clusters[[length(clust.pler.st$clusters)]]) inout.brho.st<-mod.in.out(clust.brho.st$clusters[[length(clust.brho.st$clusters)]]) inout.laca.st<-mod.in.out(clust.laca.st$clusters[[length(clust.laca.st$clusters)]]) inout.mido.st<-mod.in.out(clust.mido.st$clusters[[length(clust.mido.st$clusters)]]) inout.vumi.st<-mod.in.out(clust.vumi.st$clusters[[length(clust.vumi.st$clusters)]]) inout.camu.st<-mod.in.out(clust.camu.st$clusters[[length(clust.camu.st$clusters)]]) inout.goph.st<-mod.in.out(clust.goph.st$clusters[[length(clust.goph.st$clusters)]]) inout.pler.lt<-mod.in.out(clust.pler.lt$clusters[[length(clust.pler.lt$clusters)]]) inout.brho.lt<-mod.in.out(clust.brho.lt$clusters[[length(clust.brho.lt$clusters)]]) inout.laca.lt<-mod.in.out(clust.laca.lt$clusters[[length(clust.laca.lt$clusters)]]) inout.mido.lt<-mod.in.out(clust.mido.lt$clusters[[length(clust.mido.lt$clusters)]]) inout.vumi.lt<-mod.in.out(clust.vumi.lt$clusters[[length(clust.vumi.lt$clusters)]]) inout.camu.lt<-mod.in.out(clust.camu.lt$clusters[[length(clust.camu.lt$clusters)]]) inout.goph.lt<-mod.in.out(clust.goph.lt$clusters[[length(clust.goph.lt$clusters)]]) ``` Significance testing modularity ```{r modularity sig tests, cache=TRUE, echo=TRUE} ClustSignif<-function(d,numreps,burnin) { #get covariance matrix, which should be the same as the correlation matrix covd<-cov(d) if (max(abs(diag(covd)-1))>1e-5) { stop("Error at loc 1: d should have a covariance matrix with 1s on the diagonal") } #get autocorrelations, one for each plot numplots<-dim(d)[2] rhos<-NA*numeric(numplots) for (counter in 1:numplots) { h<-d[,counter] x<-unname(h[1:(length(h)-1)]) y<-unname(h[2:length(h)]) modl<-lm(y~x) rhos[counter]<-unname(coef(modl)[2]) } #put autocorrelations in a convenient form for later #rhos<-matrix(rep(rhos,times=numreps),numplots,numreps) #makes sims have the same autocorrelations, in each plot, as the data rhos<-matrix(mean(rhos),numplots,numreps) #makes sims have the average autocorrelation in all plots #do the simulations lendat<-dim(d)[1] surrs<-array(NA,c(burnin+lendat,numplots,numreps)) surrs[1,,]<-0 mncov<-(sum(covd)-sum(diag(covd)))/(numplots^2-numplots) sig<-matrix(mncov,numplots,numplots) diag(sig)<-1 noise<-MASS::mvrnorm((burnin+lendat-1)*numreps,mu=rep(0,numplots),Sigma=sig) noise<-array(noise,c(burnin+lendat-1,numreps,numplots)) noise<-aperm(noise,c(1,3,2)) for (counter in 2:(burnin+lendat)) { surrs[counter,,]<-rhos*surrs[counter-1,,]+sqrt(1-rhos^2)*noise[counter-1,,] } surrs<-surrs[(burnin+1):(dim(surrs)[1]),,] #get adjacency matrices, perform clustering, and get modularity for the real data fakecoords<-data.frame(X=1:numplots,Y=1:numplots) clustres_short<-wsyn::clust(dat=t(d),times=1:lendat,coords=fakecoords,method="ReXWT",tsrange=c(2,6)) datmod_short<-clustres_short$modres[[length(clustres_short$modres)]]$totQ clustres_long<-wsyn::clust(dat=t(d),times=1:lendat,coords=fakecoords,method="ReXWT",tsrange=c(6,Inf)) datmod_long<-clustres_long$modres[[length(clustres_long$modres)]]$totQ #get adjacency matrices and perform clustering for each surrogate dataset surrmod_short<-NA*numeric(numreps) surrmod_long<-NA*numeric(numreps) nclust_short<-rep(NA, numreps) nclust_long<-rep(NA, numreps) for (counter in 1:numreps) { #print(paste0("Working on surrogate ",counter," of ",numreps)) surrdat<-t(surrs[,,counter]) surrdat<-wsyn::cleandat(surrdat,times=1:lendat,clev=2)$cdat h_short<-wsyn::clust(dat=surrdat,times=1:lendat,coords=fakecoords,method="ReXWT",tsrange=c(2,6)) surrmod_short[counter]<-h_short$modres[[length(h_short$modres)]]$totQ nclust_short[counter]<-length(h_short$clusters) h_long<-wsyn::clust(dat=surrdat,times=1:lendat,coords=fakecoords,method="ReXWT",tsrange=c(6,Inf)) surrmod_long[counter]<-h_long$modres[[length(h_long$modres)]]$totQ nclust_long[counter]<-length(h_long$clusters) } #get info for p values fracmore_short<-sum(surrmod_short>datmod_short)/numreps fracmore_long<-sum(surrmod_long>datmod_long)/numreps return(list(d=d,surrs=surrs, datmod_short=datmod_short,surrmod_short=surrmod_short,fracmore_short=fracmore_short, datmod_long=datmod_long,surrmod_long=surrmod_long,fracmore_long=fracmore_long, nclust_short=nclust_short, nclust_long=nclust_long)) } d_plantago<-t(pler.cln) d_bromus<-t(brho.cln) d_lasthenia<-t(laca.cln) d_calycadenia<-t(camu.cln) d_microseris<-t(mido.cln) d_vulpia<-t(vumi.cln) #now run everything alldat<-list(platago=d_plantago,bromus=d_bromus,calycadenia=d_calycadenia, lasthenia=d_lasthenia,microseris=d_microseris,vulpia=d_vulpia) allres<-list() for (counter in 1:length(alldat)) { allres[[counter]]<-ClustSignif(alldat[[counter]],numreps=1000,burnin=500) } names(allres)<-names(alldat) #plantago ind<-1 allres[[ind]]$datmod_short allres[[ind]]$fracmore_short #p-value results for short timescale significance of clustering allres[[ind]]$datmod_long allres[[ind]]$fracmore_long #p-value results for long timescale significance of clustering summary(allres[[ind]]$nclust_short) summary(allres[[ind]]$nclust_long) #bromus ind<-2 allres[[ind]]$datmod_short allres[[ind]]$fracmore_short allres[[ind]]$datmod_long allres[[ind]]$fracmore_long summary(allres[[ind]]$nclust_short) summary(allres[[ind]]$nclust_long) #calycadenia ind<-3 allres[[ind]]$datmod_short allres[[ind]]$fracmore_short allres[[ind]]$datmod_long allres[[ind]]$fracmore_long summary(allres[[ind]]$nclust_short) summary(allres[[ind]]$nclust_long) #lasthenia ind<-4 allres[[ind]]$datmod_short allres[[ind]]$fracmore_short allres[[ind]]$datmod_long allres[[ind]]$fracmore_long summary(allres[[ind]]$nclust_short) summary(allres[[ind]]$nclust_long) #microseris ind<-5 allres[[ind]]$datmod_short allres[[ind]]$fracmore_short allres[[ind]]$datmod_long allres[[ind]]$fracmore_long summary(allres[[ind]]$nclust_short) summary(allres[[ind]]$nclust_long) #vulpia ind<-6 allres[[ind]]$datmod_short allres[[ind]]$fracmore_short allres[[ind]]$datmod_long allres[[ind]]$fracmore_long summary(allres[[ind]]$nclust_short) summary(allres[[ind]]$nclust_long) ``` ```{r geog synch figs, cache=FALSE, echo=FALSE} # Distance decay, short timescales (plants) tiff("./Fig4_synchrony_distdecay_reXWT_short.tif", units="in", res=300, width = 6.5, height=4) par(mfrow=c(2,3), mar=c(2,2,1.5,1), mgp=c(2,0.8,0), tcl=-0.2, oma=c(2,2,0,0), xpd=F) plot(spline.pler.st$emp.spline$x, spline.pler.st$emp.spline$y, type="l", main="Plantago", xlab="", ylab="", ylim=c(-0.5,1), cex.main=1, lwd=2, xlim=c(1.5,31)) #abline(h=mean(lower(clust.pler.st$adj)),col="grey") abline(h=0, col="grey") lines(spline.pler.st$emp.spline$x, spline.pler.st$emp.spline$y,lwd=2) lines(spline.pler.st$emp.spline$x, spline.pler.st$spline.quantiles[1,], lty=2) lines(spline.pler.st$emp.spline$x, spline.pler.st$spline.quantiles[2,], lty=2) points(geog.dist, lower(clust.pler.st$adj), pch=20, col=colors[lower(inout.pler.st)],cex=0.4) points(31,mean(lower(clust.pler.st$adj*ifelse(inout.pler.st==2,1,NA)),na.rm=T),col="blue",pch=18,cex=2) points(31,mean(lower(clust.pler.st$adj*ifelse(inout.pler.st==1,1,NA)),na.rm=T),col="red",pch=18,cex=2) text(2,0.99,"a)") plot(spline.brho.st$emp.spline$x, spline.brho.st$emp.spline$y, type="l", main="Bromus", xlab="", ylab="", ylim=c(-0.5,1), cex.main=1, lwd=2, xlim=c(1.5,31)) #abline(h=mean(lower(clust.brho.st$adj)),col="grey") abline(h=0, col="grey") lines(spline.brho.st$emp.spline$x, spline.brho.st$emp.spline$y,lwd=2) lines(spline.brho.st$emp.spline$x, spline.brho.st$spline.quantiles[1,], lty=2) lines(spline.brho.st$emp.spline$x, spline.brho.st$spline.quantiles[2,], lty=2) points(geog.dist, lower(clust.brho.st$adj), pch=20, col=colors[lower(inout.brho.st)],cex=0.4) points(31,mean(lower(clust.brho.st$adj*ifelse(inout.brho.st==2,1,NA)),na.rm=T),col="blue",pch=18,cex=2) points(31,mean(lower(clust.brho.st$adj*ifelse(inout.brho.st==1,1,NA)),na.rm=T),col="red",pch=18,cex=2) text(2,0.99,"b)") plot(spline.laca.st$emp.spline$x, spline.laca.st$emp.spline$y, type="l", main="Lasthenia", xlab="", ylab="", ylim=c(-0.5,1), cex.main=1, lwd=2, xlim=c(1.5,31)) #abline(h=mean(lower(clust.laca.st$adj)),col="grey") abline(h=0, col="grey") lines(spline.laca.st$emp.spline$x, spline.laca.st$emp.spline$y,lwd=2) lines(spline.laca.st$emp.spline$x, spline.laca.st$spline.quantiles[1,], lty=2) lines(spline.laca.st$emp.spline$x, spline.laca.st$spline.quantiles[2,], lty=2) points(geog.dist, lower(clust.laca.st$adj), pch=20, col=colors[inout.laca.st],cex=0.4) points(31,mean(lower(clust.laca.st$adj*ifelse(inout.laca.st==2,1,NA)),na.rm=T),col="blue",pch=18,cex=2) points(31,mean(lower(clust.laca.st$adj*ifelse(inout.laca.st==1,1,NA)),na.rm=T),col="red",pch=18,cex=2) text(2,0.99,"c)") plot(spline.mido.st$emp.spline$x, spline.mido.st$emp.spline$y, type="l", main="Microseris", xlab="", ylab="", ylim=c(-0.5,1), cex.main=1, lwd=2, xlim=c(1.5,31)) abline(h=0, col="grey") #abline(h=mean(lower(clust.mido.st$adj)),col="grey") lines(spline.mido.st$emp.spline$x, spline.mido.st$emp.spline$y,lwd=2) lines(spline.mido.st$emp.spline$x, spline.mido.st$spline.quantiles[1,], lty=2) lines(spline.mido.st$emp.spline$x, spline.mido.st$spline.quantiles[2,], lty=2) points(geog.dist, lower(clust.mido.st$adj), pch=20, col=colors[lower(inout.mido.st)],cex=0.4) points(31,mean(lower(clust.mido.st$adj*ifelse(inout.mido.st==2,1,NA)),na.rm=T),col="blue",pch=18,cex=2) points(31,mean(lower(clust.mido.st$adj*ifelse(inout.mido.st==1,1,NA)),na.rm=T),col="red",pch=18,cex=2) text(2,0.99,"d)") plot(spline.vumi.st$emp.spline$x, spline.vumi.st$emp.spline$y, type="l", main="Vulpia", xlab="", ylab="", ylim=c(-0.5,1), cex.main=1, lwd=2, xlim=c(1.5,31)) #abline(h=mean(lower(clust.vumi.st$adj)),col="grey") abline(h=0, col="grey") lines(spline.vumi.st$emp.spline$x, spline.vumi.st$emp.spline$y,lwd=2) lines(spline.vumi.st$emp.spline$x, spline.vumi.st$spline.quantiles[1,], lty=2) lines(spline.vumi.st$emp.spline$x, spline.vumi.st$spline.quantiles[2,], lty=2) points(geog.dist, lower(clust.vumi.st$adj), pch=20, col=colors[lower(inout.vumi.st)],cex=0.4) points(31,mean(lower(clust.vumi.st$adj*ifelse(inout.vumi.st==2,1,NA)),na.rm=T),col="blue",pch=18,cex=2) points(31,mean(lower(clust.vumi.st$adj*ifelse(inout.vumi.st==1,1,NA)),na.rm=T),col="red",pch=18,cex=2) text(2,0.99,"e)") plot(spline.camu.st$emp.spline$x, spline.camu.st$emp.spline$y, type="l", main="Calycadenia", xlab="", ylab="", ylim=c(-0.5,1), cex.main=1, lwd=2, xlim=c(1.5,31)) abline(h=0, col="grey") #abline(h=mean(lower(clust.camu.st$adj)),col="grey") lines(spline.camu.st$emp.spline$x, spline.camu.st$emp.spline$y,lwd=2) lines(spline.camu.st$emp.spline$x, spline.camu.st$spline.quantiles[1,], lty=2) lines(spline.camu.st$emp.spline$x, spline.camu.st$spline.quantiles[2,], lty=2) points(geog.dist, lower(clust.camu.st$adj), pch=20, col=colors[lower(inout.camu.st)],cex=0.4) points(31,mean(lower(clust.camu.st$adj*ifelse(inout.camu.st==2,1,NA)),na.rm=T),col="blue",pch=18,cex=2) points(31,mean(lower(clust.camu.st$adj*ifelse(inout.camu.st==1,1,NA)),na.rm=T),col="red",pch=18,cex=2) text(2,0.99,"f)") mtext("Synchrony",2,outer=T,line=0.5, cex=1) mtext("Distance between plots (m)",1, outer=T, line=0.5, cex=1) dev.off() #Distance decay, long timescales, plants tiff("./FigS2_synchrony_distdecay_reXWT_long.tif", units="in", res=300, width = 6.5, height=4) par(mfrow=c(2,3), mar=c(2,2,1.5,1), mgp=c(2,0.8,0), tcl=-0.2, oma=c(2,2,0,0), xpd=F) plot(spline.pler.lt$emp.spline$x, spline.pler.lt$emp.spline$y, type="l", main="Plantago", xlab="", ylab="", ylim=c(-0.5,1), cex.main=1, lwd=2, xlim=c(1.5,31)) #abline(h=mean(lower(clust.pler.lt$adj)),col="grey") abline(h=0, col="grey") lines(spline.pler.lt$emp.spline$x, spline.pler.lt$emp.spline$y,lwd=2) lines(spline.pler.lt$emp.spline$x, spline.pler.lt$spline.quantiles[1,], lty=2) lines(spline.pler.lt$emp.spline$x, spline.pler.lt$spline.quantiles[2,], lty=2) points(geog.dist, lower(clust.pler.lt$adj), pch=20, col=colors[lower(inout.pler.lt)],cex=0.4) points(31,mean(lower(clust.pler.lt$adj*ifelse(inout.pler.lt==2,1,NA)),na.rm=T),col="blue",pch=18,cex=2) points(31,mean(lower(clust.pler.lt$adj*ifelse(inout.pler.lt==1,1,NA)),na.rm=T),col="red",pch=18,cex=2) text(2,0.99,"a)") plot(spline.brho.lt$emp.spline$x, spline.brho.lt$emp.spline$y, type="l", main="Bromus", xlab="", ylab="", ylim=c(-0.5,1), cex.main=1, lwd=2, xlim=c(1.5,31)) #abline(h=mean(lower(clust.brho.lt$adj)),col="grey") abline(h=0, col="grey") lines(spline.brho.lt$emp.spline$x, spline.brho.lt$emp.spline$y,lwd=2) lines(spline.brho.lt$emp.spline$x, spline.brho.lt$spline.quantiles[1,], lty=2) lines(spline.brho.lt$emp.spline$x, spline.brho.lt$spline.quantiles[2,], lty=2) points(geog.dist, lower(clust.brho.lt$adj), pch=20, col=colors[lower(inout.brho.lt)],cex=0.4) points(31,mean(lower(clust.brho.lt$adj*ifelse(inout.brho.lt==2,1,NA)),na.rm=T),col="blue",pch=18,cex=2) points(31,mean(lower(clust.brho.lt$adj*ifelse(inout.brho.lt==1,1,NA)),na.rm=T),col="red",pch=18,cex=2) text(2,0.99,"b)") plot(spline.laca.lt$emp.spline$x, spline.laca.lt$emp.spline$y, type="l", main="Lasthenia", xlab="", ylab="", ylim=c(-0.5,1), cex.main=1, lwd=2, xlim=c(1.5,31)) #abline(h=mean(lower(clust.laca.lt$adj)),col="grey") abline(h=0, col="grey") lines(spline.laca.lt$emp.spline$x, spline.laca.lt$emp.spline$y,lwd=2) lines(spline.laca.lt$emp.spline$x, spline.laca.lt$spline.quantiles[1,], lty=2) lines(spline.laca.lt$emp.spline$x, spline.laca.lt$spline.quantiles[2,], lty=2) points(geog.dist, lower(clust.laca.lt$adj), pch=20, col=colors[inout.laca.lt],cex=0.4) points(31,mean(lower(clust.laca.lt$adj*ifelse(inout.laca.lt==2,1,NA)),na.rm=T),col="blue",pch=18,cex=2) points(31,mean(lower(clust.laca.lt$adj*ifelse(inout.laca.lt==1,1,NA)),na.rm=T),col="red",pch=18,cex=2) text(2,0.99,"c)") plot(spline.mido.lt$emp.spline$x, spline.mido.lt$emp.spline$y, type="l", main="Microseris", xlab="", ylab="", ylim=c(-0.5,1), cex.main=1, lwd=2, xlim=c(1.5,31)) abline(h=0, col="grey") #abline(h=mean(lower(clust.mido.lt$adj)),col="grey") lines(spline.mido.lt$emp.spline$x, spline.mido.lt$emp.spline$y,lwd=2) lines(spline.mido.lt$emp.spline$x, spline.mido.lt$spline.quantiles[1,], lty=2) lines(spline.mido.lt$emp.spline$x, spline.mido.lt$spline.quantiles[2,], lty=2) points(geog.dist, lower(clust.mido.lt$adj), pch=20, col=colors[lower(inout.mido.lt)],cex=0.4) points(31,mean(lower(clust.mido.lt$adj*ifelse(inout.mido.lt==2,1,NA)),na.rm=T),col="blue",pch=18,cex=2) points(31,mean(lower(clust.mido.lt$adj*ifelse(inout.mido.lt==1,1,NA)),na.rm=T),col="red",pch=18,cex=2) text(2,0.99,"d)") plot(spline.vumi.lt$emp.spline$x, spline.vumi.lt$emp.spline$y, type="l", main="Vulpia", xlab="", ylab="", ylim=c(-1,1), cex.main=1, lwd=2, xlim=c(1.5,31)) #abline(h=mean(lower(clust.vumi.lt$adj)),col="grey") abline(h=0, col="grey") lines(spline.vumi.lt$emp.spline$x, spline.vumi.lt$emp.spline$y,lwd=2) lines(spline.vumi.lt$emp.spline$x, spline.vumi.lt$spline.quantiles[1,], lty=2) lines(spline.vumi.lt$emp.spline$x, spline.vumi.lt$spline.quantiles[2,], lty=2) points(geog.dist, lower(clust.vumi.lt$adj), pch=20, col=colors[lower(inout.vumi.lt)],cex=0.4) points(31,mean(lower(clust.vumi.lt$adj*ifelse(inout.vumi.lt==2,1,NA)),na.rm=T),col="blue",pch=18,cex=2) points(31,mean(lower(clust.vumi.lt$adj*ifelse(inout.vumi.lt==1,1,NA)),na.rm=T),col="red",pch=18,cex=2) text(2,0.99,"e)") plot(spline.camu.lt$emp.spline$x, spline.camu.lt$emp.spline$y, type="l", main="Calycadenia", xlab="", ylab="", ylim=c(-0.5,1), cex.main=1, lwd=2, xlim=c(1.5,31)) abline(h=0, col="grey") #abline(h=mean(lower(clust.camu.lt$adj)),col="grey") lines(spline.camu.lt$emp.spline$x, spline.camu.lt$emp.spline$y,lwd=2) lines(spline.camu.lt$emp.spline$x, spline.camu.lt$spline.quantiles[1,], lty=2) lines(spline.camu.lt$emp.spline$x, spline.camu.lt$spline.quantiles[2,], lty=2) points(geog.dist, lower(clust.camu.lt$adj), pch=20, col=colors[lower(inout.camu.lt)],cex=0.4) points(31,mean(lower(clust.camu.lt$adj*ifelse(inout.camu.lt==2,1,NA)),na.rm=T),col="blue",pch=18,cex=2) points(31,mean(lower(clust.camu.lt$adj*ifelse(inout.camu.lt==1,1,NA)),na.rm=T),col="red",pch=18,cex=2) text(2,0.99,"f)") mtext("Synchrony",2,outer=T, cex=1) mtext("Distance between plots (m)",1, outer=T, line=0.5, cex=1) dev.off() # ---------------------------------------------------------------------- # 1b. Networks and modularity source("./plotClusterMap_JRG.R") treatment=rep(c(1,1,1,2,2,2,3,3,3),times=4) pchar=15:17 #Networks, short timescales, plants tiff("./Fig5_networks_ReXWT_short.tif", units="in", res=300, width = 6.5, height=6) layout(matrix(c(1:6,rep(7,3)),nrow=3,ncol=3, byrow = T),heights=c(0.5,0.5,0.05)) par(mar=c(1,1,1.5,1)) plotClusterMap(clust.pler.st, nodewgt=NULL, edgewgt=0.9, title="Plantago", nodesize=c(2,2), edgesize = c(0.5,0.5), addlegend = F, nodechar=pchar[treatment]) text(1.5,28.5,"a)",cex=1.5) plotClusterMap(clust.brho.st, nodewgt=NULL, edgewgt=0.9, title="Bromus", nodesize=c(2,2), edgesize = c(0.5,0.5), addlegend = F, nodechar=pchar[treatment]) text(1.5,28.5,"b)",cex=1.5) plotClusterMap(clust.laca.st, nodewgt=NULL, edgewgt=0.9, title="Lasthenia", nodesize=c(2,2), edgesize = c(0.5,0.5), addlegend = F, nodechar=pchar[treatment]) text(1.5,28.5,"c)",cex=1.5) plotClusterMap(clust.mido.st, nodewgt=NULL, edgewgt=0.9, title="Microseris", nodesize=c(2,2), edgesize = c(0.5,0.5), addlegend = F, nodechar=pchar[treatment]) text(1.5,28.5,"d)",cex=1.5) plotClusterMap(clust.vumi.st, nodewgt=NULL, edgewgt=0.9, title="Vulpia", nodesize=c(2,2), edgesize = c(0.5,0.5), addlegend = F, nodechar=pchar[treatment]) text(1.5,28.5,"e)",cex=1.5) plotClusterMap(clust.camu.st, nodewgt=NULL, edgewgt=0.9, title="Calycadenia", nodesize=c(2,2), edgesize = c(0.5,0.5), addlegend = F, nodechar=pchar[treatment]) text(1.5,28.5,"f)",cex=1.5) par(mar=c(0.5,0.1,0.1,0.1)) plot(NA,NA,bty="n",xlim=c(0,1),ylim=c(0,1),xaxt="n",yaxt="n") legend("center", c("Module 1","Module 2","Module 3"), pch=20, col=brewer.pal(3,"Set1"),ncol=3, cex=1.5, bty="n", pt.cex=2, inset=c(0,0.5)) dev.off() #Networks, long timescales, plants tiff("./FigS3_networks_ReXWT_long.tif", units="in", res=300, width = 6.5, height=6) layout(matrix(c(1:6,rep(7,3)),nrow=3,ncol=3, byrow = T),heights=c(0.5,0.5,0.05)) par(mar=c(1,1,1.5,1)) plotClusterMap(clust.pler.lt, nodewgt=NULL, edgewgt=0.9, title="Plantago", nodesize=c(2,2), edgesize = c(0.5,0.5), addlegend = F, nodechar=pchar[treatment]) text(1.5,28.5,"a)",cex=1.5) plotClusterMap(clust.brho.lt, nodewgt=NULL, edgewgt=0.9, title="Bromus", nodesize=c(2,2), edgesize = c(0.5,0.5), addlegend = F, nodechar=pchar[treatment]) text(1.5,28.5,"b)",cex=1.5) plotClusterMap(clust.laca.lt, nodewgt=NULL, edgewgt=0.9, title="Lasthenia", nodesize=c(2,2), edgesize = c(0.5,0.5), addlegend = F, nodechar=pchar[treatment]) text(1.5,28.5,"c)",cex=1.5) plotClusterMap(clust.mido.lt, nodewgt=NULL, edgewgt=0.9, title="Microseris", nodesize=c(2,2), edgesize = c(0.5,0.5), addlegend = F, nodechar=pchar[treatment]) text(1.5,28.5,"d)",cex=1.5) plotClusterMap(clust.vumi.lt, nodewgt=NULL, edgewgt=0.9, title="Vulpia", nodesize=c(2,2), edgesize = c(0.5,0.5), addlegend = F, nodechar=pchar[treatment]) text(1.5,28.5,"e)",cex=1.5) plotClusterMap(clust.camu.lt, nodewgt=NULL, edgewgt=0.9, title="Calycadenia", nodesize=c(2,2), edgesize = c(0.5,0.5), addlegend = F, nodechar=pchar[treatment]) text(1.5,28.5,"f)",cex=1.5) par(mar=c(0.5,0.1,0.1,0.1)) plot(NA,NA,bty="n",xlim=c(0,1),ylim=c(0,1),xaxt="n",yaxt="n") legend("center", c("Module 1","Module 2"), pch=20, col=brewer.pal(3,"Set1")[1:2],ncol=2, cex=1.5, bty="n", pt.cex=2, inset=c(0,0.5)) dev.off() #Gophers, short timescales coordx<-c(plot.coords$X) coordy<-c(plot.coords$Y) treatment=rep(c(1,1,1,2,2,2,3,3,3),times=4) pchar=15:17 greypal<-colorRampPalette(c("grey90","grey10")) sdlevs<-levels(cut(soil.depth$soilDepth, breaks=5)) sdtxt<-paste(substr(sdlevs,2,5),"-",substr(sdlevs,7,10)) sdtxt<-gsub("]",".0",sdtxt) tiff("./Fig3_gophers_short.tif", units="in", res=300, width=6.5, height=2) par(mfrow=c(1,3)) par(mar=c(3,3,1.5,1), mgp=c(1.5,0.2,0), tcl=-0.2) plot(spline.goph.st$emp.spline$x, spline.goph.st$emp.spline$y, type="l", main="", xlab="Distance (m)", ylab="Synchrony", ylim=c(-.5,1), cex.main=1, lwd=2, xlim=c(1.5,31)) #abline(h=mean(lower(clust.goph.st$adj)),col="grey") abline(h=0,col="grey") lines(spline.goph.st$emp.spline$x, spline.goph.st$emp.spline$y,lwd=2) lines(spline.goph.st$emp.spline$x, spline.goph.st$spline.quantiles[1,], lty=2) lines(spline.goph.st$emp.spline$x, spline.goph.st$spline.quantiles[2,], lty=2) points(geog.dist, lower(clust.goph.st$adj), pch=20, col=colors[lower(inout.goph.st)],cex=0.4) points(31,mean(lower(clust.goph.st$adj*ifelse(inout.goph.st==2,1,NA)),na.rm=T),col="blue",pch=18,cex=2) points(31,mean(lower(clust.goph.st$adj*ifelse(inout.goph.st==1,1,NA)),na.rm=T),col="red",pch=18,cex=2) text(2,0.99,"a)") par(mar=c(1,1,1.5,6.1), par(xpd=T)) plot(coordx, coordy, pch=pchar[treatment], col=greypal(5)[as.numeric(cut(soil.depth$soilDepth, breaks=5))], bty="o", asp=1, main="", xaxt="n",yaxt="n", ylab="", xlab="", cex=1.4) lines(x=c(min(coordx),min(coordx)+5),y=c(1.5,1.5)) text(min(coordx)+2.5,2.4,"5m") legx<-par('usr')[2] + 0.005*abs(diff(par('usr')[1:2])) legy<-par('usr')[4] legend(legx,legy,fill=greypal(5)[1:5],legend=sdtxt, bty="n", title="Soil depth (cm)", cex=0.8) text(0.5,28.5,"b)") par(mar=c(1,1,1.5,6.1),xpd=T) plotClusterMap(clust.goph.st, nodewgt=NULL, edgewgt=0.9, title="", nodesize=c(1.5,1.5), edgesize = c(0.5,0.5), nodechar=pchar[treatment]) text(0.2,28.5,"c)") dev.off() #Gophers, long timescales tiff("./FigS1_gophers_long.tif", units="in", res=300, width=6.5, height=2) par(mfrow=c(1,3)) par(mar=c(3,3,1.5,1), mgp=c(1.5,0.2,0), tcl=-0.2) plot(spline.goph.lt$emp.spline$x, spline.goph.lt$emp.spline$y, type="l", main="", xlab="Distance (m)", ylab="Synchrony", ylim=c(-1,1), cex.main=1, lwd=2, xlim=c(1.5,31)) #abline(h=mean(lower(clust.goph.lt$adj)),col="grey") abline(h=0,col="grey") lines(spline.goph.lt$emp.spline$x, spline.goph.lt$emp.spline$y,lwd=2) lines(spline.goph.lt$emp.spline$x, spline.goph.lt$spline.quantiles[1,], lty=2) lines(spline.goph.lt$emp.spline$x, spline.goph.lt$spline.quantiles[2,], lty=2) points(geog.dist, lower(clust.goph.lt$adj), pch=20, col=colors[lower(inout.goph.lt)],cex=0.4) points(31,mean(lower(clust.goph.lt$adj*ifelse(inout.goph.lt==2,1,NA)),na.rm=T),col="blue",pch=18,cex=2) points(31,mean(lower(clust.goph.lt$adj*ifelse(inout.goph.lt==1,1,NA)),na.rm=T),col="red",pch=18,cex=2) text(2,0.99,"a)") par(mar=c(1,1,1.5,6.1), par(xpd=T)) plot(coordx, coordy, pch=pchar[treatment], col=greypal(5)[as.numeric(cut(soil.depth$soilDepth, breaks=5))], bty="o", asp=1, main="", xaxt="n",yaxt="n", ylab="", xlab="", cex=1.4) lines(x=c(min(coordx),min(coordx)+5),y=c(1.5,1.5)) text(min(coordx)+2.5,2.4,"5m") legx<-par('usr')[2] + 0.005*abs(diff(par('usr')[1:2])) legy<-par('usr')[4] legend(legx,legy,fill=greypal(5)[1:5],legend=sdtxt, bty="n", title="Soil depth (cm)", cex=0.8) text(0.5,28.5,"b)") par(mar=c(1,1,1.5,6.1),xpd=T) plotClusterMap(clust.goph.lt, nodewgt=NULL, edgewgt=0.9, title="", nodesize=c(1.5,1.5), edgesize = c(0.5,0.5), nodechar=pchar[treatment]) text(0.2,28.5,"c)") dev.off() ``` MRM Analyses ```{r MRM, echo=TRUE, cache=TRUE} # ---------------------------------------------------------------------- # 1c. Matrix regression soild.sim<-1-(dist(soil.depth$soilDepth)/max(dist(soil.depth$soilDepth))) geog.prox<-1-(geog.dist/max(geog.dist)) ## Short Timescales ################################## #Plantago MRM(as.dist(clust.pler.st$adj)~as.dist(geog.prox),nperm=10000) #** MRM(as.dist(clust.pler.st$adj)~as.dist(clust.goph.st$adj),nperm=10000) #*** MRM(as.dist(clust.pler.st$adj)~as.dist(soild.sim),nperm=10000) #* MRM(as.dist(clust.pler.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(soild.sim)+as.dist(clust.brho.st$adj),nperm=10000) MRM(as.dist(clust.pler.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(soild.sim)+as.dist(clust.laca.st$adj),nperm=10000) MRM(as.dist(clust.pler.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(soild.sim)+as.dist(clust.mido.st$adj),nperm=10000) MRM(as.dist(clust.pler.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(soild.sim)+as.dist(clust.vumi.st$adj),nperm=10000) MRM(as.dist(clust.pler.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(soild.sim)+as.dist(clust.camu.st$adj),nperm=10000) #Bromus MRM(as.dist(clust.brho.st$adj)~as.dist(geog.prox),nperm=10000) #* MRM(as.dist(clust.brho.st$adj)~as.dist(clust.goph.st$adj),nperm=10000) MRM(as.dist(clust.brho.st$adj)~as.dist(soild.sim),nperm=10000) #* MRM(as.dist(clust.brho.st$adj)~as.dist(geog.prox)+as.dist(soild.sim)+as.dist(clust.pler.st$adj),nperm=10000) MRM(as.dist(clust.brho.st$adj)~as.dist(geog.prox)+as.dist(soild.sim)+as.dist(clust.laca.st$adj),nperm=10000) MRM(as.dist(clust.brho.st$adj)~as.dist(geog.prox)+as.dist(soild.sim)+as.dist(clust.mido.st$adj),nperm=10000) MRM(as.dist(clust.brho.st$adj)~as.dist(geog.prox)+as.dist(soild.sim)+as.dist(clust.vumi.st$adj),nperm=10000) MRM(as.dist(clust.brho.st$adj)~as.dist(geog.prox)+as.dist(soild.sim)+as.dist(clust.camu.st$adj),nperm=10000) #Lasthenia MRM(as.dist(clust.laca.st$adj)~as.dist(geog.prox),nperm=10000) #*** MRM(as.dist(clust.laca.st$adj)~as.dist(clust.goph.st$adj),nperm=10000) #*** MRM(as.dist(clust.laca.st$adj)~as.dist(soild.sim),nperm=10000) MRM(as.dist(clust.laca.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(clust.pler.st$adj),nperm=10000) MRM(as.dist(clust.laca.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(clust.brho.st$adj),nperm=10000) MRM(as.dist(clust.laca.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(clust.mido.st$adj),nperm=10000) MRM(as.dist(clust.laca.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(clust.vumi.st$adj),nperm=10000) MRM(as.dist(clust.laca.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(clust.camu.st$adj),nperm=10000) #Microseris MRM(as.dist(clust.mido.st$adj)~as.dist(geog.prox),nperm=10000) #*** MRM(as.dist(clust.mido.st$adj)~as.dist(clust.goph.st$adj),nperm=10000) #*** MRM(as.dist(clust.mido.st$adj)~as.dist(soild.sim),nperm=10000) MRM(as.dist(clust.mido.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(clust.pler.st$adj),nperm=10000) MRM(as.dist(clust.mido.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(clust.brho.st$adj),nperm=10000) MRM(as.dist(clust.mido.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(clust.laca.st$adj),nperm=10000) MRM(as.dist(clust.mido.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(clust.vumi.st$adj),nperm=10000) MRM(as.dist(clust.mido.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(clust.camu.st$adj),nperm=10000) #Vulpia MRM(as.dist(clust.vumi.st$adj)~as.dist(geog.prox),nperm=10000) #*** MRM(as.dist(clust.vumi.st$adj)~as.dist(clust.goph.st$adj),nperm=10000) #*** MRM(as.dist(clust.vumi.st$adj)~as.dist(soild.sim),nperm=10000) MRM(as.dist(clust.vumi.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(clust.pler.st$adj),nperm=10000) MRM(as.dist(clust.vumi.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(clust.brho.st$adj),nperm=10000) MRM(as.dist(clust.vumi.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(clust.laca.st$adj),nperm=10000) MRM(as.dist(clust.vumi.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(clust.mido.st$adj),nperm=10000) MRM(as.dist(clust.vumi.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(clust.camu.st$adj),nperm=10000) #Calycadenia MRM(as.dist(clust.camu.st$adj)~as.dist(geog.prox),nperm=10000) #*** MRM(as.dist(clust.camu.st$adj)~as.dist(clust.goph.st$adj),nperm=10000) #* MRM(as.dist(clust.camu.st$adj)~as.dist(soild.sim),nperm=10000) MRM(as.dist(clust.camu.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(soild.sim)+as.dist(clust.pler.st$adj),nperm=10000) MRM(as.dist(clust.camu.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(soild.sim)+as.dist(clust.brho.st$adj),nperm=10000) MRM(as.dist(clust.camu.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(soild.sim)+as.dist(clust.laca.st$adj),nperm=10000) MRM(as.dist(clust.camu.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(soild.sim)+as.dist(clust.mido.st$adj),nperm=10000) MRM(as.dist(clust.camu.st$adj)~as.dist(geog.prox)+as.dist(clust.goph.st$adj)+as.dist(soild.sim)+as.dist(clust.vumi.st$adj),nperm=10000) ## Long Timescales ################### #Plantago MRM(as.dist(clust.pler.lt$adj)~as.dist(geog.prox),nperm=10000) #* MRM(as.dist(clust.pler.lt$adj)~as.dist(clust.goph.lt$adj),nperm=10000) #** MRM(as.dist(clust.pler.lt$adj)~as.dist(soild.sim),nperm=10000) MRM(as.dist(clust.pler.lt$adj)~as.dist(geog.prox)+as.dist(clust.goph.lt$adj)+as.dist(clust.brho.lt$adj),nperm=10000) MRM(as.dist(clust.pler.lt$adj)~as.dist(geog.prox)+as.dist(clust.goph.lt$adj)+as.dist(clust.laca.lt$adj),nperm=10000) MRM(as.dist(clust.pler.lt$adj)~as.dist(geog.prox)+as.dist(clust.goph.lt$adj)+as.dist(clust.mido.lt$adj),nperm=10000) MRM(as.dist(clust.pler.lt$adj)~as.dist(geog.prox)+as.dist(clust.goph.lt$adj)+as.dist(clust.vumi.lt$adj),nperm=10000) MRM(as.dist(clust.pler.lt$adj)~as.dist(geog.prox)+as.dist(clust.goph.lt$adj)+as.dist(clust.camu.lt$adj),nperm=10000) #Bromus MRM(as.dist(clust.brho.lt$adj)~as.dist(geog.prox),nperm=10000) MRM(as.dist(clust.brho.lt$adj)~as.dist(clust.goph.lt$adj),nperm=10000) MRM(as.dist(clust.brho.lt$adj)~as.dist(soild.sim),nperm=10000) MRM(as.dist(clust.brho.lt$adj)~as.dist(clust.pler.lt$adj),nperm=10000) MRM(as.dist(clust.brho.lt$adj)~as.dist(clust.laca.lt$adj),nperm=10000) MRM(as.dist(clust.brho.lt$adj)~as.dist(clust.mido.lt$adj),nperm=10000) MRM(as.dist(clust.brho.lt$adj)~as.dist(clust.vumi.lt$adj),nperm=10000) MRM(as.dist(clust.brho.lt$adj)~as.dist(clust.camu.lt$adj),nperm=10000) #Lasthenia MRM(as.dist(clust.laca.lt$adj)~as.dist(geog.prox),nperm=10000) MRM(as.dist(clust.laca.lt$adj)~as.dist(clust.goph.lt$adj),nperm=10000) #** MRM(as.dist(clust.laca.lt$adj)~as.dist(soild.sim),nperm=10000) MRM(as.dist(clust.laca.lt$adj)~as.dist(clust.goph.lt$adj)+as.dist(clust.pler.lt$adj),nperm=10000) MRM(as.dist(clust.laca.lt$adj)~as.dist(clust.goph.lt$adj)+as.dist(clust.brho.lt$adj),nperm=10000) MRM(as.dist(clust.laca.lt$adj)~as.dist(clust.goph.lt$adj)+as.dist(clust.mido.lt$adj),nperm=10000) MRM(as.dist(clust.laca.lt$adj)~as.dist(clust.goph.lt$adj)+as.dist(clust.vumi.lt$adj),nperm=10000) MRM(as.dist(clust.laca.lt$adj)~as.dist(clust.goph.lt$adj)+as.dist(clust.camu.lt$adj),nperm=10000) #Microseris MRM(as.dist(clust.mido.lt$adj)~as.dist(geog.prox),nperm=10000) #*** MRM(as.dist(clust.mido.lt$adj)~as.dist(clust.goph.lt$adj),nperm=10000) #** MRM(as.dist(clust.mido.lt$adj)~as.dist(soild.sim),nperm=10000) #*** MRM(as.dist(clust.mido.lt$adj)~as.dist(geog.prox)+as.dist(clust.goph.lt$adj)+as.dist(soild.sim)+as.dist(clust.pler.lt$adj),nperm=10000) MRM(as.dist(clust.mido.lt$adj)~as.dist(geog.prox)+as.dist(clust.goph.lt$adj)+as.dist(soild.sim)+as.dist(clust.brho.lt$adj),nperm=10000) MRM(as.dist(clust.mido.lt$adj)~as.dist(geog.prox)+as.dist(clust.goph.lt$adj)+as.dist(soild.sim)+as.dist(clust.laca.lt$adj),nperm=10000) MRM(as.dist(clust.mido.lt$adj)~as.dist(geog.prox)+as.dist(clust.goph.lt$adj)+as.dist(soild.sim)+as.dist(clust.vumi.lt$adj),nperm=10000) MRM(as.dist(clust.mido.lt$adj)~as.dist(geog.prox)+as.dist(clust.goph.lt$adj)+as.dist(soild.sim)+as.dist(clust.camu.lt$adj),nperm=10000) #Vulpia MRM(as.dist(clust.vumi.lt$adj)~as.dist(geog.prox),nperm=10000) #* MRM(as.dist(clust.vumi.lt$adj)~as.dist(clust.goph.lt$adj),nperm=10000) MRM(as.dist(clust.vumi.lt$adj)~as.dist(soild.sim),nperm=10000) MRM(as.dist(clust.vumi.lt$adj)~as.dist(geog.prox)+as.dist(clust.pler.lt$adj),nperm=10000) MRM(as.dist(clust.vumi.lt$adj)~as.dist(geog.prox)+as.dist(clust.brho.lt$adj),nperm=10000) MRM(as.dist(clust.vumi.lt$adj)~as.dist(geog.prox)+as.dist(clust.laca.lt$adj),nperm=10000) MRM(as.dist(clust.vumi.lt$adj)~as.dist(geog.prox)+as.dist(clust.mido.lt$adj),nperm=10000) MRM(as.dist(clust.vumi.lt$adj)~as.dist(geog.prox)+as.dist(clust.camu.lt$adj),nperm=10000) #Calycadenia MRM(as.dist(clust.camu.lt$adj)~as.dist(geog.prox),nperm=10000) MRM(as.dist(clust.camu.lt$adj)~as.dist(clust.goph.lt$adj),nperm=10000) MRM(as.dist(clust.camu.lt$adj)~as.dist(soild.sim),nperm=10000) MRM(as.dist(clust.camu.lt$adj)~as.dist(clust.pler.lt$adj),nperm=10000) MRM(as.dist(clust.camu.lt$adj)~as.dist(clust.brho.lt$adj),nperm=10000) MRM(as.dist(clust.camu.lt$adj)~as.dist(clust.laca.lt$adj),nperm=10000) MRM(as.dist(clust.camu.lt$adj)~as.dist(clust.mido.lt$adj),nperm=10000) MRM(as.dist(clust.camu.lt$adj)~as.dist(clust.vumi.lt$adj),nperm=10000) ``` ## Multivariate Wavelet Coherence Models Test coherence in multivariate models and compute fractions of synchrony explained and cross-terms. First, test bivariate coherences with control variables ```{r coh, echo=TRUE, cache=TRUE} st<-c(2,6) lt<-c(6,Inf) plerXgoph<-coh(pler.cln, goph.cln, years, norm="powall", sigmethod="aaftsurrog2", nrand=10000) plerXgoph<-bandtest(plerXgoph, st); plerXgoph<-bandtest(plerXgoph, lt) print(plerXgoph$bandp) plerXpdsi<-coh(pler.cln, PDSI.cln, years, norm="powall", sigmethod="fast", nrand=10000) plerXpdsi<-bandtest(plerXpdsi, st); plerXpdsi<-bandtest(plerXpdsi, lt) print(plerXpdsi$bandp) plerXfallp<-coh(pler.cln, fallprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) plerXfallp<-bandtest(plerXfallp, st); plerXfallp<-bandtest(plerXfallp, lt) print(plerXfallp$bandp) plerXwintp<-coh(pler.cln, wintprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) plerXwintp<-bandtest(plerXwintp, st); plerXwintp<-bandtest(plerXwintp, lt) print(plerXwintp$bandp) brhoXgoph<-coh(brho.cln, goph.cln, years, norm="powall", sigmethod="aaftsurrog2", nrand=10000) brhoXgoph<-bandtest(brhoXgoph, st); brhoXgoph<-bandtest(brhoXgoph, lt) print(brhoXgoph$bandp) brhoXpdsi<-coh(brho.cln, PDSI.cln, years, norm="powall", sigmethod="fast", nrand=10000) brhoXpdsi<-bandtest(brhoXpdsi, st); brhoXpdsi<-bandtest(brhoXpdsi, lt) print(brhoXpdsi$bandp) brhoXfallp<-coh(brho.cln, fallprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) brhoXfallp<-bandtest(brhoXfallp, st); brhoXfallp<-bandtest(brhoXfallp, lt) print(brhoXfallp$bandp) brhoXwintp<-coh(brho.cln, wintprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) brhoXwintp<-bandtest(brhoXwintp, st); brhoXwintp<-bandtest(brhoXwintp, lt) print(brhoXwintp$bandp) lacaXgoph<-coh(laca.cln, goph.cln, years, norm="powall", sigmethod="aaftsurrog2", nrand=10000) lacaXgoph<-bandtest(lacaXgoph, st); lacaXgoph<-bandtest(lacaXgoph, lt) print(lacaXgoph$bandp) lacaXpdsi<-coh(laca.cln, PDSI.cln, years, norm="powall", sigmethod="fast", nrand=10000) lacaXpdsi<-bandtest(lacaXpdsi, st); lacaXpdsi<-bandtest(lacaXpdsi, lt) print(lacaXpdsi$bandp) lacaXfallp<-coh(laca.cln, fallprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) lacaXfallp<-bandtest(lacaXfallp, st); lacaXfallp<-bandtest(lacaXfallp, lt) print(lacaXfallp$bandp) lacaXwintp<-coh(laca.cln, wintprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) lacaXwintp<-bandtest(lacaXwintp, st); lacaXwintp<-bandtest(lacaXwintp, lt) print(lacaXwintp$bandp) midoXgoph<-coh(mido.cln, goph.cln, years, norm="powall", sigmethod="aaftsurrog2", nrand=10000) midoXgoph<-bandtest(midoXgoph, st); midoXgoph<-bandtest(midoXgoph, lt) print(midoXgoph$bandp) midoXpdsi<-coh(mido.cln, PDSI.cln, years, norm="powall", sigmethod="fast", nrand=10000) midoXpdsi<-bandtest(midoXpdsi, st); midoXpdsi<-bandtest(midoXpdsi, lt) print(midoXpdsi$bandp) midoXfallp<-coh(mido.cln, fallprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) midoXfallp<-bandtest(midoXfallp, st); midoXfallp<-bandtest(midoXfallp, lt) print(midoXfallp$bandp) midoXwintp<-coh(mido.cln, wintprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) midoXwintp<-bandtest(midoXwintp, st); midoXwintp<-bandtest(midoXwintp, lt) print(midoXwintp$bandp) vumiXgoph<-coh(vumi.cln, goph.cln, years, norm="powall", sigmethod="aaftsurrog2", nrand=10000) vumiXgoph<-bandtest(vumiXgoph, st); vumiXgoph<-bandtest(vumiXgoph, lt) print(vumiXgoph$bandp) vumiXpdsi<-coh(vumi.cln, PDSI.cln, years, norm="powall", sigmethod="fast", nrand=10000) vumiXpdsi<-bandtest(vumiXpdsi, st); vumiXpdsi<-bandtest(vumiXpdsi, lt) print(vumiXpdsi$bandp) vumiXfallp<-coh(vumi.cln, fallprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) vumiXfallp<-bandtest(vumiXfallp, st); vumiXfallp<-bandtest(vumiXfallp, lt) print(vumiXfallp$bandp) vumiXwintp<-coh(vumi.cln, wintprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) vumiXwintp<-bandtest(vumiXwintp, st); vumiXwintp<-bandtest(vumiXwintp, lt) print(vumiXwintp$bandp) camuXgoph<-coh(camu.cln, goph.cln, years, norm="powall", sigmethod="aaftsurrog2", nrand=10000) camuXgoph<-bandtest(camuXgoph, st); camuXgoph<-bandtest(camuXgoph, lt) print(camuXgoph$bandp) camuXpdsi<-coh(camu.cln, PDSI.cln, years, norm="powall", sigmethod="fast", nrand=10000) camuXpdsi<-bandtest(camuXpdsi, st); camuXpdsi<-bandtest(camuXpdsi, lt) print(camuXpdsi$bandp) camuXfallp<-coh(camu.cln, fallprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) camuXfallp<-bandtest(camuXfallp, st); camuXfallp<-bandtest(camuXfallp, lt) print(camuXfallp$bandp) camuXwintp<-coh(camu.cln, wintprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) camuXwintp<-bandtest(camuXwintp, st); camuXwintp<-bandtest(camuXwintp, lt) print(camuXwintp$bandp) ``` Now get bivariate coherences (for phases) between plants ```{r coh_plants, echo=TRUE, cache=TRUE} #Plantago plerXbrho<-coh(pler.cln, brho.cln, years, norm="powall", sigmethod="fast", nrand=10000) plerXbrho<-bandtest(plerXbrho, st); plerXbrho<-bandtest(plerXbrho, lt) print(plerXbrho$bandp) plerXlaca<-coh(pler.cln, laca.cln, years, norm="powall", sigmethod="fast", nrand=10000) plerXlaca<-bandtest(plerXlaca, st); plerXlaca<-bandtest(plerXlaca, lt) print(plerXlaca$bandp) plerXmido<-coh(pler.cln, mido.cln, years, norm="powall", sigmethod="fast", nrand=10000) plerXmido<-bandtest(plerXmido, st); plerXmido<-bandtest(plerXmido, lt) print(plerXmido$bandp) plerXvumi<-coh(pler.cln, vumi.cln, years, norm="powall", sigmethod="fast", nrand=10000) plerXvumi<-bandtest(plerXvumi, st); plerXvumi<-bandtest(plerXvumi, lt) print(plerXvumi$bandp) plerXcamu<-coh(pler.cln, camu.cln, years, norm="powall", sigmethod="fast", nrand=10000) plerXcamu<-bandtest(plerXcamu, st); plerXcamu<-bandtest(plerXcamu, lt) print(plerXcamu$bandp) #Bromus brhoXpler<-coh(brho.cln, pler.cln, years, norm="powall", sigmethod="fast", nrand=10000) brhoXpler<-bandtest(brhoXpler, st); brhoXpler<-bandtest(brhoXpler, lt) print(brhoXpler$bandp) brhoXlaca<-coh(brho.cln, laca.cln, years, norm="powall", sigmethod="fast", nrand=10000) brhoXlaca<-bandtest(brhoXlaca, st); brhoXlaca<-bandtest(brhoXlaca, lt) print(brhoXlaca$bandp) brhoXmido<-coh(brho.cln, mido.cln, years, norm="powall", sigmethod="fast", nrand=10000) brhoXmido<-bandtest(brhoXmido, st); brhoXmido<-bandtest(brhoXmido, lt) print(brhoXmido$bandp) brhoXvumi<-coh(brho.cln, vumi.cln, years, norm="powall", sigmethod="fast", nrand=10000) brhoXvumi<-bandtest(brhoXvumi, st); brhoXvumi<-bandtest(brhoXvumi, lt) print(brhoXvumi$bandp) brhoXcamu<-coh(brho.cln, camu.cln, years, norm="powall", sigmethod="fast", nrand=10000) brhoXcamu<-bandtest(brhoXcamu, st); brhoXcamu<-bandtest(brhoXcamu, lt) print(brhoXcamu$bandp) #Lasthenia lacaXpler<-coh(laca.cln, pler.cln, years, norm="powall", sigmethod="fast", nrand=10000) lacaXpler<-bandtest(lacaXpler, st); lacaXpler<-bandtest(lacaXpler, lt) print(lacaXpler$bandp) lacaXbrho<-coh(laca.cln, brho.cln, years, norm="powall", sigmethod="fast", nrand=10000) lacaXbrho<-bandtest(lacaXbrho, st); lacaXbrho<-bandtest(lacaXbrho, lt) print(lacaXbrho$bandp) lacaXmido<-coh(laca.cln, mido.cln, years, norm="powall", sigmethod="fast", nrand=10000) lacaXmido<-bandtest(lacaXmido, st); lacaXmido<-bandtest(lacaXmido, lt) print(lacaXmido$bandp) lacaXvumi<-coh(laca.cln, vumi.cln, years, norm="powall", sigmethod="fast", nrand=10000) lacaXvumi<-bandtest(lacaXvumi, st); lacaXvumi<-bandtest(lacaXvumi, lt) print(lacaXvumi$bandp) lacaXcamu<-coh(laca.cln, camu.cln, years, norm="powall", sigmethod="fast", nrand=10000) lacaXcamu<-bandtest(lacaXcamu, st); lacaXcamu<-bandtest(lacaXcamu, lt) print(lacaXcamu$bandp) #Microseris midoXpler<-coh(mido.cln, pler.cln, years, norm="powall", sigmethod="fast", nrand=10000) midoXpler<-bandtest(midoXpler, st); midoXpler<-bandtest(midoXpler, lt) print(midoXpler$bandp) midoXbrho<-coh(mido.cln, brho.cln, years, norm="powall", sigmethod="fast", nrand=10000) midoXbrho<-bandtest(midoXbrho, st); midoXbrho<-bandtest(midoXbrho, lt) print(midoXbrho$bandp) midoXlaca<-coh(mido.cln, laca.cln, years, norm="powall", sigmethod="fast", nrand=10000) midoXlaca<-bandtest(midoXlaca, st); midoXlaca<-bandtest(midoXlaca, lt) print(midoXlaca$bandp) midoXvumi<-coh(mido.cln, vumi.cln, years, norm="powall", sigmethod="fast", nrand=10000) midoXvumi<-bandtest(midoXvumi, st); midoXvumi<-bandtest(midoXvumi, lt) print(midoXvumi$bandp) midoXcamu<-coh(mido.cln, camu.cln, years, norm="powall", sigmethod="fast", nrand=10000) midoXcamu<-bandtest(midoXcamu, st); midoXcamu<-bandtest(midoXcamu, lt) print(midoXcamu$bandp) #Vulpia vumiXpler<-coh(vumi.cln, pler.cln, years, norm="powall", sigmethod="fast", nrand=10000) vumiXpler<-bandtest(vumiXpler, st); vumiXpler<-bandtest(vumiXpler, lt) print(vumiXpler$bandp) vumiXbrho<-coh(vumi.cln, brho.cln, years, norm="powall", sigmethod="fast", nrand=10000) vumiXbrho<-bandtest(vumiXbrho, st); vumiXbrho<-bandtest(vumiXbrho, lt) print(vumiXbrho$bandp) vumiXlaca<-coh(vumi.cln, laca.cln, years, norm="powall", sigmethod="fast", nrand=10000) vumiXlaca<-bandtest(vumiXlaca, st); vumiXlaca<-bandtest(vumiXlaca, lt) print(vumiXlaca$bandp) vumiXmido<-coh(vumi.cln, mido.cln, years, norm="powall", sigmethod="fast", nrand=10000) vumiXmido<-bandtest(vumiXmido, st); vumiXmido<-bandtest(vumiXmido, lt) print(vumiXmido$bandp) vumiXcamu<-coh(vumi.cln, camu.cln, years, norm="powall", sigmethod="fast", nrand=10000) vumiXcamu<-bandtest(vumiXcamu, st); vumiXcamu<-bandtest(vumiXcamu, lt) print(vumiXcamu$bandp) #Calycadenia camuXpler<-coh(camu.cln, pler.cln, years, norm="powall", sigmethod="fast", nrand=10000) camuXpler<-bandtest(camuXpler, st); camuXpler<-bandtest(camuXpler, lt) print(camuXpler$bandp) camuXbrho<-coh(camu.cln, brho.cln, years, norm="powall", sigmethod="fast", nrand=10000) camuXbrho<-bandtest(camuXbrho, st); camuXbrho<-bandtest(camuXbrho, lt) print(camuXbrho$bandp) camuXlaca<-coh(camu.cln, laca.cln, years, norm="powall", sigmethod="fast", nrand=10000) camuXlaca<-bandtest(camuXlaca, st); camuXlaca<-bandtest(camuXlaca, lt) print(camuXlaca$bandp) camuXmido<-coh(camu.cln, mido.cln, years, norm="powall", sigmethod="fast", nrand=10000) camuXmido<-bandtest(camuXmido, st); camuXmido<-bandtest(camuXmido, lt) print(camuXmido$bandp) camuXvumi<-coh(camu.cln, vumi.cln, years, norm="powall", sigmethod="fast", nrand=10000) camuXvumi<-bandtest(camuXvumi, st); camuXvumi<-bandtest(camuXvumi, lt) print(camuXvumi$bandp) ``` ```{r check predictor coherences, echo=TRUE, cache=TRUE} pdsiXfallp<-coh(PDSI.cln, fallprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) pdsiXfallp<-bandtest(pdsiXfallp, st); pdsiXfallp<-bandtest(pdsiXfallp, lt) print(pdsiXfallp$bandp) plotmag(pdsiXfallp) pdsiXwintp<-coh(PDSI.cln, wintprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) pdsiXwintp<-bandtest(pdsiXwintp, st); pdsiXwintp<-bandtest(pdsiXwintp, lt) print(pdsiXwintp$bandp) plotmag(pdsiXwintp) fallpXwintp<-coh(fallprcp.cln, wintprcp.cln, years, norm="powall", sigmethod="fast", nrand=1000) fallpXwintp<-bandtest(fallpXwintp, st); fallpXwintp<-bandtest(fallpXwintp, lt) print(fallpXwintp$bandp) plotmag(fallpXwintp) gophXpdsi<-coh(goph.cln, PDSI.cln, years, norm="powall", sigmethod = "fftsurrog2", nrand=10000) gophXpdsi<-bandtest(gophXpdsi, st); gophXpdsi<-bandtest(gophXpdsi, lt) print(gophXpdsi$bandp) gophXwintp<-coh(goph.cln, wintprcp.cln, years, norm="powall", sigmethod = "fftsurrog2", nrand=10000) gophXwintp<-bandtest(gophXwintp, st); gophXwintp<-bandtest(gophXwintp, lt) print(gophXwintp$bandp) gophXfallp<-coh(goph.cln, fallprcp.cln, years, norm="powall", sigmethod = "fftsurrog2", nrand=10000) gophXfallp<-bandtest(gophXfallp, st); gophXfallp<-bandtest(gophXfallp, lt) print(gophXfallp$bandp) plotmag(gophXfallp) ``` Next build multivariate wavelet models ````{r wlm, echo=TRUE, cache=TRUE} dlist<-list(pler.cln, brho.cln, laca.cln, mido.cln, vumi.cln, camu.cln, goph.cln, PDSI.cln, fallprcp.cln, wintprcp.cln) bandavg<-function(syncexpl.obj, tsrange){ return(colMeans(syncexpl.obj[syncexpl.obj$timescales > min(tsrange) & syncexpl.obj$timescales < max(tsrange),], na.rm=T)) } ## Short Timescales ############################## #Plantago wlm.plerXbrho.st<-wlm(dlist, years, resp=1, pred=c(2,7,8,9), norm="powall") wlmtest.plerXbrho.st<-wlmtest(wlm.plerXbrho.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.plerXbrho.st<-bandtest(wlmtest.plerXbrho.st, st) print(wlmtest.plerXbrho.st$bandp) print(bandavg(syncexpl(wlm.plerXbrho.st), st)) wlm.plerXlaca.st<-wlm(dlist, years, resp=1, pred=c(3,7,8,9), norm="powall") wlmtest.plerXlaca.st<-wlmtest(wlm.plerXlaca.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.plerXlaca.st<-bandtest(wlmtest.plerXlaca.st, st) print(wlmtest.plerXlaca.st$bandp) print(bandavg(syncexpl(wlm.plerXlaca.st), st)) wlm.plerXmido.st<-wlm(dlist, years, resp=1, pred=c(4,7,8,9), norm="powall") wlmtest.plerXmido.st<-wlmtest(wlm.plerXmido.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.plerXmido.st<-bandtest(wlmtest.plerXmido.st, st) print(wlmtest.plerXmido.st$bandp) print(bandavg(syncexpl(wlm.plerXmido.st), st)) wlm.plerXvumi.st<-wlm(dlist, years, resp=1, pred=c(5,7,8,9), norm="powall") wlmtest.plerXvumi.st<-wlmtest(wlm.plerXvumi.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.plerXvumi.st<-bandtest(wlmtest.plerXvumi.st, st) print(wlmtest.plerXvumi.st$bandp) print(bandavg(syncexpl(wlm.plerXvumi.st), st)) wlm.plerXcamu.st<-wlm(dlist, years, resp=1, pred=c(6,7,8,9), norm="powall") wlmtest.plerXcamu.st<-wlmtest(wlm.plerXcamu.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.plerXcamu.st<-bandtest(wlmtest.plerXcamu.st, st) print(wlmtest.plerXcamu.st$bandp) print(bandavg(syncexpl(wlm.plerXcamu.st), st)) #Bromus--no control variables were significant wlm.brhoXpler.st<-wlm(dlist, years, resp=2, pred=c(1), norm="powall") wlmtest.brhoXpler.st<-wlmtest(wlm.brhoXpler.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.brhoXpler.st<-bandtest(wlmtest.brhoXpler.st, st) print(wlmtest.brhoXpler.st$bandp) print(bandavg(syncexpl(wlm.brhoXpler.st), st)) wlm.brhoXlaca.st<-wlm(dlist, years, resp=2, pred=c(3), norm="powall") wlmtest.brhoXlaca.st<-wlmtest(wlm.brhoXlaca.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.brhoXlaca.st<-bandtest(wlmtest.brhoXlaca.st, st) print(wlmtest.brhoXlaca.st$bandp) print(bandavg(syncexpl(wlm.brhoXlaca.st), st)) wlm.brhoXmido.st<-wlm(dlist, years, resp=2, pred=c(4), norm="powall") wlmtest.brhoXmido.st<-wlmtest(wlm.brhoXmido.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.brhoXmido.st<-bandtest(wlmtest.brhoXmido.st, st) print(wlmtest.brhoXmido.st$bandp) print(bandavg(syncexpl(wlm.brhoXmido.st), st)) wlm.brhoXvumi.st<-wlm(dlist, years, resp=2, pred=c(5), norm="powall") wlmtest.brhoXvumi.st<-wlmtest(wlm.brhoXvumi.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.brhoXvumi.st<-bandtest(wlmtest.brhoXvumi.st, st) print(wlmtest.brhoXvumi.st$bandp) print(bandavg(syncexpl(wlm.brhoXvumi.st), st)) wlm.brhoXcamu.st<-wlm(dlist, years, resp=2, pred=c(6), norm="powall") wlmtest.brhoXcamu.st<-wlmtest(wlm.brhoXcamu.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.brhoXcamu.st<-bandtest(wlmtest.brhoXcamu.st, st) print(wlmtest.brhoXcamu.st$bandp) print(bandavg(syncexpl(wlm.brhoXcamu.st), st)) #Lasthenia wlm.lacaXpler.st<-wlm(dlist, years, resp=3, pred=c(1,7), norm="powall") wlmtest.lacaXpler.st<-wlmtest(wlm.lacaXpler.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.lacaXpler.st<-bandtest(wlmtest.lacaXpler.st, st) print(wlmtest.lacaXpler.st$bandp) print(bandavg(syncexpl(wlm.lacaXpler.st), st)) wlm.lacaXbrho.st<-wlm(dlist, years, resp=3, pred=c(2,7), norm="powall") wlmtest.lacaXbrho.st<-wlmtest(wlm.lacaXbrho.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.lacaXbrho.st<-bandtest(wlmtest.lacaXbrho.st, st) print(wlmtest.lacaXbrho.st$bandp) print(bandavg(syncexpl(wlm.lacaXbrho.st), st)) wlm.lacaXmido.st<-wlm(dlist, years, resp=3, pred=c(4,7), norm="powall") wlmtest.lacaXmido.st<-wlmtest(wlm.lacaXmido.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.lacaXmido.st<-bandtest(wlmtest.lacaXmido.st, st) print(wlmtest.lacaXmido.st$bandp) print(bandavg(syncexpl(wlm.lacaXmido.st), st)) wlm.lacaXvumi.st<-wlm(dlist, years, resp=3, pred=c(5,7), norm="powall") wlmtest.lacaXvumi.st<-wlmtest(wlm.lacaXvumi.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.lacaXvumi.st<-bandtest(wlmtest.lacaXvumi.st, st) print(wlmtest.lacaXvumi.st$bandp) print(bandavg(syncexpl(wlm.lacaXvumi.st), st)) wlm.lacaXcamu.st<-wlm(dlist, years, resp=3, pred=c(6,7), norm="powall") wlmtest.lacaXcamu.st<-wlmtest(wlm.lacaXcamu.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.lacaXcamu.st<-bandtest(wlmtest.lacaXcamu.st, st) print(wlmtest.lacaXcamu.st$bandp) print(bandavg(syncexpl(wlm.lacaXcamu.st), st)) #Microseris wlm.midoXpler.st<-wlm(dlist, years, resp=4, pred=c(1,8), norm="powall") wlmtest.midoXpler.st<-wlmtest(wlm.midoXpler.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.midoXpler.st<-bandtest(wlmtest.midoXpler.st, st) print(wlmtest.midoXpler.st$bandp) print(bandavg(syncexpl(wlm.midoXpler.st), st)) wlm.midoXbrho.st<-wlm(dlist, years, resp=4, pred=c(2,8), norm="powall") wlmtest.midoXbrho.st<-wlmtest(wlm.midoXbrho.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.midoXbrho.st<-bandtest(wlmtest.midoXbrho.st, st) print(wlmtest.midoXbrho.st$bandp) print(bandavg(syncexpl(wlm.midoXbrho.st), st)) wlm.midoXlaca.st<-wlm(dlist, years, resp=4, pred=c(3,8), norm="powall") wlmtest.midoXlaca.st<-wlmtest(wlm.midoXlaca.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.midoXlaca.st<-bandtest(wlmtest.midoXlaca.st, st) print(wlmtest.midoXlaca.st$bandp) print(bandavg(syncexpl(wlm.midoXlaca.st), st)) wlm.midoXvumi.st<-wlm(dlist, years, resp=4, pred=c(5,8), norm="powall") wlmtest.midoXvumi.st<-wlmtest(wlm.midoXvumi.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.midoXvumi.st<-bandtest(wlmtest.midoXvumi.st, st) print(wlmtest.midoXvumi.st$bandp) print(bandavg(syncexpl(wlm.midoXvumi.st), st)) wlm.midoXcamu.st<-wlm(dlist, years, resp=4, pred=c(6,8), norm="powall") wlmtest.midoXcamu.st<-wlmtest(wlm.midoXcamu.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.midoXcamu.st<-bandtest(wlmtest.midoXcamu.st, st) print(wlmtest.midoXcamu.st$bandp) print(bandavg(syncexpl(wlm.midoXcamu.st), st)) #Vulpia--no control variables were significant wlm.vumiXpler.st<-wlm(dlist, years, resp=5, pred=c(1), norm="powall") wlmtest.vumiXpler.st<-wlmtest(wlm.vumiXpler.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.vumiXpler.st<-bandtest(wlmtest.vumiXpler.st, st) print(wlmtest.vumiXpler.st$bandp) print(bandavg(syncexpl(wlm.vumiXpler.st), st)) wlm.vumiXbrho.st<-wlm(dlist, years, resp=5, pred=c(2), norm="powall") wlmtest.vumiXbrho.st<-wlmtest(wlm.vumiXbrho.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.vumiXbrho.st<-bandtest(wlmtest.vumiXbrho.st, st) print(wlmtest.vumiXbrho.st$bandp) print(bandavg(syncexpl(wlm.vumiXbrho.st), st)) wlm.vumiXlaca.st<-wlm(dlist, years, resp=5, pred=c(3), norm="powall") wlmtest.vumiXlaca.st<-wlmtest(wlm.vumiXlaca.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.vumiXlaca.st<-bandtest(wlmtest.vumiXlaca.st, st) print(wlmtest.vumiXlaca.st$bandp) print(bandavg(syncexpl(wlm.vumiXlaca.st), st)) wlm.vumiXmido.st<-wlm(dlist, years, resp=5, pred=c(4), norm="powall") wlmtest.vumiXmido.st<-wlmtest(wlm.vumiXmido.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.vumiXmido.st<-bandtest(wlmtest.vumiXmido.st, st) print(wlmtest.vumiXmido.st$bandp) print(bandavg(syncexpl(wlm.vumiXmido.st), st)) wlm.vumiXcamu.st<-wlm(dlist, years, resp=5, pred=c(6), norm="powall") wlmtest.vumiXcamu.st<-wlmtest(wlm.vumiXcamu.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.vumiXcamu.st<-bandtest(wlmtest.vumiXcamu.st, st) print(wlmtest.vumiXcamu.st$bandp) print(bandavg(syncexpl(wlm.vumiXcamu.st), st)) #Calycadenia wlm.camuXpler.st<-wlm(dlist, years, resp=6, pred=c(1,7,8), norm="powall") wlmtest.camuXpler.st<-wlmtest(wlm.camuXpler.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.camuXpler.st<-bandtest(wlmtest.camuXpler.st, st) print(wlmtest.camuXpler.st$bandp) print(bandavg(syncexpl(wlm.camuXpler.st), st)) wlm.camuXbrho.st<-wlm(dlist, years, resp=6, pred=c(2,7,8), norm="powall") wlmtest.camuXbrho.st<-wlmtest(wlm.camuXbrho.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.camuXbrho.st<-bandtest(wlmtest.camuXbrho.st, st) print(wlmtest.camuXbrho.st$bandp) print(bandavg(syncexpl(wlm.camuXbrho.st), st)) wlm.camuXlaca.st<-wlm(dlist, years, resp=6, pred=c(3,7,8), norm="powall") wlmtest.camuXlaca.st<-wlmtest(wlm.camuXlaca.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.camuXlaca.st<-bandtest(wlmtest.camuXlaca.st, st) print(wlmtest.camuXlaca.st$bandp) print(bandavg(syncexpl(wlm.camuXlaca.st), st)) wlm.camuXmido.st<-wlm(dlist, years, resp=6, pred=c(4,7,8), norm="powall") wlmtest.camuXmido.st<-wlmtest(wlm.camuXmido.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.camuXmido.st<-bandtest(wlmtest.camuXmido.st, st) print(wlmtest.camuXmido.st$bandp) print(bandavg(syncexpl(wlm.camuXmido.st), st)) wlm.camuXvumi.st<-wlm(dlist, years, resp=6, pred=c(5,7,8), norm="powall") wlmtest.camuXvumi.st<-wlmtest(wlm.camuXvumi.st, drop=2, sigmethod="fft", nrand=10000) wlmtest.camuXvumi.st<-bandtest(wlmtest.camuXvumi.st, st) print(wlmtest.camuXvumi.st$bandp) print(bandavg(syncexpl(wlm.camuXvumi.st), st)) ## Long Timescales #################################### #Plantago wlm.plerXbrho.lt<-wlm(dlist, years, resp=2, pred=c(1), norm="powall") wlmtest.plerXbrho.lt<-wlmtest(wlm.plerXbrho.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.plerXbrho.lt<-bandtest(wlmtest.plerXbrho.lt, lt) print(wlmtest.plerXbrho.lt$bandp) print(bandavg(syncexpl(wlm.plerXbrho.lt), lt)) wlm.plerXlaca.lt<-wlm(dlist, years, resp=2, pred=c(3), norm="powall") wlmtest.plerXlaca.lt<-wlmtest(wlm.plerXlaca.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.plerXlaca.lt<-bandtest(wlmtest.plerXlaca.lt, lt) print(wlmtest.plerXlaca.lt$bandp) print(bandavg(syncexpl(wlm.plerXlaca.lt), lt)) wlm.plerXmido.lt<-wlm(dlist, years, resp=2, pred=c(4), norm="powall") wlmtest.plerXmido.lt<-wlmtest(wlm.plerXmido.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.plerXmido.lt<-bandtest(wlmtest.plerXmido.lt, lt) print(wlmtest.plerXmido.lt$bandp) print(bandavg(syncexpl(wlm.plerXmido.lt), lt)) wlm.plerXvumi.lt<-wlm(dlist, years, resp=2, pred=c(5), norm="powall") wlmtest.plerXvumi.lt<-wlmtest(wlm.plerXvumi.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.plerXvumi.lt<-bandtest(wlmtest.plerXvumi.lt, lt) print(wlmtest.plerXvumi.lt$bandp) print(bandavg(syncexpl(wlm.plerXvumi.lt), lt)) wlm.plerXcamu.lt<-wlm(dlist, years, resp=2, pred=c(6), norm="powall") wlmtest.plerXcamu.lt<-wlmtest(wlm.plerXcamu.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.plerXcamu.lt<-bandtest(wlmtest.plerXcamu.lt, lt) print(wlmtest.plerXcamu.lt$bandp) print(bandavg(syncexpl(wlm.plerXcamu.lt), lt)) #Bromus wlmtest.brhoXpler.lt<-bandtest(wlmtest.brhoXpler.st, lt) print(wlmtest.brhoXpler.lt$bandp) print(bandavg(syncexpl(wlm.brhoXpler.st), lt)) wlmtest.brhoXlaca.lt<-bandtest(wlmtest.brhoXlaca.st, lt) print(wlmtest.brhoXlaca.lt$bandp) print(bandavg(syncexpl(wlm.brhoXlaca.st), lt)) wlmtest.brhoXmido.lt<-bandtest(wlmtest.brhoXmido.st, lt) print(wlmtest.brhoXmido.lt$bandp) print(bandavg(syncexpl(wlm.brhoXmido.st), lt)) wlmtest.brhoXvumi.lt<-bandtest(wlmtest.brhoXvumi.st, lt) print(wlmtest.brhoXvumi.lt$bandp) print(bandavg(syncexpl(wlm.brhoXvumi.st), lt)) wlmtest.brhoXcamu.lt<-bandtest(wlmtest.brhoXcamu.st, lt) print(wlmtest.brhoXcamu.lt$bandp) print(bandavg(syncexpl(wlm.brhoXcamu.st), lt)) #Lasthenia wlm.lacaXpler.lt<-wlm(dlist, years, resp=3, pred=c(1), norm="powall") wlmtest.lacaXpler.lt<-wlmtest(wlm.lacaXpler.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.lacaXpler.lt<-bandtest(wlmtest.lacaXpler.lt, lt) print(wlmtest.lacaXpler.lt$bandp) print(bandavg(syncexpl(wlm.lacaXpler.lt), lt)) wlm.lacaXbrho.lt<-wlm(dlist, years, resp=3, pred=c(2), norm="powall") wlmtest.lacaXbrho.lt<-wlmtest(wlm.lacaXbrho.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.lacaXbrho.lt<-bandtest(wlmtest.lacaXbrho.lt, lt) print(wlmtest.lacaXbrho.lt$bandp) print(bandavg(syncexpl(wlm.lacaXbrho.lt), lt)) wlm.lacaXmido.lt<-wlm(dlist, years, resp=3, pred=c(4), norm="powall") wlmtest.lacaXmido.lt<-wlmtest(wlm.lacaXmido.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.lacaXmido.lt<-bandtest(wlmtest.lacaXmido.lt, lt) print(wlmtest.lacaXmido.lt$bandp) print(bandavg(syncexpl(wlm.lacaXmido.lt), lt)) wlm.lacaXvumi.lt<-wlm(dlist, years, resp=3, pred=c(5), norm="powall") wlmtest.lacaXvumi.lt<-wlmtest(wlm.lacaXvumi.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.lacaXvumi.lt<-bandtest(wlmtest.lacaXvumi.lt, lt) print(wlmtest.lacaXvumi.lt$bandp) print(bandavg(syncexpl(wlm.lacaXvumi.lt), lt)) wlm.lacaXcamu.lt<-wlm(dlist, years, resp=3, pred=c(6), norm="powall") wlmtest.lacaXcamu.lt<-wlmtest(wlm.lacaXcamu.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.lacaXcamu.lt<-bandtest(wlmtest.lacaXcamu.lt, lt) print(wlmtest.lacaXcamu.lt$bandp) print(bandavg(syncexpl(wlm.lacaXcamu.lt), lt)) #Microseris wlm.midoXpler.lt<-wlm(dlist, years, resp=4, pred=c(1), norm="powall") wlmtest.midoXpler.lt<-wlmtest(wlm.midoXpler.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.midoXpler.lt<-bandtest(wlmtest.midoXpler.lt, lt) print(wlmtest.midoXpler.lt$bandp) print(bandavg(syncexpl(wlm.midoXpler.lt), lt)) wlm.midoXbrho.lt<-wlm(dlist, years, resp=4, pred=c(2), norm="powall") wlmtest.midoXbrho.lt<-wlmtest(wlm.midoXbrho.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.midoXbrho.lt<-bandtest(wlmtest.midoXbrho.lt, lt) print(wlmtest.midoXbrho.lt$bandp) print(bandavg(syncexpl(wlm.midoXbrho.lt), lt)) wlm.midoXlaca.lt<-wlm(dlist, years, resp=4, pred=c(3), norm="powall") wlmtest.midoXlaca.lt<-wlmtest(wlm.midoXlaca.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.midoXlaca.lt<-bandtest(wlmtest.midoXlaca.lt, lt) print(wlmtest.midoXlaca.lt$bandp) print(bandavg(syncexpl(wlm.midoXlaca.lt), lt)) wlm.midoXvumi.lt<-wlm(dlist, years, resp=4, pred=c(5), norm="powall") wlmtest.midoXvumi.lt<-wlmtest(wlm.midoXvumi.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.midoXvumi.lt<-bandtest(wlmtest.midoXvumi.lt, lt) print(wlmtest.midoXvumi.lt$bandp) print(bandavg(syncexpl(wlm.midoXvumi.lt), lt)) wlm.midoXcamu.lt<-wlm(dlist, years, resp=4, pred=c(6), norm="powall") wlmtest.midoXcamu.lt<-wlmtest(wlm.midoXcamu.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.midoXcamu.lt<-bandtest(wlmtest.midoXcamu.lt, lt) print(wlmtest.midoXcamu.lt$bandp) print(bandavg(syncexpl(wlm.midoXcamu.lt), lt)) #Vulpia wlmtest.vumiXpler.lt<-bandtest(wlmtest.vumiXpler.st, lt) print(wlmtest.vumiXpler.lt$bandp) print(bandavg(syncexpl(wlm.vumiXpler.st), lt)) wlmtest.vumiXbrho.lt<-bandtest(wlmtest.vumiXbrho.st, lt) print(wlmtest.vumiXbrho.lt$bandp) print(bandavg(syncexpl(wlm.vumiXbrho.st), lt)) wlmtest.vumiXlaca.lt<-bandtest(wlmtest.vumiXlaca.st, lt) print(wlmtest.vumiXlaca.lt$bandp) print(bandavg(syncexpl(wlm.vumiXlaca.st), lt)) wlmtest.vumiXmido.lt<-bandtest(wlmtest.vumiXmido.st, lt) print(wlmtest.vumiXmido.lt$bandp) print(bandavg(syncexpl(wlm.vumiXmido.st), lt)) wlmtest.vumiXcamu.lt<-bandtest(wlmtest.vumiXcamu.st, lt) print(wlmtest.vumiXcamu.lt$bandp) print(bandavg(syncexpl(wlm.vumiXcamu.st), lt)) #Calycadenia wlm.camuXpler.lt<-wlm(dlist, years, resp=6, pred=c(1), norm="powall") wlmtest.camuXpler.lt<-wlmtest(wlm.camuXpler.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.camuXpler.lt<-bandtest(wlmtest.camuXpler.lt, lt) print(wlmtest.camuXpler.lt$bandp) print(bandavg(syncexpl(wlm.camuXpler.lt), lt)) wlm.camuXbrho.lt<-wlm(dlist, years, resp=6, pred=c(2), norm="powall") wlmtest.camuXbrho.lt<-wlmtest(wlm.camuXbrho.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.camuXbrho.lt<-bandtest(wlmtest.camuXbrho.lt, lt) print(wlmtest.camuXbrho.lt$bandp) print(bandavg(syncexpl(wlm.camuXbrho.lt), lt)) wlm.camuXlaca.lt<-wlm(dlist, years, resp=6, pred=c(3), norm="powall") wlmtest.camuXlaca.lt<-wlmtest(wlm.camuXlaca.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.camuXlaca.lt<-bandtest(wlmtest.camuXlaca.lt, lt) print(wlmtest.camuXlaca.lt$bandp) print(bandavg(syncexpl(wlm.camuXlaca.lt), lt)) wlm.camuXmido.lt<-wlm(dlist, years, resp=6, pred=c(4), norm="powall") wlmtest.camuXmido.lt<-wlmtest(wlm.camuXmido.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.camuXmido.lt<-bandtest(wlmtest.camuXmido.lt, lt) print(wlmtest.camuXmido.lt$bandp) print(bandavg(syncexpl(wlm.camuXmido.lt), lt)) wlm.camuXvumi.lt<-wlm(dlist, years, resp=6, pred=c(5), norm="powall") wlmtest.camuXvumi.lt<-wlmtest(wlm.camuXvumi.lt, drop=2, sigmethod="fft", nrand=10000) wlmtest.camuXvumi.lt<-bandtest(wlmtest.camuXvumi.lt, lt) print(wlmtest.camuXvumi.lt$bandp) print(bandavg(syncexpl(wlm.camuXvumi.lt), lt)) ```` And finally, coherence and wavelet multiple regression for aggregate plant cover ````{r aggcover_coh, echo=TRUE, cache=TRUE} ## coherences paggXgoph<-coh(pagg.cln, goph.cln, years, norm="powall", sigmethod="aaftsurrog2", nrand=10000) paggXgoph<-bandtest(paggXgoph, st); paggXgoph<-bandtest(paggXgoph, lt) print(paggXgoph$bandp) paggXpdsi<-coh(pagg.cln, PDSI.cln, years, norm="powall", sigmethod="fast", nrand=10000) paggXpdsi<-bandtest(paggXpdsi, st); paggXpdsi<-bandtest(paggXpdsi, lt) print(paggXpdsi$bandp) paggXfallp<-coh(pagg.cln, fallprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) paggXfallp<-bandtest(paggXfallp, st); paggXfallp<-bandtest(paggXfallp, lt) print(paggXfallp$bandp) paggXwintp<-coh(pagg.cln, wintprcp.cln, years, norm="powall", sigmethod="fast", nrand=10000) paggXwintp<-bandtest(paggXwintp, st); paggXwintp<-bandtest(paggXwintp, lt) print(paggXwintp$bandp) ##wsynchrony explained dlist<-list(pagg.cln, goph.cln, PDSI.cln, fallprcp.cln, wintprcp.cln) wlm.paggXgoph<-wlm(dlist, years, resp=1, pred=2, norm="powall") print(round(bandavg(syncexpl(wlm.paggXgoph), st),4)) print(round(bandavg(syncexpl(wlm.paggXgoph), lt),4)) wlm.paggXpdsi<-wlm(dlist, years, resp=1, pred=3, norm="powall") print(round(bandavg(syncexpl(wlm.paggXpdsi), st),4)) print(round(bandavg(syncexpl(wlm.paggXpdsi), lt),4)) wlm.paggXfallp<-wlm(dlist, years, resp=1, pred=4, norm="powall") print(round(bandavg(syncexpl(wlm.paggXfallp), st),4)) print(round(bandavg(syncexpl(wlm.paggXfallp), lt),4)) wlm.paggXwintp<-wlm(dlist, years, resp=1, pred=5, norm="powall") print(round(bandavg(syncexpl(wlm.paggXwintp), st),4)) print(round(bandavg(syncexpl(wlm.paggXwintp), lt),4)) ```` Make pedagogical figure for supplement ```{r pedagogical figure, echo=FALSE, cache=FALSE} rm(list=ls()) library(mvtnorm) set.seed(101) #***A function to be used for the plot #A function which can simulate the dynamics of three habitat patches, and can be used #to illustrate the four mechanisms of the geography of synchrony of Walter et al (2017), #and the difference between a cause of synchrony and a cause of the geography of #synchrony. # #Args #r A length 4 vector. Growth rates in each patch. #K A length 4 vector. Carrying capacities in each patch. #alpha A length 4 vector. Strength of influence of the first noise in each patch. #beta A length 4 vector. Strength of influence of the second noise in each patch. #Sigdelta Covariance matrix of first noise. 4 by 4. #Sigepsilon Covariance matrix of second noise. 4 by 4. #D Dispersal matrix. 4 by 4. #burnin Length of burnin to discard #lensims Length of simulations to retain after discarding burnin #numsims Number of simulations to do # #Output - a named list consisting of the following entries #delta The noise delta used #epsilon The noise epsilon used #pops The populations, an array, lensims by 4 by numsims # #Details #K is used for an initial condition # simtor<-function(r,K,alpha,beta,Sigdelta,Sigepsilon,D,burnin,lensims,numsims) { res<-array(NA,c(lensims+burnin,4,numsims)) res[1,,]<-rep(K,times=numsims) delta<-rmvnorm(numsims*(lensims+burnin-1),mean=rep(0,4),sigma=Sigdelta) delta<-array(delta,c(lensims+burnin-1,numsims,4)) delta<-aperm(delta,c(1,3,2)) epsilon<-rmvnorm(numsims*(lensims+burnin-1),mean=rep(0,4),sigma=Sigepsilon) epsilon<-array(epsilon,c(lensims+burnin-1,numsims,4)) epsilon<-aperm(epsilon,c(1,3,2)) for (counter in 2:(lensims+burnin)) { rm<-matrix(rep(r,times=numsims),4,numsims) Km<-matrix(rep(K,times=numsims),4,numsims) alpham<-matrix(rep(alpha,times=numsims),4,numsims) betam<-matrix(rep(beta,times=numsims),4,numsims) res[counter,,]<-D%*%(res[counter-1,,]*exp(rm*(1-res[counter-1,,]/Km)+alpham*delta[counter-1,,]+betam*epsilon[counter-1,,])) } res<-res[(burnin+1):(lensims+burnin),,] delta<-delta[(burnin+1):(lensims+burnin-1),,] epsilon<-epsilon[(burnin+1):(lensims+burnin-1),,] return(list(delta=delta,epsilon=epsilon,pops=res)) } #***Tests of the function, comment out after passing # r<-c(.1,.5,.9,1.1) # K<-c(1,10,100,1000) # alpha<-c(.1,.2,.3,.4) # beta<-c(.4,.3,.2,.1) # Sigdelta<-matrix(c(1,.2,.2,.1, # .2,1,.4,.4, # .2,.4,1,.3, # .1,.4,.3,1),4,4) # eigen(Sigdelta,only.values = TRUE)$values # Sigepsilon<-matrix(c(2,.2,.2,.1, # .2,1,.4,.4, # .2,.4,2,.3, # .1,.4,.3,1),4,4) # eigen(Sigepsilon,only.values = TRUE)$values # D<-diag(4) # burnin<-100 # lensims<-100 # numsims<-3 # res<-simtor(r,K,alpha,beta,Sigdelta,Sigepsilon,D,burnin,lensims,numsims) # pops10<-res$pops[10,,] # pops11<-res$pops[11,,] # delta10<-res$delta[10,,] # epsilon10<-res$epsilon[10,,] # #the following should all be 4 # sum(pops11[,1]==D%*%(pops10[,1]*exp(r*(1-pops10[,1]/K)+alpha*delta10[,1]+beta*epsilon10[,1]))) # sum(pops11[,2]==D%*%(pops10[,2]*exp(r*(1-pops10[,2]/K)+alpha*delta10[,2]+beta*epsilon10[,2]))) # sum(pops11[,3]==D%*%(pops10[,3]*exp(r*(1-pops10[,3]/K)+alpha*delta10[,3]+beta*epsilon10[,3]))) #***Now do a case where you get geography of synchrony from geography of synchrony in the driver. #Mechanism A from Walter et al 2017 #The same deterministic model params every (red noise), to make sure geography of synchrony #is not coming from different dd dynamics in the locations r<-rep(.5,4) K<-rep(100,4) r_panC<-r K_panC<-K #Make the populations only sensitive to delta, and all equally so alpha<-rep(1,4) beta<-rep(0,4) #Make delta well correlated across the A locations, and across the B locations, but less #between the A and B locations Sigdelta<-matrix(c(1,.8,.4,.4, .8,1,.4,.4, .4,.4,1,.8, .4,.4,.8,1),4,4,byrow = TRUE) eigen(Sigdelta,only.values = TRUE)$values #Make epsilon uncorrelated between all pairs of locations, though it does not matter because #the populations are not sensitive to epsilon Sigepsilon<-matrix(0,4,4) diag(Sigepsilon)<-1 eigen(Sigepsilon,only.values = TRUE)$values #No dispersal D<-diag(4) burnin<-250 lensims<-100000 numsims<-3 res1<-simtor(r,K,alpha,beta,Sigdelta,Sigepsilon,D,burnin,lensims,numsims) cor(res1$pops[,,1]) #So you see a geography of synchrony. The Moran effect causes the synchrony, and geographic #variation in the synchrony of the environmental driver causes the geography of synchrony. #***Now do a case where synchrony is from Moran, but geography is from variation in dd. Mechanism #B from Walter et al 2017 #Different r in the A and B locations r<-c(.25,.25,1.75,1.75) K<-rep(100,4) r_panH<-r K_panH<-K #Make the populations only sensitive to delta, and all equally so alpha<-rep(1,4) beta<-rep(0,4) #Make delta equally correlated across all locations Sigdelta<-matrix(.6,4,4) diag(Sigdelta)<-1 eigen(Sigdelta,only.values = TRUE)$values #Make epsilon uncorrelated between all pairs of locations, though it does not matter because #the populations are not sensitive to epsilon Sigepsilon<-matrix(0,4,4) diag(Sigepsilon)<-1 eigen(Sigepsilon,only.values = TRUE)$values #No dispersal D<-diag(4) burnin<-250 numsims<-3 res2<-simtor(r,K,alpha,beta,Sigdelta,Sigepsilon,D,burnin,lensims,numsims) cor(res2$pops[,,1]) #So you again see a geography of synchrony. The Moran effect causes the synchrony, but geographic #variation in dd causes the geography of synchrony. #***Now do a case where limiting factors of population growth vary spatially. Mechanism C from #Walter et al 2017 #The same deterministic model params every (red noise), to make sure geography of synchrony #is not coming from different dd dynamics in the locations r<-rep(.5,4) K<-rep(100,4) #Make the A populations more sensitive to delta and the B populations more sensitive to epsilon alpha<-c(1,1,.1,.1) beta<-c(.1,.1,1,1) #Make delta equally correlated across all locations Sigdelta<-matrix(.6,4,4) diag(Sigdelta)<-1 eigen(Sigdelta,only.values = TRUE)$values #Make epsilon the same way, but independent from delta #Make delta equally correlated across all locations Sigepsilon<-matrix(.6,4,4) diag(Sigepsilon)<-1 eigen(Sigepsilon,only.values = TRUE)$values #No dispersal D<-diag(4) burnin<-250 numsims<-3 res3<-simtor(r,K,alpha,beta,Sigdelta,Sigepsilon,D,burnin,lensims,numsims) cor(res3$pops[,,1]) #So you again see a geography of synchrony. The Moran effect causes the synchrony, but geographic #variation which populations are responding more to which synchronous environmental driver causes #the geography of synchrony #***Now do a case where dispersal has a geography. Mechanism D from Walter et al 2017. #The same deterministic model params everywhere, to make sure geography of synchrony #is not coming from different dd dynamics in the locations r<-rep(2.2,4) K<-rep(100,4) #Make the populations only sensitive to delta, and all equally so alpha<-rep(1,4) beta<-rep(0,4) #Make delta uncorrelated between all pairs of locations Sigdelta<-matrix(0,4,4) diag(Sigdelta)<-1 eigen(Sigdelta,only.values = TRUE)$values #Make epsilon uncorrelated between all pairs of locations, though it does not matter because #the populations are not sensitive to epsilon Sigepsilon<-matrix(0,4,4) diag(Sigepsilon)<-1 eigen(Sigepsilon,only.values = TRUE)$values #Dispersal is more between A locations and between B locations, but less between A and B locations D<-matrix(c(.7,.3,0,0, .3,.7,0,0, 0,0,.7,.3, 0,0,.3,.7),4,4,byrow = TRUE) burnin<-250 numsims<-3 res4<-simtor(r,K,alpha,beta,Sigdelta,Sigepsilon,D,burnin,lensims,numsims) cor(res4$pops[,,1]) #So you again see a geography of synchrony. Dispersal causes the synchrony, and geographic #variation in dispersal causes the geography of synchrony. #***Make a figure out of mechanisms A and B. It will only cover two of the mechanisms of #synchrony, but one case will be synchrony and its geography kind of coming "from the same #cause", and the other will be synchrony and its geography coming from completely unrelated #causes. #plot dimensions, units inches xmarght<-.5 ymargwd<-.5 totwd<-3.5 gap<-0.1 panwd<-(totwd-2*gap-ymargwd)/2 panht<-panwd headht<-.3 totht<-2*xmarght+5*panht+4*gap+headht textsz<-.9 smalltextsz<-.5 midtextsz<-.75 #pdf(file="NewPedagogFig.pdf",width=totwd,height=totht) png(file="./NewPedagogFig.png",res=300,width=totwd,height=totht,units="in") #other common plotting parameters tsplotlen<-25 simtoplot<-1 A1col<-"red" A2col<-"pink" B1col<-"blue" B2col<-"light blue" #example 1 below corresponds to mechansim A above #example 1, panel A pancolnum<-0 par(fig=c((ymargwd+pancolnum*(gap+panwd))/totwd, (ymargwd+pancolnum*(gap+panwd)+panwd)/totwd, (3*gap+4*panht+2*xmarght)/totht, (3*gap+5*panht+2*xmarght)/totht), mai=c(0,0,0,0),mgp=c(3,.15,0),tcl=-.25) xlimits_sq<-c(-1.25,1) ylimits_sq<-c(-1,1.25) plot(c(-1,1),c(0,0),type="n",lty="dashed",xlim=xlimits_sq,ylim=ylimits_sq, xaxt="n",yaxt="n",bty="o") lines(c(0,0),c(-1,1),type="n") text(-.5,1,expression(alpha),adj=c(.5,0)) text(.5,1,expression(beta),adj=c(.5,0)) text(-1,.5,"1",adj=c(1,.5)) text(-1,-.5,"2",adj=c(1,.5)) mtext("A",3,line=-1,adj=.025) points(-.5,.5,type="p",pch=20,cex=3,col=A1col) points(-.5,-.5,type="p",pch=20,cex=3,col=A2col) points(.5,.5,type="p",pch=20,cex=3,col=B1col) points(.5,-.5,type="p",pch=20,cex=3,col=B2col) cors<-c(cor(res1$delta[,1:2,simtoplot])[1,2], cor(res1$delta[,3:4,simtoplot])[1,2], cor(res1$delta[,c(1,3),simtoplot])[1,2], cor(res1$delta[,c(2,4),simtoplot])[1,2], cor(res1$delta[,c(1,4),simtoplot])[1,2], cor(res1$delta[,c(2,3),simtoplot])[1,2]) lwds<-5*cors lines(c(-.5,-.5),c(-.5,.5),lwd=lwds[1]) #connect the A locs text(-.5,0,round(cors[1],2),adj=c(1.2,.5),cex=smalltextsz) lines(c(.5,.5),c(-.5,.5),lwd=lwds[2]) #then the B locs text(.5,0,round(cors[2],2),adj=c(-.2,.5),cex=smalltextsz) lines(c(-.5,.5),c(.5,.5),lwd=lwds[3]) #then the 1 locs text(0,.5,round(cors[3],2),adj=c(.5,-.3),cex=smalltextsz) lines(c(-.5,.5),c(-.5,-.5),lwd=lwds[4]) #then the 2 locs text(0,-.5,round(cors[4],2),adj=c(.5,1.3),cex=smalltextsz) lines(c(-.5,.5),c(.5,-.5),lwd=lwds[5]) #then A1 to B2 text(-.25,.25,round(cors[5],2),adj=c(-0.05,-0.05),cex=smalltextsz) lines(c(-.5,.5),c(-.5,.5),lwd=lwds[6]) #then A2 to B1 text(-.25,-.25,round(cors[6],2),adj=c(-0.05,1.05),cex=smalltextsz) mtext("Sim. 1: Moran driver",3,line=1,cex=textsz) mtext("causing geog.",3,line=0.1,cex=textsz) mtext("Noise correlations",2,line=1.2,cex=textsz) mtext("across locations",2,line=0.4,cex=textsz) #example 1, panel B par(fig=c((ymargwd+pancolnum*(gap+panwd))/totwd, (ymargwd+pancolnum*(gap+panwd)+panwd)/totwd, (2*gap+3*panht+2*xmarght)/totht, (2*gap+4*panht+2*xmarght)/totht), mai=c(0,0,0,0),mgp=c(3,.15,0),tcl=-.25,new=TRUE) Adelt_panB<-res1$delta[1:(tsplotlen+1),1:2,simtoplot]+2 Bdelt_panB<-res1$delta[1:(tsplotlen+1),3:4,simtoplot]-2 Adelt_panG<-res2$delta[1:(tsplotlen+1),1:2,simtoplot]+2 Bdelt_panG<-res2$delta[1:(tsplotlen+1),3:4,simtoplot]-2 ylimits_panBG<-range(Adelt_panB,Bdelt_panB,Adelt_panG,Bdelt_panG) plot(0:tsplotlen,Adelt_panB[,1],type="l",col=A1col, ylim=ylimits_panBG,xaxt="n",yaxt="n") lines(0:tsplotlen,Adelt_panB[,2],type="l",col=A2col) lines(0:tsplotlen,Bdelt_panB[,1],type="l",col=B1col) lines(0:tsplotlen,Bdelt_panB[,2],type="l",col=B2col) axis(side=1,labels=TRUE,cex.axis=textsz) mtext("Time step",side=1,line=1.2,cex=textsz) mtext("Noise (with",2,line=1.2,cex=textsz) mtext("vertical shifts)",2,line=0.4,cex=textsz) mtext("B",3,line=-1,adj=.025) #example 1, panel C par(fig=c((ymargwd+pancolnum*(gap+panwd))/totwd, (ymargwd+pancolnum*(gap+panwd)+panwd)/totwd, (2*gap+2*panht+xmarght)/totht, (2*gap+3*panht+xmarght)/totht), mai=c(0,0,0,0),mgp=c(3,.15,0),tcl=-.25,new=TRUE) plot(c(-1,1),c(0,0),type="n",lty="dashed",xlim=xlimits_sq,ylim=ylimits_sq, xaxt="n",yaxt="n",bty="o") lines(c(0,0),c(-1,1),type="n") text(-.5,1,expression(alpha),adj=c(.5,0)) text(.5,1,expression(beta),adj=c(.5,0)) text(-1,.5,"1",adj=c(1,.5)) text(-1,-.5,"2",adj=c(1,.5)) mtext("C",3,line=-1,adj=.025) text(-.5,.5,paste0("r=",round(r_panC[1],2)),adj=c(.5,-.2),cex=midtextsz) text(-.5,.5,paste0("K=",round(K_panC[1],2)),adj=c(.5,1.2),cex=midtextsz) text(-.5,-.5,paste0("r=",round(r_panC[2],2)),adj=c(.5,-.2),cex=midtextsz) text(-.5,-.5,paste0("K=",round(K_panC[2],2)),adj=c(.5,1.2),cex=midtextsz) text(.5,.5,paste0("r=",round(r_panC[3],2)),adj=c(.5,-.2),cex=midtextsz) text(.5,.5,paste0("K=",round(K_panC[3],2)),adj=c(.5,1.2),cex=midtextsz) text(.5,-.5,paste0("r=",round(r_panC[4],2)),adj=c(.5,-.2),cex=midtextsz) text(.5,-.5,paste0("K=",round(K_panC[4],2)),adj=c(.5,1.2),cex=midtextsz) mtext("Model parameters",2,line=1.2,cex=textsz) mtext("in each patch",2,line=0.4,cex=textsz) #example 1, panel D par(fig=c((ymargwd+pancolnum*(gap+panwd))/totwd, (ymargwd+pancolnum*(gap+panwd)+panwd)/totwd, (gap+panht+xmarght)/totht, (gap+2*panht+xmarght)/totht), mai=c(0,0,0,0),mgp=c(3,.15,0),tcl=-.25,new=TRUE) Apops_panD<-sqrt(res1$pops[1:(tsplotlen+1),1:2,simtoplot])+30 Bpops_panD<-sqrt(res1$pops[1:(tsplotlen+1),3:4,simtoplot]) Apops_panI<-sqrt(res2$pops[1:(tsplotlen+1),1:2,simtoplot])+30 Bpops_panI<-sqrt(res2$pops[1:(tsplotlen+1),3:4,simtoplot]) ylimits_panDI<-range(Apops_panD,Bpops_panD,Apops_panI,Bpops_panI) plot(0:tsplotlen,Apops_panD[,1],type="l",col=A1col, ylim=ylimits_panDI,xaxt="n",yaxt="n") lines(0:tsplotlen,Apops_panD[,2],type="l",col=A2col) lines(0:tsplotlen,Bpops_panD[,1],type="l",col=B1col) lines(0:tsplotlen,Bpops_panD[,2],type="l",col=B2col) axis(side=1,labels=TRUE,cex.axis=textsz) mtext("Time step",side=1,line=1.2,cex=textsz) mtext("Sqrt. pops. (with",2,line=1.2,cex=textsz) mtext("vertical shifts)",2,line=0.4,cex=textsz) mtext("D",3,line=-1,adj=.025) #example 1, panel E par(fig=c((ymargwd+pancolnum*(gap+panwd))/totwd, (ymargwd+pancolnum*(gap+panwd)+panwd)/totwd, (gap)/totht, (gap+panht)/totht), mai=c(0,0,0,0),mgp=c(3,.15,0),tcl=-.25,new=TRUE) plot(c(-1,1),c(0,0),type="n",lty="dashed",xlim=xlimits_sq,ylim=ylimits_sq, xaxt="n",yaxt="n",bty="o") lines(c(0,0),c(-1,1),type="n") text(-.5,1,expression(alpha),adj=c(.5,0)) text(.5,1,expression(beta),adj=c(.5,0)) text(-1,.5,"1",adj=c(1,.5)) text(-1,-.5,"2",adj=c(1,.5)) mtext("E",3,line=-1,adj=.025) points(-.5,.5,type="p",pch=20,cex=3,col=A1col) points(-.5,-.5,type="p",pch=20,cex=3,col=A2col) points(.5,.5,type="p",pch=20,cex=3,col=B1col) points(.5,-.5,type="p",pch=20,cex=3,col=B2col) cors<-c(cor(res1$pops[,1:2,simtoplot])[1,2], cor(res1$pops[,3:4,simtoplot])[1,2], cor(res1$pops[,c(1,3),simtoplot])[1,2], cor(res1$pops[,c(2,4),simtoplot])[1,2], cor(res1$pops[,c(1,4),simtoplot])[1,2], cor(res1$pops[,c(2,3),simtoplot])[1,2]) lwds<-5*cors lines(c(-.5,-.5),c(-.5,.5),lwd=lwds[1]) #connect the A locs text(-.5,0,round(cors[1],2),adj=c(1.2,.5),cex=smalltextsz) lines(c(.5,.5),c(-.5,.5),lwd=lwds[2]) #then the B locs text(.5,0,round(cors[2],2),adj=c(-.2,.5),cex=smalltextsz) lines(c(-.5,.5),c(.5,.5),lwd=lwds[3]) #then the 1 locs text(0,.5,round(cors[3],2),adj=c(.5,-.3),cex=smalltextsz) lines(c(-.5,.5),c(-.5,-.5),lwd=lwds[4]) #then the 2 locs text(0,-.5,round(cors[4],2),adj=c(.5,1.3),cex=smalltextsz) lines(c(-.5,.5),c(.5,-.5),lwd=lwds[5]) #then A1 to B2 text(-.25,.25,round(cors[5],2),adj=c(-0.05,-0.05),cex=smalltextsz) lines(c(-.5,.5),c(-.5,.5),lwd=lwds[6]) #then A2 to B1 text(-.25,-.25,round(cors[6],2),adj=c(-0.05,1.05),cex=smalltextsz) mtext("Pop. correlations",2,line=1.2,cex=textsz) mtext("across locations",2,line=0.4,cex=textsz) #example 2, panel F pancolnum<-1 par(fig=c((ymargwd+pancolnum*(gap+panwd))/totwd, (ymargwd+pancolnum*(gap+panwd)+panwd)/totwd, (3*gap+4*panht+2*xmarght)/totht, (3*gap+5*panht+2*xmarght)/totht), mai=c(0,0,0,0),mgp=c(3,.15,0),tcl=-.25,new=TRUE) plot(c(-1,1),c(0,0),type="n",lty="dashed",xlim=xlimits_sq,ylim=ylimits_sq, xaxt="n",yaxt="n",bty="o") lines(c(0,0),c(-1,1),type="n") text(-.5,1,expression(alpha),adj=c(.5,0)) text(.5,1,expression(beta),adj=c(.5,0)) text(-1,.5,"1",adj=c(1,.5)) text(-1,-.5,"2",adj=c(1,.5)) mtext("F",3,line=-1,adj=.025) points(-.5,.5,type="p",pch=20,cex=3,col=A1col) points(-.5,-.5,type="p",pch=20,cex=3,col=A2col) points(.5,.5,type="p",pch=20,cex=3,col=B1col) points(.5,-.5,type="p",pch=20,cex=3,col=B2col) cors<-c(cor(res2$delta[,1:2,simtoplot])[1,2], cor(res2$delta[,3:4,simtoplot])[1,2], cor(res2$delta[,c(1,3),simtoplot])[1,2], cor(res2$delta[,c(2,4),simtoplot])[1,2], cor(res2$delta[,c(1,4),simtoplot])[1,2], cor(res2$delta[,c(2,3),simtoplot])[1,2]) lwds<-5*cors lines(c(-.5,-.5),c(-.5,.5),lwd=lwds[1]) #connect the A locs text(-.5,0,round(cors[1],2),adj=c(1.2,.5),cex=smalltextsz) lines(c(.5,.5),c(-.5,.5),lwd=lwds[2]) #then the B locs text(.5,0,round(cors[2],2),adj=c(-.2,.5),cex=smalltextsz) lines(c(-.5,.5),c(.5,.5),lwd=lwds[3]) #then the 1 locs text(0,.5,round(cors[3],2),adj=c(.5,-.3),cex=smalltextsz) lines(c(-.5,.5),c(-.5,-.5),lwd=lwds[4]) #then the 2 locs text(0,-.5,round(cors[4],2),adj=c(.5,1.3),cex=smalltextsz) lines(c(-.5,.5),c(.5,-.5),lwd=lwds[5]) #then A1 to B2 text(-.25,.25,round(cors[5],2),adj=c(-0.05,-0.05),cex=smalltextsz) lines(c(-.5,.5),c(-.5,.5),lwd=lwds[6]) #then A2 to B1 text(-.25,-.25,round(cors[6],2),adj=c(-0.05,1.05),cex=smalltextsz) mtext("Sim. 2: Moran driver",3,line=1,cex=textsz) mtext("not causing geog.",3,line=0.1,cex=textsz) #example 2, panel G par(fig=c((ymargwd+pancolnum*(gap+panwd))/totwd, (ymargwd+pancolnum*(gap+panwd)+panwd)/totwd, (2*gap+3*panht+2*xmarght)/totht, (2*gap+4*panht+2*xmarght)/totht), mai=c(0,0,0,0),mgp=c(3,.15,0),tcl=-.25,new=TRUE) plot(0:tsplotlen,Adelt_panG[,1],type="l",col=A1col, ylim=ylimits_panBG,xaxt="n",yaxt="n") lines(0:tsplotlen,Adelt_panG[,2],type="l",col=A2col) lines(0:tsplotlen,Bdelt_panG[,1],type="l",col=B1col) lines(0:tsplotlen,Bdelt_panG[,2],type="l",col=B2col) axis(side=1,labels=TRUE,cex.axis=textsz) mtext("Time step",side=1,line=1.2,cex=textsz) mtext("G",3,line=-1,adj=.025) #example 2, panel H par(fig=c((ymargwd+pancolnum*(gap+panwd))/totwd, (ymargwd+pancolnum*(gap+panwd)+panwd)/totwd, (2*gap+2*panht+xmarght)/totht, (2*gap+3*panht+xmarght)/totht), mai=c(0,0,0,0),mgp=c(3,.15,0),tcl=-.25,new=TRUE) plot(c(-1,1),c(0,0),type="n",lty="dashed",xlim=xlimits_sq,ylim=ylimits_sq, xaxt="n",yaxt="n",bty="o") lines(c(0,0),c(-1,1),type="n") text(-.5,1,expression(alpha),adj=c(.5,0)) text(.5,1,expression(beta),adj=c(.5,0)) text(-1,.5,"1",adj=c(1,.5)) text(-1,-.5,"2",adj=c(1,.5)) mtext("H",3,line=-1,adj=.025) text(-.5,.5,paste0("r=",round(r_panH[1],2)),adj=c(.5,-.2),cex=midtextsz) text(-.5,.5,paste0("K=",round(K_panH[1],2)),adj=c(.5,1.2),cex=midtextsz) text(-.5,-.5,paste0("r=",round(r_panH[2],2)),adj=c(.5,-.2),cex=midtextsz) text(-.5,-.5,paste0("K=",round(K_panH[2],2)),adj=c(.5,1.2),cex=midtextsz) text(.5,.5,paste0("r=",round(r_panH[3],2)),adj=c(.5,-.2),cex=midtextsz) text(.5,.5,paste0("K=",round(K_panH[3],2)),adj=c(.5,1.2),cex=midtextsz) text(.5,-.5,paste0("r=",round(r_panH[4],2)),adj=c(.5,-.2),cex=midtextsz) text(.5,-.5,paste0("K=",round(K_panH[4],2)),adj=c(.5,1.2),cex=midtextsz) #example 2, panel I par(fig=c((ymargwd+pancolnum*(gap+panwd))/totwd, (ymargwd+pancolnum*(gap+panwd)+panwd)/totwd, (gap+panht+xmarght)/totht, (gap+2*panht+xmarght)/totht), mai=c(0,0,0,0),mgp=c(3,.15,0),tcl=-.25,new=TRUE) plot(0:tsplotlen,Apops_panI[,1],type="l",col=A1col, ylim=ylimits_panDI,xaxt="n",yaxt="n") lines(0:tsplotlen,Apops_panI[,2],type="l",col=A2col) lines(0:tsplotlen,Bpops_panI[,1],type="l",col=B1col) lines(0:tsplotlen,Bpops_panI[,2],type="l",col=B2col) axis(side=1,labels=TRUE,cex.axis=textsz) mtext("Time step",side=1,line=1.2,cex=textsz) mtext("I",3,line=-1,adj=.025) #example 2, panel J par(fig=c((ymargwd+pancolnum*(gap+panwd))/totwd, (ymargwd+pancolnum*(gap+panwd)+panwd)/totwd, (gap)/totht, (gap+panht)/totht), mai=c(0,0,0,0),mgp=c(3,.15,0),tcl=-.25,new=TRUE) plot(c(-1,1),c(0,0),type="n",lty="dashed",xlim=xlimits_sq,ylim=ylimits_sq, xaxt="n",yaxt="n",bty="o") lines(c(0,0),c(-1,1),type="n") text(-.5,1,expression(alpha),adj=c(.5,0)) text(.5,1,expression(beta),adj=c(.5,0)) text(-1,.5,"1",adj=c(1,.5)) text(-1,-.5,"2",adj=c(1,.5)) mtext("J",3,line=-1,adj=.025) points(-.5,.5,type="p",pch=20,cex=3,col=A1col) points(-.5,-.5,type="p",pch=20,cex=3,col=A2col) points(.5,.5,type="p",pch=20,cex=3,col=B1col) points(.5,-.5,type="p",pch=20,cex=3,col=B2col) cors<-c(cor(res2$pops[,1:2,simtoplot])[1,2], cor(res2$pops[,3:4,simtoplot])[1,2], cor(res2$pops[,c(1,3),simtoplot])[1,2], cor(res2$pops[,c(2,4),simtoplot])[1,2], cor(res2$pops[,c(1,4),simtoplot])[1,2], cor(res2$pops[,c(2,3),simtoplot])[1,2]) lwds<-5*cors lines(c(-.5,-.5),c(-.5,.5),lwd=lwds[1]) #connect the A locs text(-.5,0,round(cors[1],2),adj=c(1.2,.5),cex=smalltextsz) lines(c(.5,.5),c(-.5,.5),lwd=lwds[2]) #then the B locs text(.5,0,round(cors[2],2),adj=c(-.2,.5),cex=smalltextsz) lines(c(-.5,.5),c(.5,.5),lwd=lwds[3]) #then the 1 locs text(0,.5,round(cors[3],2),adj=c(.5,-.3),cex=smalltextsz) lines(c(-.5,.5),c(-.5,-.5),lwd=lwds[4]) #then the 2 locs text(0,-.5,round(cors[4],2),adj=c(.5,1.3),cex=smalltextsz) lines(c(-.5,.5),c(.5,-.5),lwd=lwds[5]) #then A1 to B2 text(-.25,.25,round(cors[5],2),adj=c(-0.05,-0.05),cex=smalltextsz) lines(c(-.5,.5),c(-.5,.5),lwd=lwds[6]) #then A2 to B1 text(-.25,-.25,round(cors[6],2),adj=c(-0.05,1.05),cex=smalltextsz) dev.off() ```