R Tutorial

In development. Not ready for teaching

State space models

• Time-series model
• Only first order autoregressive component
• Separately model
  • the process model - how the system evolves in time or space
  • the observation model - observation error or indirect observations
• Estimates the true value of the underlying latent state variables

Model

Draw on board while walking through models

y_t-1    y_t    y_t+1
  |       |       |
x_t-1 -> x_t -> x_t+1   Process model

Process model

• What is actually happening in the system
• First order autoregressive component

x_t+1 = f(x_t) + e_t

• Simple linear model is AR1:

x_t+1 = b0 + b1 * x_t + e_t

Observation model

• Counts of rodents in traps aren’t perfect measures of the number of rodents (which are what should be changing in the process model and what we care about)
• So model this imperfect observation

y_t = Pois(x_t)

• Can be much more complicated

…

• State space model of AR1 + rain w/Poisson error

state_space_model = mvgam(abundance ~ 1,
                          trend_formula = ~ mintemp,
                          trend_model = "AR1",
                          family = poisson(link = "log"),
                          data = data_train,
                          newdata = data_test)
plot(state_space_model, type = "forecast")

output = state_space_model$model_output

• Normally would want several chains with different starting positions to avoid local minima

• Send to JAGS

j.model   <- jags.model (file = textConnection(RandomWalk),
                         data = data,
                         inits = init,
                         n.chains = 1)

• Burn in

jags.out   <- coda.samples (model = j.model,
                            variable.names = c("tau_proc","tau_obs"),
                            n.iter = 10000)
plot(jags.out)

• Sample from MCMC with full vector of X’s
• This starts sampling from the point were the previous run of coda.samples ends so it gets rid of the burn-in samples

jags.out   <- coda.samples (model = j.model,
                            variable.names = c("x","tau_proc","tau_obs"),
                            n.iter = 10000)

• Visualize
• Convert the output into a matrix & drop parameters

out <- as.matrix(jags.out)
xs <- out[,3:ncol(out)]

• Point predictions are averages across MCMC samples

predictions <- colMeans(xs)
plot(time, predictions, type = "l")

• And this looks very similar to the observed dynamics of y

• Add prediction intervals as range containing 95% of MCMC samples

ci <- apply(xs, 2, quantile, c(0.025, 0.975))
lines(time, ci[1,], lty = "dashed", col = "blue")
lines(time, ci[2,], lty = "dashed", col = "blue")

• These are very narrow prediction intervals, so the model appears to be very confident
• But it’s important to keep in mind that when fitting the value of x at time t, the model has access to the value of y at time t
• And the y is present it isn’t being estimated, it’s just the observed value
• So, will this model forecast well?

Forecasting

• To make forecasts using a JAGS model we include data for y that is NA
• This tells the model that we don’t know the values and therefore the model estimates them as part of the fitting process
• To make a true forecast we would add one NA to the end of y for each time step we wanted to forecast
• To hindcast or backcast like we replace the values for y that are part of the test set with NA
• We’ll hindcast, so to do this we’ll replace the last year of y values with NA and then compare the final year of data to our predictions

Make these changes at top of script and rerun

data$y[(length(y)-51):length(y)] = NA
jags.out   <- coda.samples (model = j.model,
                            variable.names = c("y","tau_proc","tau_obs"),
                            n.iter = 10000)

• We can see from plotting the predictions that the forecast doesn’t look promising
• Without the observed data to influence the estimates of x[t] the model predicts little change over the forecast year
• We can directly compare this to the empirical data by adding it to the plot

lines(time, y)

• So the point estimates don’t perform well
• This raises the question of whether the model accurately predicts that it is uncertain when making forecasts
• Plotting the prediction intervals suggests that it does
• They very quickly expand towards zero and the upper limits of the data