Example 1: a two-species system with nonlinear trends
First we will use a simple simulation in which multiple replicate
observations are taken at each timepoint for two different species. The
simulation produces observations at a single site over six years, with
five replicate surveys per year. Each species is simulated to have
different nonlinear temporal trends and different detection
probabilities. For now, detection probability is fixed (i.e. it does not
change over time or in association with any covariates). Notice that we
add the cap variable, which does not need to be static, to
define the maximum possible value that we think the latent abundance
could be for each timepoint. This simply needs to be large enough that
we get a reasonable idea of which latent N values are most likely,
without adding too much computational cost:
set.seed(999)
# Simulate observations for species 1, which shows a declining trend and 0.7 detection probability
data.frame(site = 1,
# five replicates per year; six years
replicate = rep(1:5, 6),
time = sort(rep(1:6, 5)),
species = 'sp_1',
# true abundance declines nonlinearly
truth = c(rep(28, 5),
rep(26, 5),
rep(23, 5),
rep(16, 5),
rep(14, 5),
rep(14, 5)),
# observations are taken with detection prob = 0.7
obs = c(rbinom(5, 28, 0.7),
rbinom(5, 26, 0.7),
rbinom(5, 23, 0.7),
rbinom(5, 15, 0.7),
rbinom(5, 14, 0.7),
rbinom(5, 14, 0.7))) %>%
# add 'series' information, which is an identifier of site, replicate and species
dplyr::mutate(series = paste0('site_', site,
'_', species,
'_rep_', replicate),
time = as.numeric(time),
# add a 'cap' variable that defines the maximum latent N to
# marginalize over when estimating latent abundance; in other words
# how large do we realistically think the true abundance could be?
cap = 100) %>%
dplyr::select(- replicate) -> testdat
# Now add another species that has a different temporal trend and a smaller
# detection probability (0.45 for this species)
testdat = testdat %>%
dplyr::bind_rows(data.frame(site = 1,
replicate = rep(1:5, 6),
time = sort(rep(1:6, 5)),
species = 'sp_2',
truth = c(rep(4, 5),
rep(7, 5),
rep(15, 5),
rep(16, 5),
rep(19, 5),
rep(18, 5)),
obs = c(rbinom(5, 4, 0.45),
rbinom(5, 7, 0.45),
rbinom(5, 15, 0.45),
rbinom(5, 16, 0.45),
rbinom(5, 19, 0.45),
rbinom(5, 18, 0.45))) %>%
dplyr::mutate(series = paste0('site_', site,
'_', species,
'_rep_', replicate),
time = as.numeric(time),
cap = 50) %>%
dplyr::select(-replicate))This data format isn’t too difficult to set up, but it does differ from the traditional multidimensional array setup that is commonly used for fitting N-mixture models in other software packages. Next we ensure that species and series IDs are included as factor variables, in case we’d like to allow certain effects to vary by species
testdat$species <- factor(testdat$species,
levels = unique(testdat$species))
testdat$series <- factor(testdat$series,
levels = unique(testdat$series))Preview the dataset to get an idea of how it is structured:
dplyr::glimpse(testdat)
#> Rows: 60
#> Columns: 7
#> $ site <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
#> $ time <dbl> 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5,…
#> $ species <fct> sp_1, sp_1, sp_1, sp_1, sp_1, sp_1, sp_1, sp_1, sp_1, sp_1, sp…
#> $ truth <dbl> 28, 28, 28, 28, 28, 26, 26, 26, 26, 26, 23, 23, 23, 23, 23, 16…
#> $ obs <int> 20, 19, 23, 17, 18, 21, 18, 21, 19, 18, 17, 16, 20, 11, 19, 9,…
#> $ series <fct> site_1_sp_1_rep_1, site_1_sp_1_rep_2, site_1_sp_1_rep_3, site_…
#> $ cap <dbl> 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 10…
head(testdat, 12)
#> site time species truth obs series cap
#> 1 1 1 sp_1 28 20 site_1_sp_1_rep_1 100
#> 2 1 1 sp_1 28 19 site_1_sp_1_rep_2 100
#> 3 1 1 sp_1 28 23 site_1_sp_1_rep_3 100
#> 4 1 1 sp_1 28 17 site_1_sp_1_rep_4 100
#> 5 1 1 sp_1 28 18 site_1_sp_1_rep_5 100
#> 6 1 2 sp_1 26 21 site_1_sp_1_rep_1 100
#> 7 1 2 sp_1 26 18 site_1_sp_1_rep_2 100
#> 8 1 2 sp_1 26 21 site_1_sp_1_rep_3 100
#> 9 1 2 sp_1 26 19 site_1_sp_1_rep_4 100
#> 10 1 2 sp_1 26 18 site_1_sp_1_rep_5 100
#> 11 1 3 sp_1 23 17 site_1_sp_1_rep_1 100
#> 12 1 3 sp_1 23 16 site_1_sp_1_rep_2 100Setting up the trend_map
Finally, we need to set up the trend_map object. This is
crucial for allowing multiple observations to be linked to the same
latent process model (see more information about this argument in the Shared latent states vignette). In this case, the
mapping operates by species and site to state that each set of replicate
observations from the same time point should all share the exact same
latent abundance model:
testdat %>%
# each unique combination of site*species is a separate process
dplyr::mutate(trend = as.numeric(factor(paste0(site, species)))) %>%
dplyr::select(trend, series) %>%
dplyr::distinct() -> trend_map
trend_map
#> trend series
#> 1 1 site_1_sp_1_rep_1
#> 2 1 site_1_sp_1_rep_2
#> 3 1 site_1_sp_1_rep_3
#> 4 1 site_1_sp_1_rep_4
#> 5 1 site_1_sp_1_rep_5
#> 6 2 site_1_sp_2_rep_1
#> 7 2 site_1_sp_2_rep_2
#> 8 2 site_1_sp_2_rep_3
#> 9 2 site_1_sp_2_rep_4
#> 10 2 site_1_sp_2_rep_5Notice how all of the replicates for species 1 in site 1 share the
same process (i.e. the same trend). This will ensure that
all replicates are Binomial draws of the same latent N.
Modelling with the nmix() family
Now we are ready to fit a model using mvgam(). This
model will allow each species to have different detection probabilities
and different temporal trends. We will use Cmdstan as the
backend, which by default will use Hamiltonian Monte Carlo for full
Bayesian inference
mod <- mvgam(
# the observation formula sets up linear predictors for
# detection probability on the logit scale
formula = obs ~ species - 1,
# the trend_formula sets up the linear predictors for
# the latent abundance processes on the log scale
trend_formula = ~ s(time, by = trend, k = 4) + species,
# the trend_map takes care of the mapping
trend_map = trend_map,
# nmix() family and data
family = nmix(),
data = testdat,
# priors can be set in the usual way
priors = c(prior(std_normal(), class = b),
prior(normal(1, 1.5), class = Intercept_trend)),
samples = 1000)View the automatically-generated Stan code to get a
sense of how the marginalization over latent N works
code(mod)
#> // Stan model code generated by package mvgam
#> data {
#> int<lower=0> total_obs; // total number of observations
#> int<lower=0> n; // number of timepoints per series
#> int<lower=0> n_sp_trend; // number of trend smoothing parameters
#> int<lower=0> n_lv; // number of dynamic factors
#> int<lower=0> n_series; // number of series
#> matrix[n_series, n_lv] Z; // matrix mapping series to latent states
#> int<lower=0> num_basis; // total number of basis coefficients
#> int<lower=0> num_basis_trend; // number of trend basis coefficients
#> vector[num_basis_trend] zero_trend; // prior locations for trend basis coefficients
#> matrix[total_obs, num_basis] X; // mgcv GAM design matrix
#> matrix[n * n_lv, num_basis_trend] X_trend; // trend model design matrix
#> array[n, n_series] int<lower=0> ytimes; // time-ordered matrix (which col in X belongs to each [time, series] observation?)
#> array[n, n_lv] int ytimes_trend;
#> int<lower=0> n_nonmissing; // number of nonmissing observations
#> array[total_obs] int<lower=0> cap; // upper limits of latent abundances
#> array[total_obs] int ytimes_array; // sorted ytimes
#> array[n, n_series] int<lower=0> ytimes_pred; // time-ordered matrix for prediction
#> int<lower=0> K_groups; // number of unique replicated observations
#> int<lower=0> K_reps; // maximum number of replicate observations
#> array[K_groups] int<lower=0> K_starts; // col of K_inds where each group starts
#> array[K_groups] int<lower=0> K_stops; // col of K_inds where each group ends
#> array[K_groups, K_reps] int<lower=0> K_inds; // indices of replicated observations
#> matrix[3, 6] S_trend1; // mgcv smooth penalty matrix S_trend1
#> matrix[3, 6] S_trend2; // mgcv smooth penalty matrix S_trend2
#> array[total_obs] int<lower=0> flat_ys; // flattened observations
#> }
#> transformed data {
#> matrix[total_obs, num_basis] X_ordered = X[ytimes_array, : ];
#> array[K_groups] int<lower=0> Y_max;
#> array[K_groups] int<lower=0> N_max;
#> for (k in 1 : K_groups) {
#> Y_max[k] = max(flat_ys[K_inds[k, K_starts[k] : K_stops[k]]]);
#> N_max[k] = max(cap[K_inds[k, K_starts[k] : K_stops[k]]]);
#> }
#> }
#> parameters {
#> // raw basis coefficients
#> vector[num_basis] b_raw;
#> vector[num_basis_trend] b_raw_trend;
#>
#> // smoothing parameters
#> vector<lower=0>[n_sp_trend] lambda_trend;
#> }
#> transformed parameters {
#> // detection probability
#> vector[total_obs] p;
#>
#> // latent states
#> matrix[n, n_lv] LV;
#>
#> // latent states and loading matrix
#> vector[n * n_lv] trend_mus;
#> matrix[n, n_series] trend;
#>
#> // basis coefficients
#> vector[num_basis] b;
#> vector[num_basis_trend] b_trend;
#>
#> // observation model basis coefficients
#> b[1 : num_basis] = b_raw[1 : num_basis];
#>
#> // process model basis coefficients
#> b_trend[1 : num_basis_trend] = b_raw_trend[1 : num_basis_trend];
#>
#> // detection probability
#> p = X_ordered * b;
#>
#> // latent process linear predictors
#> trend_mus = X_trend * b_trend;
#> for (j in 1 : n_lv) {
#> LV[1 : n, j] = trend_mus[ytimes_trend[1 : n, j]];
#> }
#>
#> // derived latent states
#> for (i in 1 : n) {
#> for (s in 1 : n_series) {
#> trend[i, s] = dot_product(Z[s, : ], LV[i, : ]);
#> }
#> }
#> }
#> model {
#> // prior for speciessp_1...
#> b_raw[1] ~ std_normal();
#>
#> // prior for speciessp_2...
#> b_raw[2] ~ std_normal();
#>
#> // dynamic process models
#>
#> // prior for (Intercept)_trend...
#> b_raw_trend[1] ~ normal(1, 1.5);
#>
#> // prior for speciessp_2_trend...
#> b_raw_trend[2] ~ std_normal();
#>
#> // prior for s(time):trendtrend1_trend...
#> b_raw_trend[3 : 5] ~ multi_normal_prec(zero_trend[3 : 5],
#> S_trend1[1 : 3, 1 : 3]
#> * lambda_trend[1]
#> + S_trend1[1 : 3, 4 : 6]
#> * lambda_trend[2]);
#>
#> // prior for s(time):trendtrend2_trend...
#> b_raw_trend[6 : 8] ~ multi_normal_prec(zero_trend[6 : 8],
#> S_trend2[1 : 3, 1 : 3]
#> * lambda_trend[3]
#> + S_trend2[1 : 3, 4 : 6]
#> * lambda_trend[4]);
#> lambda_trend ~ normal(5, 30);
#> {
#> // likelihood functions
#> array[total_obs] real flat_trends;
#> array[total_obs] real flat_ps;
#> flat_trends = to_array_1d(trend);
#> flat_ps = to_array_1d(p);
#>
#> // loop over replicate sampling window (each site*time*species combination)
#> for (k in 1 : K_groups) {
#> // all log_lambdas are identical because they represent site*time
#> // covariates; so just use the first measurement
#> real log_lambda = flat_trends[K_inds[k, 1]];
#> vector[N_max[k] - Y_max[k] + 1] terms;
#> int l = 0;
#>
#> // marginalize over latent abundance
#> for (Ni in Y_max[k] : N_max[k]) {
#> l = l + 1;
#> // factor for poisson prob of latent Ni; compute
#>
#> // only once per sampling window
#> terms[l] = poisson_log_lpmf(Ni | log_lambda)
#> + // for each replicate observation, binomial prob observed is
#> // computed in a vectorized statement
#> binomial_logit_lpmf(flat_ys[K_inds[k, K_starts[k] : K_stops[k]]] | Ni, flat_ps[K_inds[k, K_starts[k] : K_stops[k]]]);
#> }
#> target += log_sum_exp(terms);
#> }
#> }
#> }
#> generated quantities {
#> vector[n_lv] penalty = rep_vector(1e12, n_lv);
#> vector[n_sp_trend] rho_trend = log(lambda_trend);
#> }The posterior summary of this model shows that it has converged nicely
summary(mod)
#> GAM observation formula:
#> obs ~ species - 1
#> <environment: 0x00000237248751d0>
#>
#> GAM process formula:
#> ~s(time, by = trend, k = 4) + species
#> <environment: 0x00000237248751d0>
#>
#> Family:
#> nmix
#>
#> Link function:
#> log
#>
#> Trend model:
#> None
#>
#> N process models:
#> 2
#>
#> N series:
#> 10
#>
#> N timepoints:
#> 6
#>
#> Status:
#> Fitted using Stan
#> 4 chains, each with iter = 1500; warmup = 500; thin = 1
#> Total post-warmup draws = 4000
#>
#>
#> GAM observation model coefficient (beta) estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> speciessp_1 -0.28 0.720 1.40 1 1630
#> speciessp_2 -1.20 0.035 0.88 1 1913
#>
#> GAM process model coefficient (beta) estimates:
#> 2.5% 50% 97.5% Rhat n_eff
#> (Intercept)_trend 2.700 3.0000 3.500 1.00 1453
#> speciessp_2_trend -1.200 -0.6300 0.140 1.00 1587
#> s(time):trendtrend1.1_trend -0.078 0.0160 0.210 1.01 931
#> s(time):trendtrend1.2_trend -0.250 0.0058 0.270 1.00 2384
#> s(time):trendtrend1.3_trend -0.450 -0.2500 -0.040 1.00 1747
#> s(time):trendtrend2.1_trend -0.210 -0.0120 0.092 1.00 849
#> s(time):trendtrend2.2_trend -0.190 0.0290 0.530 1.00 627
#> s(time):trendtrend2.3_trend 0.052 0.3300 0.640 1.00 2547
#>
#> Approximate significance of GAM process smooths:
#> edf Ref.df F p-value
#> s(time):seriestrend1 1.16 3 1.04 0.64
#> s(time):seriestrend2 1.12 3 1.58 0.62
#>
#> Stan MCMC diagnostics:
#> n_eff / iter looks reasonable for all parameters
#> Rhat looks reasonable for all parameters
#> 0 of 4000 iterations ended with a divergence (0%)
#> 0 of 4000 iterations saturated the maximum tree depth of 12 (0%)
#> E-FMI indicated no pathological behavior
#>
#> Samples were drawn using NUTS(diag_e) at Thu May 09 9:46:13 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)loo() functionality works just as it does for all
mvgam models to aid in model comparison / selection (though
note that Pareto K values often give warnings for mixture models so
these may not be too helpful)
loo(mod)
#> Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
#>
#> Computed from 4000 by 60 log-likelihood matrix
#>
#> Estimate SE
#> elpd_loo -452.5 77.5
#> p_loo 254.0 64.5
#> looic 905.0 155.0
#> ------
#> Monte Carlo SE of elpd_loo is NA.
#>
#> Pareto k diagnostic values:
#> Count Pct. Min. n_eff
#> (-Inf, 0.5] (good) 29 48.3% 823
#> (0.5, 0.7] (ok) 1 1.7% 1302
#> (0.7, 1] (bad) 9 15.0% 49
#> (1, Inf) (very bad) 21 35.0% 1
#> See help('pareto-k-diagnostic') for details.Plot the estimated smooths of time from each species’ latent abundance process (on the log scale)
marginaleffects support allows for more useful
prediction-based interrogations on different scales (though note that at
the time of writing this Vignette, you must have the development version
of marginaleffects installed for nmix() models
to be supported; use
remotes::install_github('vincentarelbundock/marginaleffects')
to install). Objects that use family nmix() have a few
additional prediction scales that can be used (i.e. link,
response, detection or latent_N).
For example, here are the estimated detection probabilities per species,
which show that the model has done a nice job of estimating these
parameters:
plot_predictions(mod, condition = 'species',
type = 'detection') +
ylab('Pr(detection)') +
ylim(c(0, 1)) +
theme_classic() +
theme(legend.position = 'none')
A common goal in N-mixture modelling is to estimate the true latent abundance. The model has automatically generated predictions for the unknown latent abundance that are conditional on the observations. We can extract these and produce decent plots using a small function
hc <- hindcast(mod, type = 'latent_N')
# Function to plot latent abundance estimates vs truth
plot_latentN = function(hindcasts, data, species = 'sp_1'){
all_series <- unique(data %>%
dplyr::filter(species == !!species) %>%
dplyr::pull(series))
# Grab the first replicate that represents this series
# so we can get the true simulated values
series <- as.numeric(all_series[1])
truths <- data %>%
dplyr::arrange(time, series) %>%
dplyr::filter(series == !!levels(data$series)[series]) %>%
dplyr::pull(truth)
# In case some replicates have missing observations,
# pull out predictions for ALL replicates and average over them
hcs <- do.call(rbind, lapply(all_series, function(x){
ind <- which(names(hindcasts$hindcasts) %in% as.character(x))
hindcasts$hindcasts[[ind]]
}))
# Calculate posterior empirical quantiles of predictions
pred_quantiles <- data.frame(t(apply(hcs, 2, function(x)
quantile(x, probs = c(0.05, 0.2, 0.3, 0.4,
0.5, 0.6, 0.7, 0.8, 0.95)))))
pred_quantiles$time <- 1:NROW(pred_quantiles)
pred_quantiles$truth <- truths
# Grab observations
data %>%
dplyr::filter(series %in% all_series) %>%
dplyr::select(time, obs) -> observations
# Plot
ggplot(pred_quantiles, aes(x = time, group = 1)) +
geom_ribbon(aes(ymin = X5., ymax = X95.), fill = "#DCBCBC") +
geom_ribbon(aes(ymin = X30., ymax = X70.), fill = "#B97C7C") +
geom_line(aes(x = time, y = truth),
colour = 'black', linewidth = 1) +
geom_point(aes(x = time, y = truth),
shape = 21, colour = 'white', fill = 'black',
size = 2.5) +
geom_jitter(data = observations, aes(x = time, y = obs),
width = 0.06,
shape = 21, fill = 'darkred', colour = 'white', size = 2.5) +
labs(y = 'Latent abundance (N)',
x = 'Time',
title = species)
}Latent abundance plots vs the simulated truths for each species are shown below. Here, the red points show the imperfect observations, the black line shows the true latent abundance, and the ribbons show credible intervals of our estimates:
We can see that estimates for both species have correctly captured the true temporal variation and magnitudes in abundance
