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 `trend`

s). 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}