R Tutorial

Complex Time-Series Models in R 1

Material heavily influenced by Nicholas J Clark’s excellent course on Ecological forecasting with mvgam and brms

Installation

• Check if your installation worked

library(cmdstanr)

cmdstan_version()

If this returns a version number (e.g., "2.37.0") then things are working properly.

Introduction

• So far we’ve developed relatively simple time-series models
• Linear time-series dependance
• Linear responses to environmental factors
• Normally distributed errors
• No model of observation error
• No interactions between species
• Can’t handle missing data
• Most of these are violated in ecological systems
• So start fitting more complex models

mvgam

• Models like this are not trivial to fit
• Use STAN
• Uses MCMC to explore parameter space to fit the model using Bayesian methods
• Typically requires learning a separate language - STAN is it’s own language
• This lets you write arbitrarily complex models, but really needs a course in Bayesian methods
• So, we’re going to use an R package called mvgam to implement our models
• We’re going to use it because it’s the simplest way to make complex time-series model in R

library(mvgam)
library(dplyr)
library(tsibble)

Data

• Data on the population dynamics of the Desert Pocket Mouse

pp_data <- read.csv("pp_abundance_by_month.csv")

• mvgam doesn’t currently work with tsibbles
• It requires an integer time variable be present in the data (like index)
• Also requires a a series variable, which needs to be a factor
• Helps when analyzing multiple time series at once (e.g., multiple species)

pp_data = read.csv("pp_abundance_by_month.csv") |>
    mutate(month = yearmonth(month)) |>
    mutate(time = as.numeric(month)) |>
    mutate(series = as.factor('PP')) |>
    tibble()
pp_data

• 124 months of data
• Reserve 24 for testing

data_train <- filter(pp_data, time <= max(time) - 24)
data_test <- filter(pp_data, time > max(time) - 24)

Simply time-series models in mvgam

• mvgam comes with a time-series visualization function like we had in fable

plot_mvgam_series(data = pp_data, y = 'abundance')

• Start by building something similar to what we’ve done before with the number of desert pocket mice
• Model has an exogenous driver - minimum temperature
• An autoregressive component
• Gaussian error

\[y_t = c + \beta_1 x_{1,t} + \beta_2 y_{t-1} + \mathcal{N}(0,\sigma^{2})\]

• Which can also be written as

\[y_t = \mathcal{N}(\mu_t,\sigma^{2})\] \[u_t = c + \beta_1 x_{1,t} + \beta_2 y_{t-1}\]

• Fit using the mvgam() function
• Follows a Base R model structure so start with the model formula
• Instead of including the AR component in the model we add a separate trend_model argument
• We’ll use AR(p = 1)
• Specify the error family = gaussian()
• Then we can specify the data for fitting the model and the data for making/evaluating forecasts

baseline_model = mvgam(abundance ~ mintemp,
                       trend_model = AR(p = 1),
                       family = gaussian(),
                       data = data_train)

Bayesian model fitting

• That’s a lot of output for fitting a model
• What’s going on?
• It is difficult to fit more complex time-series models
• One way to fit them is using Bayesian methods
• These methods iteratively search parameter space for the best parameter values
• Using something called Markov Chain Monte Carlo (MCMC)
• Draw 2D parameter space owith \(\beta_1\) and \(\beta_2\) on axes
• MCMC starts with some initial point
• Draw a point
• And the searches through the parameter space
• Draw iterative steps with connected lines
• Using an algorithm that eventually leads to points being sampled around the best fitting values
• Draw converged sampling with connected lines
• The "Iteration" lines are telling us that the model is working it’s way through this process
• The "Warmup" lines are iterations that are used to get in the right parameter space but then thrown away
• The "Sampling" lines are iterations that serve as samples for each of the parameters we are fitting
• Having multiple samples gives us the uncertainty in the parameters by looking at how variable those values are
• The different "chains" are because we typically go through this process multiple times with different starting points to check than we are converging to the right values
• We can look at the result of this fitting process using mcmc_plot

mcmc_plot(baseline_model, type = "trace", variable = c("mintemp", "ar1[1]", "sigma[1]"))

• The red lines at the bottom something isn’t quite right
• The model isn’t converged yet
• We could try running it longer
• But part of what’s going on here is that we’re using a poorly specified model given the data because we’ve assumed Gaussian errors
• So let’s look at the model & forecast and then keep going

• We can look at the model structure using summary()

summary(baseline_model)

• If we plot the model we’ll get some diagnostic plots

plot(baseline_model)

• The residuals are normally distributed
• And the residual autocorrelation isn’t significant

Bayesian model forecasting

• To make forecasts in mvgam we use the forecast() function (just like fable)
• Use newdata (not new_data as in fable) to include the test data for the driver forecasts

baseline_forecast = forecast(baseline_model, newdata = data_test)

• We can then plot the forecast

plot(baseline_forecast)

• These forecasts look similar to those we generated using fable
• Prediction intervals are regularly negative because we’ve assumed normally distributed error
• Actual counts can only be non-negative integers: 0, 1, 2…

Modeling count data

• To model count data explicitly we need observations that are integers
• We can use use Poisson error to accomplish this
• Instead of \(N(\mu, \sigma^2)\)

\[y_t = \mathrm{Pois}(\lambda_t)\]

• The Poisson distribution has one parameter \(\lambda\)
• Which is both the mean and the variance
• It generates only integer draws based on a mean

rpois(n = 10, lambda = 4.5)
hist(rpois(n = 1000, lambda = 4.5))

• If the mean is 4.5 sometimes you’ll draw a 4, sometimes a 5, sometimes a 0 or a 10

• \(\lambda_t\) can be a decimal
• We could expect an average of 4.5 rodents based on the environment even though we can only observe an integer number
• But \(\lambda_t\) does have to be positive because we can’t reasonably expect to see negative rodents

• To handle this we use a log link function to give us only positive values of \(\lambda_t\)
• The log link means that instead of modeling \(\lambda_t\) directly we model \(log(\lambda_t)\)

\[y_t = \mathrm{Pois}(\lambda_t)\] \[\mathrm{log(\lambda_t)} = c + \beta_1 x_{1,t} + \beta_2 y_{t-1}\]

• \(\lambda_t = e^{log(\lambda_t)}\) so even if the linear fit yields a negative value \(\lambda_t\) will always be positive

exp(-10)

• Change the error model and the link function in family

poisson_model = mvgam(abundance ~ mintemp,
                      trend_model = AR(p = 1),
                      family = poisson(link = "log"),
                      data = data_train)

• No more warnings at the end
• Look at the model

summary(poisson_model)

• *Discuss summary: sometimes rhat will still flag and increasing burnin may fix it

plot(poisson_model)

• Model structure looks OK, be we still have some residual seasonal lag we haven’t captured

• Now let’s look at the forecast

poisson_forecast = forecast(poisson_model, newdata = data_test)
plot(poisson_forecast)

• Now all of our predictions are positive!
• But while the point estimates seem reasonable the uncertainties seem really large