Time-varying effects in mvgam

Simulating time-varying effects

Simulate a time-varying coefficient using a squared exponential Gaussian Process function with length scale ρ=10. We will do this using an internal function from mvgam (the sim_gp function):

set.seed(1111)
N <- 200
beta_temp <- mvgam:::sim_gp(rnorm(1),
                            alpha_gp = 0.75,
                            rho_gp = 10,
                            h = N) + 0.5

A plot of the time-varying coefficient shows that it changes smoothly through time:

plot(beta_temp, type = 'l', lwd = 3, 
     bty = 'l', xlab = 'Time', ylab = 'Coefficient',
     col = 'darkred')
box(bty = 'l', lwd = 2)

Next we need to simulate the values of the covariate, which we will call temp (to represent temperature). In this case we just use a standard normal distribution to simulate this covariate:

temp <- rnorm(N, sd = 1)

Finally, simulate the outcome variable, which is a Gaussian observation process (with observation error) over the time-varying effect of temperature

out <- rnorm(N, mean = 4 + beta_temp * temp,
             sd = 0.25)
time <- seq_along(temp)
plot(out,  type = 'l', lwd = 3, 
     bty = 'l', xlab = 'Time', ylab = 'Outcome',
     col = 'darkred')
box(bty = 'l', lwd = 2)

Gather the data into a data.frame for fitting models, and split the data into training and testing folds.

data <- data.frame(out, temp, time)
data_train <- data[1:190,]
data_test <- data[191:200,]

The dynamic() function

Time-varying coefficients can be fairly easily set up using the s() or gp() wrapper functions in mvgam formulae by fitting a nonlinear effect of time and using the covariate of interest as the numeric by variable (see ?mgcv::s or ?brms::gp for more details). The dynamic() formula wrapper offers a way to automate this process, and will eventually allow for a broader variety of time-varying effects (such as random walk or AR processes). Depending on the arguments that are specified to dynamic, it will either set up a low-rank GP smooth function using s() with bs = 'gp' and a fixed value of the length scale parameter ρ, or it will set up a Hilbert space approximate GP using the gp() function with c=5/4 so that ρ is estimated (see ?dynamic for more details). In this first example we will use the s() option, and will mis-specify the ρ parameter here as, in practice, it is never known. This call to dynamic() will set up the following smooth: s(time, by = temp, bs = "gp", m = c(-2, 8, 2), k = 40)

mod <- mvgam(out ~ dynamic(temp, rho = 8, stationary = TRUE, k = 40),
             family = gaussian(),
             data = data_train,
             silent = 2)

Inspect the model summary, which shows how the dynamic() wrapper was used to construct a low-rank Gaussian Process smooth function:

summary(mod, include_betas = FALSE)

Because this model used a spline with a gp basis, it’s smooths can be visualised just like any other gam. We can plot the estimates for the in-sample and out-of-sample periods to see how the Gaussian Process function produces sensible smooth forecasts. Here we supply the full dataset to the newdata argument in plot_mvgam_smooth to inspect posterior forecasts of the time-varying smooth function. Overlay the true simulated function to see that the model has adequately estimated it’s dynamics in both the training and testing data partitions

plot_mvgam_smooth(mod, smooth = 1, newdata = data)
abline(v = 190, lty = 'dashed', lwd = 2)
lines(beta_temp, lwd = 2.5, col = 'white')
lines(beta_temp, lwd = 2)

We can also use plot_predictions from the marginaleffects package to visualise the time-varying coefficient for what the effect would be estimated to be at different values of temperature:

require(marginaleffects)
range_round = function(x){
  round(range(x, na.rm = TRUE), 2)
}
plot_predictions(mod, 
                 newdata = datagrid(time = unique,
                                    temp = range_round),
                 by = c('time', 'temp', 'temp'),
                 type = 'link')

This results in sensible forecasts of the observations as well

fc <- forecast(mod, newdata = data_test)
plot(fc)

The syntax is very similar if we wish to estimate the parameters of the underlying Gaussian Process, this time using a Hilbert space approximation. We simply omit the rho argument in dynamic to make this happen. This will set up a call similar to gp(time, by = 'temp', c = 5/4, k = 40).

mod <- mvgam(out ~ dynamic(temp, k = 40),
             family = gaussian(),
             data = data_train,
             silent = 2)

This model summary now contains estimates for the marginal deviation and length scale parameters of the underlying Gaussian Process function:

summary(mod, include_betas = FALSE)

Effects for gp() terms can also be plotted as smooths:

plot_mvgam_smooth(mod, smooth = 1, newdata = data)
abline(v = 190, lty = 'dashed', lwd = 2)
lines(beta_temp, lwd = 2.5, col = 'white')
lines(beta_temp, lwd = 2)

A State-Space Beta regression

mvgam can easily handle data that are bounded at 0 and 1 with a Beta observation model (using the mgcv function betar(), see ?mgcv::betar for details). First we will fit a simple State-Space model that uses an AR1 dynamic process model with no predictors and a Beta observation model:

mod0 <- mvgam(formula = survival ~ 1,
              trend_model = AR(),
              noncentred = TRUE,
              priors = prior(normal(-3.5, 0.5), class = Intercept),
              family = betar(),
              data = model_data,
              silent = 2)

The summary of this model shows good behaviour of the Hamiltonian Monte Carlo sampler and provides useful summaries on the Beta observation model parameters:

summary(mod0)

A plot of the underlying dynamic component shows how it has easily handled the temporal evolution of the time series:

plot(mod0, type = 'trend')

Including time-varying upwelling effects

Now we can increase the complexity of our model by constructing and fitting a State-Space model with a time-varying effect of the coastal upwelling index in addition to the autoregressive dynamics. We again use a Beta observation model to capture the restrictions of our proportional observations, but this time will include a dynamic() effect of CUI.apr in the latent process model. We do not specify the ρ parameter, instead opting to estimate it using a Hilbert space approximate GP:

mod1 <- mvgam(formula = survival ~ 1,
              trend_formula = ~ dynamic(CUI.apr, k = 25, scale = FALSE) - 1,
              trend_model = AR(),
              noncentred = TRUE,
              priors = prior(normal(-3.5, 0.5), class = Intercept),
              family = betar(),
              data = model_data,
              silent = 2)

The summary for this model now includes estimates for the time-varying GP parameters:

summary(mod1, include_betas = FALSE)

The estimates for the underlying dynamic process, and for the hindcasts, haven’t changed much:

plot(mod1, type = 'trend')

plot(mod1, type = 'forecast')

But the process error parameter σ is slightly smaller for this model than for the first model:

# Extract estimates of the process error 'sigma' for each model
mod0_sigma <- as.data.frame(mod0, variable = 'sigma', regex = TRUE) %>%
  dplyr::mutate(model = 'Mod0')
mod1_sigma <- as.data.frame(mod1, variable = 'sigma', regex = TRUE) %>%
  dplyr::mutate(model = 'Mod1')
sigmas <- rbind(mod0_sigma, mod1_sigma)

# Plot using ggplot2
require(ggplot2)
ggplot(sigmas, aes(y = `sigma[1]`, fill = model)) +
  geom_density(alpha = 0.3, colour = NA) +
  coord_flip()

Why does the process error not need to be as flexible in the second model? Because the estimates of this dynamic process are now informed partly by the time-varying effect of upwelling, which we can visualise on the link scale using plot():

plot(mod1, type = 'smooths', trend_effects = TRUE)

Comparing model predictive performances

A key question when fitting multiple time series models is whether one of them provides better predictions than the other. There are several options in mvgam for exploring this quantitatively. First, we can compare models based on in-sample approximate leave-one-out cross-validation as implemented in the popular loo package:

loo_compare(mod0, mod1)

The second model has the larger Expected Log Predictive Density (ELPD), meaning that it is slightly favoured over the simpler model that did not include the time-varying upwelling effect. However, the two models certainly do not differ by much. But this metric only compares in-sample performance, and we are hoping to use our models to produce reasonable forecasts. Luckily, mvgam also has routines for comparing models using approximate leave-future-out cross-validation. Here we refit both models to a reduced training set (starting at time point 30) and produce approximate 1-step ahead forecasts. These forecasts are used to estimate forecast ELPD before expanding the training set one time point at a time. We use Pareto-smoothed importance sampling to reweight posterior predictions, acting as a kind of particle filter so that we don’t need to refit the model too often (you can read more about how this process works in Bürkner et al. 2020).

lfo_mod0 <- lfo_cv(mod0, min_t = 30)
lfo_mod1 <- lfo_cv(mod1, min_t = 30)

The model with the time-varying upwelling effect tends to provides better 1-step ahead forecasts, with a higher total forecast ELPD

sum(lfo_mod0$elpds)
sum(lfo_mod1$elpds)

We can also plot the ELPDs for each model as a contrast. Here, values less than zero suggest the time-varying predictor model (Mod1) gives better 1-step ahead forecasts:

plot(x = 1:length(lfo_mod0$elpds) + 30,
     y = lfo_mod0$elpds - lfo_mod1$elpds,
     ylab = 'ELPDmod0 - ELPDmod1',
     xlab = 'Evaluation time point',
     pch = 16,
     col = 'darkred',
     bty = 'l')
abline(h = 0, lty = 'dashed')

A useful exercise to further expand this model would be to think about what kinds of predictors might impact measurement error, which could easily be implemented into the observation formula in mvgam(). But for now, we will leave the model as-is.