SuppMat.Rmd

---
title: "Supporting material: A new approach to interspecific synchrony in population ecology using tail association"
author: "Shyamolina Ghosh, Lawrence W. Sheppard, Philip C. Reid, Daniel Reuman"
fontsize: 12pt
geometry: "left=1.5cm,right=1.5cm,top=1.5cm,bottom=1.5cm"

output: 
  pdf_document:
    number_sections: yes
    keep_tex: yes
    fig_caption: yes

header-includes:
      - \usepackage{xr} \externaldocument[MT-]{MainText}
      - \input{head_SuppMat.sty}

mainfont: Times New Roman         
tables: True
link-citations: True
urlcolor : blue
indent : True
csl: TheAmericanNaturalist.csl
bibliography: REF_CIS.bib
---

\tableofcontents
\listoftables
\listoffigures

<!--Basic setup-->
```{r setup, echo=F, message=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.pos = "H")
seed<-101
library(dplyr)
library(kableExtra)
# families<-c(1,3:10,13,14,16:20)
source("mtime.R") #A function needed for caching
#options(scipen = 1, digits = 5) #This option round all numbers appeared in the inline r code upto 5th digit
```

<!--introductory fig, also used in methods-->
```{r fig_intro,echo=F,results="hide"}
set.seed(seed)
source("./pedagog_fig.R")
```

<!--read plankton data in our format with all 22 species and growing season temperature-->
```{r read_plankton_data,echo=F,results="hide"}
set.seed(seed)
source("./plankton_css_datacall.R")
d0_plankton<-as.matrix(read.csv("./Data/Plankton_UK_Sea_data/planktontimeseries201117_4.csv",header=F))
data_plankton<-plankton_css_datacall(d0_plankton)

# create res folder for aphids
if(!dir.exists("./Results/plankton_results/")){
  dir.create("./Results/plankton_results/")
}

# Now, include growing season temp. data for planktons 
temp_gr<-as.matrix(read.csv("./Data/Plankton_UK_Sea_data/northseaweather/timeseriestemp0202153to9.csv",header=F))
colnames(temp_gr)<-c(1945:2014)
ys<-which(colnames(temp_gr)==1958)
ye<-which(colnames(temp_gr)==2013)
temp_gr<-temp_gr[,c(ys:ye)]
#-------------------------------------------------------------------
#plot for growing season avg. yearly temperature for all location
pdf("./Results/plankton_results/temp_gr_all_loc.pdf",height=15,width = 15)
op<-par(mfrow=c(6,5),mar=c(3,3,3,3), mgp=c(1.5,0.5,0))
for(loc in 1:26){
  plot(x=c(1958:2013),y=temp_gr[loc,],col=rgb(1,0,0,0.3),pch=19,
       xlab="Year",
       ylab="Temp_growing_season_yearly_avg.",
       type="b",
       xlim=c(1958,2013),
       ylim=c(min(temp_gr),max(temp_gr)))
  grid(col="black")
  mtext(paste0("location=",loc))
}
par(op)
dev.off()

#-------------------------------------------
# include Temp. data
for(i in 1:26){
  t0<-data.frame(Year=c(1958:2013),Dat=temp_gr[i,])
  data_plankton[[i]]<-append(data_plankton[[i]],list(t0))
  ind<-length(data_plankton[[i]])
  names(data_plankton[[i]])[ind]<-"Temp_gr"
}

#-----------------------

source("vivj_matrix.R")

nsp<-length(data_plankton[[1]]) # at any location it should be 4+1(zooplankton)+1(temp)
n_each_loc<-data.frame(loc=1:26,data_pt=NA*numeric(26))
    
pdf("./Results/plankton_results/plankton_all_cop.pdf",width=2*nsp, height=2*nsp)
op<-par(mfrow=c(nsp,nsp),mar=c(4,4,4,4), mgp=c(2.5,0.5,0))
  
  for(loc in 1:length(data_plankton)){
    
    for(i in 1:nsp){
      for(j in 1:nsp){
        
        ms<-vivj_matrix(d_allsp = data_plankton,loc=loc,i=i,j=j,level=0.05,ploton=T,onbounds=T,lb=1/3,ub=2/3)
        
      }
    }
  }
par(op)
dev.off()

#============ NOTES : How to interprete ./Results/plankton_results/plankton_all_cop.pdf ? =========================

# Each of 26 pages of this pdf indicates (6 X 6 panel) copula plots for each location. significantly +ve correlated copulas are 
# represented in grey scale, whereas significantly -ve correlated copulas are represented in green. [sp_i,sp_j] copula
# should be symmetric with [sp_j,sp_i] copula unless they are significantly -ve correlated. If they are -vely correlated,
# then consider upper tringular part of 6 X 6 matrix plot arrangements; where sp_j is reversed.
```

<!--screen plankton data for the locations which had ceratium sp. and CALFIN data for more than 35 years-->
```{r screening_plankton_data,echo=F,results="hide"}
set.seed(seed)
#---------------------------
# This function to screen locations
call_summary_each_loc_plankton<-function(rawdata,ploton){
  
  summary_each_loc<-data.frame(sp=c(1:5),n=NA,nT=NA,n0=NA,avg.Abun.=NA)

  for(sp in c(1:5)){
    temp0<-rawdata[[sp]]
    temp<-na.omit(temp0)
    summary_each_loc$n[sp]<-nrow(temp) # number of data_pt. available
    summary_each_loc$nT[sp]<-length(temp$Dat)-length(unique(temp$Dat)) #no. of ties
    n0<-sum(temp$Dat==0) # number of years when a species was absent at a given location 
    summary_each_loc$n0[sp]<-n0
    #if(n0>0){
    #  summary_each_loc$nrep0[sp]<-n0-1# no. of repeat of 0
    #}else{
    #  summary_each_loc$nrep0[sp]<-NA
    #}
    summary_each_loc$avg.Abun.[sp]<-mean(temp$Dat) # avg. abundance
    
    
    if(ploton==T){
     
      plot(temp0$Year,temp0$Dat,col=rgb(0,0,1,0.3),xlab="Year",ylab="Abundance",xlim=c(1958,2013),pch=19)
      mtext(paste0("sp=",sp,
                  ", n=",summary_each_loc$n[sp],
                 ", nT=",summary_each_loc$nT[sp],
                 ", n0=",summary_each_loc$n0[sp],
                 ", avg.Ab.=",round(summary_each_loc$avg.Abun.[sp],3),sep=""))
    }
  }
 
 return(summary_each_loc)
  
}

#----------------------------
summary_all_loc_raw_plankton<-vector("list",26)
names(summary_all_loc_raw_plankton)<-paste0("loc",c(1:26),sep="")

locwise_n<-data.frame(loc=c(1:26),n=NA)


for(loc in 1:26){
  pdf(paste0("./Results/plankton_results/rawdataplot_loc_",loc,".pdf",sep=""),width=12,height=4)
  op<-par(mfrow=c(2,3),mar=c(3,3,3,3), mgp=c(1.5,0.5,0))
  rawdata<-data_plankton[[loc]]
  summary_all_loc_raw_plankton[[loc]]<-call_summary_each_loc_plankton(rawdata=rawdata,ploton=T)
  locwise_n$n[loc]<-summary_all_loc_raw_plankton[[loc]]$n[1]
  par(op)
  dev.off() 
}
#------------------------------

loc_ts_n<-which(locwise_n$n>35)  # screen location based on 1st threshold : check on n

#min_n<-min(locwise_n$n[loc_ts_n]) # minimum number of data points among all locations after 1st threshold

thrs<-0.1 # percentage threshold
splist_updated<-vector("list",26)
names(splist_updated)<-paste0("loc",1:26,sep="")

good_loc_plankton<-c()
for(i in 1:length(loc_ts_n)){
  splist<-c(1:5)
  loc<-loc_ts_n[i]
  #badsp<-which(summary_all_loc_raw_plankton[[loc]]$n0>thrs*min_n) # sp absent for more than (10% of 36 years) data were not allowed
  
  # screen species : sp absent for more than (10% of vailable years of data for a given selected location) was not allowed
  badsp<-which(summary_all_loc_raw_plankton[[loc]]$n0>thrs*summary_all_loc_raw_plankton[[loc]]$n)
  temp<-setdiff(splist,badsp)
  splist_upd<-temp
  splist_updated[[loc]]<-splist_upd
  if(length(splist_upd)>=3 & 5%in% splist_upd){  # atleast 2 ceratium species present in that loc 
                                     #including sp5(zoopl.) : 3rd threshold
    good_loc_plankton<-c(good_loc_plankton,loc)
  }
}

splist_updated<-splist_updated[good_loc_plankton]
saveRDS(splist_updated,"./Results/plankton_results/splist_updated.RDS")
saveRDS(good_loc_plankton,"./Results/plankton_results/good_loc_plankton.RDS")
```

<!--A map showing selected locations among 26 sampling sites for planktons-->
```{r map_plot_css_plankton, echo=F, results="hide"}
#source("good_splist.R")
library(maps)

longs<-read.csv("./Data/Plankton_UK_Sea_data/boxcornerlongs201117.csv",header=F)
lats<-read.csv("./Data/Plankton_UK_Sea_data/boxcornerlats201117.csv",header=F)
longs<-longs[[1]]
lats<-lats[[1]]

locno<-length(longs)
df_loc_info_org<-data.frame(loc_no=c(1:length(lats)),longs=longs,lats=lats)
df_loc_info<-df_loc_info_org[good_loc_plankton,]
xlim=c(-12,10)
ylim=c(48, 60)

pdf("./Results/plankton_results/map_plot_plankton_selected_loc.pdf",height=8,width=8)
map("world", fill=TRUE, col="grey84", bg="white",border=0, xlim=xlim, ylim=ylim)
abline(h=seq(from=ylim[1],to=ylim[2],by=2), v=seq(from=xlim[1],to=xlim[2],by=2), col="black", lty=2)

map.axes(cex.axis=2.5)
points(x = df_loc_info_org$longs+1, y = df_loc_info_org$lats+1, col="red",pch=16)
points(x = df_loc_info$longs+1, y = df_loc_info$lats+1, col="blue", pch=0, cex=1.5)
text(df_loc_info_org$loc_no, x = df_loc_info_org$longs+1, y = df_loc_info_org$lats+1, col="red",  pos=1)
dev.off()
```

<!--Computing non-parametric statistics at two tails (nbin=3) of plankton copula-->
```{r plankton_npa, echo=F, cache=T, cache.extra=list(seed,data_plankton,splist_updated,good_loc_plankton,mtime("NonParamStat.R"),mtime("vivj_matrix.R"),mtime("CopulaFunctions.R"), mtime("CopulaFunctions_flexible.R"))}

source("NonParamStat.R")

set.seed(seed) 
resloc<-"./Results/plankton_results/npa_stat_results/"

if(!dir.exists(resloc)){
  dir.create(resloc)
}

splist_updated<-readRDS("./Results/plankton_results/splist_updated.RDS")
good_loc<-readRDS("./Results/plankton_results/good_loc_plankton.RDS")

numloc<-length(good_loc)
cor_stat_plankton_ln_all<-vector("list",numloc)
names(cor_stat_plankton_ln_all)<-paste("loc",good_loc,sep="")

P_stat_plankton_ln_all<-cor_stat_plankton_ln_all
D2_stat_plankton_ln_all<-cor_stat_plankton_ln_all


for(loc in c(1:numloc)){
  nm<-names(cor_stat_plankton_ln_all)[loc]
  if (!dir.exists(paste(resloc,nm,sep=""))){
  dir.create(paste(resloc,nm,sep=""))
  }
  
  good_sp <- splist_updated[[loc]]
  good_sp<-c(good_sp,6:length(data_plankton[[good_loc[loc]]])) # to include temp. 
                                                            # or other factors as extra species
  #cat(paste("start-time: ",Sys.time(),"\n"))
  #cat("good_loc = ",good_loc[loc],"-------good_sp = ",good_sp) 
  resloc2<-paste(resloc,names(cor_stat_plankton_ln_all)[loc],"/",sep="")
  
  #----------------------- cor npa stat -----------------------------------
  cor_stat_plankton<-multcall(d_allsp=data_plankton,
                          loc=good_loc[loc],
                          pfname=paste(resloc2,"plankton",sep=''),
                          good_sp=good_sp,
                          nbin=3,npa_stats = "cor")
  cor_stat_plankton_ln_all[[loc]]<-cor_stat_plankton
  
  #----------------------- P npa stat -----------------------------------
  P_stat_plankton<-multcall(d_allsp=data_plankton,
                          loc=good_loc[loc],
                          pfname=paste(resloc2,"plankton",sep=''),
                          good_sp=good_sp,
                          nbin=3,npa_stats = "P")
  P_stat_plankton_ln_all[[loc]]<-P_stat_plankton
  
  #----------------------- D2 npa stat -----------------------------------
  D2_stat_plankton<-multcall(d_allsp=data_plankton,
                          loc=good_loc[loc],
                          pfname=paste(resloc2,"plankton",sep=''),
                          good_sp=good_sp,
                          nbin=3,npa_stats = "D2")
  
  D2_stat_plankton_ln_all[[loc]]<-D2_stat_plankton
  
  #cat(paste("stop-time: ",Sys.time(),"\n"))
}

saveRDS(cor_stat_plankton_ln_all,paste(resloc,file="cor_stat_plankton_ln_all.RDS",sep=''))
saveRDS(P_stat_plankton_ln_all,paste(resloc,file="P_stat_plankton_ln_all.RDS",sep=''))
saveRDS(D2_stat_plankton_ln_all,paste(resloc,file="D2_stat_plankton_ln_all.RDS",sep=''))
```

<!--Genarating plots with results from non-parametric statistics computed on plankton copula-->
```{r plot_res_plankton_npa, echo=F, results="hide",message=F,cache=T,warning=F, cache.extra=list(seed,cor_stat_plankton_ln_all,P_stat_plankton_ln_all,D2_stat_plankton_ln_all,mtime("NonParamStat_matrixplot.R"),mtime("mycorrplot.R"))}
set.seed(seed)
source("NonParamStat_matrixplot.R")
resloc<-"./Results/plankton_results/npa_stat_results/"

#----------------------- cor npa stat ---------------------------------------------------------
cor_npa_diff_plankton_ln_all<-NonParamStat_matrixplot(data_ln_all = cor_stat_plankton_ln_all,
                        resloc=resloc,
                        nvar=2,nvar_names=c("C. fin.","Temp."),tagon=T,npa_stats = "cor",wd=22,ht=16)
saveRDS(cor_npa_diff_plankton_ln_all,paste(resloc,file="cor_npa_diff_plankton_ln_all.RDS",sep=''))

#----------------------- P npa stat ---------------------------------------------------------
P_npa_diff_plankton_ln_all<-NonParamStat_matrixplot(data_ln_all = P_stat_plankton_ln_all,
                        resloc=resloc,
                        nvar=2,nvar_names=c("C. fin.","Temp."),tagon=T,npa_stats = "P",wd=22,ht=16)
saveRDS(P_npa_diff_plankton_ln_all,paste(resloc,file="P_npa_diff_plankton_ln_all.RDS",sep=''))

#----------------------- D2 npa stat ---------------------------------------------------------
D2_npa_diff_plankton_ln_all<-NonParamStat_matrixplot(data_ln_all = D2_stat_plankton_ln_all,
                        resloc=resloc,
                        nvar=2,nvar_names=c("C. fin.","Temp."),tagon=T,npa_stats = "D2",wd=22,ht=16)
saveRDS(D2_npa_diff_plankton_ln_all,paste(resloc,file="D2_npa_diff_plankton_ln_all.RDS",sep=''))
```

<!--plots with plankton results: 1) summary of tail dep. results on map answering Q1, 2) scatter plot answering Q2-->
```{r some_plots_plankton,echo=F,results="hide"}

#-------------- First generate a map for planktons with Corl-Coru stats summary for sp-sp interaction matrix ------------

summary_LU_plankton_sp_only<-cor_npa_diff_plankton_ln_all$summary_LU_CorlmCoru

# on map
pdf("./Results/plankton_results/npa_stat_results/Corstat_LmU_values_on_map_sp_only.pdf",height=7,width=8)
op<-par(mar=c(1,8,0,0), mgp=c(1.5,1,0))
#op<-par(mfrow=c(1,2),mar=c(1,3,1,3), mgp=c(1.5,1,0))

library(maps)
longs<-read.csv("./Data/Plankton_UK_Sea_data/boxcornerlongs201117.csv",header=F)
lats<-read.csv("./Data/Plankton_UK_Sea_data/boxcornerlats201117.csv",header=F)
longs<-longs[[1]]
lats<-lats[[1]]

locno<-length(longs)
df_loc_info_org<-data.frame(loc_no=c(1:length(lats)),longs=longs,lats=lats)
df_loc_info<-df_loc_info_org[good_loc_plankton,]
xlim=c(-12,10)
ylim=c(48, 60)

map("world", fill=T, col="grey84", bg="white",border=0, xlim=xlim, ylim=ylim)
abline(h=seq(from=ylim[1],to=ylim[2],by=2), v=seq(from=xlim[1],to=xlim[2],by=2), col="black", lty=2)

#axis(side = 1, at = seq(from=xlim[1],to=xlim[2],by=2), cex.axis=2)
map.axes(cex.axis=2.5)

frcval_L<-summary_LU_plankton_sp_only[1,]/(summary_LU_plankton_sp_only[1,]+abs(summary_LU_plankton_sp_only[2,]))
frcval_U<-abs(summary_LU_plankton_sp_only[2,])/(summary_LU_plankton_sp_only[1,]+abs(summary_LU_plankton_sp_only[2,]))

moreL<-which(frcval_L>frcval_U)
moreU<-which(frcval_L<frcval_U)
id_indep<-which(!is.finite(frcval_L))

frcval_L[is.na(frcval_L)]<-0
frcval_U[is.na(frcval_U)]<-0

points(x = df_loc_info$longs+1, y = df_loc_info$lats+1-0.4, col=rgb(1,0,0,frcval_L),pch=25,bg=rgb(1,0,0,frcval_L),cex=1.5)
points(x = df_loc_info$longs+1, y = df_loc_info$lats+1+0.4, col=rgb(0,0,1,frcval_U),pch=24,bg=rgb(0,0,1,frcval_U),cex=1.5)

text(df_loc_info$loc_no[moreU], x = df_loc_info$longs[moreU]+1, y = df_loc_info$lats[moreU]+1, 
     col="blue",cex=1.5)

text(df_loc_info$loc_no[moreL], x = df_loc_info$longs[moreL]+1, y = df_loc_info$lats[moreL]+1, 
     col="red",cex=1.5)

text(df_loc_info$loc_no[id_indep], x = df_loc_info$longs[id_indep]+1, y = df_loc_info$lats[id_indep]+1, 
     col="black",cex=1.5)

legend("bottomright", legend=c(expression(F[paste(C,",",L)]),0.2,0.4,0.6,0.8,1),
       col=c(rgb(0,0,0,0),rgb(1,0,0,0.2),rgb(1,0,0,0.4),rgb(1,0,0,0.6),rgb(1,0,0,0.8),rgb(1,0,0,1)),
       pch=25,
       pt.bg=c(rgb(0,0,0,0),rgb(1,0,0,0.2),rgb(1,0,0,0.4),rgb(1,0,0,0.6),rgb(1,0,0,0.8),rgb(1,0,0,1)),
       pt.cex=1.5,horiz=F,bty="o",bg="white",cex=2,y.intersp = 0.75,x.intersp = 0.6)

legend("topleft", legend=c(expression("|"*F[paste(C,",",R)]*"|"),0.2,0.4,0.6,0.8,1),
       col=c(rgb(0,0,0,0),rgb(0,0,1,0.2),rgb(0,0,1,0.4),rgb(0,0,1,0.6),rgb(0,0,1,0.8),rgb(0,0,1,1)),
       pch=24,
       pt.bg=c(rgb(0,0,0,0),rgb(0,0,1,0.2),rgb(0,0,1,0.4),rgb(0,0,1,0.6),rgb(0,0,1,0.8),rgb(0,0,1,1)),
       pt.cex=1.5,horiz=F,bty="o",bg="white",cex=2,y.intersp = 0.75,x.intersp = 0.6)

#mtext("Tail-dep for Ceratium sp-sp matrix",cex=2,line=1.05)

par(op)
dev.off()

#-------------- Now generate a scatter plot for planktons with Corl-Coru stats summary from 
#                total values of sp-sp CorlmCoru matrix and total values from sp-nvar, (nvar=CALFIN,Temp.) values------------
nvar<-2 # CALFIN, Temp.
xx<-c() # to store values of total taildep. 20 sp-nvar corl-coru contribution from each location
yy<-c() # to store values of total taildep. 20X20 sp-sp corl-coru contribution from each location

x<-cor_npa_diff_plankton_ln_all$CorlmCoru_all_ln_list
x_detailed<-cor_stat_plankton_ln_all
for(loc in c(1:length(x))){
  tempo<-x[[loc]]
  dl2<-ncol(tempo)-nvar
  
  tempo_sp_poscor<-tempo # this considers sp-sp interaction matrix with only positive correlated cells
  ind_N<-x_detailed[[loc]]$posnN
  tempo_sp_poscor[ind_N]<-NA
  tempo_spmat<-tempo_sp_poscor[1:dl2,1:dl2]
  
  dl<-ncol(tempo)-nvar+1
  tempo_spnvar<-tempo[c(1:dl2),c(dl:ncol(tempo))]
  
  xx<-c(xx,sum(tempo_spnvar,na.rm=T))
  yy<-c(yy,sum(tempo_spmat,na.rm=T))
}
xxyy<-as.data.frame(cbind(xx,yy))
id_indep<-which(rowSums(xxyy)==0) # locations which were independent
xxyy<-xxyy[-id_indep,]
pdf("./Results/plankton_results/npa_stat_results/Corstat_scatter_LmU_values.pdf",height=8,width=8)
op<-par(mar=c(7,8, 4.1, 2), mgp=c(4, 1, 0))
plot(xxyy$xx,xxyy$yy,
     xlab=expression(paste("Total community-driver tail association, ",A[D]^n)),
     ylab=expression(paste("Total community tail association, ",A[C]^n)),
     col=rgb(0,0,0,0.2),pch=19,cex=2,cex.lab=2,cex.axis=2)
text(xxyy$xx,xxyy$yy,good_loc_plankton[-id_indep],pos=1,col="black",cex=1.2)
lines(c(-100,100),c(0,0),type='l',lty='dashed',col="grey")
lines(c(0,0),c(-100,100),type='l',lty='dashed',col="grey")
dat<-data.frame(x=xxyy$xx,y=xxyy$yy)
mylm<-lm(y~x,data=dat)
abline(mylm)
c<-cor.test(dat$x,dat$y,method = "pearson",alternative = "t")
mtext(paste0("Pearson correlation = ",round(unname(c$estimate),3),", p = ",round(c$p.value,4),sep=""),cex=2,line=1)
dev.off()
```


<!--read aphid data in our format-->
```{r read_aphid_data, echo=F}
source("aphid_css_datacall.R")
firstflight<-as.matrix(read.csv("./Data/Aphid_data/APHID_DATA_RAW/aphidtsfirstflight141117.csv",header=F))
data_ff<-aphid_css_datacall(d0=firstflight)

# create res folder for aphids
if(!dir.exists("./Results/aphid_results/")){
  dir.create("./Results/aphid_results/")
}
```

<!--read avg. winter temperature DEC-MAR for each site and make some basic plots with that-->
```{r winter_avg_temp,echo=F,results="hide"}
# Time series analysis aphids
temp_aphid_raw<-as.matrix(read.csv("./Data/Aphid_data/APHID_DATA_RAW/aphidmonthlytempdata.csv",header = F))

temp_aphid<-vector("list",dim(temp_aphid_raw)[1])
names(temp_aphid)<-paste("loc",c(1:dim(temp_aphid_raw)[1]),sep="")

for(il in 1:length(temp_aphid)){
  
  s<-matrix(temp_aphid_raw[il,],36,12,byrow = T)
  rownames(s)<-c(1975:2010)
  colnames(s)<-c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
  
  temp_aphid[[il]]<-s
}

winter_avg_temp<-matrix(NA,11,35)
rownames(winter_avg_temp)<-paste("loc",c(1:dim(temp_aphid_raw)[1]),sep="")
colnames(winter_avg_temp)<-c(1976:2010)

for(il in 1:length(temp_aphid)){
  temploc<-temp_aphid[[il]]
  for(i in 2:dim(temploc)[1]){
    t<-c(temploc[i-1,12],temploc[i,1:3]) # previous yr dec to present year march averaging temp.
    winter_avg_temp[il,i-1]<-mean(t,na.rm = T)
  }
}
#------------------------------------------------------------------
pdf("./Results/aphid_results/all_loc_winter_avg_temp_timeseries_multipanel_plot.pdf",height=20,width=20)
op<-par(mfrow=c(4,3),mar=c(3,3,3,3), mgp=c(1.5,0.5,0))
for(loc in 1:nrow(winter_avg_temp)){
  y<-winter_avg_temp[loc,]
  x<-c(1976:2010)
  plot(x,y,type="b",xlab="Year",ylab="Avg. winter temp.",ylim=c(0,8),
       panel.first = grid(lwd=2,col="red"))
  mtext(paste0("loc =",loc),side = 3, line=0.15, adj=0.5, col="red")
}
par(op)
dev.off()
#----------------------
pdf("./Results/aphid_results/all_loc_winter_avg_temp_timeseries_in_single_plot.pdf",height=10,width=10)
plot(NA,xlim=c(1976,2010),ylim=c(1,8),xlab="Year",ylab="winter_avg_temp",panel.first = grid(lwd=2,col="grey"))
cl<-rainbow(11)
for(loc in 1:nrow(winter_avg_temp)){
  y<-winter_avg_temp[loc,]
  x<-c(1976:2010)
  lines(x,y,col=cl[loc],type="l")
  text(x,y,labels = loc,col=cl[loc])
}
par(op)
dev.off()
#-----------------------------------------------------------------------
pdf("./Results/aphid_results/all_loc_winter_avg_temp_multipanel_hist.pdf",height=20,width=20)
op<-par(mfrow=c(4,3),mar=c(3,3,3,3), mgp=c(1.5,0.5,0))
for(loc in 1:nrow(winter_avg_temp)){
  y<-winter_avg_temp[loc,]
  hist(y,breaks = 50,main = paste0("loc =",loc),cex.main=2,xlim=c(0,8))
}
par(op)
dev.off()
```

<!--formatted data with aphid first flight and avg. winter temp. for 11 sites-->
```{r aphid_data_with_winter_avg_temp,echo=F}

data_ff_with_temp<-vector("list",length(data_ff)) #initializing
names(data_ff_with_temp)<-names(data_ff)

for(loc in 1:length(data_ff_with_temp)){
  s<-data.frame(Year=c(1976:2010),Dat=winter_avg_temp[loc,])
  s<-list(winter.avg.temp=s)
  data_ff_with_temp[[loc]]<-append(data_ff[[loc]],s)
}

saveRDS(data_ff_with_temp,"./Results/aphid_results/data_aphid_ff_winter_avg_temp.RDS")
```

<!--showing map for 10 selected sites with atleast 30 years of aphid first flight data-->
```{r map_plot_css_ff, echo=F, results="hide"}
source("good_splist.R")
library(maps)

longs<-c(-4.567,0.57,-3.069,-3.312,-2.637,-1.682,-2.763,-0.356,-3.454,0.427,0.939)
lats<-c(55.477,52.26,56.457,55.949,52.125,55.213,53.854,51.807,50.628,51.733,51.185)

ff_selected_loc<-c()
for(loc in c(1:11)){
  temp<-good_splist(d_allsp=data_ff_with_temp,loc=loc,data_pt_thrs=30) # at least 30 datapts for each pairwise species taken
  nsp<-length(data_ff_with_temp[[loc]])
  if(length(temp)==nsp){ # considering that is true for all 20+1 species
    ff_selected_loc<-c(ff_selected_loc,loc)
  }
}

saveRDS(ff_selected_loc,"./Results/aphid_results/ff_selected_loc_data_pt_thrs_30.RDS")

locno<-length(longs)
df_loc_info_org<-data.frame(loc_no=c(1:length(lats)),longs=longs,lats=lats)
df_loc_info<-df_loc_info_org[ff_selected_loc,]
xlim=c(-10,2)
ylim=c(50, 60)

pdf("./Results/aphid_results/map_plot_css_ff_selected_loc.pdf",height=8,width=8)
map("world", fill=TRUE, col="grey84", bg="white",border=0, xlim=xlim, ylim=ylim)
#grid(col="black")
map.axes(cex.axis=2.5)
points(x = df_loc_info_org$longs, y = df_loc_info_org$lats, col="red",pch=16)
points(x = df_loc_info$longs, y = df_loc_info$lats, col="blue", pch=0, cex=1.5)
text(df_loc_info_org$loc_no, x = df_loc_info_org$longs, y = df_loc_info_org$lats, col="red",  pos=1)
dev.off()
```

<!--21 by 21 (20 aphid sp, 1 winter avg temp.) copula plot at 10 selected locations-->
```{r save_plot_css_aphid_with_winter_avg_temp, echo=F, results="hide", cache=T, cache.extra=list(seed,ff_selected_loc,data_ff_with_temp,mtime("vivj_matrix.R"))}
source("./vivj_matrix.R")
set.seed(seed)
 
nsp<-length(data_ff_with_temp[[1]]) # at any location it should be 20+1(temp) sp.
numloc<-length(ff_selected_loc)
level<-0.05 # p value threshold to check independence

pdf("./Results/aphid_results/copulaplot_ff_with_temp_all_loc.pdf",width=nsp*2, height=nsp*2)
op<-par(mfrow=c(nsp,nsp),mar=c(4,4,4,4), mgp=c(3,0.5,0))
  
for(loc in 1:numloc){

    for(i in 1:nsp){
      for(j in 1:nsp){
        
        ms<-vivj_matrix(d_allsp = data_ff_with_temp,loc=ff_selected_loc[loc],
                        i=i,j=j,level=0.05,ploton=T,onbounds=T,lb=0,ub=0.5)
      }
    }
}
par(op)
dev.off()
#============ NOTES : How to interprete ./Results/aphid_results/copulaplot_ff_with_temp_all_loc.pdf ? =========================

# Each of 11 pages of this pdf indicates (21 X 21 panel) copula plots for each location. significantly +ve correlated copulas are 
# represented in grey scale, whereas significantly -ve correlated copulas are represented in green. [sp_i,sp_j] copula
# should be symmetric with [sp_j,sp_i] copula unless they are significantly -ve correlated. If they are -vely correlated,
# then consider upper tringular part of 21 X 21 matrix plot arrangements; where sp_j is reversed.
```

<!--Computing non-parametric statistics at two tails (nbin=2) of aphid copula-->
```{r ff_npa, echo=F, cache=T, cache.extra=list(seed,data_ff_with_temp,ff_selected_loc,mtime("good_splist.R"),mtime("NonParamStat.R"),mtime("vivj_matrix.R"),mtime("CopulaFunctions.R"), mtime("CopulaFunctions_flexible.R"))}

source("good_splist.R")
source("NonParamStat.R")

set.seed(seed) 

resloc<-"./Results/aphid_results/ff_npa_stat_results/"

if(!dir.exists(resloc)){
  dir.create(resloc)
}

numloc<-length(ff_selected_loc)
cor_stat_ff_ln_all<-vector("list",numloc)
names(cor_stat_ff_ln_all)<-paste("loc",ff_selected_loc,sep="")

P_stat_ff_ln_all<-cor_stat_ff_ln_all
D2_stat_ff_ln_all<-cor_stat_ff_ln_all

for(loc in c(1:numloc)){
  nm<-names(cor_stat_ff_ln_all)[loc]
  if (!dir.exists(paste(resloc,nm,sep=""))){
  dir.create(paste(resloc,nm,sep=""))
  }
  
  good_sp <- good_splist(d_allsp=data_ff_with_temp,loc=ff_selected_loc[loc],data_pt_thrs=30)
  
  #cat(paste("start-time: ",Sys.time(),"\n"))
  resloc2<-paste(resloc,names(cor_stat_ff_ln_all)[loc],"/",sep="")
  
  #------------------------- cor npa stat ------------------------------
  cor_stat_ff<-multcall(d_allsp=data_ff_with_temp,
                          loc=ff_selected_loc[loc],
                          pfname=paste(resloc2,"ff",sep=''),
                          good_sp=good_sp,
                          nbin=2,npa_stats = "cor")
  cor_stat_ff_ln_all[[loc]]<-cor_stat_ff
  
  #------------------------- P npa stat ------------------------------
  P_stat_ff<-multcall(d_allsp=data_ff_with_temp,
                          loc=ff_selected_loc[loc],
                          pfname=paste(resloc2,"ff",sep=''),
                          good_sp=good_sp,
                          nbin=2,npa_stats = "P")
  P_stat_ff_ln_all[[loc]]<-P_stat_ff
  
  #------------------------- D2 npa stat ------------------------------
  D2_stat_ff<-multcall(d_allsp=data_ff_with_temp,
                          loc=ff_selected_loc[loc],
                          pfname=paste(resloc2,"ff",sep=''),
                          good_sp=good_sp,
                          nbin=2,npa_stats = "D2")
  D2_stat_ff_ln_all[[loc]]<-D2_stat_ff
  
  #cat(paste("stop-time: ",Sys.time(),"\n"))
}

saveRDS(cor_stat_ff_ln_all,paste(resloc,file="cor_stat_ff_ln_all.RDS",sep=''))
saveRDS(P_stat_ff_ln_all,paste(resloc,file="P_stat_ff_ln_all.RDS",sep=''))
saveRDS(D2_stat_ff_ln_all,paste(resloc,file="D2_stat_ff_ln_all.RDS",sep=''))
```

<!--Genarating plots with results from non-parametric statistics computed on aphid copula-->
```{r plot_res_css_ff_npa, echo=F, results="hide",message=F,cache=T, cache.extra=list(seed,cor_stat_ff_ln_all,P_stat_ff_ln_all,D2_stat_ff_ln_all,mtime("NonParamStat_matrixplot.R"),mtime("mycorrplot.R"))}

set.seed(seed)
source("NonParamStat_matrixplot.R")
resloc<-"./Results/aphid_results/ff_npa_stat_results/"

#------------------------- cor npa stat ------------------------------
cor_npa_diff_ff_ln_all<-NonParamStat_matrixplot(data_ln_all = cor_stat_ff_ln_all,
                        resloc=resloc,
                        nvar=1,nvar_names=c("Temp."),tagon=T,npa_stats="cor",wd=24,ht=24)
saveRDS(cor_npa_diff_ff_ln_all,paste(resloc,file="cor_npa_diff_ff_ln_all.RDS",sep=''))

#------------------------- P npa stat ------------------------------
P_npa_diff_ff_ln_all<-NonParamStat_matrixplot(data_ln_all = P_stat_ff_ln_all,
                        resloc=resloc,
                        nvar=1,nvar_names=c("Temp."),tagon=T,npa_stats="P",wd=24,ht=24)
saveRDS(P_npa_diff_ff_ln_all,paste(resloc,file="P_npa_diff_ff_ln_all.RDS",sep=''))

#------------------------- D2 npa stat ------------------------------
D2_npa_diff_ff_ln_all<-NonParamStat_matrixplot(data_ln_all = D2_stat_ff_ln_all,
                        resloc=resloc,
                        nvar=1,nvar_names=c("Temp."),tagon=T,npa_stats="D2",wd=24,ht=24)
saveRDS(D2_npa_diff_ff_ln_all,paste(resloc,file="D2_npa_diff_ff_ln_all.RDS",sep=''))
```

<!--plots with aphid results: 1) summary of tail dep. results on map answering Q1, 2) scatter plot, 3) multipanel plot answering Q2-->
```{r some_plots_aphid_ff,echo=F,results="hide"}

#-------------- First generate a map for aphids ff with Corl-Coru stats summary for 20X20 sp-sp interaction matrix ------------

summary_LU_ff_sp_only<-cor_npa_diff_ff_ln_all$summary_LU_CorlmCoru

# on map
pdf("./Results/aphid_results/ff_npa_stat_results/Corstat_LmU_values_on_map_sp_only.pdf",height=9,width=9)
op<-par(mar=c(1,8,1,3), mgp=c(1.5,1,0))

library(maps)
longs<-c(-4.567,0.57,-3.069,-3.312,-2.637,-1.682,-2.763,-0.356,-3.454,0.427,0.939)
lats<-c(55.477,52.26,56.457,55.949,52.125,55.213,53.854,51.807,50.628,51.733,51.185)

locno<-length(longs)
df_loc_info<-data.frame(loc_no=c(1:length(lats)),longs=longs,lats=lats)

df_loc_info<-df_loc_info[-1,] # removing 1st location based on 30 years of data

xlim=c(-10,2)
ylim=c(50, 60)

map("world", fill=T, col="grey84", bg="white",border=0, xlim=xlim, ylim=ylim)
map.axes(cex.axis=2.5)

frcval_L<-summary_LU_ff_sp_only[1,]/(summary_LU_ff_sp_only[1,]+abs(summary_LU_ff_sp_only[2,]))
frcval_U<-abs(summary_LU_ff_sp_only[2,])/(summary_LU_ff_sp_only[1,]+abs(summary_LU_ff_sp_only[2,]))
frcval_L[is.na(frcval_L)]<-0
frcval_U[is.na(frcval_U)]<-0

points(x = df_loc_info$longs+0.4, y = df_loc_info$lats-0.1, col=rgb(1,0,0,frcval_L),pch=25,bg=rgb(1,0,0,frcval_L),cex=1.5)
points(x = df_loc_info$longs+0.4, y = df_loc_info$lats+0.1, col=rgb(0,0,1,frcval_U),pch=24,bg=rgb(0,0,1,frcval_U),cex=1.5)

moreL<-which(frcval_L>frcval_U)
moreU<-which(frcval_L<frcval_U)

text(df_loc_info$loc_no[moreU], x = df_loc_info$longs[moreU], y = df_loc_info$lats[moreU], 
     col="blue",cex=1.5)

text(df_loc_info$loc_no[moreL], x = df_loc_info$longs[moreL], y = df_loc_info$lats[moreL], 
     col="red",cex=1.5)


legend("topleft", legend=c(expression(F[paste(C,",",L)]),0.2,0.4,0.6,0.8,1),cex=2,
       col=c(rgb(0,0,0,0),rgb(1,0,0,0.2),rgb(1,0,0,0.4),rgb(1,0,0,0.6),rgb(1,0,0,0.8),rgb(1,0,0,1)),
       pch=25,pt.bg=c(rgb(0,0,0,0),rgb(1,0,0,0.2),rgb(1,0,0,0.4),rgb(1,0,0,0.6),rgb(1,0,0,0.8),rgb(1,0,0,1)),
       pt.cex=2,horiz=F,bty="o") 
legend("topright", legend=c(expression("|"*F[paste(C,",",R)]*"|"),0.2,0.4,0.6,0.8,1), cex=2,
       col=c(rgb(0,0,0,0),rgb(0,0,1,0.2),rgb(0,0,1,0.4),rgb(0,0,1,0.6),rgb(0,0,1,0.8),rgb(0,0,1,1)),
       pch=24,pt.bg=c(rgb(0,0,0,0),rgb(0,0,1,0.2),rgb(0,0,1,0.4),rgb(0,0,1,0.6),rgb(0,0,1,0.8),rgb(0,0,1,1)),
       pt.cex=2,horiz=F,bty="o")
#mtext("Summary of tail-dep. for 20X20 sp-sp interaction matrix",cex=1.5)

par(op)
dev.off()

#-------------- Now generate a scatter plot for aphids ff with Corl-Coru stats summary from 
#                total values of 20X20 sp-sp CorlmCoru matrix and total values from 20 sp-temp values------------

xx<-c() # to store values of total taildep. 20 sp-nvar corl-coru contribution from each location
yy<-c() # to store values of total taildep. 20X20 sp-sp corl-coru contribution from each location

x<-cor_npa_diff_ff_ln_all$CorlmCoru_all_ln_list
for(loc in c(1:length(x))){
  temp<-x[[loc]]
  xx<-c(xx,sum(temp[1:20,21],na.rm=T)) # for aphids we know all aphid sp are +vely correlated always with each other
  yy<-c(yy,sum(temp[1:20,1:20],na.rm=T))
}

pdf("./Results/aphid_results/ff_npa_stat_results/Corstat_scatter_LmU_values.pdf",height=8,width=8)
op<-par(mar=c(7,8, 4.1, 2), mgp=c(4, 1, 0))
plot(xx,yy,
     xlab=expression(paste("Total community-driver tail association, ",A[D]^n)),
     ylab=expression(paste("Total community tail association, ",A[C]^n)),
     col=rgb(0,0,0,0.2),pch=19,cex=2,cex.lab=2,cex.axis=2,ylim=c(min(yy),max(yy)+1))
text(xx,yy,ff_selected_loc,pos=1,col="black",cex=1.2)
lines(c(-100,100),c(0,0),type='l',lty='dashed',col="grey")
lines(c(0,0),c(-100,100),type='l',lty='dashed',col="grey")
dat<-data.frame(x=xx,y=yy)
mylm<-lm(y~x,data=dat)
abline(mylm)
c<-cor.test(dat$x,dat$y,method = "pearson",alternative = "t")
mtext(paste0("Pearson correlation = ",round(unname(c$estimate),3),", p = ",round(c$p.value,4),sep=""),cex=2,line=1)
dev.off()

#----- Now generate another version of above approach to get multipanel plots for each location --------------------
LTmUTdep_matlist_ff<-cor_npa_diff_ff_ln_all$CorlmCoru_all_ln_list
pdf("./Results/aphid_results/ff_npa_stat_results/Corstat_LmU_avg_sp_temp_multipanel_plot.pdf",height = 33,width = 48)
op<-par(mfrow=c(4,3),mar=c(20,25,15,1.5), mgp=c(8,1,0),oma=c(3,3,0,0))

for(loc in 1:length(LTmUTdep_matlist_ff)){
  st<-LTmUTdep_matlist_ff[[loc]]
  s<-st[1:20,1:20]
  avg_s<-apply(X=s,MARGIN = 1,FUN = mean,na.rm=T)
  s_temp<-st[1:20,21]
  ylim<-range(avg_s,s_temp,na.rm = T)
  ylim[2]<-ylim[2]+0.1
  t<-cbind(avg_s,s_temp)
  t<-na.omit(t)
  ct<-cor.test(t[,1],t[,2],alternative = "greater",method = "pearson",conf.level = 0.95)
  pval<-ct$p.value
  tau<-unname(ct$estimate)
  nt<-nrow(t)
  t<-t[,1]/t[,2]
  t<-sum(t>0,na.rm=T)
  plot(c(1:20),avg_s,type="b",col="black",cex=4,pch=16,ylim = ylim, lty="longdash",
       ylab="",
       #ylab=expression(alpha[X]^n*(i)),
       cex.lab=12,cex.axis=8,xlab="",xaxt="n")
  axis(side = 1, at = seq(from=1,to=20,by=1), cex.axis=8, line=NA,padj=0.8)
  mtext((expression(i^th * species)), side=1, line=18, col="black", cex=7)
  #mtext(paste("loc=",loc,", s.sign=",t," out of",nt,", tau=",round(tau,2),", p=",round(pval,3)),cex=3.5)
  mtext(paste("P =",round(tau,2),", p=",formatC(pval,3,format="f")),cex=6,side=3, line=1.5) #pval rounded upto 3digits to show
  mtext(paste("n = ",ff_selected_loc[loc],sep=""),cex=6,side=3,line=-7,adj=0.06)
  points(c(1:20),s_temp,type="b",col="grey",cex=4,pch=17)
  lines(c(1,20),c(0,0),type='l',lty='dotted',lwd=0.5)
  #legend("topleft", legend=c("avg_sp", "sp-temp"),
  #       col=c("red", "blue"), pch=c(16,16),cex=3.5,horiz=T,bty="n")
}
mtext(expression(paste("Tail association, ", alpha[X]^n*(i))),cex=6,side=2,line=14,adj=-1.5)
#op2<-par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 1, 0, 0), new = TRUE)

par(op)
dev.off()

pdf("./Results/aphid_results/ff_npa_stat_results/Corstat_LmU_avg_sp_temp_multipanel_plot_legend.pdf", height=2, width=26)
op<-par(mar=c(0,1,0,0), mgp=c(0,0,0),oma=c(0,0,0,0))
plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n", ann=F)
legend("center", 
       c(expression(paste("X = C, species-community tail association, ", alpha[C]^n*(i))), 
         expression(paste("X = D, species-driver tail association, ", alpha[D]^n*(i)))), 
       col = c("black","grey"),bty="n", text.width = c(0.95,0.95),
       cex = 3, pch = c(16, 17), lty=c(5,1), x.intersp = c(0,0),y.intersp = c(1,1),
       xpd = TRUE, horiz = T, inset = c(0.06,0.05)) 
par(op)
dev.off()
```

\newpage
# The bound $b$ \label{boundb}
In order to mitigate multiple-testing problems, we selected, prior to any analyses, one 
value of b for each dataset and used it. This is equivalent to selecting one statistic 
(the one which corresponds to the value if b we selected) for each dataset and using it.
We do not present and did not compute results for other values of $b$, in order to
reduce the chance that some values of $b$ may produce "significant" results by chance 
alone. The values of $b$ we selected for each dataset ($b=1/2$ for the aphid data and 
$b=1/3$ for the plankton data) were based on the lengths of the datasets, and past
experience on the performance of the tail association statistics for datasets of a given
length.  

<!--plankton results : sampling sites in map-->
\begin{figure}[!h]
\begin{center}
\includegraphics[width=13 cm]{./Results/plankton_results/map_plot_plankton_selected_loc.pdf} 
\caption[Sampling sites on map for Plankton abundance data]{Plankton abundance data were obtained from 26 locations (2$^\circ$ $\times$ 2$^\circ$ grid cells) for the years 1958-2013. We analyzed data from 15 locations (blue boxes) based on criteria listed in the Data section of the main text.\label{fig_plankton_map}}
\end{center}
\end{figure}

<!--aphid results : sampling sites in map-->
\begin{figure}[!h]
\begin{center}
\includegraphics[width=13 cm]{./Results/aphid_results/map_plot_css_ff_selected_loc.pdf} 
\caption[Position of 11 suction traps on map that sample flying aphids]{Aphid first flight
data were obtained from 11 suction traps located across the UK for the years 1976-2010. 
We analyzed data from 10 of the traps (blue boxes) based on criteria listed in the Data
section of the main text.\label{fig_aphid_map}}
\end{center}
\end{figure}

<!--plankton results : corl-coru plot for all location-->
\begin{figure}[!h]
\begin{center}
\hspace{-2 cm}
\includegraphics[width=6 cm]{./Results/plankton_results/npa_stat_results/loc2/loc2_Corl-Coru.pdf} 
\includegraphics[width=6 cm]{./Results/plankton_results/npa_stat_results/loc7/loc7_Corl-Coru.pdf}
\includegraphics[width=6 cm]{./Results/plankton_results/npa_stat_results/loc10/loc10_Corl-Coru.pdf}\\
\hspace{-2 cm}
\includegraphics[width=6 cm]{./Results/plankton_results/npa_stat_results/loc11/loc11_Corl-Coru.pdf} 
\includegraphics[width=6 cm]{./Results/plankton_results/npa_stat_results/loc14/loc14_Corl-Coru.pdf}
\includegraphics[width=6 cm]{./Results/plankton_results/npa_stat_results/loc17/loc17_Corl-Coru.pdf}\\
\hspace{-2 cm}
\includegraphics[width=6 cm]{./Results/plankton_results/npa_stat_results/loc18/loc18_Corl-Coru.pdf} 
\includegraphics[width=6 cm]{./Results/plankton_results/npa_stat_results/loc20/loc20_Corl-Coru.pdf}
\includegraphics[width=6 cm]{./Results/plankton_results/npa_stat_results/loc21/loc21_Corl-Coru.pdf}\\
\hspace{-2 cm}
\includegraphics[width=6 cm]{./Results/plankton_results/npa_stat_results/loc23/loc23_Corl-Coru.pdf}
\includegraphics[width=6 cm]{./Results/plankton_results/npa_stat_results/loc24/loc24_Corl-Coru.pdf}
\includegraphics[width=6 cm]{./Results/plankton_results/npa_stat_results/loc25/loc25_Corl-Coru.pdf}
\caption[Planktons: $C^n$ matrix plot for selected locations]{As Fig. 
\ref{MT-fig_CorlmCoru_plankton_map_loc12_26}A, B, but for other plankton sampling locations.
Location 3 is not shown because the hypothesis that \emph{Ceratium} times series were independent
could not be rejected for that location (Methods). The two green squares for location 18
were the only two \emph{Ceratium} time series for which the association was significantly negative.
This association was not considered for our analyses, which focussed on synchronous dynamics
(see Introduction).\label{fig_CorlmCoru_plankton_all_loc}}
\end{center}
\end{figure}


<!--aphid results : npa : corl - coru all location plot-->
\begin{figure}[!h]
\begin{center}
\includegraphics[width=6 cm]{./Results/aphid_results/ff_npa_stat_results/loc3/loc3_Corl-Coru.pdf}
\includegraphics[width=6 cm]{./Results/aphid_results/ff_npa_stat_results/loc4/loc4_Corl-Coru.pdf}
\includegraphics[width=6 cm]{./Results/aphid_results/ff_npa_stat_results/loc6/loc6_Corl-Coru.pdf}\\ 
\includegraphics[width=6 cm]{./Results/aphid_results/ff_npa_stat_results/loc7/loc7_Corl-Coru.pdf} 
\includegraphics[width=6 cm]{./Results/aphid_results/ff_npa_stat_results/loc8/loc8_Corl-Coru.pdf} 
\includegraphics[width=6 cm]{./Results/aphid_results/ff_npa_stat_results/loc9/loc9_Corl-Coru.pdf}\\
\includegraphics[width=6 cm]{./Results/aphid_results/ff_npa_stat_results/loc10/loc10_Corl-Coru.pdf} 
\includegraphics[width=6 cm]{./Results/aphid_results/ff_npa_stat_results/loc11/loc11_Corl-Coru.pdf} 
\caption[Aphids: $C^n$ matrix plot for selected locations]{As Fig. 
\ref{MT-fig_CorlmCoru_ff_map_loc2_5}A, B, but for other aphid sampling locations.\label{fig_CorlmCoru_ff_all_loc}}
\end{center}
\end{figure}

```{r tab_notation, echo=F, results='asis',message=F}
library(tinytex)
library(tibble)
library(kableExtra)
library(dplyr)

dt1<-tibble("Symbols" = c("$s_{i}^n(t)$",
                          "$\\text{cor}_{l_b,u_b}$",
                          "$C^n$",
                          "$D^n$",
                          "$N_{L}^n$",
                          "$N_{R}^n$",
                          "$A_{C,L}^n$",
                          "$A_{C,R}^n$",
                          "$A_{C}^n$",
                          "$F_{C,L}^n$",
                          "$F_{C,R}^n$",
                          "$A_{D}^n$",
                          "$\\alpha_{C}^n(i)$",
                          "$\\alpha_{D}^n(i)$"
))

dt2<-tibble("Definition"=c("Data (population or phenological) for species $i$ in site $n$ at time $t$.",
                           "Partial Spearman correlation measuring association of two variables between the bounds $l_b$ and $u_b$",
                           "Community tail association matrix for site $n$",
                           "Community-driver tail association matrix for site $n$",
                           "Number of species pairs at site $\\textit{n}$ that showed stronger left- than right-tail associations. Equals the number of positive entries of $C^n$.",
                           "Number of species pairs at site $\\textit{n}$ that showed stronger right- than left-tail associations. Equals the number of negative entries of $C^n$.",
                           "Sum of all positive entries of $C^n$",
                           "Sum of all negative entries of $C^n$",
                           "Total community tail association. Equals $A_{C,L}^n + A_{C,R}^n$, i.e., the sum of all non-NA entries of $C^n$.",
                           "A normalized analogue of $A_{C,L}^n$. Equal to $A_{C,L}^n /(A_{C,L}^n + |A_{C,R}^n|)$.",
                           "A normalized analogue of $A_{C,R}^n$. Equal to $A_{C,R}^n /(A_{C,L}^n + |A_{C,R}^n|)$.",
                           "Total community-driver tail association. Equals the sum of all non-NA entries of $D^n$.",
                           "Species-community tail association for species $i$. Equals the sum of non-NA entries of $C^n(i,j)$ over all $j$ for which data were available.",
                           "Species-driver tail association. Equals the sum of non-NA entries of $D^n(i,k)$ over all $k$ for which data were available."
                           ))
  
dt<-cbind(dt1,dt2)

knitr::kable(dt, "latex", booktabs = T, linesep = "\\addlinespace",align="l", escape=F,
      caption = "Notation used in the main text. \\label{tab_notation}")%>%column_spec(2, width = "9 cm")
#%>%add_header_above(header=c("Symbols" = 1, "Definition" = 1))
```

<!--
# References
\setlength{\parindent}{-0.2in}
\setlength{\leftskip}{0.2in}
\setlength{\parskip}{1pt}
\noindent-->