This is an R Markdown document containing the results for the analyses of recruitment timeseries for the NSF CNH project on CCLME.

library(wsyn, warn.conflicts = F, quietly = T)#for wavelet analyses
library(reshape2)#to use melt function
library(ggplot2)#to plot times series
library(dplyr, warn.conflicts = F, quietly = T)
library(corrplot)#to plot correlation matrix

library(Hmisc)#for pearsons and spearmans correlation matrices

library(colorednoise)#to figure out autocorrelation values of timeseries

Data for analyses

Recruitment timeseries for multiple groundfish species (recruitment deviations), dungeness crabs and pink shrimp (log-recruitment) were obtained from stock assessments. Survival estimates were used for Chinook salmon from NCA and SFB (logit transformation), resource availability indices were used for albacore tuna from Yutaro's model (raw data from fish ticket data). All data from all sites (if available) were compiled and a box-cox transformation done on each separate timeseries using the function in the package wsyn.

sp13cd<-read.csv("D:/Rutgers_postdoc/data/recruitment timeseries/csv files for wavelet analyses/sp13cd_20191021.csv")
rownames(sp13cd)<-sp13cd[,1]
sp13cd<-sp13cd[,-1]
colnames(sp13cd)<-seq(1986, 2006, 1)
sp13cd<-as.matrix(sp13cd)

Exploration of timeseries

par(mar=c(2,2,1,1), mfrow=c(5,4), oma=c(0.5,0.5,0.1, 0.1))
for(i in 1:19){
  acf(sp13cd[i,], lag.max=8)
  mtext(dimnames(sp13cd)[[1]][i], side=3, line=-1.5)
}

Plots of transformed data in ggplot, first with all groundfish, then non-groundfish

sp13df<-as.data.frame(t(sp13cd))
sp13df$year<-seq(1986,2006,1)
sp13df2<-melt(sp13df, id.vars="year", value.name = "cleandat")
sp13df2<-sp13df2 %>% mutate(sp=substr(variable, 1,4), loc=substr(variable, 6,8))
sp13df2$sp<-as.factor(sp13df2$sp)
sp13df2$loc<-as.factor(sp13df2$loc)
grdfishsp<-c("DSOL", "HAKE", "LSTH", "POPF", "PSOL", "SABF", "SSTH", "WDRF", "YTRF")
grdfishdf<- filter(sp13df2, sp %in% grdfishsp)
grdfish<-ggplot(grdfishdf, aes(x=year, y=cleandat, color=sp))
grdfish + geom_line() + theme_bw() + labs(x="Year", y="Transformed index", title="Groundfish species")

nongrdfishdf<-filter(sp13df2, !sp %in% grdfishsp)
nongrdfishdf$loc<-factor(nongrdfishdf$loc, levels=c("WA", "OR","NCA", "SFB", "CA", "CCA"))
oth4<-ggplot(nongrdfishdf, aes(x=year, y=cleandat, color=sp))
oth4 + geom_line() + theme_bw() + facet_wrap(~loc, scales="free", nrow=2) + 
  labs(x="Year", y="Transformed index", title="Non-groundfish species")

Potentially important timescales for wavelet coherence

sp13wmf<-wmf(sp13cd, times=1986:2006, scale.min = 2, scale.max.input = 7)
plotmag(sp13wmf, colorbar=T, title="13 species 1986-2006 wmf")

sp13wpmf<-wpmf(sp13cd, times=1986:2006, scale.min = 2, scale.max.input = 7, sigmethod = "quick", nrand=1000)
plotmag(sp13wpmf, colorbar=T, title="13 species 1986-2006 wpmf", sigthresh=0.95)

Wmf plot: 6yr ts strong from 1992-2000 that seems to be decreasing, strong 2 yr ts 1993-1996

Wavelet coherence for all 19 variables (13 species) from 1986-2006

Note that wavelet coherence analyses were run and results saved in allmat list. Coherence values were obtained and p-values were adjusted with false detection rate

load(file="allmat.RData")
spcoh2.4<-allmat[[1]]
rownames(spcoh2.4)<-dimnames(sp13cd)[[1]]
colnames(spcoh2.4)<-dimnames(sp13cd)[[1]]
spcoh2.4fdr<-allmat[[2]]
colbwr<-colorRampPalette(c("blue", "white", "red"))#to specify colour palette
corrplot(spcoh2.4, method="number", type="lower", tl.pos="ld", tl.srt=20, tl.offset=0.2, 
         col=colbwr(10), is.corr=TRUE, diag=F, tl.col="black", tl.cex=0.6, number.cex=0.6, p.mat=spcoh2.4fdr, sig.level=0.1, insig="blank", mar=c(0.5,0.5,1.5,0.5))
mtext("13 sp coh, timescale 2-4 years (fdr<10%)", side=3, line=1)

Result: With fdr<10%, only LSTH and CRAB_WA, and POPF and SABF were coherent at 2-4 year timescales.

spcoh5.7<-allmat[[3]]
rownames(spcoh5.7)<-dimnames(sp13cd)[[1]]
colnames(spcoh5.7)<-dimnames(sp13cd)[[1]]
spcoh5.7fdr<-allmat[[4]]
length(which(spcoh5.7fdr<0.10))

Other correlation analyses. Pearsons correlation is shown here. Spearman's rank correlation was also done, in which only POPF and SABF were correlated at fdr<10%.

#sp13cdmat<-t(sp13cd)
#pcor.all.sp<-rcorr(sp13cdmat, type="pearson")
pcor.r.sp<-allmat[[7]]
pcor.qv.sp<-allmat[[8]]
corrplot(pcor.r.sp, method="number", type="lower", tl.pos="ld", tl.srt=20, tl.offset=0.2, 
         col=colbwr(10), is.corr=TRUE, diag=F, tl.col="black", p.mat=pcor.qv.sp, 
         tl.cex=0.6, number.cex=0.6, sig.level=0.1, insig="blank", mar=c(0.5,0.5,1.5,0.5))
mtext("13 sp coh, Pearson's corr (fdr<10%)", side=3, line=1)

Result: With fdr<10%, only POPF & SABF, and WDRF & YTRF were positively correlated (r>0.72).

Plots of timeseries with significant coherence and phase analyses

Crab_WA and longspine thornyhead timeseries plot

crab.lsth<-ggplot(subset(sp13df2, variable=="CRAB_WA"|variable=="LSTH_ALL"), aes(x=year, y=cleandat, color=variable))
crab.lsth + geom_line() + theme_bw() + labs(x="Year", y="Transformed index", title="Crab and longspine thornyhead", color=NULL) +
  theme(legend.position = c(0.70, 0.15))

Phase analyses for crab_WA and LSTH

crab.lsth.coh<-allmat[[5]]
plotmag(crab.lsth.coh)

Plot shows that coherence between these two timeseries were significant from 2-5 yr timescales, where the dotted red line is above the two black lines.

plotphase(crab.lsth.coh)

Plot shows that at the timescales that were significant from the previous plots, the two timeseries were in-phase at 2-3 yr timescales, and that crab was lagging by about pi/2 cycle (~2 yrs) at ~4yr timescales. This can also be seen in the timeseries plots of the transformed data, with crab lagging before 2000 and both are in-phase after 2000.

Sablefish and Pacific ocean perch timeseries plot

sab.pop<-ggplot(subset(sp13df2, variable=="SABF_ALL"|variable=="POPF_ALL"), aes(x=year, y=cleandat, color=variable))
sab.pop + geom_line() + theme_bw() + labs(x="Year", y="Transformed index", title="Sablefish and Pacific ocean perch", color=NULL) +
  theme(legend.position = c(0.70, 0.15))

Phase analyses for Sablefish and POP

sab.pop.coh<-allmat[[6]]
plotmag(sab.pop.coh)

Plot shows that coherence between these two timeseries were significant from 2-4 and 5-7 yr timescales, where the dotted red line is above the two black lines.

plotphase(sab.pop.coh)

Plot shows that at the timescales that were significant from the previous plots, the two timeseries were in-phase at ~2 yr timescales and at 5-7 yr timescales. This can also be seen in the timeseries plots of the transformed data.

Null models with appropriate autocorrelation

Null models were run using random red noise timeseries with autocorrelation = 0.4. Mean autocorrelation of actual timeseries was 0.15, with 80% of the 19 timeseries having autocorrelation <0.4, so to be conservative I chose an autocorrelation value of 0.4.

Coherence analyses for 2-4 year timescales were run on the null models 500 times, using the same 10,000 surrogates for each pairwise coherence test. Plot below shows a histogram of the number of coherent pairs for each null model. Red line is the number of coherent pairs obtained from the real data (2 pairs), p-value =3+1/500+1=0.008.

redqv24<-read.csv("nullredfdr10.csv")
red04p1<-ggplot(redqv24, aes(nfdr10))
red04p1 + geom_histogram(binwidth=1) + theme_classic() + labs(title="Null model for 2-4yr timescales", x="Number of coherent pairs", y="Count") +
  theme(text=element_text(size=14)) + scale_y_continuous(expand=c(0,0), limits=c(0,500)) +
  geom_vline(xintercept=2, col="red")

Conclusion from appropriately autocorrelated null models is that the 2 coherent pairs we see from the data are unlikely to occur from random timeseries, hence we have greater confidence that these 2 pairs of coherent relationships are not by chance.

Recruitment

Introduction