jsdgam: Fit Joint Species Distribution Models in 'mvgam'

In mvgam: Multivariate (Dynamic) Generalized Additive Models

Fit Joint Species Distribution Models in mvgam

Description

This function sets up a Joint Species Distribution Model whereby the residual associations among species can be modelled in a reduced-rank format using a set of latent factors. The factor specification is extremely flexible, allowing users to include spatial, temporal or any other type of predictor effects to more efficiently capture unmodelled residual associations, while the observation model can also be highly flexible (including all smooth, GP and other effects that mvgam can handle)

Usage

jsdgam(
  formula,
  factor_formula = ~-1,
  knots,
  factor_knots,
  data,
  newdata,
  family = poisson(),
  unit = time,
  species = series,
  share_obs_params = FALSE,
  priors,
  n_lv = 2,
  backend = getOption("brms.backend", "cmdstanr"),
  algorithm = getOption("brms.algorithm", "sampling"),
  control = list(max_treedepth = 10, adapt_delta = 0.8),
  chains = 4,
  burnin = 500,
  samples = 500,
  thin = 1,
  parallel = TRUE,
  threads = 1,
  silent = 1,
  run_model = TRUE,
  return_model_data = FALSE,
  residuals = TRUE,
  ...
)

Arguments

formula	A formula object specifying the GAM observation model formula. These are exactly like the formula for a GLM except that smooth terms, s(), te(), ti(), t2(), as well as time-varying dynamic() terms, nonparametric gp() terms and offsets using offset(), can be added to the right hand side to specify that the linear predictor depends on smooth functions of predictors (or linear functionals of these). Details of the formula syntax used by mvgam can be found in mvgam_formulae
factor_formula	A formula object specifying the linear predictor effects for the latent factors. Use by = trend within calls to functional terms (i.e. s(), te(), ti(), t2(), dynamic(), or gp()) to ensure that each factor captures a different axis of variation. See the example below as an illustration
knots	An optional list containing user specified knot values for basis construction. For most bases the user simply supplies the knots to be used, which must match up with the k value supplied. Different terms can use different numbers of knots, unless they share a covariate.
factor_knots	An optional list containing user specified knot values to be used for basis construction of any smooth terms in factor_formula. For most bases the user simply supplies the knots to be used, which must match up with the k value supplied (note that the number of knots is not always just k). Different terms can use different numbers of knots, unless they share a covariate
data	A dataframe or list containing the model response variable and covariates required by the GAM formula and factor_formula objects
newdata	Optional dataframe or list of test data containing the same variables as in data. If included, observations in variable y will be set to NA when fitting the model so that posterior simulations can be obtained.
family	family specifying the observation family for the outcomes. Currently supported families are: gaussian() for real-valued data betar() for proportional data on ⁠(0,1)⁠ lognormal() for non-negative real-valued data student_t() for real-valued data Gamma() for non-negative real-valued data bernoulli() for binary data poisson() for count data nb() for overdispersed count data binomial() for count data with imperfect detection when the number of trials is known; note that the cbind() function must be used to bind the discrete observations and the discrete number of trials beta_binomial() as for binomial() but allows for overdispersion Default is poisson(). See mvgam_families for more details
unit	The unquoted name of the variable that represents the unit of analysis in data over which latent residuals should be correlated. This variable should be either a numeric or integer variable in the supplied data. Defaults to time to be consistent with other functionalities in mvgam, though note that the data need not be time series in this case. See examples below for further details and explanations
species	The unquoted name of the factor variable that indexes the different response units in data (usually 'species' in a JSDM). Defaults to series to be consistent with other mvgam models
share_obs_params	logical. If TRUE and the family has additional family-specific observation parameters (e.g., variance components, dispersion parameters), these will be shared across all outcome variables. Useful when multiple outcomes share properties. Default is FALSE.
priors	An optional data.frame with prior definitions (in Stan syntax) or, preferentially, a vector containing objects of class brmsprior (see. prior for details). See get_mvgam_priors and for more information on changing default prior distributions
n_lv	integer the number of latent factors to use for modelling residual associations. Cannot be ⁠> n_species⁠. Defaults arbitrarily to 2
backend	Character string naming the package for Stan model fitting. Options are "cmdstanr" (default) or "rstan". Can be set globally via "brms.backend" option. See https://mc-stan.org/rstan/ and https://mc-stan.org/cmdstanr/ for details.
algorithm	Character string naming the estimation approach: "sampling": MCMC (default) "meanfield": Variational inference with factorized normal distributions "fullrank": Variational inference with multivariate normal distribution "laplace": Laplace approximation (cmdstanr only) "pathfinder": Pathfinder algorithm (cmdstanr only) Can be set globally via "brms.algorithm" option. Limited testing suggests "meanfield" performs best among non-MCMC approximations for dynamic GAMs.
control	Named list for controlling sampler behaviour. Valid elements include max_treedepth, adapt_delta and init.
chains	integer specifying the number of parallel chains for the model. Ignored for variational inference algorithms.
burnin	integer specifying the number of warmup iterations to tune sampling algorithms. Ignored for variational inference algorithms.
samples	integer specifying the number of post-warmup iterations for sampling the posterior distribution.
thin	Thinning interval for monitors. Ignored for variational inference algorithms.
parallel	logical specifying whether to use multiple cores for parallel MCMC simulation. If TRUE, uses min(c(chains, parallel::detectCores() - 1)) cores.
threads	integer Experimental option to use multithreading for within-chain parallelisation in Stan. We recommend its use only if you are experienced with Stan's reduce_sum function and have a slow running model that cannot be sped up by any other means. Currently works for all families when using cmdstanr as the backend
silent	Verbosity level between 0 and 2. If 1 (default), most informational messages are suppressed. If 2, even more messages are suppressed. Sampling progress is still printed - set refresh = 0 to disable. For backend = "rstan", also set open_progress = FALSE to prevent additional progress bars.
run_model	logical. If FALSE, the model is not fitted but instead the function returns the model file and the data/initial values needed to fit the model outside of mvgam.
return_model_data	logical. If TRUE, the list of data needed to fit the model is returned, along with initial values for smooth and AR parameters, once the model is fitted. Helpful for users who wish to modify the model file to add other stochastic elements. Default is FALSE unless run_model == FALSE.
residuals	logical. Whether to compute series-level randomized quantile residuals. Default is TRUE. Set to FALSE to save time and reduce object size (can add later using add_residuals).
...	Other arguments to pass to mvgam

Details

Joint Species Distribution Models allow for responses of multiple species to be learned hierarchically, whereby responses to environmental variables in formula can be partially pooled and any latent, unmodelled residual associations can also be learned. In mvgam, both of these effects can be modelled with the full power of latent factor Hierarchical GAMs, providing unmatched flexibility to model full communities of species. When calling jsdgam, an initial State-Space model using trend = 'None' is set up and then modified to include the latent factors and their linear predictors. Consequently, you can inspect priors for these models using get_mvgam_priors by supplying the relevant formula, factor_formula, data and family arguments and keeping the default trend = 'None'.

In a JSDGAM, the expectation of response Y_{ij} is modelled with

g(\mu_{ij}) = X_i\beta + u_i\theta_j,

where g(.) is a known link function, X is a design matrix of linear predictors (with associated \beta coefficients), u are n_{lv}-variate latent factors (n_{lv}<<n_{species}) and \theta_j are species-specific loadings on the latent factors, respectively. The design matrix X and \beta coefficients are constructed and modelled using formula and can contain any of mvgam's predictor effects, including random intercepts and slopes, multidimensional penalized smooths, GP effects etc... The factor loadings \theta_j are constrained for identifiability but can be used to reconstruct an estimate of the species' residual variance-covariance matrix using \Theta \Theta' (see the example below and residual_cor() for details). The latent factors are further modelled using:

u_i \sim \text{Normal}(Q_i\beta_{factor}, 1)

where the second design matrix Q and associated \beta_{factor} coefficients are constructed and modelled using factor_formula. Again, the effects that make up this linear predictor can contain any of mvgam's allowed predictor effects, providing enormous flexibility for modelling species' communities.

Value

A list object of class mvgam containing model output, the text representation of the model file, the mgcv model output (for easily generating simulations at unsampled covariate values), Dunn-Smyth residuals for each species and key information needed for other functions in the package. See mvgam-class for details. Use methods(class = "mvgam") for an overview on available methods

Author(s)

Nicholas J Clark

References

Nicholas J Clark & Konstans Wells (2023). Dynamic generalised additive models (DGAMs) for forecasting discrete ecological time series. Methods in Ecology and Evolution. 14:3, 771-784.

David I Warton, F Guillaume Blanchet, Robert B O'Hara, Otso Ovaskainen, Sara Taskinen, Steven C Walker & Francis KC Hui (2015). So many variables: joint modeling in community ecology. Trends in Ecology & Evolution 30:12, 766-779.

See Also

mvgam(), residual_cor()

Examples

## Not run: 
# ========================================================================
# Example 1: Basic JSDGAM with Portal Data
# ========================================================================

# Fit a JSDGAM to the portal_data captures
mod <- jsdgam(
  formula = captures ~
    # Fixed effects of NDVI and mintemp, row effect as a GP of time
    ndvi_ma12:series + mintemp:series + gp(time, k = 15),
  factor_formula = ~ -1,
  data = portal_data,
  unit = time,
  species = series,
  family = poisson(),
  n_lv = 2,
  silent = 2,
  chains = 2
)

# Plot covariate effects
library(ggplot2); theme_set(theme_bw())
plot_predictions(
  mod,
  condition = c('ndvi_ma12', 'series', 'series')
)

plot_predictions(
  mod,
  condition = c('mintemp', 'series', 'series')
)

# A residual correlation plot
plot(residual_cor(mod))

# An ordination biplot can also be constructed
# from the factor scores and their loadings
if(requireNamespace('ggrepel', quietly = TRUE)){
  ordinate(mod, alpha = 0.7)
}

# ========================================================================
# Example 2: Advanced JSDGAM with Spatial Predictors
# ========================================================================

# Simulate latent count data for 500 spatial locations and 10 species
set.seed(0)
N_points <- 500
N_species <- 10

# Species-level intercepts (on the log scale)
alphas <- runif(N_species, 2, 2.25)

# Simulate a covariate and species-level responses to it
temperature <- rnorm(N_points)
betas <- runif(N_species, -0.5, 0.5)

# Simulate points uniformly over a space
lon <- runif(N_points, min = 150, max = 155)
lat <- runif(N_points, min = -20, max = -19)

# Set up spatial basis functions as a tensor product of lat and lon
sm <- mgcv::smoothCon(
  mgcv::te(lon, lat, k = 5),
  data = data.frame(lon, lat),
  knots = NULL
)[[1]]

# The design matrix for this smooth is in the 'X' slot
des_mat <- sm$X
dim(des_mat)

# Function to generate a random covariance matrix where all variables
# have unit variance (i.e. diagonals are all 1)
random_Sigma = function(N){
  L_Omega <- matrix(0, N, N);
  L_Omega[1, 1] <- 1;
  for (i in 2 : N) {
    bound <- 1;
    for (j in 1 : (i - 1)) {
      L_Omega[i, j] <- runif(1, -sqrt(bound), sqrt(bound));
      bound <- bound - L_Omega[i, j] ^ 2;
    }
    L_Omega[i, i] <- sqrt(bound);
  }
  Sigma <- L_Omega %*% t(L_Omega);
  return(Sigma)
}

# Simulate a variance-covariance matrix for the correlations among
# basis coefficients
Sigma <- random_Sigma(N = NCOL(des_mat))

# Now simulate the species-level basis coefficients hierarchically, where
# spatial basis function correlations are a convex sum of a base correlation
# matrix and a species-level correlation matrix
basis_coefs <- matrix(NA, nrow = N_species, ncol = NCOL(Sigma))
base_field <- mgcv::rmvn(1, mu = rep(0, NCOL(Sigma)), V = Sigma)
for(t in 1:N_species){
  corOmega <- (cov2cor(Sigma) * 0.7) +
    (0.3 * cov2cor(random_Sigma(N = NCOL(des_mat))))
  basis_coefs[t, ] <- mgcv::rmvn(1, mu = rep(0, NCOL(Sigma)), V = corOmega)
}

# Simulate the latent spatial processes
st_process <- do.call(rbind, lapply(seq_len(N_species), function(t){
  data.frame(
    lat = lat,
    lon = lon,
    species = paste0('species_', t),
    temperature = temperature,
    process = alphas[t] +
      betas[t] * temperature +
      des_mat %*% basis_coefs[t,]
  )
}))

# Now take noisy observations at some of the points (60)
obs_points <- sample(1:N_points, size = 60, replace = FALSE)
obs_points <- data.frame(
  lat = lat[obs_points],
  lon = lon[obs_points],
  site = 1:60
)

# Keep only the process data at these points
st_process %>%
  dplyr::inner_join(obs_points, by = c('lat', 'lon')) %>%
  # now take noisy Poisson observations of the process
  dplyr::mutate(count = rpois(NROW(.), lambda = exp(process))) %>%
  dplyr::mutate(species = factor(
    species,
    levels = paste0('species_', 1:N_species)
  )) %>%
  dplyr::group_by(lat, lon) -> dat

# View the count distributions for each species
ggplot(dat, aes(x = count)) +
  geom_histogram() +
  facet_wrap(~ species, scales = 'free')

ggplot(dat, aes(x = lon, y = lat, col = log(count + 1))) +
  geom_point(size = 2.25) +
  facet_wrap(~ species, scales = 'free') +
  scale_color_viridis_c()

# ------------------------------------------------------------------------
# Model Fitting with Custom Priors
# ------------------------------------------------------------------------

# Inspect default priors for a joint species model with three spatial factors
priors <- get_mvgam_priors(
  formula = count ~
    # Environmental model includes random slopes for
    # a linear effect of temperature
    s(species, bs = 're', by = temperature),

  # Each factor estimates a different nonlinear spatial process, using
  # 'by = trend' as in other mvgam State-Space models
  factor_formula = ~ gp(lon, lat, k = 6, by = trend) - 1,
  n_lv = 3,

  # The data and grouping variables
  data = dat,
  unit = site,
  species = species,

  # Poisson observations
  family = poisson()
)
head(priors)

# Fit a JSDM that estimates hierarchical temperature responses
# and that uses three latent spatial factors
mod <- jsdgam(
  formula = count ~
    # Environmental model includes random slopes for a
    # linear effect of temperature
    s(species, bs = 're', by = temperature),

  # Each factor estimates a different nonlinear spatial process, using
  # 'by = trend' as in other mvgam State-Space models
  factor_formula = ~ gp(lon, lat, k = 6, by = trend) - 1,
  n_lv = 3,

  # Change default priors for fixed random effect variances and
  # factor GP marginal deviations to standard normal
  priors = c(
    prior(std_normal(), class = sigma_raw),
    prior(std_normal(), class = `alpha_gp_trend(lon, lat):trendtrend1`),
    prior(std_normal(), class = `alpha_gp_trend(lon, lat):trendtrend2`),
    prior(std_normal(), class = `alpha_gp_trend(lon, lat):trendtrend3`)
  ),

  # The data and the grouping variables
  data = dat,
  unit = site,
  species = species,

  # Poisson observations
  family = poisson(),
  chains = 2,
  silent = 2
)

# ------------------------------------------------------------------------
# Model Visualization and Diagnostics
# ------------------------------------------------------------------------

# Plot the implicit species-level intercept estimates
plot_predictions(mod, condition = 'species', type = 'link')

# Plot species' hierarchical responses to temperature
plot_predictions(
  mod,
  condition = c('temperature', 'species', 'species'),
  type = 'link'
)

# Plot posterior median estimates of the latent spatial factors
plot(mod, type = 'smooths', trend_effects = TRUE)

# Or using gratia, if you have it installed
if(requireNamespace('gratia', quietly = TRUE)){
  gratia::draw(mod, trend_effects = TRUE, dist = 0)
}

# Plot species' randomized quantile residual distributions
# as a function of latitude
pp_check(
  mod,
  type = 'resid_ribbon_grouped',
  group = 'species',
  x = 'lat',
  ndraws = 200
)

# ------------------------------------------------------------------------
# Residual Correlation Analysis
# ------------------------------------------------------------------------

# Calculate residual spatial correlations
post_cors <- residual_cor(mod)
names(post_cors)

# Look at lower and upper credible interval estimates for
# some of the estimated correlations
post_cors$cor[1:5, 1:5]
post_cors$cor_upper[1:5, 1:5]
post_cors$cor_lower[1:5, 1:5]

# Plot of the posterior median correlations for those estimated
# to be non-zero
plot(post_cors, cluster = TRUE)

# An ordination biplot can also be constructed
# from the factor scores and their loadings
if(requireNamespace('ggrepel', quietly = TRUE)){
  ordinate(mod)
}

# ------------------------------------------------------------------------
# Model Validation and Prediction
# ------------------------------------------------------------------------

# Posterior predictive checks and ELPD-LOO can ascertain model fit
pp_check(
  mod,
  type = "pit_ecdf_grouped",
  group = "species",
  ndraws = 200
)
loo(mod)

# Forecast log(counts) for entire region (site value doesn't matter as long
# as each spatial location has a different and unique site identifier);
# note this calculation takes a few minutes because of the need to calculate
# draws from the stochastic latent factors
newdata <- st_process %>%
  dplyr::mutate(species = factor(
    species,
    levels = paste0('species_', 1:N_species)
  )) %>%
  dplyr::group_by(lat, lon) %>%
  dplyr::mutate(site = dplyr::cur_group_id()) %>%
  dplyr::ungroup()
preds <- predict(mod, newdata = newdata)

# Plot the median log(count) predictions on a grid
newdata$log_count <- preds[,1]
ggplot(newdata, aes(x = lon, y = lat, col = log_count)) +
  geom_point(size = 1.5) +
  facet_wrap(~ species, scales = 'free') +
  scale_color_viridis_c() +
  theme_classic()

# Not needed for general use; cleans up connections for automated testing
closeAllConnections()

## End(Not run)

mvgam documentation built on Jan. 21, 2026, 9:07 a.m.