Formating data for use with mvgam - defining series

Interfaces
1

This feels like a very naive question probably due to my lack of experience with time series data. @nicholasjclark, thank you for mvgam and its well-documented use-cases. I am having issues defining the ‘series’ and ‘time’ variables for my specific data structure. I have data on Diptera captures on sticky traps along agricultural field margins. Some fields were sampled only some years of the study but as I understand, this is not an issue. For each given year of sampling for a given field, anywhere between 1-6 sticky traps were used. The precise location of each trap along the field margin is not known such that the trap ID within a given field is not representative of anything meaningful. Thus, for a given time point within a field, I have multiple measures of Diptera abundance. My issue arises when trying to define the ‘series’ variable. Biologically, it should be field ID but when I try to run get_mvgam_priors() on a model with field ID as the ‘series’ and multiple distinct ‘time’ values for a given series, I get the message Error: One or more series in data is missing observations for one or more timepoints which I don’t get when trying to fit a model with distinct ‘time’ values per series.

Reproducible example:

expand.grid(series = factor(1:3), time = 1:10) %>%
  mutate(y = rpois(n(), 1)) -> dat

get_mvgam_priors(y ~ 1 , family = poisson, data = dat)
# gives the expected output

dat %>%
expand_grid(trap_ID = 1:5) -> dat2

get_mvgam_priors(y ~ 1 , family = poisson, data = dat2)
# Error: One or more series in data is missing observations for one or more timepoints

I would not have any conceptual issue with defining trap_ID as a ‘series’ IF its location along a field margin would have been constant through time but it is not. Also, some years had a different number of traps for a given field.

I guess my question is: does it make sense to define trap_ID as a ‘series’ in this specific use case? I would thus specify the hierarchical structure of the data by using field_ID as random intercepts/slopes.

Thank you for your time.

2

Thanks for the question @allenbushbeaupre. My suggestion would be to look through the shared latent states vignette and work out how you might need to map the observations to the shared latent states. mvgam() can set up and fit State-Space models where multiple observations are taken of the same latent, unobserved dynamic processes. This is advantageous as some of your observations may have unique biases that you might like to account for in the observation model formula (i.e. one of your traps may have malfunctioned, for example) but you can also include covariates that you think impact the “true” latent abundances (by placing covariates in the trend_formula. The key is to supply a trend_map object that contains the instructions for mapping from the observations (in the series column) to the latent states (the unique trends). In the below example, I’ve assumed trap_ID indicates the top level of the hierarchy, but you may need to use something like field instead.

library(mvgam)
library(tidyverse)

expand.grid(series = factor(1:3), time = 1:10) %>%
  mutate(y = rpois(n(), 1)) -> dat

get_mvgam_priors(y ~ 1 , family = poisson, data = dat)
#>    param_name param_length  param_info                               prior
#> 1 (Intercept)            1 (Intercept) (Intercept) ~ student_t(3, 0, 2.5);
#>                example_change new_lowerbound new_upperbound
#> 1 (Intercept) ~ normal(0, 1);             NA             NA
# gives the expected output

dat %>%
  expand_grid(trap_ID = 1:5) -> dat2

# Each series must only have one observation per timepoint; but you can set up
# a State-Space model that allows multiple series to contain observations of the 
# same latent dynamic process; here I assume that trap_ID contains the mapping
dat2 %>%
  mutate(orig_series = series,
         series = as.factor(paste0(series, '_', trap_ID))) -> dat3

head(dat3, 16)
#> # A tibble: 16 × 5
#>    series  time     y trap_ID orig_series
#>    <fct>  <int> <int>   <int> <fct>      
#>  1 1_1        1     2       1 1          
#>  2 1_2        1     2       2 1          
#>  3 1_3        1     2       3 1          
#>  4 1_4        1     2       4 1          
#>  5 1_5        1     2       5 1          
#>  6 2_1        1     0       1 2          
#>  7 2_2        1     0       2 2          
#>  8 2_3        1     0       3 2          
#>  9 2_4        1     0       4 2          
#> 10 2_5        1     0       5 2          
#> 11 3_1        1     1       1 3          
#> 12 3_2        1     1       2 3          
#> 13 3_3        1     1       3 3          
#> 14 3_4        1     1       4 3          
#> 15 3_5        1     1       5 3          
#> 16 1_1        2     0       1 1

# Add a trend_map object that contains the mapping from observations to latent
# temporal dynamics (one latent temporal series per trap_ID)
dat3 %>%
  select(series, orig_series) %>%
  distinct() %>%
  mutate(trend = as.numeric(orig_series)) %>%
  select(-orig_series) -> trend_map

trend_map
#> # A tibble: 15 × 2
#>    series trend
#>    <fct>  <dbl>
#>  1 1_1        1
#>  2 1_2        1
#>  3 1_3        1
#>  4 1_4        1
#>  5 1_5        1
#>  6 2_1        2
#>  7 2_2        2
#>  8 2_3        2
#>  9 2_4        2
#> 10 2_5        2
#> 11 3_1        3
#> 12 3_2        3
#> 13 3_3        3
#> 14 3_4        3
#> 15 3_5        3

# Now confirm that the setup works, here assuming a random walk for each
# of the five latent time series
get_mvgam_priors(y ~ 1,
                 trend_formula = ~ 1,
                 trend_map = trend_map,
                 trend_model = RW(),
                 family = poisson(), 
                 data = dat3)
#>                     param_name param_length                param_info
#> 1                  (Intercept)            1               (Intercept)
#> 2 vector<lower=0>[n_lv] sigma;            3          process error sd
#> 3            (Intercept)_trend            1 (Intercept) for the trend
#>                                     prior                    example_change
#> 1     (Intercept) ~ student_t(3, 0, 2.5);       (Intercept) ~ normal(0, 1);
#> 2           sigma ~ student_t(3, 0, 2.5);        sigma ~ exponential(0.66);
#> 3 (Intercept)_trend ~ student_t(3, 0, 2); (Intercept)_trend ~ normal(0, 1);
#>   new_lowerbound new_upperbound
#> 1           <NA>           <NA>
#> 2           <NA>           <NA>
#> 3           <NA>           <NA>

# Fit the model and inspect the latent trend estimates
mod <- mvgam(y ~ 1,
             trend_formula = ~ 1,
             trend_map = trend_map,
             trend_model = RW(),
             family = poisson(), 
             data = dat3)
#> Your model may benefit from using "noncentred = TRUE"
#> Compiling Stan program using cmdstanr
#> 
#> Start sampling
#> Running MCMC with 4 parallel chains...
#> 
#> Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
#> Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
#> Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
#> Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
#> Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
#> Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
#> Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
#> Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
#> Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
#> Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
#> Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
#> Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
#> Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
#> Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
#> Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
#> Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
#> Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
#> Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
#> Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
#> Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
#> Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
#> Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
#> Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
#> Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
#> Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
#> Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
#> Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
#> Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
#> Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
#> Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
#> Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
#> Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
#> Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
#> Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
#> Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
#> Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
#> Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
#> Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
#> Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
#> Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
#> Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
#> Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
#> Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
#> Chain 1 finished in 2.4 seconds.
#> Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
#> Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
#> Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
#> Chain 4 finished in 2.7 seconds.
#> Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
#> Chain 3 finished in 2.9 seconds.
#> Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
#> Chain 2 finished in 3.1 seconds.
#> 
#> All 4 chains finished successfully.
#> Mean chain execution time: 2.8 seconds.
#> Total execution time: 3.3 seconds.

# Observations that share the same trap_ID will have identical
# latent dynamics
plot(mod, type = 'trend', series = 1)


plot(mod, type = 'trend', series = 2)



plot(mod, type = 'trend', series = 6)


plot(mod, type = 'trend', series = 7)

Created on 2024-08-27 with reprex v2.0.2

1 Like
3

Thank you @nicholasjclark for the detailed answer. This seems to be the solution I was looking for. I will spend some time reading up on shared latent states. Great package by the way, thank you again.

1 Like