R Tutorial

More complex time-series models: Part 2

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

Starting Material

• Start where we left off last time

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

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

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

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


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

Visualizing environmental drivers

• As our models get more complicate it’s important to make sure we understand what they are doing
• We can visualize the environmental component of the model using type = pterms
• pterms is short for “parametric terms”
• Which is what we have since we modeled a linear relationship with a fixed slope

plot(poisson_model, type = 'pterms')

• This shows a linear relationship
• But that is the relationship between log(abundance) and mintemp
• So let’s look at the actual relationship

plot_predictions(poisson_model, condition = "mintemp")

• Since log(abundance) is linearly related to temperature
• This makes the response to untransformed abundance exponential
• That doesn’t feel right
• Could be part of the reason for our really large upper prediction intervals

Non-linear responses

• Linear relationships between abundance and the environment are unlikely to hold over reasonable environmental gradients
• Typically we think of species as having some optimal environmental value
• With decreasing performance as you move away from that environmental value in both directions
• Draw standard performance response curve
• We can model using in mvgam by using a Generalized Additive Model that fits a smoother to the relationship

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

• Visualizing the environmental component of the model shows that it is now non-linear

plot(poisson_gam_model, type = 'smooths')

• It doesn’t have an optimal value, but it does saturate
• This means that once it’s warm enough then temperature stops influencing abundance
• This is the smooth on log(abundance)
• Look at the relationship with raw abundance

plot_predictions(poisson_gam_model, condition = "mintemp")

• The GAM now gives us a more realistic response to temperature
• When it’s cold abundance is always near zero
• Even if gets a little less cold
• Abundance then increases with temperature
• But starts to asymptote at higher temperatures

• How does this influence the forecasts

poisson_gam_forecast = forecast(poisson_gam_model, newdata = data_test)
plot(poisson_gam_forecast)

• The extremely wide prediction intervals have been reduced substantially
• If we look at the diagnostics

plot(poisson_gam_model)

• We can also see that the residual season autocorrelation is improved

State space models (optional)

• To add an explicit observation model we use the

state_space_model = mvgam(abundance ~ 1,
                          trend_formula = ~ s(mintemp),
                          trend_model = AR(p = 1),
                          family = poisson(link = "log"),
                          data = data_train)
state_space_forecast = forecast(state_space_model, newdata = data_test)
plot(state_space_forecast)

Some notes on Bayesian forecasting

• Since forecasting and model fitting are all part of the same process in Bayesian models
• We can also just incorporate the test data directly in the model fitting step

poisson_gam_model = mvgam(abundance ~ s(mintemp),
                      trend_model = AR(p = 1),
                      family = poisson(link = "log"),
                      data = data_train,
                      newdat = data_test)

• And then plot forecasts directly from the model object

plot(poisson_gam_model, type = "forecast")

• If you want to know why we’re using mvgam
• Here’s what the code to build this model directly in STAN looks like

code(poisson_gam_model)

Evaluation and comparison

• To evaluate the forecasts from mvgam we can use the score function

scores <- score(poisson_gam_forecast)

• The output shows us two things we’ve seen before
• Whether or not each observed value falls within the prediction interval
• With a default prediction interval of 90%
• And the continuous rank probability score
• Other scores are available by changing the optional score argument
• We can also change the prediction interval using the interval_width argument

scores <- score(poisson_gam_forecast, interval_width = 0.5)

• If we want to calculate the coverage we can sum

in_interval = scores$PP$in_interval
length(in_interval[in_interval == 1]) / length(in_interval)