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


# MARSS models  {#chap-mss}
\chaptermark{Multivariate state-space models}

This lab will show you how to fit multivariate state-space (MARSS) models using the **MARSS** package. This class of time-series model is also called vector autoregressive state-space (VARSS) models.  This chapter works through an example which uses model selection to test different population structures in west coast harbor seals.  See \citet{Holmesetal2014} for a fuller version of this example. 

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

### Data and packages  {-}

All the data used in the chapter are in the **MARSS** package. For most examples, we will use the ``MARSS()`` function to fit models via maximum-likelihood.  We also show how to fit a Bayesian model using JAGS and Stan.  For these sectiosn you will need the **R2jags**, **coda** and **rstan**  packages.  To run the JAGS code, you will also need [JAGS](http://mcmc-jags.sourceforge.net/) installed.  See Chapter \@ref(chap-jags) for more details on JAGS and Chapter \@ref(chap-stan) for more details on Stan.  
```{r mss-loadpackages, results='hide', message=FALSE, warning=FALSE}
library(MARSS)
library(R2jags)
library(coda)
library(rstan)
```

## Overview {#sec-mss-overview}

As discussed in Chapter \@ref(chap-univariate-state-space), the **MARSS** package fits multivariate state-space models in this form:
\begin{equation}
\begin{gathered}
\mathbf{x}_t = \mathbf{B} \mathbf{x}_{t-1}+\mathbf{u}+\mathbf{w}_t \text{ where } \mathbf{w}_t \sim \,\text{N}(0,\mathbf{Q}) \\
\mathbf{y}_t = \mathbf{Z}\mathbf{x}_t+\mathbf{a}+\mathbf{v}_t \text{ where } \mathbf{v}_t \sim \,\text{N}(0,\mathbf{R}) \\
\mathbf{x}_0 = \boldsymbol{\mu}
\end{gathered}   
(\#eq:mss-marss)
\end{equation}
where each of the bolded terms are matrices.  Those that are bolded and small (not capitalized) have one column only, so are column matrices. 

To fit a multivariate time series model with the **MARSS** package, you need to first determine the size and structure of each of the parameter matrices: $\mathbf{B}$, $\mathbf{u}$, $\mathbf{Q}$, $\mathbf{Z}$, $\mathbf{a}$, $\mathbf{R}$ and $\boldsymbol{\mu}$.  This requires first writing down your model in matrix form.  We will illustarte this with a series of models for the temporal population dynamics of West coast harbor seals.

## West coast harbor seals counts {#sec-mss-west-coast-harbor-seals-counts}

In this example, we will use multivariate state-space models to combine surveys from four survey regions to estimate the average long-term population growth rate and the year-to-year variability in that population growth rate.  

We have five regions (or sites) where harbor seals were censused from 1978-1999 while hauled out of land\footnote{Jeffries et al. 2003.  Trends and status of harbor seals in Washington State: 1978-1999. Journal of Wildlife Management 67(1):208--219 }.  During the period of this dataset, harbor seals were recovering steadily after having been reduced to low levels by hunting prior to protection.  We will assume that the underlying population process is a stochastic exponential growth process with mean rates of increase that were not changing through 1978-1999.  

The survey methodologies were consistent throughout the 20 years of the data but we do not know what fraction of the population that each region represents nor do we know the observation-error variance for each region.  Given differences between the numbers of haul-outs in each region, the observation errors may be quite different.  The regions have had different levels of sampling; the best sampled region has only 4 years missing while the worst has over half the years missing (Figure \@ref(fig:mss-fig1)).  

```{r mss-noshowlegend, echo=FALSE, results='hide'}
d <- MARSS::harborSealWA
legendnames <- (unlist(dimnames(d)[2]))[2:ncol(d)]
for (i in 1:length(legendnames)) cat(paste(i, legendnames[i], "\n", sep = " "))
```

(ref:mss-fig1) Plot of the of the count data from the five harbor seal regions (Jeffries et al. 2003). The numbers on each line denote the different regions: 1) Strait of Juan de Fuca (SJF), 2) San Juan Islands (SJI), 2) Eastern Bays (EBays), 4) Puget Sound (PSnd), and 5) Hood Canal (HC).  Each region is an index of the total harbor seal population in each region. 

```{r mss-fig1, fig=TRUE, echo=FALSE, fig.width=5, fig.height=5, fig.cap='(ref:mss-fig1)', warning=FALSE}
temp <- as.data.frame(MARSS::harborSealWA)
pdat <- reshape2::melt(temp, id.vars = "Year", variable.name = "region")
p <- ggplot(pdat, aes(x = Year, y = value, col = region)) +
  geom_point() +
  geom_line()
p + ggtitle("Puget Sound Harbor Seal Surveys")
```
 


### Load the harbor seal data {#sec-mss-load-the-harbor-seal-data}

The harbor seal data are included in the **MARSS** package as matrix with years in column 1 and the logged counts in the other columns. Let's look at the first few years of data:
```{r mss-Cs2-showdata}
data(harborSealWA, package = "MARSS")
print(harborSealWA[1:8, ], digits = 3)
```
We are going to leave out Hood Canal (HC) since that region is somewhat isolated from the others and experiencing very different conditions due to hypoxic events and periodic intense killer whale predation.  We will set up the data as follows:

```{r mss-Cs2-readindata}
dat <- MARSS::harborSealWA
years <- dat[, "Year"]
dat <- dat[, !(colnames(dat) %in% c("Year", "HC"))]
dat <- t(dat) # transpose to have years across columns
colnames(dat) <- years
n <- nrow(dat) - 1
```



## A single well-mixed population {#sec-mss-a-single-well-mixed-population}

When we are looking at data over a large geographic region, we might make the assumption that the different census regions are measuring a single population if we think animals are moving sufficiently such that the whole area (multiple regions together) is "well-mixed".  We write a model of the total  population abundance for this case as:
\begin{equation}
n_t = \,\text{exp}(u + w_t) n_{t-1},
(\#eq:mss-expstoc)
\end{equation}
where $n_t$ is the total count in year $t$, $u$ is the mean population growth rate, and $w_t$ is the deviation from that average in year $t$. 
We then take the log of both sides and write the model in log space:
\begin{equation}
x_t = x_{t-1} + u + w_t, \textrm{ where } w_t \sim \,\text{N}(0,q)
(\#eq:mss-seg)
\end{equation}
$x_t=\log{n_t}$. When there is one effective population, there is one $x$, therefore $\mathbf{x}_t$ is a $1 \times 1$ matrix.  This is our **state** model and $x$ is called the "state".  This is just the jargon used in this type of model (state-space model) for the hidden state that you are estimating from the data.  "Hidden" means that you observe this state with error.



### The observation process {#sec-mss-the-observation-process}

We assume that all four regional time series are observations of this one population trajectory but they are scaled up or down relative to that trajectory.   In effect, we think of each regional survey as an index of the total population.  With this model, we do not think the regions represent independent subpopulations but rather independent observations of one population.
Our model for the data, $\mathbf{y}_t = \mathbf{Z} \mathbf{x}_t + \mathbf{a} + \mathbf{v}_t$, is written as:
\begin{equation}
 \left[ \begin{array}{c}
    y_{1} \\
    y_{2} \\
    y_{3} \\
    y_{4}  \end{array} \right]_t = 
    \left[ \begin{array}{c}
    1\\
    1\\
    1\\
    1\end{array} \right] x_t +  
    \left[ \begin{array}{c}
    0 \\
    a_2 \\
    a_3 \\
    a_4  \end{array} \right] + 
    \left[ \begin{array}{c}
    v_{1} \\
    v_{2} \\
    v_{3} \\
    v_{4}  \end{array} \right]_t 
(\#eq:mss-meas)
\end{equation}
Each $y_{i}$ is the observed time series of counts for a different region.  The $a$'s are the bias between the regional sample and the total population.  $\mathbf{Z}$ specifies which observation time series, $y_i$, is associated with which population trajectory, $x_j$.  In this case, $\mathbf{Z}$ is a matrix with 1 column since each region is an observation of the one population trajectory.

We allow that each region could have a unique observation variance and that the observation errors are independent between regions.  We assume that the observations errors on log(counts) are normal and thus the errors on (counts) are log-normal. The assumption of normality is not unreasonable since these regional counts are the sum of counts across multiple haul-outs.  We specify independent observation errors with different variances by specifying  that $\mathbf{v} \sim \,\text{MVN}(0,\mathbf{R})$, where
\begin{equation}
\mathbf{R} = \begin{bmatrix}
    r_1 & 0 & 0 & 0 \\
    0 & r_2 & 0 & 0\\
    0 & 0 & r_3 & 0 \\
    0 & 0 & 0 & r_4 \end{bmatrix}
(\#eq:mss-Rdiag)
\end{equation}
This is a diagonal matrix with unequal variances.  The shortcut for this structure in ```MARSS()``` is ```"diagonal and unequal"```.


### Fitting the model {#sec-mss-fitting-the-model}

We need to write the model in the form of Equation \@ref(eq:mss-marss) with each parameter written as a matrix.  The observation model (Equation \@ref(eq:mss-meas)) is already in matrix form.  Let's write the state model in matrix form too:
\begin{equation}
[x]_t = [1][x]_{t-1} + [u] + [w]_t, \textrm{ where } [w]_t \sim \,\text{N}(0,[q])
(\#eq:mss-seg-mat)
\end{equation}
It is very simple since all terms are $1 \times 1$ matrices.

To fit our model with ```MARSS()```, we set up a list which precisely describes the size and structure of each parameter matrix.  Fixed values in a matrix are designated with their numeric value and estimated values are given a character name and put in quotes.  Our model list for a single well-mixed population is:
```{r mss-fit.0.model}
mod.list.0 <- list(
  B = matrix(1),
  U = matrix("u"),
  Q = matrix("q"),
  Z = matrix(1, 4, 1),
  A = "scaling",
  R = "diagonal and unequal",
  x0 = matrix("mu"),
  tinitx = 0
)
```
and fit:
```{r mss-fit.0.fit}
fit.0 <- MARSS(dat, model = mod.list.0)
```
We already discussed that the short-cut ```"diagonal and unequal"``` means a diagonal matrix with each diagonal element having a different value.  The short-cut ```"scaling"``` means the form of $\mathbf{a}$ in Equation \@ref(eq:mss-meas) with one value set to 0 and the rest estimated.  You should run the code in the list to make sure you see that each parameter in the list has the same form as in our mathematical equation for the model.


### Model residuals {#sec-mss-model-residuals}

The model fits fine but look at the model residuals (Figure \@ref(fig:mss-model-resids-plot)).  They have problems.
```{r mss-model-resids, fig.show='hide'}
par(mfrow = c(2, 2))
resids <- MARSSresiduals(fit.0, type="tt1")
for (i in 1:4) {
  plot(resids$model.residuals[i, ], ylab = "model residuals", xlab = "")
  abline(h = 0)
  title(rownames(dat)[i])
}
```


(ref:mss-model-resids-plot) The model residuals for the first model.  SJI and EBays do not look good.

```{r mss-model-resids-plot, echo=FALSE, fig=TRUE, fig.cap='(ref:mss-model-resids-plot)'}
par(mfrow = c(2, 2))
resids <- MARSSresiduals(fit.0, type="tt1")
for (i in 1:4) {
  plot(resids$model.residuals[i, ], ylab = "model residuals", xlab = "")
  abline(h = 0)
  title(rownames(dat)[i])
}
```
 

\clearpage


## Four subpopulations with temporally uncorrelated errors {#sec-mss-segind}

The model for one well-mixed population was not very good.  Another reasonable assumption is that the different census regions are measuring four different temporally independent subpopulations.  We write a model of the log subpopulation abundances for this case as:
\begin{equation}
\begin{gathered}
\begin{bmatrix}x_1\\x_2\\x_3\\x_4\end{bmatrix}_t = 
\begin{bmatrix}
    1 & 0 & 0 & 0 \\
    0 & 1 & 0 & 0 \\
    0 & 0 & 1 & 0 \\
    0 & 0 & 0 & 1 
\end{bmatrix}
\begin{bmatrix}x_1\\x_2\\x_3\\x_4\end{bmatrix}_{t-1} +
\begin{bmatrix}u\\u\\u\\u\end{bmatrix} + 
\begin{bmatrix}w_1\\w_2\\w_3\\w_4\end{bmatrix}_t \\
\textrm{ where } \mathbf{w}_t \sim \,\text{MVN}\begin{pmatrix}0,
\begin{bmatrix}
    q & 0 & 0 & 0 \\
    0 & q & 0 & 0\\
    0 & 0 & q & 0 \\
    0 & 0 & 0 & q \end{bmatrix}\end{pmatrix}\\
\begin{bmatrix}x_1\\x_2\\x_3\\x_4\end{bmatrix}_0 = \begin{bmatrix}\mu_1\\\mu_2\\\mu_3\\\mu_4\end{bmatrix}_t 
\end{gathered}
(\#eq:mss-seg-mod1)
\end{equation}
The $\mathbf{Q}$ matrix is diagonal with one variance value.  This means that the process variance (variance in year-to-year population growth rates) is independent (good and bad years are not correlated) but the level of variability is the same across regions.  We made the $\mathbf{u}$ matrix with one $u$ value.  This means that we assume the population growth rates are the same across regions.

Notice that we set the $\mathbf{B}$ matrix equal to a diagonal matrix with 1 on the diagonal.  This is the "identity" matrix and it is like a 1 but for matrices.  We do not need $\mathbf{B}$ for our model, but ```MARSS()``` requires a value.


### The observation process {#sec-mss-the-observation-process-2}

In this model, each survey is an observation of a different $x$:
\begin{equation}
 \left[ \begin{array}{c}
    y_{1} \\
    y_{2} \\
    y_{3} \\
    y_{4}  \end{array} \right]_t = 
\begin{bmatrix}
    1 & 0 & 0 & 0 \\
    0 & 1 & 0 & 0\\
    0 & 0 & 1 & 0 \\
    0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}x_1\\x_2\\x_3\\x_4\end{bmatrix}_t +  
    \left[ \begin{array}{c}
    0 \\
    0 \\
    0 \\
    0  \end{array} \right] + 
    \left[ \begin{array}{c}
    v_{1} \\
    v_{2} \\
    v_{3} \\
    v_{4}  \end{array} \right]_t 
(\#eq:mss-meas-mod1)
\end{equation}
No $a$'s can be estimated since we do not have multiple observations of a given $x$ time series. Our $\mathbf{R}$ matrix doesn't change; the observation errors are still assumed to the independent with different variances.

Notice that our $\mathbf{Z}$ matrix changed.  $\mathbf{Z}$ is specifying which $y_i$ goes to which $x_j$.  The one we have specified means that $y_1$ is observing $x_1$, $y_2$ observes $x_2$, etc. We could have set up $\mathbf{Z}$ like so
\begin{equation}
\begin{bmatrix}
    0 & 1 & 0 & 0 \\
    1 & 0 & 0 & 0 \\
    0 & 0 & 0 & 1 \\
    0 & 0 & 1 & 0 
\end{bmatrix}
\end{equation}

This would mean that $y_1$ observes $x_2$, $y_2$ observes $x_1$, $y_3$ observes $x_4$, and $y_4$ observes $x_3$.  Which $x$ goes to which $y$ is arbitrary; we need to make sure it is one-to-one.  We will stay with $\mathbf{Z}$ as an identity matrix since $y_i$ observing $x_i$ makes it easier to remember which $x$ goes with which $y$.  


### Fitting the model {#sec-mss-fitting-mod1}

We set up the model list for `MARSS()` as:
```{r mss-fit-1-model}
mod.list.1 <- list(
  B = "identity",
  U = "equal",
  Q = "diagonal and equal",
  Z = "identity",
  A = "scaling",
  R = "diagonal and unequal",
  x0 = "unequal",
  tinitx = 0
)
```
We introduced a few more short-cuts.  ```"equal"``` means all the values in the matrix are the same.   ```"diagonal and equal"``` means that the matrix is diagonal with one value on the diagonal.  ```"unequal"``` means that all values in the matrix are different.

We can then fit our model for 4 subpopulations as:
```{r mss-fit.1.fit, results='hide'}
fit.1 <- MARSS::MARSS(dat, model = mod.list.1)
```


## Four subpopulations with temporally correlated errors {#sec-mss-four-subpopulations-with-temporally-correlated-errors}

Another reasonable assumption is that the different census regions are measuring different subpopulations but that the year-to-year population growth rates are correlated (good and bad year coincide).  The only parameter that changes is the $\mathbf{Q}$ matrix:
\begin{equation}
\mathbf{Q}=\begin{bmatrix}
    q & c & c & c \\
    c & q & c & c\\
    c & c & q & c \\
    c & c & c & q \end{bmatrix}
(\#eq:mss-qseg-mod2)
\end{equation}
This $\mathbf{Q}$ matrix structure means that the process variance (variance in year-to-year population growth rates) is the same across regions and the covariance in year-to-year population growth rates is also the same across regions.


### Fitting the model {#sec-mss-fitting-mod2}

Set up the model list for `MARSS()` as:
```{r mss-fit-2-model}
mod.list.2 <- mod.list.1
mod.list.2$Q <- "equalvarcov"
```
```"equalvarcov"``` is a shortcut for the matrix form in Equation \@ref(eq:mss-qseg-mod2).

Fit the model with:
```{r mss-fit-1-fit, results='hide'}
fit.2 <- MARSS::MARSS(dat, model = mod.list.2)
```
Results are not shown, but here are the AICc.  This last model is much better:
```{r mss-fits-aicc}
c(fit.0$AICc, fit.1$AICc, fit.2$AICc)
```


### Model residuals {#sec-mss-residuals-mod2}

Look at the model residuals (Figure \@ref(fig:mss-model-resids-2)).  They are also much better.

(ref:mss-model-resids-2) The model residuals for the model with four temporally correlated subpopulations.

```{r mss-model-resids-2, echo=FALSE, fig=TRUE, fig.cap='(ref:mss-model-resids-2)'}
par(mfrow = c(2, 2))
resids <- MARSSresiduals(fit.2, type="tt1")
for (i in 1:4) {
  plot(resids$model.residuals[i, ], ylab = "model residuals", xlab = "")
  abline(h = 0)
  title(rownames(dat)[i])
}
```
 

Figure \@ref(fig:mss-fig2-plot) shows the estimated states for each region using this code:
```{r mss-fig2, fig.show='hide'}
par(mfrow = c(2, 2))
for (i in 1:4) {
  plot(years, fit.2$states[i, ], ylab = "log subpopulation estimate", xlab = "", type = "l")
  lines(years, fit.2$states[i, ] - 1.96 * fit.2$states.se[i, ], type = "l", lwd = 1, lty = 2, col = "red")
  lines(years, fit.2$states[i, ] + 1.96 * fit.2$states.se[i, ], type = "l", lwd = 1, lty = 2, col = "red")
  title(rownames(dat)[i])
}
```


(ref:mss-fig2-plot) Plot of the estimate of log harbor seals in each region. The 95\% confidence intervals on the population estimates are the dashed lines.  These are not the confidence intervals on the observations, and the observations (the numbers) will not fall between the confidence interval lines.

```{r mss-fig2-plot, fig=TRUE, echo=FALSE, fig.width=6, fig.height=6, fig.cap='(ref:mss-fig2-plot)'}
par(mfrow = c(2, 2))
for (i in 1:4) {
  plot(years, fit.2$states[i, ], ylab = "log subpopulation estimate", xlab = "", type = "l")
  lines(years, fit.2$states[i, ] - 1.96 * fit.2$states.se[i, ], type = "l", lwd = 1, lty = 2, col = "red")
  lines(years, fit.2$states[i, ] + 1.96 * fit.2$states.se[i, ], type = "l", lwd = 1, lty = 2, col = "red")
  title(rownames(dat)[i])
}
```
 


## Using MARSS models to study spatial structure {#sec-mss-using-marss-models-to-study-spatial-structure}

For our next example, we will use MARSS models to test hypotheses about the population structure of harbor seals on the west coast.   For this example, we will evaluate the support for different population structures (numbers of subpopulations) using different $\mathbf{Z}$s to specify how survey regions map onto subpopulations.  We will assume correlated process errors with the same magnitude of process variance and covariance.  We will assume independent observations errors with equal variances at each site. We could do unequal variances but it takes a long time to fit so for this example, the observation variances are set equal.

The dataset we will use is ```harborSeal```, a 29-year dataset of abundance indices for 12 regions along the U.S. west coast between 1975-2004 (Figure \@ref(fig:mss-Cs02-fig1)). 

We start by setting up our data matrix.  We will leave off Hood Canal.
```{r mss-Cs01-setup-data}
dat <- MARSS::harborSeal
years <- dat[, "Year"]
good <- !(colnames(dat) %in% c("Year", "HoodCanal"))
sealData <- t(dat[, good])
```


(ref:mss-Cs02-fig1) Plot of log counts at each survey region in the harborSeal dataset. Each region is an index of the harbor seal abundance in that region. 

```{r mss-Cs02-fig1, fig=TRUE, echo=FALSE, fig.width=6, fig.height=6, fig.cap='(ref:mss-Cs02-fig1)', warning=FALSE}
# par(mfrow=c(4,3),mar=c(2,2,2,2))
# years = MARSS::harborSeal[,"Year"]
# for(i in 2:dim(MARSS::harborSeal)[2]) {
#     plot(years, MARSS::harborSeal[,i], xlab="", ylab="", main=colnames(MARSS::harborSeal)[i])
# }

temp <- as.data.frame(MARSS::harborSeal)
pdat <- reshape2::melt(temp, id.vars = "Year", variable.name = "region")
p <- ggplot(pdat, aes(Year, value)) +
  geom_point()
p + facet_wrap(~region)
```
 



## Hypotheses regarding spatial structure {#sec-mss-hypotheses-regarding-spatial-structure}

We will evaluate the data support for the following hypotheses about the population structure: 


*   H1: ```stock```  3 subpopulations defined by management units
*   H2: ```coast+PS```  2 subpopulations defined by coastal versus WA inland
*   H3: ```N+S```  2 subpopulations defined by north and south split in the middle of Oregon
*   H4:```NC+strait+PS+SC```  4 subpopulations defined by N coastal, S coastal, SJF+Georgia Strait, and Puget Sound
*   H5: ```panmictic```  All regions are part of the same panmictic population
*   H6: ```site```  Each of the 11 regions is a subpopulation


These hypotheses translate to these $\mathbf{Z}$ matrices (H6 not shown; it is an identity matrix):
\begin{equation*}
\begin{array}{rcccc}
&H1&H2&H4&H5\\
&\text{pnw ps ca}&\text{coast pc}&\text{nc is ps sc}&\text{pan}\\
\hline
\begin{array}{r}\text{Coastal Estuaries}\\ \text{Olympic Peninsula} \\ \text{Str. Juan de Fuca} \\ \text{San Juan Islands} \\ 
\text{Eastern Bays} \\ \text{Puget Sound} \\ \text{CA Mainland} \\ \text{CA Channel Islands} \\ \text{OR North Coast} \\ 
\text{OR South Coast} \\ \text{Georgia Strait} \end{array}&
\begin{bmatrix}
1 & 0 & 0 \\
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 1 & 0 \\
0 & 1 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1 \\
0 & 0 & 1 \\
1 & 0 & 0 \\
1 & 0 & 0 \\
0 & 1 & 0 
\end{bmatrix}&
\begin{bmatrix}
1 & 0  \\
1 & 0  \\
0 & 1  \\
0 & 1  \\
0 & 1  \\
0 & 1  \\
1 & 0  \\
1 & 0  \\
1 & 0  \\
1 & 0  \\
0 & 1  
\end{bmatrix}&
\begin{bmatrix}
1 & 0 & 0 & 0\\
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & 1 & 0\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1\\
0 & 0 & 0 & 1\\
1 & 0 & 0 & 0\\
0 & 0 & 0 & 1\\
0 & 1 & 0 & 0
\end{bmatrix}&
\begin{bmatrix}
1 \\
1 \\
1 \\
1 \\
1 \\
1 \\
1 \\
1 \\
1 \\
1 \\
1 
\end{bmatrix}
\end{array}
\end{equation*}

To tell ```MARSS()``` the form of $\mathbf{Z}$, we construct the same matrix in R.  For example, for hypotheses 1, we can write:
```{r mss-Zmodel, tidy=FALSE}
Z.model <- matrix(0, 11, 3)
Z.model[c(1, 2, 9, 10), 1] <- 1 # which elements in col 1 are 1
Z.model[c(3:6, 11), 2] <- 1 # which elements in col 2 are 1
Z.model[7:8, 3] <- 1 # which elements in col 3 are 1
```

Or we can use a short-cut by specifying $\mathbf{Z}$ as a factor that has the name of the subpopulation associated with each row in $\mathbf{y}$.  For hypothesis 1, this is
```{r mss-Zmodel1}
Z1 <- factor(c("pnw", "pnw", rep("ps", 4), "ca", "ca", "pnw", "pnw", "ps"))
```
Notice it is 11 elements in length; one element for each row of data. 


## Set up the hypotheses as different models {#sec-mss-set-up-the-hypotheses-as-different-models}

Only the $\mathbf{Z}$ matrices change for our model.  We will set up a base model list used for all models.
```{r mss-model-list, tidy=FALSE}
mod.list <- list(
  B = "identity",
  U = "unequal",
  Q = "equalvarcov",
  Z = "placeholder",
  A = "scaling",
  R = "diagonal and equal",
  x0 = "unequal",
  tinitx = 0
)
```

Then we set up the $\mathbf{Z}$ matrices using the factor short-cut.
```{r mss-set-up-Zs, tidy=FALSE}
Z.models <- list(
  H1 = factor(c("pnw", "pnw", rep("ps", 4), "ca", "ca", "pnw", "pnw", "ps")),
  H2 = factor(c(rep("coast", 2), rep("ps", 4), rep("coast", 4), "ps")),
  H3 = factor(c(rep("N", 6), "S", "S", "N", "S", "N")),
  H4 = factor(c("nc", "nc", "is", "is", "ps", "ps", "sc", "sc", "nc", "sc", "is")),
  H5 = factor(rep("pan", 11)),
  H6 = factor(1:11) # site
)
names(Z.models) <-
  c("stock", "coast+PS", "N+S", "NC+strait+PS+SC", "panmictic", "site")
```


### Fit the models  {#sec-fit-models-hyp}

We loop through the models, fit and store the results:
```{r mss-Cs05-run-models, cache=TRUE}
out.tab <- NULL
fits <- list()
for (i in 1:length(Z.models)) {
  mod.list$Z <- Z.models[[i]]
  fit <- MARSS::MARSS(sealData,
    model = mod.list,
    silent = TRUE, control = list(maxit = 1000)
  )
  out <- data.frame(
    H = names(Z.models)[i],
    logLik = fit$logLik, AICc = fit$AICc,
    num.param = fit$num.params,
    m = length(unique(Z.models[[i]])),
    num.iter = fit$numIter,
    converged = !fit$convergence
  )
  out.tab <- rbind(out.tab, out)
  fits <- c(fits, list(fit))
}
```


We will use AICc and AIC weights to summarize the data support for the different hypotheses.  First we will sort the fits based on AICc:
```{r mss-Cs06-sort-results}
min.AICc <- order(out.tab$AICc)
out.tab.1 <- out.tab[min.AICc,]
```
Next we add the $\Delta$AICc values by subtracting the lowest AICc:
```{r mss-Cs07-add-delta-aicc}
out.tab.1 <- cbind(out.tab.1,
  delta.AICc = out.tab.1$AICc - out.tab.1$AICc[1]
)
```
Relative likelihood is defined as $\,\text{exp}(-\Delta \mathrm{AICc}/2)$.
```{r mss-Cs08-add-delta-aicc}
out.tab.1 <- cbind(out.tab.1,
  rel.like = exp(-1 * out.tab.1$delta.AICc / 2)
)
```
The AIC weight for a model is its relative likelihood divided by the sum of all the relative likelihoods.  
```{r mss-Cs09-aic-weight}
out.tab.1 <- cbind(out.tab.1,
  AIC.weight = out.tab.1$rel.like / sum(out.tab.1$rel.like)
)
```

Let's look at the model weights (```out.tab.1```):
```{r mss-Cs10-print-table, echo=FALSE}
out.tab.1$delta.AICc <- round(out.tab.1$delta.AICc, digits = 2)
out.tab.1$AIC.weight <- round(out.tab.1$AIC.weight, digits = 3)
print(out.tab.1[, c("H", "delta.AICc", "AIC.weight", "converged")], row.names = FALSE)
```


## Fitting a MARSS model with JAGS {#sec-mss-multivariate-state-space-models-with-jags}

Here we show you how to fit a MARSS model for the harbor seal data using JAGS. We will focus on four time series from inland Washington and set up the data as follows:

```{r mss-set-up-seal-data-jags}
data(harborSealWA, package = "MARSS")
sites <- c("SJF", "SJI", "EBays", "PSnd")
Y <- harborSealWA[, sites]
Y <- t(Y) # time across columns
```

We will fit the model with four temporally independent subpopulations with the same population growth rate ($u$) and year-to-year variance ($q$).  This is the model in Section \@ref(sec-mss-segind).


### Writing the model in JAGS {#sec-mss-writing-the-model-in-jags}

The first step is to write this model in JAGS.  See Chapter \@ref(chap-jags) for more information on and examples of JAGS models.

```{r mss-jagsscript}
jagsscript <- cat("
model {  
   U ~ dnorm(0, 0.01);
   tauQ~dgamma(0.001,0.001);
   Q <- 1/tauQ;

   # Estimate the initial state vector of population abundances
   for(i in 1:nSites) {
      X[i,1] ~ dnorm(3,0.01); # vague normal prior 
   }

   # Autoregressive process for remaining years
   for(t in 2:nYears) {
      for(i in 1:nSites) {
         predX[i,t] <- X[i,t-1] + U;
         X[i,t] ~ dnorm(predX[i,t], tauQ);
      }
   }

   # Observation model
   # The Rs are different in each site
   for(i in 1:nSites) {
     tauR[i]~dgamma(0.001,0.001);
     R[i] <- 1/tauR[i];
   }
   for(t in 1:nYears) {
     for(i in 1:nSites) {
       Y[i,t] ~ dnorm(X[i,t],tauR[i]);
     }
   }
}  

",file="marss-jags.txt")
```

### Fit the JAGS model
{#sec-mss-fit-jags}

Then we write the data list, parameter list, and pass the model to the `jags()` function:

```{r mss-marss-jags, results='hide', message=FALSE, cache=TRUE}
jags.data <- list("Y" = Y, nSites = nrow(Y), nYears = ncol(Y)) # named list
jags.params <- c("X", "U", "Q", "R")
model.loc <- "marss-jags.txt" # name of the txt file
mod_1 <- jags(jags.data,
  parameters.to.save = jags.params,
  model.file = model.loc, n.chains = 3,
  n.burnin = 5000, n.thin = 1, n.iter = 10000, DIC = TRUE
)
```


### Plot the posteriors for the estimated states {#sec-mss-plot-the-posteriors-for-the-estimated-states}

We can plot any of the variables we chose to return to R in the `jags.params` list. Let's focus on the `X`. When we look at the dimension of the `X`, we can use the `apply()` function to calculate the means and 95 percent CIs of the estimated states.


(ref:NA) Plot of the posterior means and credible intervals for the estimated states.

```{r mss-plot-jags-states, fig.cap='(ref:NA)'}
#attach.jags attaches the jags.params to our workspace
attach.jags(mod_1)
means <- apply(X, c(2, 3), mean)
upperCI <- apply(X, c(2, 3), quantile, 0.975)
lowerCI <- apply(X, c(2, 3), quantile, 0.025)
par(mfrow = c(2, 2))
nYears <- ncol(Y)
for (i in 1:nrow(means)) {
  plot(means[i, ],
    lwd = 3, ylim = range(c(lowerCI[i, ], upperCI[i, ])),
    type = "n", main = colnames(Y)[i], ylab = "log abundance", xlab = "time step"
  )
  polygon(c(1:nYears, nYears:1, 1),
    c(upperCI[i, ], rev(lowerCI[i, ]), upperCI[i, 1]),
    col = "skyblue", lty = 0
  )
  lines(means[i, ], lwd = 3)
  title(rownames(Y)[i])
}
detach.jags()
```

```{r marss-jags-reset, include=FALSE, purl=FALSE}
file.remove("marss-jags.txt")
```

## Fitting a MARSS model with Stan {#sec-marss-fitting-with-stan}

Let's fit the same model as in Section \@ref(sec-mss-multivariate-state-space-models-with-jags) with Stan using the **rstan** package. If you have not already, you will need to install the **rstan** package.  This package depends on a number of other packages which should install automatically when you install **rstan**.

First we write the model.  We could write this to a file (recommended), but for this example, we write as a character object. Though the syntax is different from the JAGS code, it has many similarities.  Note that Stan does not allow missing values in the data, thus we need to pass in only the non-missing values along with the row and column indices of those values.  The latter is so we can match them to the appropriate state ($x$) values.

```{r marss-stan-model}
scode <- "
data {
  int<lower=0> TT; // length of ts
  int<lower=0> N; // num of ts; rows of y
  int<lower=0> n_pos; // number of non-NA values in y
  int<lower=0> col_indx_pos[n_pos]; // col index of non-NA vals
  int<lower=0> row_indx_pos[n_pos]; // row index of non-NA vals
  vector[n_pos] y;
}
parameters {
  vector[N] x0; // initial states
  real u;
  vector[N] pro_dev[TT]; // refed as pro_dev[TT,N]
  real<lower=0> sd_q;
  real<lower=0> sd_r[N]; // obs variances are different
}
transformed parameters {
  vector[N] x[TT]; // refed as x[TT,N]
  for(i in 1:N){
    x[1,i] = x0[i] + u + pro_dev[1,i];
    for(t in 2:TT) {
      x[t,i] = x[t-1,i] + u + pro_dev[t,i];
    }
  }
}
model {
  sd_q ~ cauchy(0,5);
  for(i in 1:N){
    x0[i] ~ normal(y[i],10); // assume no missing y[1]
    sd_r[i] ~ cauchy(0,5);
    for(t in 1:TT){
    pro_dev[t,i] ~ normal(0, sd_q);
    }
  }
  u ~ normal(0,2);
  for(i in 1:n_pos){
    y[i] ~ normal(x[col_indx_pos[i], row_indx_pos[i]], sd_r[row_indx_pos[i]]);
  }
}
generated quantities {
  vector[n_pos] log_lik;
  for (n in 1:n_pos) log_lik[n] = normal_lpdf(y[n] | x[col_indx_pos[n], row_indx_pos[n]], sd_r[row_indx_pos[n]]);
}
"
```

Then we call `stan()` and pass in the data, names of parameter we wish to have returned, and information on number of chains, samples (iter), and thinning. The output is verbose (hidden here) and may have some warnings.

```{r marss-stan-fit-model, message=FALSE, warning=FALSE, results='hide', cache=TRUE}
ypos <- Y[!is.na(Y)]
n_pos <- length(ypos) # number on non-NA ys
indx_pos <- which(!is.na(Y), arr.ind = TRUE) # index on the non-NAs
col_indx_pos <- as.vector(indx_pos[, "col"])
row_indx_pos <- as.vector(indx_pos[, "row"])
mod <- rstan::stan(
  model_code = scode,
  data = list(
    y = ypos, TT = ncol(Y), N = nrow(Y), n_pos = n_pos,
    col_indx_pos = col_indx_pos, row_indx_pos = row_indx_pos
  ),
  pars = c("sd_q", "x", "sd_r", "u", "x0"),
  chains = 3,
  iter = 1000,
  thin = 1
)
```

We use `extract()` to extract the parameters from the fitted model and then the means and 95\% credible intervals.
```{r marss-stan-extract, message=FALSE}
pars <- rstan::extract(mod)
means <- apply(pars$x, c(2,3), mean)
upperCI <- apply(pars$x, c(2,3), quantile, 0.975)
lowerCI <- apply(pars$x, c(2,3), quantile, 0.025)
colnames(means) <- colnames(upperCI) <- colnames(lowerCI) <- rownames(Y)
```

```{r marss-stan-plot, fig.cap="Estimated level and 95 percent credible intervals.", echo=FALSE}
temp <- as.data.frame(means)
pdat1 <- reshape2::melt(temp, variable.name = "region", value.name="mean")
temp <- as.data.frame(upperCI)
pdat2 <- reshape2::melt(temp, variable.name = "region", value.name="upperCI")
temp <- as.data.frame(lowerCI)
pdat3 <- reshape2::melt(temp, variable.name = "region", value.name="lowerCI")
pdat <- cbind(year=MARSS::harborSealWA[,"Year"],pdat1, high=pdat2[,2], low=pdat3[,2])
ggplot(pdat , aes(x = year , y = mean)) +
  facet_wrap(~region) +
  geom_line() +
  geom_ribbon(aes(x=year, ymin=low, ymax=high, group=region), alpha=0.2)+
  theme_bw()
```

\clearpage

## Problems {#sec-mss-problems}

For these questions, use the `harborSealWA` data set in **MARSS**.  The data are already logged, but you will need to remove the year column and have time going across the columns not down the rows.
```{r mss-problems-data}
require(MARSS)
data(harborSealWA, package="MARSS")
dat <- t(harborSealWA[,2:6])
```
The sites are San Juan de Fuca (SJF 3), San Juan Islands (SJI 4), Eastern Bays (EBays 5), Puget Sound (PSnd 6) and Hood Canal (HC 7).

![Regions in the harbor seal surveys](images/harborSealWA_regions.png)

1. Plot the harbor seal data.  Use whatever plotting functions you wish (e.g. `ggplot()`, `plot(); points(); lines()`, `matplot()`).

2. Fit a panmictic population model that assumes that each of the 5 sites is observing one "Inland WA" harbor seal population with trend $u$. Assume the observation errors are independent and identical. This means 1 variance on diagonal and 0s on off-diagonal.  This is the default assumption for `MARSS()`.

    a. Write the $\mathbf{Z}$ for this model.
The code to use for making a matrix in Rmarkdown is 
```
$$\begin{bmatrix}a & b & 0\\d & e & f\\0 & h & i\end{bmatrix}$$
```

    b. Write the $\mathbf{Z}$ matrix in R using `Z=matrix(...)` and using the factor short-cut for specifying $\mathbf{Z}$.  `Z=factor(c(...)`.
    
    c. Fit the model using `MARSS()`.  What is the estimated trend ($u$)? How fast was the population increasing (percent per year) based on this estimated $u$?
    
    d. Compute the confidence intervals for the parameter estimates. Compare the intervals using the Hessian approximation and using a parametric bootstrap. What differences do you see between the two approaches? Use this code:
```
library(broom)
tidy(fit)
# set nboot low so it doesn't take forever
tidy(fit, method="parametric",nboot=100)
```

    e. What does an estimate of $\mathbf{Q}=0$ mean? What would the estimated state ($x$) look like when $\mathbf{Q}=0$?

3. Using the same panmictic population model, compare 3 assumptions about the observation error structure.  
    * The observation errors are independent with different variances.
    * The observation errors are independent with the same variance.
    * The observation errors are correlated with the same variance and same correlation.
    
    a. Write the $\mathbf{R}$ variance-covariance matrices for each assumption.  
    
    b. Create each R matrix in R.  To combine, numbers and characters in a matrix use a list matrix like so:
```
A <- matrix(list(0),3,3)
A[1,1] <- "sigma2"
```

    c. Fit each model using `MARSS()` and compute the confidence intervals (CIs) for the estimated parameters.  Compare the estimated $u$ (the population long-term trend) along with their CIs.  Does the assumption about the observation errors change the $u$ estimate?
    
    d. Plot the state residuals, the ACF of the state residuals, and the histogram of the state residuals for each fit.  Are there any issues that you see?  Use this code to get your state residuals:
```
MARSSresiduals(fit)$state.residuals[1,]
```
You need the `[1,]` since the residuals are returned as a matrix.

4. Fit a model with 3 subpopulations. 1=SJF,SJI; 2=PS,EBays; 3=HC.  The $x$ part of the model is the population structure.  Assume that the observation errors are identical and independent (`R="diagonal and equal"`).  Assume that the process errors are unique and independent (`Q="diagonal and unequal"`). Assume that the $u$ are unique among the 3 subpopulation.

    a. Write the $\mathbf{x}$ equation. Make sure each matrix in the equation has the right number of rows and columns.
    
    b. Write the $\mathbf{Z}$ matrix.
    
    c. Write the $\mathbf{Z}$ in R using `Z=matrix(...)` and using the factor shortcut `Z=factor(c(...))`.
    
    d. Fit the model with `MARSS()`.
    
    e. What do the estimated $u$ and $\mathbf{Q}$ imply about the population dynamics in the 3 subpopulations?
    
5. Repeat the fit from Question 4 but assume that the 3 subpopulations covary. Use `Q="unconstrained"`.  

    a. What does the estimated  $\mathbf{Q}$ matrix tell you about how the 3 subpopulation covary?
    
    b. Compare the AICc from the model in Question 4 and the one with `Q="unconstrained"`.  Which is more supported?
    
    c. Fit the model with `Q="equalvarcov"`.  Is this more supported based on AICc?

6. Develop the following alternative models for the structure of the inland harbor seal population.  For each model assume that the observation errors are identical and independent (`R="diagonal and equal"`).  Assume that the process errors covary with equal variance and covariances (`Q="equalvarcov"`).

    * 5 subpopulations with unique $u$.
    * 5 subpopulations with shared (equal) $u$.
    * 5 subpopulations but with $u$ shared in some regions: SJF+SJI shared, PS+EBays shared, HC unique.
    * 1 panmictic population.
    * 3 subpopulations, 1=SJF,SJI, 2=PS,EBays, 3=HC, with unique $u$
    * 2 subpopulations, 1=SJF,SJI,PS,EBays, 2=HC, with unique $u$
    
    a. Fit each model using `MARSS()`.
    
    b. Prepare a table of each model with a column for the AICc values. And a column for $\Delta AICc$ (AICc minus the lowest AICc in the group).  What is the most supported model?

7. Do diagnostics on the model innovations residuals for the 3 subpopulation model from question 4. Use the following code to get your model residuals.  This will put NAs in the model residuals where there is missing data. Then do the tests on each row of `resids`.
    ```{r mss-resids, eval=FALSE}
    resids <- MARSSresiduals(fit, type="tt1")$model.residuals
    resids[is.na(dat)] <- NA
    ```

    a. Plot the model residuals.
   
    b. Plot the ACF of the model residuals. Use `acf(..., na.action=na.pass)`.
   
    c. Plot the histogram of the model residuals.  
    
    d. Fit an ARIMA() model to your model residuals using `forecast::auto.arima()`.  Are the best fit models what you want?  Note, we cannot use the Augmented Dickey-Fuller or KPSS tests when there are missing values in our residuals time series.

<!--
8.  We cannot use the Augmented Dickey-Fuller test when there are missing values in our residuals time series.  Test the model residuals for autocorrelation using the Augmented Dickey-Fuller test using lag 1 and `urca::ur.df(..., type="trend", lags=1)`. Note this is the same as `tseries::adf.test(..., k=1)`.
    
    a. Compare the test using `type="none"` and `type="trend"`.  Which one do you think is appropriate for your residuals in this case?
    
         b. Plot the ACF of the model residuals. Use `acf(..., na.action=na.pass)`.
   
     c. Plot the histogram of the model residuals.  
    
     d. Fit an ARIMA() model to your model residuals using `forecast::auto.arima()`.  Are the best fit models what you want?  Note, we cannot use the Augmented Dickey-Fuller or KPSS tests when there are missing values in our residuals time series.


-->

