---
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 "Micro-scale geography of synchrony in a serpentine plant community," by Jonathan Walter, Lauren Hallett, Lawrence Sheppard, Thomas Anderson, Lei Zhao, Richard Hobbs, Katherine Suding, and Daniel 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("plantago.csv", header=TRUE)
bromus<-read.csv("bromus.csv", header=TRUE)
lasthenia<-read.csv("lasthenia.csv", header=TRUE)
microseris<-read.csv("microseris.csv", header=TRUE)
vulpia<-read.csv("vulpia.csv", header=TRUE)
calycadenia<-read.csv("calycadenia.csv", header=TRUE)
agg.all<-read.csv("all_plants.csv",header=TRUE)
gophers<-read.csv("gophers.csv", header=TRUE)
plot.coords<-read.csv("plot_coords.csv", header=TRUE)
soil.depth<-read.csv("soil_depth.csv", header=TRUE)
PDSI.raw<-read.csv("pdsi.csv")
precip.raw<-read.csv("precip_monthly.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 <- precip.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)]])

# Distance decay, short timescales (plants)
tiff("./Fig3_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))
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)
abline(h=mean(lower(clust.pler.st$adj)),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=".", col=colors[lower(inout.pler.st)],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)
abline(h=mean(lower(clust.brho.st$adj)),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=".", col=colors[lower(inout.brho.st)],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)
abline(h=mean(lower(clust.laca.st$adj)),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=".", col=colors[inout.laca.st],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)
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=".", col=colors[lower(inout.mido.st)],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)
abline(h=mean(lower(clust.vumi.st$adj)),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=".", col=colors[lower(inout.vumi.st)],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)
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=".", col=colors[lower(inout.camu.st)],cex=2)
text(2,0.99,"f)")

mtext("Synchrony",2,outer=T,line=0.5, cex=1)
mtext("Distance (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))
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)
abline(h=mean(lower(clust.pler.lt$adj)),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=".", col=colors[lower(inout.pler.lt)],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)
abline(h=mean(lower(clust.brho.lt$adj)),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=".", col=colors[lower(inout.brho.lt)],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)
abline(h=mean(lower(clust.laca.lt$adj)),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=".", col=colors[lower(inout.laca.lt)],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)
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=".", col=colors[lower(inout.mido.lt)],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(-0.5,1), cex.main=1, lwd=2)
abline(h=mean(lower(clust.vumi.lt$adj)),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=".", col=colors[lower(inout.vumi.lt)],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)
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=".", col=colors[lower(inout.camu.lt)],cex=2)
text(2,0.99,"f)")

mtext("Synchrony",2,outer=T, cex=1)
mtext("Distance (m)",1, outer=T, line=0.5, cex=1)

dev.off()

# ----------------------------------------------------------------------
# 1b. Networks and modularity

source("./plotClusterMap_JRG.R")

#Networks, short timescales, plants
tiff("./Fig5_networks_ReXWT_short.tif", units="in", res=300, width = 6.5, height=4)
par(mfrow=c(2,3))

plotClusterMap(clust.pler.st, nodewgt="mod.decomp", edgewgt=0.9, title="Plantago", nodesize=c(0.7,2), edgesize = c(0.5,0.5))
text(1.5,28.5,"a)")
plotClusterMap(clust.brho.st, nodewgt="mod.decomp", edgewgt=0.9, title="Bromus", nodesize=c(0.7,2), edgesize = c(0.5,0.5))
text(1.5,28.5,"b)")
plotClusterMap(clust.laca.st, nodewgt="mod.decomp", edgewgt=0.9, title="Lasthenia", nodesize=c(0.7,2), edgesize = c(0.5,0.5))
text(1.5,28.5,"c)")
plotClusterMap(clust.mido.st, nodewgt="mod.decomp", edgewgt=0.9, title="Microseris", nodesize=c(0.7,2), edgesize = c(0.5,0.5))
text(1.5,28.5,"d)")
plotClusterMap(clust.vumi.st, nodewgt="mod.decomp", edgewgt=0.9, title="Vulpia", nodesize=c(0.7,2), edgesize = c(0.5,0.5))
text(1.5,28.5,"e)")
plotClusterMap(clust.camu.st, nodewgt="mod.decomp", edgewgt=0.9, title="Calycadenia", nodesize=c(0.7,2), edgesize = c(0.5,0.5))
text(1.5,28.5,"f)")

dev.off()

#Networks, long timescales, plants
tiff("./FigS2_networks_ReXWT_long.tif", units="in", res=300, width = 6.5, height=4)
par(mfrow=c(2,3), mar=c(1,1,1.5,6.1))

plotClusterMap(clust.pler.lt, nodewgt="mod.decomp", edgewgt=0.9, title="Plantago", nodesize=c(0.7,2), edgesize = c(0.5,0.5))
text(1.5,28.5,"a)")
plotClusterMap(clust.brho.lt, nodewgt="mod.decomp", edgewgt=0.9, title="Bromus", nodesize=c(0.7,2), edgesize = c(0.5,0.5), addlegend = T)
text(1.5,28.5,"b)")
plotClusterMap(clust.laca.lt, nodewgt="mod.decomp", edgewgt=0.9, title="Lasthenia", nodesize=c(0.7,2), edgesize = c(0.5,0.5), addlegend = T)
text(1.5,28.5,"c)")
plotClusterMap(clust.mido.lt, nodewgt="mod.decomp", edgewgt=0.9, title="Microseris", nodesize=c(0.7,2), edgesize = c(0.5,0.5))
text(1.5,28.5,"d)")
plotClusterMap(clust.vumi.lt, nodewgt="mod.decomp", edgewgt=0.9, title="Vulpia", nodesize=c(0.7,2), edgesize = c(0.5,0.5), addlegend = T)
text(1.5,28.5,"e)")
plotClusterMap(clust.camu.lt, nodewgt="mod.decomp", edgewgt=0.9, title="Calycadenia", nodesize=c(0.7,2), edgesize = c(0.5,0.5), addlegend = T)
text(1.5,28.5,"f)")

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("./Fig4_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)
abline(h=mean(lower(clust.goph.st$adj)),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=".", col=colors[lower(inout.goph.st)],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="")
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,28.5,"b)")

par(mar=c(1,1,1.5,6.1),xpd=T)
plotClusterMap(clust.goph.st, nodewgt="mod.decomp", edgewgt=0.9, title="", nodesize=c(0.7,2), edgesize = c(0.5,0.5))
text(1.2,28.5,"c)")

dev.off()

#Gophers, long timescales

tiff("./FigS2_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)
abline(h=mean(lower(clust.goph.lt$adj)),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=".", col=colors[lower(inout.goph.lt)],cex=2)
text(2,0.99,"a)")

par(xpd=T, mar=c(1,1,1.5,6.1))
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")
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,28.5,"b)")

par(mar=c(1,1,1.5,6.1))
plotClusterMap(clust.goph.lt, nodewgt="mod.decomp", edgewgt=0.9, title="", nodesize=c(0.7,2),
               edgesize = c(0.5,0.5))
text(1.2,28.5,"c)")

dev.off()

# ----------------------------------------------------------------------
# 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)

````

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

````