Shared latent states in mvgam

Checking trend_map with run_model = FALSE

Supplying this trend_map to the mvgam function for a simple model, but setting run_model = FALSE, allows us to inspect the constructed Stan code and the data objects that would be used to condition the model. Here we will set up a model in which each series has a different observation process (with only a different intercept per series in this case), and the two latent dynamic process models evolve as independent AR1 processes that also contain a shared nonlinear smooth function to capture repeated seasonality. This model is not too complicated but it does show how we can learn shared and independent effects for collections of time series in the mvgam framework:

fake_mod <- mvgam(y ~ 
                    # observation model formula, which has a 
                    # different intercept per series
                    series - 1,
                  
                  # process model formula, which has a shared seasonal smooth
                  # (each latent process model shares the SAME smooth)
                  trend_formula = ~ s(season, bs = 'cc', k = 6),
                  
                  # AR1 dynamics (each latent process model has DIFFERENT)
                  # dynamics; processes are estimated using the noncentred
                  # parameterisation for improved efficiency
                  trend_model = AR(),
                  noncentred = TRUE,
                  
                  # supplied trend_map
                  trend_map = trend_map,
                  
                  # data and observation family
                  family = poisson(),
                  data = simdat$data_train,
                  run_model = FALSE)

Inspecting the Stan code shows how this model is a dynamic factor model in which the loadings are constructed to reflect the supplied trend_map:

code(fake_mod)

Notice the line that states “lv_coefs = Z;”. This uses the supplied Z matrix to construct the loading coefficients. The supplied matrix now looks exactly like what you’d use if you were to create a similar model in the MARSS package:

fake_mod$model_data$Z

Fitting and inspecting the model

Though this model doesn’t perfectly match the data-generating process (which allowed each series to have different underlying dynamics), we can still fit it to show what the resulting inferences look like:

full_mod <- mvgam(y ~ series - 1,
                  trend_formula = ~ s(season, bs = 'cc', k = 6),
                  trend_model = AR(),
                  noncentred = TRUE,
                  trend_map = trend_map,
                  family = poisson(),
                  data = simdat$data_train,
                  silent = 2)

The summary of this model is informative as it shows that only two latent process models have been estimated, even though we have three observed time series. The model converges well

summary(full_mod)

Both series 1 and 2 share the exact same latent process estimates, while the estimates for series 3 are different:

plot(full_mod, type = 'trend', series = 1)
plot(full_mod, type = 'trend', series = 2)
plot(full_mod, type = 'trend', series = 3)

However, forecasts for series’ 1 and 2 will differ because they have different intercepts in the observation model

The shared signal model

Now we can formulate and fit a model that allows each sensor’s observation error to depend nonlinearly on temperature while allowing the true signal to depend nonlinearly on productivity. By fixing all of the values in the trend column to 1 in the trend_map, we are assuming that all observation sensors are tracking the same latent signal. We use informative priors on the two variance components (process error and observation error), which reflect our prior belief that the observation error is smaller overall than the true process error

mod <- mvgam(formula =
               # formula for observations, allowing for different
               # intercepts and hierarchical smooth effects of temperature
               observed ~ series + 
               s(temperature, k = 10) +
               s(series, temperature, bs = 'sz', k = 8),
             
             trend_formula =
               # formula for the latent signal, which can depend
               # nonlinearly on productivity
               ~ s(productivity, k = 8) - 1,
             
             trend_model =
               # in addition to productivity effects, the signal is
               # assumed to exhibit temporal autocorrelation
               AR(),
             noncentred = TRUE,
             
             trend_map =
               # trend_map forces all sensors to track the same
               # latent signal
               data.frame(series = unique(model_dat$series),
                          trend = c(1, 1, 1)),
             
             # informative priors on process error
             # and observation error will help with convergence
             priors = c(prior(normal(2, 0.5), class = sigma),
                        prior(normal(1, 0.5), class = sigma_obs)),
             
             # Gaussian observations
             family = gaussian(),
             data = model_dat,
             silent = 2)

View a reduced version of the model summary because there will be many spline coefficients in this model

summary(mod, include_betas = FALSE)

Inspecting effects on both process and observation models

Don’t pay much attention to the approximate p-values of the smooth terms. The calculation for these values is incredibly sensitive to the estimates for the smoothing parameters so I don’t tend to find them to be very meaningful. What are meaningful, however, are prediction-based plots of the smooth functions. All main effects can be quickly plotted with conditional_effects:

conditional_effects(mod, type = 'link')

conditional_effects is simply a wrapper to the more flexible plot_predictions function from the marginaleffects package. We can get more useful plots of these effects using this function for further customisation:

require(marginaleffects)
plot_predictions(mod, 
                 condition = c('temperature', 'series', 'series'),
                 points = 0.5) +
  theme(legend.position = 'none')

We have successfully estimated effects, some of them nonlinear, that impact the hidden process AND the observations. All in a single joint model. But there can always be challenges with these models, particularly when estimating both process and observation error at the same time.

Recovering the hidden signal

A final but very key question is whether we can successfully recover the true hidden signal. The trend slot in the returned model parameters has the estimates for this signal, which we can easily plot using the mvgam S3 method for plot. We can also overlay the true values for the hidden signal, which shows that our model has done a good job of recovering it:

plot(mod, type = 'trend')

# Overlay the true simulated signal
points(true_signal, pch = 16, cex = 1, col = 'white')
points(true_signal, pch = 16, cex = 0.8)