---
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()
```
