```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
set.seed(123)
```

## Topics

Multivariate state-space models

1. state (process) model

2. observation model


## Simple model for 2+ time series | Random walk observed with error

$$
x_{i,t} = x_{i,t-1} + w_{i,t} \\
y_{i,t} = x_{i,t} + a_i + v_{i,t}
$$

with

$w_{i,t} \sim \text{N}(0, q)$

$v_{i,t} \sim \text{N}(0, r)$


## Random walk observed with error

$$
x_{1,t} = x_{1,t-1} + w_{1,t} \\
x_{2,t} = x_{2,t-1} + w_{2,t} \\
\vdots \\
x_{n,t} = x_{n,t-1} + w_{n,t}
$$

$$
y_{1,t} = x_{1,t} + a_1 + v_{1,t} \\
y_{2,t} = x_{2,t} + a_2 + v_{2,t} \\
\vdots \\
y_{n,t} = x_{n,t} + a_2 + v_{n,t} \\
$$


## Random walk observed with error | In matrix form

$$
\begin{bmatrix}
x_{1,t} \\
x_{2,t} \\
\vdots \\
x_{n,t} \\
\end{bmatrix} =
\begin{bmatrix}
x_{1,t-1} \\
x_{2,t-1} \\
\vdots \\
x_{n,t-1} \\
\end{bmatrix} + 
\begin{bmatrix}
w_{1,t} \\
w_{2,t} \\
\vdots \\
w_{n,t} \\
\end{bmatrix}
$$

$$
\begin{bmatrix}
y_{1,t} \\
y_{2,t} \\
\vdots \\
y_{n,t} \\
\end{bmatrix} =
\begin{bmatrix}
x_{1,t} \\
x_{2,t} \\
\vdots \\
x_{n,t} \\
\end{bmatrix} + 
\begin{bmatrix}
a_1 \\
a_2 \\
\vdots \\
a_n \\
\end{bmatrix} + 
\begin{bmatrix}
v_{1,t} \\
v_{2,t} \\
\vdots \\
v_{n,t} \\
\end{bmatrix}
$$


## Random walk observed with error | In matrix form

$$
\mathbf{x}_t = \mathbf{x}_{t-1} + \mathbf{w}_t \\
\mathbf{y}_t = \mathbf{x}_t + \mathbf{a} + \mathbf{v}_t
$$

with 

$\mathbf{w}_t \sim \text{MVN}(\mathbf{0}, \mathbf{Q})$

$\mathbf{v}_t \sim \text{MVN}(\mathbf{0}, \mathbf{R})$


## Forms of covariances matrices $\mathbf{Q}$

No covariance

$$
\mathbf{Q} \stackrel{?}{=}
\begin{bmatrix}
 \sigma & 0 & 0 & 0 \\
 0 & \sigma & 0 & 0 \\
 0 & 0 & \sigma & 0 \\
 0 & 0 & 0 & \sigma
\end{bmatrix}
~\text{or}~~
\mathbf{Q} \stackrel{?}{=}
\begin{bmatrix}
 \sigma_1 & 0 & 0 & 0 \\
 0 & \sigma_2 & 0 & 0 \\
 0 & 0 & \sigma_3 & 0 \\
 0 & 0 & 0 & \sigma_4
\end{bmatrix}
$$


## Forms of covariances matrices $\mathbf{Q}$

With covariance

$$
\mathbf{Q} \stackrel{?}{=}
\begin{bmatrix}
 \sigma & \gamma & \gamma & \gamma \\
 \gamma & \sigma & \gamma & \gamma \\
 \gamma & \gamma & \sigma & \gamma \\
 \gamma & \gamma & \gamma & \sigma
\end{bmatrix}
~\text{or}~~
\mathbf{Q} \stackrel{?}{=}
\begin{bmatrix}
 \sigma_1 & \gamma_{1,2} & \gamma_{1,3} & \gamma_{1,4} \\
 \gamma_{1,2} & \sigma_2 & \gamma_{2,3} & \gamma_{2,4} \\
 \gamma_{1,3} & \gamma_{2,3} & \sigma_3 & \gamma_{3,4} \\
 \gamma_{1,4} & \gamma_{2,4} & \gamma_{3,4} & \sigma_4
\end{bmatrix}
$$


## Forms of covariances matrices $\mathbf{R}$

No covariance

$$
\mathbf{R} \stackrel{?}{=}
\begin{bmatrix}
 \sigma & 0 & 0 & 0 \\
 0 & \sigma & 0 & 0 \\
 0 & 0 & \sigma & 0 \\
 0 & 0 & 0 & \sigma
\end{bmatrix}
~\text{or}~~
\mathbf{R} \stackrel{?}{=}
\begin{bmatrix}
 \sigma_1 & 0 & 0 & 0 \\
 0 & \sigma_2 & 0 & 0 \\
 0 & 0 & \sigma_3 & 0 \\
 0 & 0 & 0 & \sigma_4
\end{bmatrix}
$$


## Forms of covariances matrices $\mathbf{R}$

With covariance

$$
\mathbf{R} \stackrel{?}{=}
\begin{bmatrix}
 \sigma & \gamma & \gamma & \gamma \\
 \gamma & \sigma & \gamma & \gamma \\
 \gamma & \gamma & \sigma & \gamma \\
 \gamma & \gamma & \gamma & \sigma
\end{bmatrix}
$$




## Animal tracking

Estimating animal locations from tagging is a classic example of applying this type of multivariate state-space model

$$
LAT_t = LAT_{t-1} + w_{LAT,t} \\
LON_t = LON_{t-1} + w_{LON,t} \\
$$

$$
y_t = LAT_t + v_{y,t} \\
x_t = LON_t + v_{x,t} \\
$$


## Random walk observed with error | In matrix form

$$
\mathbf{x}_t = \mathbf{x}_{t-1} + \mathbf{w}_t \\
\mathbf{y}_t = \mathbf{x}_t + \mathbf{a} + \mathbf{v}_t \\
\Downarrow \\
\mathbf{x}_t = \mathbf{x}_{t-1} + \mathbf{w}_t \\
\mathbf{y}_t = \mathbf{Z} \mathbf{x}_t + \mathbf{a} + \mathbf{v}_t
$$

with $\mathbf{Z}$ equal to the $n \times n$ Identity matrix $\mathbf{I}_n$


## Random walk observed with error | In matrix form

$$
\mathbf{x}_t = \mathbf{x}_{t-1} + \mathbf{w}_t \\
\mathbf{y}_t = \mathbf{Z} \mathbf{x}_t + \mathbf{a} + \mathbf{v}_t
$$

Note that both $\mathbf{x}_t$ and $\mathbf{y}_t$ are $n \times 1$ column vectors




## Multiple observations of one state

What if we had multiple observations of only 1 state?

If so, $\mathbf{x}_t$ would be a $1 \times 1$ vector (ie, a scalar) and $\mathbf{y}_t$ would be an $n \times 1$ column vector

$$
[x_t] = [x_{t-1}] + [w_t]
$$

$$
\begin{bmatrix}
y_{1,t} \\
y_{2,t} \\
\vdots \\
y_{n,t} \\
\end{bmatrix} =
\begin{bmatrix}
x_{1,t} \\
x_{2,t} \\
\vdots \\
x_{n,t} \\
\end{bmatrix} + 
\begin{bmatrix}
a_1 \\
a_2 \\
\vdots \\
a_n \\
\end{bmatrix} + 
\begin{bmatrix}
v_{1,t} \\
v_{2,t} \\
\vdots \\
v_{n,t} \\
\end{bmatrix}
$$


## Multiple observations of one state

$$
[x_t] = [x_{t-1}] + [w_t]
$$

$$
\begin{bmatrix}
y_{1,t} \\
y_{2,t} \\
\vdots \\
y_{n,t} \\
\end{bmatrix} =
\begin{bmatrix}
x_{1,t} \\
x_{2,t} \\
\vdots \\
x_{n,t} \\
\end{bmatrix} + 
\begin{bmatrix}
a_1 \\
a_2 \\
\vdots \\
a_n \\
\end{bmatrix} + 
\begin{bmatrix}
v_{1,t} \\
v_{2,t} \\
\vdots \\
v_{n,t} \\
\end{bmatrix}
$$

We need a way to map the $1 \times 1$ state onto the $n \times 1$ vector of observations

$$
\mathbf{y}_t = ~?~ x_t
$$


## Multiple observations of one state

We'll use the matrix $\mathbf{Z}$ we saw earlier and treat $x_t$ like a matrix

$$
\mathbf{y}_t = \mathbf{Z} \mathbf{x}_t
$$

Recall that the inner dimensions must match when multiplying matrices, such that

$$
[n \times m] [m \times p] = [n \times p]
$$

If $\mathbf{y}_t$ is $n \times 1$ and $\mathbf{x}_t$ is $1 \times 1$ then $\mathbf{Z}$ must be $n \times 1$ 

$$
[n \times 1] = [n \times 1] [1 \times 1]
$$


## Multiple observations of one state

For example

$$
\begin{bmatrix}
y_{1,t} \\
y_{2,t} \\
y_{3,t}
\end{bmatrix} =
\begin{bmatrix}
1 \\
1 \\
1
\end{bmatrix}
\begin{bmatrix}
x_t
\end{bmatrix} + 
\begin{bmatrix}
a_1 \\
a_2 \\
a_3
\end{bmatrix} + 
\begin{bmatrix}
v_{1,t} \\
v_{2,t} \\
v_{3,t}
\end{bmatrix}
$$


## Multiple obs of multiple states

What if we had $n$ observations of $m$ different states?

We just need to follow the same logic for matrix multiplication


## Multiple obs of multiple states

For example

$$
\begin{bmatrix}
 y_1 \\
 y_2 \\
 y_3 \\
 y_4 \\
 y_5 
\end{bmatrix}_t =
\begin{bmatrix}
 1 & 0 & 0 \\
 0 & 1 & 0 \\
 0 & 1 & 0 \\
 0 & 0 & 1 \\
 0 & 0 & 1 \\
\end{bmatrix} \times
\begin{bmatrix}
 x_1 \\
 x_2 \\
 x_3 
\end{bmatrix}_t +
\begin{bmatrix}
 a_1 \\
 a_2 \\
 a_3 \\
 a_4 \\
 a_5 
\end{bmatrix} +
\begin{bmatrix}
 v_1 \\
 v_2 \\
 v_3 \\
 v_4 \\
 v_5 
\end{bmatrix}_t
$$

$y_1$ observes $x_1$  
$y_2$ & $y_3$ observe $x_2$  
$y_4$ & $y_5$ observe $x_3$  


## Identification of population structure

We can use the same approach to systematically evaluate the data support for different hypotheses about the underlying population structure

$$
\mathbf{Z} =
\begin{bmatrix}
1 \\
1 \\
\vdots \\
1
\end{bmatrix}
~~ \text{vs} ~~
\mathbf{Z} =
\begin{bmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
\vdots & \vdots & \vdots \\
0 & 0 & 1 
\end{bmatrix}
~~ \text{vs} ~~
\mathbf{Z} =
\begin{bmatrix}
1 & 0 & \dots & 0 \\
1 & 0 & \dots & 0 \\
\vdots & \vdots & \ddots & 0 \\
0 & 0 & \dots & 1 
\end{bmatrix}
$$


## Popn structure of harbor seals

Here is an example of using different forms of $\mathbf{Z}$ to evaluate the possible metapopulation structure of Harbor Seals

[Harbor seal example](https://nwfsc-timeseries.github.io/atsa-labs/sec-mss-west-coast-harbor-seals-counts.html)


## Estimating community interactions

Number today is a function of the number yesterday

```{r MAR_diag_1, dpi=300, fig.height=4, fig.width=8, out.height="100%", out.width="100%", fig.align='center', echo=FALSE, warning=FALSE}
par(mai=c(0.8,0.8,0,0), omi=rep(0,4))
## boundaries
ss <- 5
nn <- 8
rr <- ss*3
cc <- ss*nn
## mid-points
xm <- seq(5,cc-ss,rr)
ymt <- rr - ss/2
ymb <- ss/2
## arrow locs
y0 <- rr - ss
y1 <- ss
## empty plot space
plot(c(0,cc), c(0,rr), type="n", xlab="", ylab="",
     xaxt="n", yaxt="n", bty="n")
## top row
symbols(x=xm[2], y=ymt, rectangles=matrix(c(2*ss,ss),1,2),
        lty="solid",  bg="#656565", fg=NA,
        inches=FALSE, add=TRUE, lwd=3)
text("Number", x=xm[2], y=ymt, cex=1.8, col="#ffffff", pos=3)
text("yesterday", x=xm[2], y=ymt, cex=1.8, col="#ffffff", pos=1)
## arrows
arrows(x0=xm[2], y0=y0, y1=y1,
       col="#656565", lwd=3, length=0.12)
## bottom row: obs
symbols(x=xm[2], y=ymb, rectangles=matrix(c(2*ss,ss),1,2),
        lty="solid",  bg="#488fdf", fg=NA,
        inches=FALSE, add=TRUE, lwd=3)
text("Number", x=xm[2], y=ymb, cex=1.8, col="#ffffff", pos=3)
text("today", x=xm[2], y=ymb, cex=1.8, col="#ffffff", pos=1)
```


## Estimating community interactions

and the number of predators, prey & competitors

```{r MAR_diag_2, dpi=300, fig.height=4, fig.width=8, out.height="100%", out.width="100%", fig.align='center', echo=FALSE, warning=FALSE}
par(mai=c(0.8,0.8,0,0), omi=rep(0,4))
## empty plot space
plot(c(0,cc), c(0,rr), type="n", xlab="", ylab="",
     xaxt="n", yaxt="n", bty="n")
## top row
symbols(x=xm, y=rep(ymt,3),
        rectangles=matrix(c(2*ss,ss),3,2,byrow=TRUE),
        lty="solid",  bg=c("#c10101","#ff8100","#844870"),
        fg=NA,
        inches=FALSE, add=TRUE, lwd=3)
text(c("Predators","Prey","Competitors"),
     x=xm, y=rep(ymt,3), cex=1.7, col="#ffffff", pos=3)
text("yesterday", x=xm, y=rep(ymt,3), cex=1.7, col="#ffffff", pos=1)
## arrows
arrows(x0=ss*seq(2,6,2), x1=ss*seq(3,5),
       y0=y0, y1=y1,
       col=c("#c10101","#ff8100","#844870"), lwd=3, length=0.12)
## bottom row: obs
symbols(x=xm[2], y=ymb, rectangles=matrix(c(2*ss,ss),1,2),
        lty="solid",  bg="#488fdf", fg=NA,
        inches=FALSE, add=TRUE, lwd=3)
text("Number", x=xm[2], y=ymb, cex=1.8, col="#ffffff", pos=3)
text("today", x=xm[2], y=ymb, cex=1.8, col="#ffffff", pos=1)
```


## State model for species interactions

and external forces at various times

```{r MAR_diag_3, dpi=300, fig.height=4, fig.width=8, out.height="100%", out.width="100%", fig.align='center', echo=FALSE, warning=FALSE}
par(mai=c(0.8,0.8,0,0), omi=rep(0,4))
## empty plot space
plot(c(0,cc), c(0,rr), type="n", xlab="", ylab="",
     xaxt="n", yaxt="n", bty="n")
## top row
symbols(x=xm[2], y=ymt, rectangles=matrix(c(2*ss,ss),1,2),
        lty="solid",  bg="#339933", fg=NA,
        inches=FALSE, add=TRUE, lwd=3)
text("External", x=xm[2], y=ymt, cex=1.8, col="#ffffff", pos=3)
text("forces", x=xm[2], y=ymt, cex=1.8, col="#ffffff", pos=1)
## arrows
arrows(x0=xm[2], y0=y0, y1=y1,
       col="#339933", lwd=3, length=0.12)
## bottom row: obs
symbols(x=xm[2], y=ymb, rectangles=matrix(c(2*ss,ss),1,2),
        lty="solid",  bg="#488fdf", fg=NA,
        inches=FALSE, add=TRUE, lwd=3)
text("Number", x=xm[2], y=ymb, cex=1.8, col="#ffffff", pos=3)
text("today", x=xm[2], y=ymb, cex=1.8, col="#ffffff", pos=1)
```

## State model for species interactions

$$
x_{i,t} = \sum^{m}_{j = 1}{b_{i,j} x_{j,t}} + w_t
$$

where $b_{i,j}$ is the effect of the $j^{th}$ species on species $i$

$b_{i,j}$ with $i = j$ is the density-dependent effect


## State model for species interactions

We can write this model in matrix notation as

$$
\mathbf{x}_t = \mathbf{B} \mathbf{x}_{t-1} + \mathbf{w}_t
$$

with

$$
\mathbf{B} =
\begin{bmatrix}
b_{1,1} & b_{1,2} & \dots & b_{1,m} \\
b_{2,1} & b_{1,2} & \dots & b_{2,m} \\
\vdots & \vdots & \ddots & \vdots \\
b_{m,1} & b_{m,2} & \dots & b_{m,m}
\end{bmatrix}
$$


## State model for species interactions

Including the effects of exogenous drivers yields

$$
\mathbf{x}_t = \mathbf{B} \mathbf{x}_{t-1} + \mathbf{C} \mathbf{c}_{t-k} + \mathbf{w}_t
$$

Note that the lag $k \geq 0$


## State model for species interactions

Including the effects of exogenous drivers

$$
\mathbf{x}_t = \mathbf{B} \mathbf{x}_{t-1} + \mathbf{C} \mathbf{c}_{t-k} + \mathbf{w}_t
$$

The $m \times p$ matrix $\mathbf{C}$ contains the effect(s) of each covariate (cols) on each state (rows)

The $p \times 1$ column vector $\mathbf{c}_{t-k}$ contains each of the $p$ covariates at time $t - k$


## Covariate effects

The effect(s) of covariates can vary by state/species/etc

$$
\mathbf{C} =
\begin{bmatrix}
C_{1, Temp} & C_{1, DO} \\ 
C_{2, Temp} & C_{2, DO} \\ 
\vdots & \vdots \\ 
C_{m, Temp} & C_{m, DO}
\end{bmatrix}
~~ \text{or} ~~
\mathbf{C} =
\begin{bmatrix}
C_{Temp} & C_{DO} \\ 
C_{Temp} & C_{DO} \\ 
\vdots & \vdots \\ 
C_{Temp} & C_{DO}
\end{bmatrix}
$$

with

$$
\mathbf{c}_{t-k} =
\begin{bmatrix}
Temp_{t-k} \\
DO_{t-k}
\end{bmatrix}
$$


## Estimating species interactions

Here is an example of estimating species interactions from time series data

[Holmes, Ward & Scheuerell (2018) Chap 13](https://cran.r-project.org/web/packages/MARSS/vignettes/UserGuide.pdf)