Exercises and associated data

The data and modelling objects created in this notebook can be downloaded directly to save computational time.


Users who wish to complete the exercises can download a small template R script. Assuming you have already downloaded the data objects above, this script will load all data objects so that the steps used to create them are not necessary to tackle the exercises.

Load libraries and time series data

This tutorial relates to content covered in Lecture 5, and relies on the following packages for manipulating data, shaping time series, fitting dynamic regression models and plotting:

library(dplyr)
library(mvgam) 
library(tidybayes)
library(bayesplot)
library(gratia)
library(ggplot2); theme_set(theme_classic())
library(marginaleffects)

This tutorial will focus on one of mvgam’s more advanced features: the ability to fit dynamic models in a State-Space representation. We have not gone into too much detail about how these models work (but see this outstanding resource from E. E. Holmes, M. D. Scheuerell, and E. J. Ward for many useful details, along with this nice lecture by E. E. Holmes).

Illustration of a basic State-Space model, which assumes that a latent dynamic process (X) can evolve independently from the way we take observations (Y) of that process
Illustration of a basic State-Space model, which assumes that a latent dynamic process (X) can evolve independently from the way we take observations (Y) of that process


Briefly, these models allow us to separately make inferences about the underlying dynamic process model that we are interested in (i.e. the evolution of a time series or a collection of time series) and the observation model (i.e. the way that we survey / measure this underlying process). This is extremely useful in ecology because our observations are always imperfect / noisy measurements of the thing we are interested in measuring. It is also helpful because we often know that some covariates will impact our ability to measure accurately (i.e. we cannot take accurate counts of rodents if there is a thunderstorm happening) while other covariate impact the underlying process (it is highly unlikely that rodent abundance responds to one storm, but instead probably responds to longer-term weather and climate variation). A State-Space model allows us to model both components in a single unified modelling framework. A major advantage of mvgam is that it can include nonlinear effects and random effects in BOTH model components while also capturing dynamic processes. I am not aware of any other packages that can easily do this, but of course there may be some.

Lake Washington plankton data

The data we will use to illustrate how we can fit State-Space models in mvgam are from a long-term monitoring study of plankton counts (cells per mL) taken from Lake Washington in Washington, USA. The data are available as part of the MARSS package and can be downloaded using the following:

load(url('https://github.com/atsa-es/MARSS/raw/master/data/lakeWAplankton.rda'))

We will work with five different groups of plankton:

outcomes <- c('Greens', 'Bluegreens', 'Diatoms', 'Unicells', 'Other.algae')

As usual, preparing the data into the correct format for mvgam modelling takes a little bit of wrangling in dplyr:

plankton_data <- xfun::cache_rds(do.call(rbind, lapply(outcomes, function(x){
  
  # create a group-specific dataframe with counts labelled 'y'
  # and the group name in the 'series' variable
  data.frame(year = lakeWAplanktonTrans[, 'Year'],
             month = lakeWAplanktonTrans[, 'Month'],
             y = lakeWAplanktonTrans[, x],
             series = x,
             temp = lakeWAplanktonTrans[, 'Temp'])})) %>%
  
  # change the 'series' label to a factor
  dplyr::mutate(series = factor(series)) %>%
  
  # filter to only include some years in the data
  dplyr::filter(year >= 1965 & year < 1975) %>%
  dplyr::arrange(year, month) %>%
  dplyr::group_by(series) %>%
  
  # z-score the counts so they are approximately standard normal
  dplyr::mutate(y = as.vector(scale(y))) %>%
  
  # add the time indicator
  dplyr::mutate(time = dplyr::row_number()) %>%
  dplyr::ungroup())
# loop across each plankton group to create the long datframe
plankton_data <- do.call(rbind, lapply(outcomes, function(x){
  
  # create a group-specific dataframe with counts labelled 'y'
  # and the group name in the 'series' variable
  data.frame(year = lakeWAplanktonTrans[, 'Year'],
             month = lakeWAplanktonTrans[, 'Month'],
             y = lakeWAplanktonTrans[, x],
             series = x,
             temp = lakeWAplanktonTrans[, 'Temp'])})) %>%
  
  # change the 'series' label to a factor
  dplyr::mutate(series = factor(series)) %>%
  
  # filter to only include some years in the data
  dplyr::filter(year >= 1965 & year < 1975) %>%
  dplyr::arrange(year, month) %>%
  dplyr::group_by(series) %>%
  
  # z-score the counts so they are approximately standard normal
  dplyr::mutate(y = as.vector(scale(y))) %>%
  
  # add the time indicator
  dplyr::mutate(time = dplyr::row_number()) %>%
  dplyr::ungroup()

Inspect the data structure

head(plankton_data)
## # A tibble: 6 × 6
##    year month       y series       temp  time
##   <dbl> <dbl>   <dbl> <fct>       <dbl> <int>
## 1  1965     1 -0.542  Greens      -1.23     1
## 2  1965     1 -0.344  Bluegreens  -1.23     1
## 3  1965     1 -0.0768 Diatoms     -1.23     1
## 4  1965     1 -1.52   Unicells    -1.23     1
## 5  1965     1 -0.491  Other.algae -1.23     1
## 6  1965     2 NA      Greens      -1.32     2
dplyr::glimpse(plankton_data)
## Rows: 600
## Columns: 6
## $ year   <dbl> 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 1965, 196…
## $ month  <dbl> 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, …
## $ y      <dbl> -0.54241769, -0.34410776, -0.07684901, -1.52243490, -0.49055442…
## $ series <fct> Greens, Bluegreens, Diatoms, Unicells, Other.algae, Greens, Blu…
## $ temp   <dbl> -1.2306562, -1.2306562, -1.2306562, -1.2306562, -1.2306562, -1.…
## $ time   <int> 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, …

Note that we have z-scored the counts in this example as that will make it easier to specify priors for this simple example (though note that this is not completely necessary; it is often better to build a model that respects the properties of the actual outcome variables)

plot_mvgam_series(data = plankton_data, series = 'all')

As usual, check the data for NAs:

image(is.na(t(plankton_data)), axes = F,
      col = c('grey80', 'darkred'))
axis(3, at = seq(0,1, len = NCOL(plankton_data)), 
     labels = colnames(plankton_data))

Manipulate data for modeling

We have some missing observations, but of course this isn’t an issue for modelling in mvgam. A useful property to understand about these counts is that they tend to be highly seasonal. Below are some plots of z-scored counts against the z-scored temperature measurements in the lake for each month:

plankton_data %>%
  dplyr::filter(series == 'Other.algae') %>%
  ggplot(aes(x = time, y = temp)) +
  geom_line(size = 1.1) +
  geom_line(aes(y = y), col = 'white',
            size = 1.2) +
  geom_line(aes(y = y), col = 'darkred',
            size = 1.1) +
  ylab('z-score') +
  xlab('Time') +
  ggtitle('Temperature (black) vs Other algae (red)')
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

plankton_data %>%
  dplyr::filter(series == 'Diatoms') %>%
  ggplot(aes(x = time, y = temp)) +
  geom_line(size = 1.1) +
  geom_line(aes(y = y), col = 'white',
            size = 1.2) +
  geom_line(aes(y = y), col = 'darkred',
            size = 1.1) +
  ylab('z-score') +
  xlab('Time') +
  ggtitle('Temperature (black) vs Diatoms (red)')

plankton_data %>%
  dplyr::filter(series == 'Greens') %>%
  ggplot(aes(x = time, y = temp)) +
  geom_line(size = 1.1) +
  geom_line(aes(y = y), col = 'white',
            size = 1.2) +
  geom_line(aes(y = y), col = 'darkred',
            size = 1.1) +
  ylab('z-score') +
  xlab('Time') +
  ggtitle('Temperature (black) vs Greens (red)')

We will have to try and capture this seasonality in our process model, which should be easy to do given the flexibility of GAMs. Next we will split the data into training and testing splits:

plankton_train <- plankton_data %>%
  dplyr::filter(time <= 112)
plankton_test <- plankton_data %>%
  dplyr::filter(time > 112)

Exercises

  1. Calculate the number of timepoints in the training data that have non-missing observations for all five time series.

Capturing seasonality

Now time to fit some models. This requires a bit of thinking about how we can best tackle the seasonal variation and the likely dependence structure in the data. These algae are interacting as part of a complex system within the same lake, so we certainly expect there to be some lagged cross-dependencies underling their dynamics. But if we do not capture the seasonal variation, our multivariate dynamic model will be forced to try and capture it, which could lead to poor convergence and unstable results (we could feasibly capture cyclic dynamics with a more complex multi-species Lotka-Volterra model, but ordinary differential equation approaches are beyond the scope of this workshop).

First we will fit a model that does not include a dynamic component, just to see if it can reproduce the seasonal variation in the observations. This model builds from our hierarchical GAMs above by introducing hierarchical multidimensional smooths. It includes a “global” tensor product of the month and temp variables, capturing our expectation that algal seasonality responds to temperature variation. But this response should depend on when in the year these temperatures are recorded (i.e. a response to warm temperatures in Spring should be different to a response to warm temperatures in Autumn). The model also fits series-specific deviation smooths (i.e. one tensor product per series) to capture how each algal group’s seasonality differs from the overall “global” seasonality. Note that we do not include series-specific intercepts in this model because each series was z-scored to have a mean of 0.

notrend_mod <- mvgam(y ~ 
                       # tensor of temp and month to capture
                       # "global" seasonality
                       te(temp, month, k = c(4, 4)) +
                       
                       # series-specific deviation tensor products
                       te(temp, month, k = c(4, 4), by = series),
                     family = gaussian(),
                     data = plankton_train,
                     newdata = plankton_test,
                     trend_model = 'None')

The “global” tensor product smooth function can be quickly visualized using gratia:

gratia::draw(notrend_mod$mgcv_model, select = 1)

We can then plot the deviation smooths for each algal group to see how they vary from the “global” pattern:

gratia::draw(notrend_mod$mgcv_model, select = 2)

gratia::draw(notrend_mod$mgcv_model, select = 3)

gratia::draw(notrend_mod$mgcv_model, select = 4)

gratia::draw(notrend_mod$mgcv_model, select = 5)

gratia::draw(notrend_mod$mgcv_model, select = 6)

These multidimensional smooths have done a good job of capturing the seasonal variation in our observations:

plot(notrend_mod, type = 'forecast', series = 1)
## Out of sample CRPS:
## 6.7267833552403

plot(notrend_mod, type = 'forecast', series = 2)
## Out of sample CRPS:
## 6.7757230198503

plot(notrend_mod, type = 'forecast', series = 3)
## Out of sample CRPS:
## 4.10374012676295

plot(notrend_mod, type = 'forecast', series = 4)
## Out of sample CRPS:
## 3.57032374064579

plot(notrend_mod, type = 'forecast', series = 5)
## Out of sample CRPS:
## 2.90199373806287

Multiseries dynamics

The basic model gives us confidence that we can capture the seasonal variation in the observations. But the model has not captured the remaining temporal dynamics, which is obvious when we inspect Dunn-Smyth residuals for each series:

plot(notrend_mod, type = 'residuals', series = 1)

plot(notrend_mod, type = 'residuals', series = 2)

plot(notrend_mod, type = 'residuals', series = 3)

plot(notrend_mod, type = 'residuals', series = 4)

plot(notrend_mod, type = 'residuals', series = 5)

Now it is time to get into multivariate State-Space models. We will fit two models that can both incorporate lagged cross-dependencies in the latent process models. The first model assumes that the process errors operate independently from one another, while the second assumes that there may be contemporaneous correlations in the process errors. Both models include a Vector Autoregressive component for the process means, and so both can model complex community dynamics. The models can be described mathematically as follows:

\[\begin{align*} \boldsymbol{count}_t & \sim \text{Normal}(\mu_{obs[t]}, \sigma_{obs}) \\ \mu_{obs[t]} & = process_t \\ process_t & \sim \text{MVNormal}(\mu_{process[t]}, \Sigma_{process}) \\ \mu_{process[t]} & = VAR * process_{t-1} + f_{global}(\boldsymbol{month},\boldsymbol{temp})_t + f_{series}(\boldsymbol{month},\boldsymbol{temp})_t \\ f_{global}(\boldsymbol{month},\boldsymbol{temp}) & = \sum_{k=1}^{K}b_{global} * \beta_{global} \\ f_{series}(\boldsymbol{month},\boldsymbol{temp}) & = \sum_{k=1}^{K}b_{series} * \beta_{series} \end{align*}\]

Here you can see that there are no terms in the observation model apart from the underlying process model. But we could easily add covariates into the observation model if we felt that they could explain some of the systematic observation errors. We also assume independent observation processes (there is no covariance structure in the observation errors \(\sigma_{obs}\)). At present, mvgam does not support multivariate observation models. But this feature will be added in future versions. However the underlying process model is multivariate, and there is a lot going on here. This component has a Vector Autoregressive part, where the process mean at time \(t\) \((\mu_{process[t]})\) is a vector that evolves as a function of where the vector-valued process model was at time \(t-1\). The \(VAR\) matrix captures these dynamics with self-dependencies on the diagonal and possibly asymmetric cross-dependencies on the off-diagonals, while also incorporating the nonlinear smooth functions that capture seasonality for each series. The contemporaneous process errors are modeled by \(\Sigma_{process}\), which can be constrained so that process errors are independent (i.e. setting the off-diagonals to 0) or can be fully parameterized using a Cholesky decomposition (using Stan’s \(LKJcorr\) distribution to place a prior on the strength of inter-species correlations). For those that are interested in the inner-workings, mvgam makes use of a recent breakthrough by Sarah Heaps to enforce stationarity of Bayesian VAR processes. This is advantageous as we often don’t expect forecast variance to increase without bound forever into the future, but many estimated VARs tend to behave this way.


Ok that was a lot to take in. Let’s fit some models to try and inspect what is going on and what they assume. But first, we need to update mvgam’s default priors for the observation and process errors. By default, mvgam uses a fairly wide Student-T prior on these parameters to avoid being overly informative. But our observations are z-scored and so we do not expect very large process or observation errors. However, we also do not expect very small observation errors either as we know these measurements are not perfect. So let’s update the priors for these parameters. In doing so, you will get to see how the formula for the latent process (i.e. trend) model is used in mvgam:

priors <- get_mvgam_priors(
  # observation formula, which has no terms in it
  y ~ -1,
  
  # process model formula, which includes the smooth functions
  trend_formula = ~ te(temp, month, k = c(4, 4)) +
    te(temp, month, k = c(4, 4), by = trend),
  
  # VAR1 model with uncorrelated process errors
  trend_model = 'VAR1',
  family = gaussian(),
  data = plankton_train)

Get names of all parameters whose priors can be modified:

priors[, 3]
##  [1] "(Intercept)"                                                                                                                                                                                                                                                           
##  [2] "process error sd"                                                                                                                                                                                                                                                      
##  [3] "diagonal autocorrelation population mean"                                                                                                                                                                                                                              
##  [4] "off-diagonal autocorrelation population mean"                                                                                                                                                                                                                          
##  [5] "diagonal autocorrelation population variance"                                                                                                                                                                                                                          
##  [6] "off-diagonal autocorrelation population variance"                                                                                                                                                                                                                      
##  [7] "shape1 for diagonal autocorrelation precision"                                                                                                                                                                                                                         
##  [8] "shape1 for off-diagonal autocorrelation precision"                                                                                                                                                                                                                     
##  [9] "shape2 for diagonal autocorrelation precision"                                                                                                                                                                                                                         
## [10] "shape2 for off-diagonal autocorrelation precision"                                                                                                                                                                                                                     
## [11] "observation error sd"                                                                                                                                                                                                                                                  
## [12] "te(temp,month) smooth parameters, te(temp,month):trendtrend1 smooth parameters, te(temp,month):trendtrend2 smooth parameters, te(temp,month):trendtrend3 smooth parameters, te(temp,month):trendtrend4 smooth parameters, te(temp,month):trendtrend5 smooth parameters"

And their default prior distributions:

priors[, 4]
##  [1] "(Intercept) ~ student_t(3, -0.1, 2.5);"
##  [2] "sigma ~ student_t(3, 0, 2.5);"         
##  [3] "es[1] = 0;"                            
##  [4] "es[2] = 0;"                            
##  [5] "fs[1] = sqrt(0.455);"                  
##  [6] "fs[2] = sqrt(0.455);"                  
##  [7] "gs[1] = 1.365;"                        
##  [8] "gs[2] = 1.365;"                        
##  [9] "hs[1] = 0.071175;"                     
## [10] "hs[2] = 0.071175;"                     
## [11] "sigma_obs ~ student_t(3, 0, 2.5);"     
## [12] "lambda_trend ~ normal(5, 30);"

Setting priors is easy in mvgam as you can use brms routines. Here we use more informative Normal priors for both error components, but we impose a lower bound of 0.2 for the observation errors:

priors <- c(prior(normal(0.5, 0.1), class = sigma_obs, lb = 0.2),
            prior(normal(0.5, 0.25), class = sigma))

You may have noticed something else unique about this model: there is no intercept term in the observation formula. This is because a shared intercept parameter can sometimes be unidentifiable with respect to the latent VAR process, particularly if our series have similar long-run averages (which they do in this case because they were z-scored). We will often get better convergence in these State-Space models if we drop this parameter. mvgam accomplishes this by fixing the coefficient for the intercept to zero. Now we can fit the first model, which assumes that process errors are contemporaneously uncorrelated

Now we can fit the first model, which assumes that process errors are contemporaneously uncorrelated

var_mod <- mvgam(  
  # observation formula, which is empty
  y ~ -1,
  
  # process model formula, which includes the smooth functions
  trend_formula = ~ te(temp, month, k = c(4, 4)) +
    te(temp, month, k = c(4, 4), by = trend),
  
  # VAR1 model with uncorrelated process errors
  trend_model = 'VAR1',
  family = gaussian(),
  data = plankton_train,
  newdata = plankton_test,
  
  # include the updated priors
  priors = priors)

Inspecting SS models

This model’s summary is a bit different to other mvgam summaries. It separates parameters based on whether they belong to the observation model or to the latent process model. This is because we may often have covariates that impact the observations but not the latent process, so we can have fairly complex models for each component. You will notice that some parameters have not fully converged, particularly for the VAR coefficients (called A in the output) and for the process errors (Sigma). Note that we set include_betas = FALSE to stop the summary from printing output for all of the spline coefficients, which can be dense and hard to interpret:

summary(var_mod, include_betas = FALSE)
## GAM observation formula:
## y ~ 1
## 
## GAM process formula:
## ~te(temp, month, k = c(4, 4)) + te(temp, month, k = c(4, 4), 
##     by = trend)
## 
## Family:
## gaussian
## 
## Link function:
## identity
## 
## Trend model:
## VAR1
## 
## N process models:
## 5 
## 
## N series:
## 5 
## 
## N timepoints:
## 120 
## 
## Status:
## Fitted using Stan 
## 4 chains, each with iter = 1500; warmup = 1000; thin = 1 
## Total post-warmup draws = 2000
## 
## 
## Observation error parameter estimates:
##              2.5%  50% 97.5% Rhat n_eff
## sigma_obs[1] 0.20 0.26  0.34 1.00   452
## sigma_obs[2] 0.26 0.41  0.54 1.03   143
## sigma_obs[3] 0.44 0.64  0.81 1.02   130
## sigma_obs[4] 0.23 0.37  0.50 1.02   199
## sigma_obs[5] 0.32 0.44  0.55 1.01   230
## 
## GAM observation model coefficient (beta) estimates:
##             2.5% 50% 97.5% Rhat n_eff
## (Intercept)    0   0     0  NaN   NaN
## 
## Process model VAR parameter estimates:
##          2.5%     50%  97.5% Rhat n_eff
## A[1,1] -0.160  0.4700 0.8300 1.05    89
## A[1,2] -0.430 -0.0460 0.2100 1.02   357
## A[1,3] -0.580 -0.0480 0.3500 1.01   339
## A[1,4] -0.330  0.0430 0.5000 1.01   446
## A[1,5] -0.086  0.1600 0.6000 1.03   211
## A[2,1] -0.180  0.0110 0.2400 1.01   380
## A[2,2]  0.620  0.7800 0.9100 1.01   506
## A[2,3] -0.410 -0.1400 0.0430 1.01   405
## A[2,4] -0.029  0.1200 0.3600 1.01   230
## A[2,5] -0.054  0.0660 0.2000 1.00   677
## A[3,1] -0.430  0.0037 0.3700 1.01   270
## A[3,2] -0.550 -0.2300 0.0068 1.03   134
## A[3,3]  0.019  0.4000 0.7000 1.03   177
## A[3,4] -0.032  0.2500 0.7000 1.02   161
## A[3,5] -0.059  0.1400 0.4300 1.02   278
## A[4,1] -0.220  0.0500 0.3600 1.02   212
## A[4,2] -0.110  0.0490 0.2500 1.01   329
## A[4,3] -0.450 -0.1300 0.1100 1.03   207
## A[4,4]  0.500  0.7500 0.9700 1.02   252
## A[4,5] -0.200 -0.0310 0.1400 1.00   536
## A[5,1] -0.300  0.0620 0.4500 1.00   327
## A[5,2] -0.460 -0.1500 0.0630 1.02   169
## A[5,3] -0.660 -0.2200 0.0970 1.03   168
## A[5,4] -0.041  0.2200 0.6900 1.02   173
## A[5,5]  0.530  0.7500 0.9900 1.01   233
## 
## Process error parameter estimates:
##             2.5%  50% 97.5% Rhat n_eff
## Sigma[1,1] 0.031 0.27  0.64 1.03    92
## Sigma[1,2] 0.000 0.00  0.00  NaN   NaN
## Sigma[1,3] 0.000 0.00  0.00  NaN   NaN
## Sigma[1,4] 0.000 0.00  0.00  NaN   NaN
## Sigma[1,5] 0.000 0.00  0.00  NaN   NaN
## Sigma[2,1] 0.000 0.00  0.00  NaN   NaN
## Sigma[2,2] 0.066 0.11  0.17 1.01   462
## Sigma[2,3] 0.000 0.00  0.00  NaN   NaN
## Sigma[2,4] 0.000 0.00  0.00  NaN   NaN
## Sigma[2,5] 0.000 0.00  0.00  NaN   NaN
## Sigma[3,1] 0.000 0.00  0.00  NaN   NaN
## Sigma[3,2] 0.000 0.00  0.00  NaN   NaN
## Sigma[3,3] 0.052 0.15  0.29 1.03    99
## Sigma[3,4] 0.000 0.00  0.00  NaN   NaN
## Sigma[3,5] 0.000 0.00  0.00  NaN   NaN
## Sigma[4,1] 0.000 0.00  0.00  NaN   NaN
## Sigma[4,2] 0.000 0.00  0.00  NaN   NaN
## Sigma[4,3] 0.000 0.00  0.00  NaN   NaN
## Sigma[4,4] 0.048 0.12  0.25 1.02   150
## Sigma[4,5] 0.000 0.00  0.00  NaN   NaN
## Sigma[5,1] 0.000 0.00  0.00  NaN   NaN
## Sigma[5,2] 0.000 0.00  0.00  NaN   NaN
## Sigma[5,3] 0.000 0.00  0.00  NaN   NaN
## Sigma[5,4] 0.000 0.00  0.00  NaN   NaN
## Sigma[5,5] 0.095 0.20  0.35 1.01   194
## 
## Approximate significance of GAM process smooths:
##                              edf Ref.df    F p-value  
## te(temp,month)              3.43     15 3.26   0.329  
## te(temp,month):seriestrend1 1.48     15 0.10   1.000  
## te(temp,month):seriestrend2 1.16     15 0.51   0.998  
## te(temp,month):seriestrend3 3.00     15 4.56   0.038 *
## te(temp,month):seriestrend4 2.71     15 0.53   0.957  
## te(temp,month):seriestrend5 2.87     15 0.21   0.998  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Stan MCMC diagnostics:
## n_eff / iter looks reasonable for all parameters
## Rhats above 1.05 found for 1 parameters
##  *Diagnose further to investigate why the chains have not mixed
## 0 of 2000 iterations ended with a divergence (0%)
## 0 of 2000 iterations saturated the maximum tree depth of 12 (0%)
## E-FMI indicated no pathological behavior
## 
## Samples were drawn using NUTS(diag_e) at Wed May 15 2:20:09 PM 2024.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split MCMC chains
## (at convergence, Rhat = 1)

We can again plot the smooth functions, which this time operate on the process model. The coefficients for this model are now accessible through the trend_mgcv_model slot in the model object:

gratia::draw(var_mod$trend_mgcv_model, select = 1)

We can see the same plot using the mvgam version by using trend_effects = TRUE in the plotting functions:

plot(var_mod, 'smooths', trend_effects = TRUE)

The VAR matrix is of particular interest here, as it captures lagged dependencies and cross-dependencies in the latent process model:

mcmc_plot(var_mod, variable = 'A', regex = TRUE, type = 'hist')

Unfortunately bayesplot doesn’t know this is a matrix of parameters so what we see is actually the transpose of the VAR matrix. A little bit of wrangling gives us these histograms in the correct order:

A_pars <- matrix(NA, nrow = 5, ncol = 5)
for(i in 1:5){
  for(j in 1:5){
    A_pars[i, j] <- paste0('A[', i, ',', j, ']')
  }
}
mcmc_plot(var_mod, 
          variable = as.vector(t(A_pars)), 
          type = 'hist')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

There is a lot happening in this matrix. Each cell captures the lagged effect of the process in the column on the process in the row in the next timestep. So for example, the effect in cell [5,2], which is quite strongly negative, means that an increase in the process for series 5 (Unicells) at time \(t\) is expected to lead to a subsequent decrease in the process for series 2 (Diatoms) at time \(t+1\). The latent process model is now capturing these effects and the smooth seasonal effects, so the trend plot shows our best estimate of what the true count should have been at each time point:

plot(var_mod, type = 'trend', series = 1)

plot(var_mod, type = 'trend', series = 3)

The process error \((\Sigma)\) captures unmodelled variation in the process models. Again, we fixed the off-diagonals to 0, so the histograms for these will look like flat boxes:

Sigma_pars <- matrix(NA, nrow = 5, ncol = 5)
for(i in 1:5){
  for(j in 1:5){
    Sigma_pars[i, j] <- paste0('Sigma[', i, ',', j, ']')
  }
}
mcmc_plot(var_mod, 
          variable = as.vector(t(Sigma_pars)), 
          type = 'hist')

The observation error estimates \((\sigma_{obs})\) represent how much the model thinks we might miss the true count when we take our imperfect measurements:

mcmc_plot(var_mod, variable = 'sigma_obs', regex = TRUE, type = 'hist')

Session Info

sessionInfo()
## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 22621)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_Australia.utf8  LC_CTYPE=English_Australia.utf8   
## [3] LC_MONETARY=English_Australia.utf8 LC_NUMERIC=C                      
## [5] LC_TIME=English_Australia.utf8    
## 
## time zone: Australia/Brisbane
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ggplot2_3.5.0               gratia_0.8.9.11            
##  [3] bayesplot_1.10.0            tidybayes_3.0.6            
##  [5] mvgam_1.1.1                 insight_0.19.7             
##  [7] marginaleffects_0.17.0.9002 brms_2.21.0                
##  [9] Rcpp_1.0.12                 mgcv_1.8-42                
## [11] nlme_3.1-162                dplyr_1.1.4                
## [13] downloadthis_0.3.3         
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1     svUnit_1.0.6         farver_2.1.1        
##  [4] scoringRules_1.1.1   loo_2.6.0            fastmap_1.1.1       
##  [7] tensorA_0.36.2.1     digest_0.6.34        timechange_0.2.0    
## [10] mime_0.12            lifecycle_1.0.4      StanHeaders_2.26.28 
## [13] magrittr_2.0.3       posterior_1.5.0      compiler_4.3.1      
## [16] rlang_1.1.3          sass_0.4.8           tools_4.3.1         
## [19] utf8_1.2.4           yaml_2.3.8           data.table_1.14.10  
## [22] knitr_1.45           labeling_0.4.3       bridgesampling_1.1-2
## [25] pkgbuild_1.4.3       plyr_1.8.9           RColorBrewer_1.1-3  
## [28] abind_1.4-5          klippy_0.0.0.9500    withr_3.0.0         
## [31] purrr_1.0.2          grid_4.3.1           stats4_4.3.1        
## [34] fansi_1.0.6          colorspace_2.1-0     inline_0.3.19       
## [37] scales_1.3.0         MASS_7.3-60          isoband_0.2.7       
## [40] cli_3.6.2            mvtnorm_1.2-4        crayon_1.5.2        
## [43] rmarkdown_2.25       generics_0.1.3       RcppParallel_5.1.7  
## [46] rstudioapi_0.15.0    reshape2_1.4.4       cachem_1.0.8        
## [49] rstan_2.32.3         stringr_1.5.1        splines_4.3.1       
## [52] assertthat_0.2.1     parallel_4.3.1       matrixStats_1.2.0   
## [55] base64enc_0.1-3      vctrs_0.6.5          Matrix_1.6-4        
## [58] jsonlite_1.8.8       patchwork_1.2.0      arrayhelpers_1.1-0  
## [61] ggdist_3.3.1         jquerylib_0.1.4      tidyr_1.3.1         
## [64] glue_1.7.0           ggokabeito_0.1.0     codetools_0.2-19    
## [67] mvnfast_0.2.8        distributional_0.3.2 lubridate_1.9.3     
## [70] stringi_1.8.3        gtable_0.3.4         QuickJSR_1.0.9      
## [73] munsell_0.5.0        tibble_3.2.1         pillar_1.9.0        
## [76] htmltools_0.5.7      Brobdingnag_1.2-9    R6_2.5.1            
## [79] evaluate_0.23        lattice_0.21-8       highr_0.10          
## [82] backports_1.4.1      bslib_0.6.1          rstantools_2.3.1.1  
## [85] bsplus_0.1.4         zip_2.3.0            coda_0.19-4.1       
## [88] gridExtra_2.3        checkmate_2.3.1      xfun_0.42           
## [91] fs_1.6.3             pkgconfig_2.0.3

  1. , https://github.com/nicholasjclark↩︎