Univariate state-space models

1. state (process) model

2. observation model

Consist of 2 parts

Describes the true state of nature over time

States of nature might be:

• animal's location

• density of organisms

• reproductive state

We can think of changes in the state of nature being driven by a combination of

• intrinsic (eg, fecundity), and

• extrinsic factors (eg, temperature)

Some of the extrinsic drivers may be unknown

In time series modeling, we often call these unknown extrinsic factors process errors

Data = Truth + Error

Consider for a moment the very act of collecting data

Information gathered depends on many factors

• Environmental conditions (eg, turbidity, currents)

• Behavior (eg, vertical migration, threat avoidance)

• Demographics (age, sex, maturity)

• Sampling design/coverage

Linear regression

State model

\[ x_t = \alpha + \beta t \]

Observation model

\[ y_t = x_t + v_t ~ \text{with} ~ v_t \sim \text{N}(0,r) \]

All of the variance in \(y_t\) is due to observation error

Errors \(v_t = y_t - x_t\)

Biased random walk

State model

\[ x_t = x_{t-1} + u + w_t ~ \text{with} ~ w_t \sim \text{N}(0,q) \]

Observation model

\[ y_t = x_t \]

All of the variance in \(y_t\) is due to process error

Errors: \(w_t = y_t - y_{t-1}\)

Biased random walk with observation error

State model

\[ x_t = x_{t-1} + u + w_t ~ \text{with} ~ w_t \sim \text{N}(0,q) \]

Observation model

\[ y_t = x_t + v_t ~ \text{with} ~ v_t \sim \text{N}(0,r) \]

The variance in \(y_t\) is due to both process & observation error

How is it possible to separate out both the process and observation errors?

They have different temporal patterns

Stochastic Gompertz model

\[ N_t = N_{t-1} \text{e}^{(r_{max} + (b - 1) \log(N_{t-1}))} \text{e}^{w_t} ~ \text{with} ~ w_t \sim \text{N}(0,q) \]

Stochastic Gompertz model

\[ N_t = N_{t-1} \text{e}^{(r_{max} + (b - 1) \log(N_{t-1}))} \text{e}^{w_t} ~ \text{with} ~ w_t \sim \text{N}(0,r) \]

Taking the log of both sides and substituting \(x_t\) for \(\log(N_t)\)

\[ \begin{align} \log(N_t) & = \log(N_{t-1}) + r_{max} + (b - 1) \log(N_{t-1}) + w_t \\ & \Downarrow \\ x_t &= x_{t-1} + r_{max} + (b - 1) x_{t-1} + w_t \\ &= x_{t-1} + r_{max} + b x_{t-1} - x_{t-1} + w_t \\ &= r_{max} + b x_{t-1} + w_t \end{align} \]

What is the mean of \(\{x_t\}\)?

\[ x_t = r_{max} + b x_{t-1} + w_t \]

What is the mean of \(\{x_t\}\)?

\[ x_t = r_{max} + b x_{t-1} + w_t \]

\[ \text{E}(x_t) = \frac{r_{max}}{1 - b} \]

\[ x_t = r_{max} + b x_{t-1} + w_t \]

It turns out that \(r_{max}\) and \(b\) are confounded and create a likelihood surface with a strong ridge

This means we need a long time series or mulitple observations of the process

In practice, we will de-mean our data and instead use a familiar AR(1) model

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

If our population censuses contain observation or sampling errors, we can use a state-space version

\[ x_t = b x_{t-1} + w_t \\ y_t = x_t + v_t \]

Simulate the AR(1) state model

## number of time steps
TT <- 50
## strength of density-dependence (0 < b < 1)
bb <- 0.5
## time series of process errors with SD = 1
ww <- rnorm(TT, 0, sqrt(1))
## initialize state & set x0 = w0
xx <- ww
## loop over time steps
for(t in 2:TT) {
  xx[t] <- bb * xx[t-1] + ww[t]
}

Add some observation error

## obs errors with var = 0.5
vv <- rnorm(TT, 0, sqrt(0.5))
## obs data
yy <- xx + vv

Plot the state and observation

plot.ts(xx, lwd = 2, type = "o", pch = 16, col = "gray",
        ylim = c(min(xx,yy), max(xx,yy)),
        ylab = expression(italic(x[t])~~or~~italic(y[t])))
lines(yy, lwd = 2, type = "o", pch = 16, col = "blue")

We can use the MARSS pkg for R to fit state-space models

To do so, we must write the model out just as we would on a white board

Here are the main function arguments

MARSS(y, model = NULL, inits = NULL, control = NULL, ...)

• y is an \(n \times T\) matrix of data (observations)

• model is a list that defines the state-space model

• inits is an optional list of initial values

• control is an optional list of control parameters

MARSS() uses a slightly expanded form of state-space model, such that

\[ \begin{align} x_t &= b x_{t-1} + w_t \\ y_t &= x_t + v_t \\ & \Downarrow \\ x_t &= b x_{t-1} + u + w_t \\ y_t &= Z x_t + a + v_t \end{align} \]

with \(u = 0\), \(Z = 1\), and \(a = 0\)

\[ x_t = b x_{t-1} + u + w_t ~ \text{with} ~ w_t \sim \text{N}(0,q) \\ y_t = Z x_t + a + v_t ~ \text{with} ~ v_t \sim \text{N}(0,r) \]

MARSS() works on matrices, so we have to define any scalar as a \(1 \times 1\) matrix

mod_list <- list(
  ## state model
  B = matrix("b"), U = matrix(0), Q = matrix("q"),
  ## obs model
  Z = matrix(1), A = matrix(0), R = matrix("r")
  )

What about the initial state \(x_0\)?

\[ x_0 \sim \text{N}(\pi,\lambda) \]

Estimating both \(\pi\) and \(\lambda\) is nearly impossible, so the default is to treat \(x_0\) as fixed, but unknown, effect rather than a random effect, such that

\[ x_0 \sim \text{N}(\pi,0) \]

MARSS() will do this for us

We can ignore the inits and control arguments for now and fit the model

## load MARSS package
library(MARSS)
## define the data as an N (rows) x T (cols) matrix
dat <- matrix(yy, nrow = 1, ncol = TT)
## fit the model
mod_fit <- MARSS(y = dat, model = mod_list)

## Success! abstol and log-log tests passed at 68 iterations.
## Alert: conv.test.slope.tol is 0.5.
## Test with smaller values (<0.1) to ensure convergence.
## 
## MARSS fit is
## Estimation method: kem 
## Convergence test: conv.test.slope.tol = 0.5, abstol = 0.001
## Estimation converged in 68 iterations. 
## Log-likelihood: -80.44073 
## AIC: 168.8815   AICc: 169.7704   
##  
##       Estimate
## R.r      0.436
## B.b      0.484
## Q.q      0.959
## x0.x0    0.260
## Initial states (x0) defined at t=0
## 
## Standard errors have not been calculated. 
## Use MARSSparamCIs to compute CIs and bias estimates.

We can include covariates (explanatory variables) in our models as well

Covariates in the state model affect the underlying process

\[ x_t = b x_{t-1} + C c_t + w_t \]

We can include covariates (explanatory variables) in our models as well

Covariates in the state model affect the underlying process

\[ x_t = b x_{t-1} + C c_t + w_t \]

Covariates in the observation model are offsets to the underlying process

\[ y_t = Z x_t + D d_t + v_t \]