```{r cylcic-sockeye-setup, include=FALSE}
knitr::opts_knit$set(unnamed.chunk.label = "cyclic-sockeye-")
knitr::opts_chunk$set(echo = TRUE, comment=NA, cache=TRUE,
                      tidy.opts=list(width.cutoff=60), tidy=TRUE,
                      fig.align='center', out.width='80%', message=FALSE,
                      warning=FALSE)
```

# Modeling cyclic sockeye {#chap-cyclic-sockeye}
\chaptermark{Cyclic sockeye}

A script with all the R code in the chapter can be downloaded  [here](./Rcode/cyclic-sockeye.R). The Rmd for this chapter can be downloaded [here](./Rmds/cyclic-sockeye.Rmd)

### Data and packages  {-}

```{r message=FALSE}
library(atsalibrary)
library(ggplot2)
library(MARSS)
```

## Overview {#sec-cyclic-sockeye-overview}

In this chapter we will use DLMs to account for cyclicity in Bristol Bay sockeye spawner escapement (returning spawners minus catch).

```{r echo=FALSE, out.width="50%"}
knitr::include_graphics("images/BB_sockeye_rivers_inset.png")
# ![](images/BB_sockeye_rivers_inset.png)
```

## Analysis goal

Our goal is to look for the underlying level (trends) in the rivers and studying whether there is evidence of synchronicity (correlated trend changes across rivers). In looking at the data, we have a couple problems. 

1. There are strong cycles in the annual data due to the life cycle of sockeye which return 4-5 years after their brood year (year they were spawned). But unfortunately the cycles do not occur regularly every 5 years. Sockeye produce cycles that are 5 years apart but the cycle might break down for a year or two and then restart after a year or two. Thus the cycle appears to shift.
2. The amplitude of the cycle changes over time.

Both of these will cause us problems when we try to estimate a stochastic level.

```{r echo=FALSE}
ggplot(sockeye, aes(x=brood_year, y=log(spawners))) + geom_line() + facet_wrap(~region, scales="free_y") + ggtitle("log spawners")
```

ACF of the spawners showing the 5-year cycle in most of the rivers.
```{r echo=FALSE}
a <- tapply(sockeye$spawners, sockeye$region, function(x){acf(x, na.action=na.pass, plot=FALSE, lag=10)$acf[,1,1]})
aa <- data.frame(acf=Reduce(c, a),
                    region=rep(names(a), each=11),
                    lag=rep(0:10, length(names(a))))
ggplot(aa, aes(x=lag, y=acf)) +
       geom_bar(stat = "identity", position = "identity") + geom_vline(xintercept=5)+
  facet_wrap(~region)+ggtitle("ACF")
```

## Modeling the cycle

As discussed in Chapter \@ref(chap-seasonal-dlm) in the covariates chapter, we can model changes in seasonality with a DLM with sine and cosine covariates. Here is a DLM model with the level and season modeled as a random walk.

$$\begin{bmatrix}x \\ \beta_1 \\ \beta_2 \end{bmatrix}_t = \begin{bmatrix}x \\ \beta_1 \\ \beta_2 \end{bmatrix}_{t-1} + \begin{bmatrix}w_1 \\ w_2 \\ w_3 \end{bmatrix}_t$$
$$y_t = \begin{bmatrix}1& \sin(2\pi t/p)&\cos(2\pi t/p)\end{bmatrix} \begin{bmatrix}x \\ \beta_1 \\ \beta_2 \end{bmatrix}_t + v_t$$

We can fit the model to the Kvichak River log spawner data and estimate the $\beta$'s and stochastic level ($x$). This is annual data so what does $p$ mean? $p$ is the time steps between peaks. For sockeye that is 5 years. So we set $p=5$. If $p$ were changing, that would cause problems, but it is not for these data (which you can confirm by looking at the ACF for different parts of the time series).

### Set up the data

```{r}
river <- "KVICHAK"
df <- subset(sockeye, region==river)
yt <- log(df$spawners)
TT <- length(yt)
p <- 5
```

### Specify the $\mathbf{Z}$ matrix

$\mathbf{Z}$ is time-varying and we set this up with an array with the 3rd dimension being time.

```{r cylcic-sockeye-Z1}
Z <- array(1, dim=c(1,3,TT))
Z[1,2,] <- sin(2*pi*(1:TT)/p)
Z[1,3,] <- cos(2*pi*(1:TT)/p)
```

### Specify the model list

Then we make our model list. We need to set $\mathbf{A}$ since `MARSS()` doesn't like the default value of `scaling` when $\mathbf{Z}$ is time-varying.
```{r cylcic-sockeye-mod-list1}
mod.list <- list(
  U = "zero",
  Q = "diagonal and unequal",
  Z = Z,
  A = "zero")
```

### Fit the model

When we fit the model we need to give `MARSS()` initial values for `x0`. It cannot come up with default ones for this model. It doesn't really matter what you pick.
```{r cyclic-sockeye-fit-1, cache=TRUE}
m <- dim(Z)[2]
fit <- MARSS(yt, model=mod.list, inits=list(x0=matrix(0,m,1)))
```

### Plot the output

The $\beta_1$ estimate is State X2 and $\beta_2$ is State X3. The estimates match what we put into the simulated data.

```{r echo=FALSE}
plot(fit, plot.type="xtT")
```

We can plot our cycle estimates and see that the peak has shifted over time. The peak has not been regularly every 5 years.

```{r echo=FALSE}
beta1s = fit$states[2,]
beta2s = fit$states[3,]
value = beta1s*sin(2*pi*(1:TT/p))+beta2s*cos(2*pi*(1:TT)/p)

plot(1:TT, value, type="l",xlab="", ylab="beta1*sin() + beta2*cos()")
abline(v=seq(0,TT,p), col="grey")
title(river)
```

Let's look at the other rivers. Write a function to do the fit.

```{r}
fitriver <- function(river, p=5){ 
df <- subset(sockeye, region==river)
yt <- log(df$spawners)
TT <- length(yt)
Z <- array(1, dim=c(1,3,TT))
Z[1,2,] <- sin(2*pi*(1:TT)/p)
Z[1,3,] <- cos(2*pi*(1:TT)/p)
mod.list <- list(
  U = "zero",
  Q = "diagonal and unequal",
  Z = Z,
  A = "zero")
fit <- MARSS(yt, model=mod.list, inits=list(x0=matrix(0,3,1)), silent=TRUE)
return(fit)
}
```

The make a list with all the fits.
```{r cyclic-sockeye-list-of-fits, cache=TRUE}
fits <- list()
for(river in names(a)){
  fits[[river]] <- fitriver(river)
}
```

Create a data frame of the amplitude of the cycle ($\sqrt{\beta_1^2+\beta_2^2}$) and the stochastic level ($x$).
```{r}
dfz <- data.frame()
for(river in names(a)){
  fit <- fits[[river]]
  tmp <- data.frame(amplitude = sqrt(fit$states[2,]^2+fit$states[3,]^2),
                    trend = fit$states[1,],
                    river=river,
                    brood_year=subset(sockeye, region==river)$brood_year)
  dfz <- rbind(dfz, tmp)
}
```

## Univariate results

Plot of the amplitude of the cycles.  All the rivers were analyzed independently. It certainly looks like there are common patterns in the amplitude of the cycles with many showing a steady decline in amplitude. Note the counts were not decreasing so this is not due to fewer spawners.
```{r}
ggplot(dfz, aes(x=brood_year, y=amplitude)) + 
  geom_line() + 
  facet_wrap(~river, scales="free_y") + 
  ggtitle("Cycle Amplitude")
```

Plot of the stochastic level. Again all the rivers were analyzed independently. It certainly looks like there are common patterns in the trends. In the next step, we can test this.
```{r}
ggplot(dfz, aes(x=brood_year, y=trend)) + 
  geom_line() + 
  facet_wrap(~river, scales="free_y") + 
  ggtitle("Stochastic Level")
```


## Multivariate DLM 1: Synchrony in levels

In the first analysis, we will look at whether the stochastic levels (underlying trends) are correlated. We will analyze all the rivers together but in the equations, we will show just two rivers to keep the equations concise.

### State model

The hidden states model will have the following components:

* Each trend $x$ will be modeled as separate but allowed to be correlated. This means either an unconstrained $\mathbf{Q}$ or an equal variance and equal covariance matrix.
* Each seasonal trend, the $\beta$'s, will also be treated as separate but independent. This means either a diagonal and equal variance $\mathbf{Q}$ or diagona and unequal variances.

The $\mathbf{x}$ equation is then:

$$\begin{bmatrix}x_a \\ x_b \\ \beta_{1a} \\ \beta_{1b} \\ \beta_{2a} \\ \beta_{2b} \end{bmatrix}_t = \begin{bmatrix}x_a \\ x_b \\ \beta_{1a} \\ \beta_{1b} \\ \beta_{2a} \\ \beta_{2b}\end{bmatrix}_{t-1} + \begin{bmatrix}w_1 \\ w_2 \\ w_3 \\ w_4 \\ w_5 \\ w_6 \end{bmatrix}_t$$

$$\begin{bmatrix}w_1 \\ w_2 \\ w_3 \\ w_4 \\ w_5 \\ w_6 \end{bmatrix}_t \sim \text{MVN}\left(0, \begin{bmatrix}
q_a & c & 0 & 0 & 0 & 0\\ 
c & q_b & 0 & 0 & 0 & 0 \\ 
0 & 0 & q_1 & 0 & 0 & 0 \\ 
0 & 0 & 0 & q_2 & 0 & 0 \\ 
0 & 0 & 0 & 0 & q_3 & 0 \\ 
0 & 0 & 0 & 0 & 0 & q_4 \end{bmatrix}\right)$$

### Observation model

The observation model will have the following components:

* Each spawner count time series will be treated as independent with independent error (equal or unequal variance).

$$\begin{bmatrix}y_a \\ y_b\end{bmatrix}_t = 
\begin{bmatrix}
1 & 0 & \sin(2\pi t/p) & 0 & \cos(2\pi t/p) & 0\\
0 & 1 & 0&\sin(2\pi t/p) & 0&\cos(2\pi t/p)
\end{bmatrix} \begin{bmatrix}
x_a \\ x_b \\ 
\beta_{1a} \\ \beta_{1b} \\
\beta_{2a} \\ \beta_{2b}
\end{bmatrix}_t + \mathbf{v}_t$$

### Fit model

Set the number of rivers.
```{r}
n <- 2
```

The following code will create the $\mathbf{Z}$ for a model with $n$ rivers. The first $\mathbf{Z}$ is shown.

```{r}
Z <- array(1, dim=c(n,n*3,TT))
Z[1:n,1:n,] <- diag(1,n)
for(t in 1:TT){
Z[,(n+1):(2*n),t] <- diag(sin(2*pi*t/p),n)
Z[,(2*n+1):(3*n),t] <- diag(cos(2*pi*t/p),n)
}
Z[,,1]
```

And this code will make the $\mathbf{Q}$ matrix:

```{r}
Q <- matrix(list(0), 3*n, 3*n)
Q[1:n,1:n] <- "c"
diag(Q) <- c(paste0("q",letters[1:n]), paste0("q",1:(2*n)))
Q
```

We will write a function to prepare the model matrices and fit. It takes the names of the rivers.
```{r}
fitriver.m <- function(river, p=5){ 
  require(tidyr)
  require(dplyr)
  require(MARSS)
df <- subset(sockeye, region %in% river)
df <- df %>% pivot_wider(id_cols=brood_year,names_from="region", values_from=spawners) %>%
  ungroup() %>% select(-brood_year)
yt <- t(log(df))
TT <- ncol(yt)
n <- nrow(yt)
Z <- array(1, dim=c(n,n*3,TT))
Z[1:n,1:n,] <- diag(1,n)
for(t in 1:TT){
Z[,(n+1):(2*n),t] <- diag(sin(2*pi*t/p),n)
Z[,(2*n+1):(3*n),t] <- diag(cos(2*pi*t/p),n)
}
Q <- matrix(list(0), 3*n, 3*n)
Q[1:n,1:n] <- paste0("c",1:(n^2))
diag(Q) <- c(paste0("q",letters[1:n]), paste0("q",1:(2*n)))
Q[lower.tri(Q)] <- t(Q)[lower.tri(Q)]
mod.list <- list(
  U = "zero",
  Q = Q,
  Z = Z,
  A = "zero")
fit <- MARSS(yt, model=mod.list, inits=list(x0=matrix(0,3*n,1)), silent=TRUE)
return(fit)
}
```

Now we can fit for two (or more) rivers. Note it didn't quite converge as some of the variances for the $\beta$'s are going to 0 (constant $\beta$ value).
```{r cyclic-sockeye-more-rivers, cache=TRUE}
river <- unique(sockeye$region)
n <- length(river)
fit <- fitriver.m(river)
```

### Look at the results

We will look at the correlation plot for the trends.
```{r cyclic-sockeye-corrplot}
require(corrplot)
Qmat <- coef(fit, type="matrix")$Q[1:n,1:n]
rownames(Qmat) <- colnames(Qmat) <- river
M <- cov2cor(Qmat)
corrplot(M, order = "hclust", addrect = 4)
```

We can compare to the locations and see that this suggests that there is small scale regional correlation in the spawner counts.

```{r echo=FALSE, out.width="50%"}
knitr::include_graphics("images/BB_sockeye_rivers_inset.png")
# ![](images/BB_sockeye_rivers_inset.png)
```
