Fitting Joint Species Distribution Models with jsdmstan#

Notebook Repository Method Repository alt text

Primary Contact: Dr. Fiona Seaton

Challenge:

Modelling ecological communities involves understanding not just how individual species respond to their environment, but also how they interact with one another. Species often co-occur or exhibit similar patterns of response to environmental conditions, meaning their distributions are not independent. Traditional species distribution models, which treat each species in isolation, miss this shared structure. The challenge is to develop models that can capture these interdependencies, enabling more accurate and ecologically realistic predictions by leveraging information across multiple species simultaneously.

Approach:

Joint Species Distribution Models (JSDMs) address this challenge by modelling all species in a community together. They incorporate environmental covariates, species-specific effects (like individual intercepts), and—critically—covariance between species to account for shared variation and co-occurrence patterns. This covariance can be specified either through a full multivariate model that directly estimates the species covariance matrix, or more efficiently via latent variable models (GLLVMs), which reduce dimensionality by capturing shared structure with a small number of latent factors.

The jsdmstan R package is an interface for fitting joint species distribution models (JSDMs) in Stan. First, we will go through an introduction to joint species distribution models in general, followed by examples of how to use the jsdmstan package to simulate data and fit JSDMs to this simulated data, and finally two examples of using jsdmstan on real ecological data.

Introduction to Joint Species Distribution Models#

Joint Species Distibution Models are models that model an entire community of species simultaneously. The idea behind these is that they allow information to be borrowed across species, such that the covariance between species can be used to inform the predictions of distributions of related or commonly co-occurring species.

In plain language (or as plain as I can manage) JSDMs involve the modelling of an entire species community as a function of some combination of intercepts, covariate data and species covariance. Therefore the change of a single species is related to not only change in the environment but also how it relates to other species. There are several decisions to be made in how to specify these models - the standard decisions on which covariates to include, whether each species should have its own intercept (generally yes) and how to represent change across sites - but also how to represent the covariance between species. There are two options for representing this species covariance in this package. First, the original way of running JSDMs was to model the entire covariance matrix between species in a multivariate generalised linear mixed model (MGLMM). However, more recently there have been methods developed that involve representing the covariance matrix with a set of linear latent variables - known as generalised linear latent variable models (GLLVM).

The jsdmstan package aims to provide an interface for fitting these models in Stan using the Stan Hamiltonian Monte Carlo sampling as a robust Bayesian methodology.

Underlying maths#

Feel free to skip this bit if you don’t want to read equations, it is largely based on Warton et al. (2015). We model the community data \(m_{ij}\) for each site \(i\) and taxon \(j\) as a function of a species intercept, environmental covariates and species covariance matrix:

\[g(m_{ij}) = \beta_{0j} + \mathbf{x}_i^\intercal\beta_j + u_{ij}\]

where \(g(\cdot)\) is the link function, \(\mathbf{x}_i^\intercal\) is the transpose of vector \(\mathbf{x}\), and for each taxon \(j\), \(\beta_{0j}\) is an intercept and \(beta_j\) is a vector of regression coefficients related to measured predictors.

A site effect \(\alpha_{i}\) can also be added to adjust for total abundance or richness:

\[g(m_{ij}) = \alpha_{i} + \beta_{0j} + \mathbf{x}_i^\intercal\beta_j + u_{ij}\]

Multivariate Generalised Linear Mixed Models#

The entire matrix of covariance between species is modelled in MGLMMs.

\[u_{ij} \sim \mathrm{N}(\mathbf{0},\mathbf{\Sigma})\]

Fitting the entire covariance matrix means that the amount of time required to fit these models scales with the number of species cubed, and the data required scales with the number of species squared. This makes these models both computationally and data intensive.

Generalised Linear Latent Variable Models#

In response to some of these issues in fitting MGLMMs, GLLVMs were developed in which \(u_{ij}\) is now specified as a linear function of a set of latent variables \(\mathbf{z_i}\):

\[y_{ij}|\mathbf{u}_i \sim \mathrm{F}(m_{ij},\phi_j)\]
\[u_{ij} = \mathbf{z}_i^\intercal \lambda_j\]

The latent variables \(\mathbf{z_i}\) are treated as random by assuming:

\[y_{ij}|\mathbf{z_i} \sim \mathrm{F}(m_{ij},\phi_j)\]
\[\mathbf{z_i} \sim \mathrm{N}(\mathbf{0},\mathbf{1})\]

Treating the species covariance as pulling from a set of latent variables greatly reduces the computational time required to fit these models.

Relationship to environmental covariates#

Within jsdmstan the response of species to environmental covariates can either be unstructured (the default) or constrained by a covariance matrix between the environmental covariates. This second option (specified by setting beta_param = "cor") assumes that if one species is strongly positively related to multiple covariates then it is more likely that other species will either also be positively related to all these covariates, or negatively related. Mathematically this corresponds to:

\[\beta_j \sim \mathrm{N}(\mathbf{0},\mathbf{\Sigma})\]

Examples with Simulated Data#

First, we’ll load in the R packages we’ll be using here. The jsdmstan package can be installed from github using the remotes package.

# remotes::install_github("NERC-CEH/jsdmstan")
library(jsdmstan)
library(dplyr)
library(ggplot2)
library(patchwork)
set.seed(3757892)

Fitting a MGLMM#

We can use the in-built functions for simulating data according to the MGLMM model - we’ll choose to simulate 15 species over 200 sites with 2 environmental covariates. The species are assumed to follow a Poisson distribution (with a log-link), and we use the defaults of including a species-specific intercept but no site-specific intercept. At the moment only default priors (standard normal distribution) are supported. We can do this using either the jsdm_sim_data() function with method = "mglmm" or with the mglmm_sim_data() function which just calls jsdm_sim_data() in the background.

nsites <- 75
nspecies <- 8
ncovar <- 2
mglmm_test_data <- mglmm_sim_data(N = nsites, S = nspecies, 
                                  K = ncovar, family = "pois")

This returns a list, which includes the Y matrix, the X matrix, plus also the exact parameters used to create the data:

names(mglmm_test_data)
[1] "Y"    "pars" "N"    "S"    "D"    "K"    "X"   

Now, to fit the model we can use the stan_jsdm() function, which interfaces to Stan through the rstan package. There are multiple ways to supply data to the stan_jsdm() function, one is to supply the data as a list with the appropriate named components (the jsdm_sim_data() functions supply data in the correct format already), the second way is to specify the Y and X matrices directly, and the third way is to use a formula for the environmental covariates and supply the environmental data to the data argument, which is what we’ll use here:

mglmm_fit <- stan_jsdm(~ V1 + V2, data = dat, Y = mglmm_test_data$Y, 
                       family = "pois", method = "mglmm", refresh = 0)

If we print the model object we will get a brief overview of the type of jSDM and the data, plus if there are any parameters with Rhat > 1.01 or effective sample size ratio (Neff/N) < 0.05 then they will be printed:

mglmm_fit
Family: poisson 
 Model type: mglmm
  Number of species: 8
  Number of sites: 75
  Number of predictors: 3

Model run on 4 chains with 4000 iterations per chain (2000 warmup).

No parameters with Rhat > 1.01 or Neff/N < 0.05

To get a summary of all the model parameters we can use summary(), there are many parameters in these models so we just include a few here:

summary(mglmm_fit, pars = "cor_species")
                   mean    sd    15%    85%  Rhat Bulk.ESS Tail.ESS
cor_species[2,1] -0.258 0.261 -0.534  0.019 1.002     1589     2574
cor_species[3,1]  0.552 0.152  0.390  0.713 1.001     2483     4145
cor_species[4,1] -0.268 0.152 -0.424 -0.111 1.002     2736     4083
cor_species[5,1]  0.450 0.190  0.251  0.651 1.001     2326     3703
cor_species[6,1] -0.096 0.164 -0.264  0.076 1.001     2650     4743
cor_species[7,1] -0.348 0.115 -0.470 -0.228 1.002     2600     4817
cor_species[8,1]  0.046 0.114 -0.072  0.162 1.000     2448     4063
cor_species[1,2] -0.258 0.261 -0.534  0.019 1.002     1589     2574
cor_species[2,2]  1.000 0.000  1.000  1.000 1.000     7833       NA
cor_species[3,2] -0.255 0.259 -0.525  0.015 1.002     1984     3235
cor_species[4,2]  0.206 0.253 -0.057  0.474 1.002     1541     2820
cor_species[5,2] -0.241 0.267 -0.520  0.041 1.002     2687     3983
cor_species[6,2]  0.207 0.265 -0.077  0.483 1.001     1641     2579
cor_species[7,2]  0.190 0.255 -0.083  0.459 1.001      957     2053
cor_species[8,2]  0.050 0.274 -0.242  0.346 1.001      831     1558
cor_species[1,3]  0.552 0.152  0.390  0.713 1.001     2483     4145
cor_species[2,3] -0.255 0.259 -0.525  0.015 1.002     1984     3235
cor_species[3,3]  1.000 0.000  1.000  1.000 1.000     8019     7569
cor_species[4,3] -0.378 0.153 -0.537 -0.219 1.000     3941     5536
cor_species[5,3]  0.380 0.199  0.170  0.592 1.003     2514     4534
cor_species[6,3] -0.194 0.176 -0.377 -0.012 1.002     2675     4889
cor_species[7,3] -0.222 0.131 -0.357 -0.088 1.001     3087     4432
cor_species[8,3] -0.121 0.133 -0.261  0.016 1.004     2096     4288
cor_species[1,4] -0.268 0.152 -0.424 -0.111 1.002     2736     4083
cor_species[2,4]  0.206 0.253 -0.057  0.474 1.002     1541     2820
cor_species[3,4] -0.378 0.153 -0.537 -0.219 1.000     3941     5536
cor_species[4,4]  1.000 0.000  1.000  1.000 1.000     8152       NA
cor_species[5,4] -0.172 0.180 -0.357  0.019 1.001     3464     5486
cor_species[6,4]  0.557 0.162  0.384  0.728 1.001     4349     5516
cor_species[7,4] -0.276 0.145 -0.425 -0.126 1.002     2578     3579
cor_species[8,4]  0.403 0.144  0.254  0.555 1.004     1868     3774
cor_species[1,5]  0.450 0.190  0.251  0.651 1.001     2326     3703
cor_species[2,5] -0.241 0.267 -0.520  0.041 1.002     2687     3983
cor_species[3,5]  0.380 0.199  0.170  0.592 1.003     2514     4534
cor_species[4,5] -0.172 0.180 -0.357  0.019 1.001     3464     5486
cor_species[5,5]  1.000 0.000  1.000  1.000 1.000     7881     7281
cor_species[6,5] -0.169 0.208 -0.382  0.049 1.001     3431     5608
cor_species[7,5] -0.380 0.193 -0.581 -0.171 1.001     1486     3017
cor_species[8,5] -0.065 0.177 -0.248  0.119 1.004     1543     2996
cor_species[1,6] -0.096 0.164 -0.264  0.076 1.001     2650     4743
cor_species[2,6]  0.207 0.265 -0.077  0.483 1.001     1641     2579
cor_species[3,6] -0.194 0.176 -0.377 -0.012 1.002     2675     4889
cor_species[4,6]  0.557 0.162  0.384  0.728 1.001     4349     5516
cor_species[5,6] -0.169 0.208 -0.382  0.049 1.001     3431     5608
cor_species[6,6]  1.000 0.000  1.000  1.000 1.000     8015       NA
cor_species[7,6] -0.139 0.151 -0.297  0.020 1.000     2866     4823
cor_species[8,6]  0.475 0.137  0.334  0.620 1.002     2467     4808
cor_species[1,7] -0.348 0.115 -0.470 -0.228 1.002     2600     4817
cor_species[2,7]  0.190 0.255 -0.083  0.459 1.001      957     2053
cor_species[3,7] -0.222 0.131 -0.357 -0.088 1.001     3087     4432
cor_species[4,7] -0.276 0.145 -0.425 -0.126 1.002     2578     3579
cor_species[5,7] -0.380 0.193 -0.581 -0.171 1.001     1486     3017
cor_species[6,7] -0.139 0.151 -0.297  0.020 1.000     2866     4823
cor_species[7,7]  1.000 0.000  1.000  1.000 1.000     7748       NA
cor_species[8,7] -0.157 0.109 -0.272 -0.041 1.002     2728     4274
cor_species[1,8]  0.046 0.114 -0.072  0.162 1.000     2448     4063
cor_species[2,8]  0.050 0.274 -0.242  0.346 1.001      831     1558
cor_species[3,8] -0.121 0.133 -0.261  0.016 1.004     2096     4288
cor_species[4,8]  0.403 0.144  0.254  0.555 1.004     1868     3774
cor_species[5,8] -0.065 0.177 -0.248  0.119 1.004     1543     2996
cor_species[6,8]  0.475 0.137  0.334  0.620 1.002     2467     4808
cor_species[7,8] -0.157 0.109 -0.272 -0.041 1.002     2728     4274
cor_species[8,8]  1.000 0.000  1.000  1.000 1.000     7797     7780

To get a better overview of the R-hat and effective sample size we can use the mcmc_plot() function to plot histograms of R-hat and ESS.

mcmc_plot(mglmm_fit, plotfun = "rhat_hist")
../../_images/7877078a88e749377894f29726364ce7b676a545ae9a63491a250f8185e0af4c.png
mcmc_plot(mglmm_fit, plotfun = "neff_hist")
../../_images/42a0e9a81746561cbb231f3d9d759d42d8e41b2cf6f58c9f0c83d8c4096f8d21.png

We can also examine the output for each parameter visually using a traceplot combined with a density plot, which is given by the default plot() command:

plot(mglmm_fit, ask = FALSE)
../../_images/ac97dd96cf9762b12f9ecb9b436dfeaba19f49176353b60d3b542d6fea8e6a2f.png ../../_images/9876d3ac854b0e554e9dd5f72a793ae3ca29e7469317c3d8a53a8da4296c0f0c.png ../../_images/128e433b152bf8c4c4fc5d94f5d349064d456a591097122bb3e99b115d9a4a78.png ../../_images/98076fbe2e38b1013c283b537ef2db44a57f902ff05c0bb50fa9fdd600041693.png

By default the plot() command plots all of the parameters with sigma or kappa in their name plus a random selection of 20 other parameters, but this can be overridden by either specifying the parameters by name (with or without regular expression matching) or changing the number of parameters to be randomly sampled. Use the get_parnames() function to get the names of parameters within a model - and the jsdm_stancode() function can also be used to see the underlying structure of the model.

All the mcmc plot types within bayesplot are supported by the mcmc_plot() function, and to see a full list either use bayesplot::available_mcmc() or run mcmc_plot() with an incorrect type and the options will be printed.

We can also view the environmental effect parameters for each species using the envplot() function.

envplot(mglmm_fit)
../../_images/b4bd9899ddc1201714896ef894c09b1d5a1bb3ee6befa5868b613affd157f712.png

Posterior predictions can be extracted from the models using either posterior_linpred() or posterior_predict(), where the linpred function extracts the linear predictor for the community composition within each draw and the predict function combines this linear predictor extraction with a random generation based on the predicted probability for the family. Both functions by default return a list of length equal to the number of draws extracted, where each element of the list is a sites by species matrix.

mglmm_pp <- posterior_predict(mglmm_fit)
length(mglmm_pp)
[1] 8000
[1] 75  8

As well as the MCMC plotting functions within bayesplot the ppc_ family of functions is also supported through the pp_check() function. This family of functions provides a graphical way to check your posterior against the data used within the model to evaluate model fit - called a posterior retrodictive check (or posterior predictive historically and when the prior only has been sampled from). To use these you need to have set save_data = TRUE within the stan_jsdm() call. Unlike in other packages by default pp_check() for jsdmStanFit objects extracts the posterior predictions then calculates summary statistics over the rows and plots those summary statistics against the same for the original data. The default behaviour is to calculate the sum of all the species per site - i.e. total abundance.

pp_check(mglmm_fit)
../../_images/1c9f45707ccb400057cf624682bcac49d85202ccc530887c4ca73b32bb850681.png

The summary statistic can be changed, as can whether it is calculated for every species or every site:

pp_check(mglmm_fit, summary_stat = "mean", calc_over = "species",
         plotfun = "ecdf_overlay")
../../_images/836690c92431c3abadde4b991bf9a7ceae11ae3afda7861874f0742fe8f138d1.png

We can examine the species-specific posterior predictive check through using multi_pp_check(), or examine how well the relationships between specific species are recovered using pp_check() with plotfun = "pairs".

As we have run the above model on simulated data and the original data list contains the parameters used to simulate the data we can use the mcmc_recover_ functions from bayesplot to see how the model did:

mcmc_plot(mglmm_fit, plotfun = "recover_hist",
          pars = paste0("sigmas_species[",1:8,"]"),
          true = mglmm_test_data$pars$sigmas_species)
../../_images/e6e6aa4b92375bd19e3e15cd8f344496388bfba2456255b076b5f4e0abcaadaa.png
mcmc_plot(mglmm_fit, plotfun = "recover_intervals",
          pars = paste0("cor_species[",rep(1:nspecies, nspecies:1),",",
                        unlist(sapply(1:8, ":",8)),"]"),
          true = c(mglmm_test_data$pars$cor_species[lower.tri(mglmm_test_data$pars$cor_species, diag = TRUE)])) +
  theme(axis.text.x = element_text(angle = 90))
../../_images/34f67d3fb7a4fe9450ff46474c3ccdc03f81c89ea2ec88719865224e53ae92e3.png

Fitting a GLLVM#

The model fitting workflow for latent variable models is very similar to that above, with the addition of specifying the number of latent variables (D) in the data simulation and model fit. Here we change the family to a Bernoulli family (i.e. the special case of the binomial where the number of trials is 1 for all observations), make the covariate effects on each species draw from a correlation matrix such that information can be shared across species, and change the prior to be a Student’s T prior on the predictor-specific sigma parameter.

set.seed(3562251)
gllvm_data <- gllvm_sim_data(N = 50, S = 12, D = 2, K = 1,
                             family = "bernoulli",
                             beta_param = "cor",
                             prior = jsdm_prior(sigmas_preds = "student_t(3,0,1)"))
gllvm_fit <- stan_jsdm(Y = gllvm_data$Y, X = gllvm_data$X,
                       D = gllvm_data$D,  
                       family = "bernoulli",
                       method = "gllvm", 
                       beta_param = "cor",
                       prior = jsdm_prior(sigmas_preds = "student_t(3,0,1)"),
                       refresh = 0)
gllvm_fit
Family: bernoulli 
 Model type: gllvm with 2 latent variables
  Number of species: 12
  Number of sites: 50
  Number of predictors: 2

Model run on 4 chains with 4000 iterations per chain (2000 warmup).

No parameters with Rhat > 1.01 or Neff/N < 0.05

Again, the diagnostic statistics seem reasonable:

mcmc_plot(gllvm_fit, plotfun = "rhat_hist")
../../_images/74c5998477bcb35c25eae190e57ef14db6ab425cdc45b6655e08fa5e0f5c8add.png
mcmc_plot(gllvm_fit, plotfun = "neff_hist")
../../_images/58cdbd3678c168c5b0a9dfdd75f743d0160f1d5f40a7d86bc3d17380b8a2edc1.png

For brevity’s sake we will not go into the detail of the different functions again here, however there is one plotting function specifically for GLLVM models - ordiplot(). This plots the species or sites scores against the latent variables from a random selection of draws:

ordiplot(gllvm_fit, errorbar_range = 0.5)
../../_images/e3c264835eeb270da2fc878bc05115874837b70722fec0e4421d01c8d1296e02.png
ordiplot(gllvm_fit, type = "sites", geom = "text", errorbar_range = 0) +
  theme(legend.position = "none")
../../_images/000703572f29e1c0afc5b4c0ddf4e19762d388920689fb0d2a0b0328e6a393b5.png

You can change the latent variables selected by specifying the choices argument, and alter the number of draws or whether you want to plot species or sites with the other arguments.

Examples with Real Data#

MGLMM with broadleaved woodland survey data#

We read in the bunce71 data from the jsdmstan package, giving a response matrix of 45 tree genera over 103 sites and an environmental information matrix that contains site location, number of plots per site as well as five predictor variables - annual rainfall, maximum temperature in summer, minimum temperature in winter, soil organic matter and soil pH. We’re also changing the plotting theme, but that’s just to make the following plots easier to parse and can be skipped if preferred.

data("bunce71")
theme_set(bayesplot::theme_default() + 
            theme(panel.grid.major.y = element_line(colour = "grey92")))

We scale and centre the environmental predictor data - this is particularly necessary as the rainfall data is centred around 1000. We also sort the columns of the data frame of abundance, but this is purely to make the graphs easier to read later on.

bunce71_env <- mutate(bunce71$env, across(pH:WinterMinTemp, scale))
bunce71_abund <- bunce71$abund[,names(sort(colSums(bunce71$abund>0),
                                           decreasing = TRUE))]

We fit the model with unstructured environmental effects, default priors, a binomial family distribution (Number of trials = number of plots per site) and the default number of iterations.

bwmod <- stan_mglmm(~Rainfall + SummerMaxTemp + WinterMinTemp + pH + SOM,
                    data = bunce71_env, Y = bunce71_abund, 
                    Ntrials = bunce71_env$Nplots, 
                    family = "binomial", 
                    cores = 4)
bwmod
Family: binomial 
 Model type: mglmm
  Number of species: 45
  Number of sites: 103
  Number of predictors: 6

Model run on 4 chains with 4000 iterations per chain (2000 warmup).

Parameters with Rhat > 1.01, or Neff/N < 0.05:
                   mean sd 15% 85% Rhat Bulk.ESS Tail.ESS
cor_species[15,15]    1  0   1   1    1     8215     7647
cor_species[17,17]    1  0   1   1    1     7959     7566

No warnings are returned, the model has run acceptably with all Rhats < 1.01 and effective sample sizes > 500. Now we can evaluate fit to data, first by looking at some summary statistics and then by looking at fit for each individual species:

pp_check(bwmod)
../../_images/b5ff01dc9e37377c63371f006d3b6f4060412941e2bb26dbe505c45b164aea90.png ../../_images/5d9a7be1097c67c27d9cc85be8b07881bf7be9a76151f2334b12567b8ce11bf0.png
multi_pp_check(bwmod, plotfun = "ecdf_overlay", species = 1:9)
../../_images/931fbe5aac6fd1b7a6efe2e788aaf6b2e7cdb2e69b0fe5e308fdae0c3c04ad42.png ../../_images/6a56f26320fe6068bb066fd5f47d21192c9c279408c82b5fc5b9df37077dbb7a.png ../../_images/b9422a37a802ee40faf164e9b20b289ad9cdc833efe415261777a5e4e06a5a2e.png ../../_images/b787bfd20d77b921810b8725476a6e13f64f3a4dc5da017f61d3ee3f539199f7.png ../../_images/c053b3b964fe581d6d185fba60267a1ad1d0caf1a09b868af2fc043fe7edf4a8.png

We can now plot the environmental effects for each species:

envplot(bwmod, widths = c(1.5,1,1))
../../_images/ddf2a663bc8424e9640634348f3e4372e2fbc716a99cc1bc0390cb6e72b5ac7d.png

We can also plot the modelled correlations between species once this environmental effect was accounted for:

corrplot(bwmod, species = c("Betula","Fraxinus","Quercus","Corylus"))
../../_images/a6d8b993cce0a778a3b602359c25a41dfa9682fad712fde958414cb02a70187d.png

A subset of these plots are combined for the paper:

e1 <- envplot(bwmod, widths = c(1.4,1), preds = c("pH","SOM"))
c1 <- corrplot(bwmod, species = c("Betula","Fraxinus"))
ggpubr::ggarrange(e1,c1, labels = c("a)","b)","c)","d)"))
../../_images/28bd5d29ef18182ce7e520e429774991d2d708812f526f197435148d8f2854b8.png

GLLVM with spider data#

theme_set(bayesplot::theme_default())

We use the spider data from the mvabund package to run a gllvm model with two latent variables and a negative binomial response. All other options are set to the default. Environmental predictors are scaled.

data("spider", package = "mvabund")
spider$x <- scale(spider$x)
spmod <- stan_gllvm(X = spider$x, Y = spider$abund, D = 2,
                    family = "neg", cores = 4)
spmod
Family: neg_binomial 
 With parameters: kappa 
Model type: gllvm with 2 latent variables
  Number of species: 12
  Number of sites: 28
  Number of predictors: 7

Model run on 4 chains with 4000 iterations per chain (2000 warmup).

No parameters with Rhat > 1.01 or Neff/N < 0.05

No warnings are returned, the model has run acceptably with all Rhats < 1.01 and effective sample sizes > 500. Now we can evaluate fit to data, first by looking at some summary statistics and then by looking at fit for each individual species:

p1 <- pp_check(spmod, plotfun = "ecdf_overlay", ndraws = 20) +
  labs(x = "Abundance") +
  coord_cartesian(xlim = c(0,750))
p2 <- pp_check(spmod,  plotfun = "ecdf_overlay", 
         summary_stat = function(x) sum(x>0), ndraws = 20) + 
  labs(x = "Richness")+
  coord_cartesian(xlim = c(0,12))
p1 + p2 + plot_annotation(tag_levels = "a")
../../_images/8a6164c699b6765eda5aba4c11b84500f098e318a2310fe9eb5f2adc8a52cf6f.png
multi_pp_check(spmod, plotfun = "ecdf_overlay")
../../_images/13968758d51eb08bd56efd63b3549d930992b482d4a0e602bdd5af136693786b.png

We can now plot the environmental effects for each species:

envplot(spmod, widths = c(1.2,1,1))
../../_images/88a299685ed466b75043e54145670305f9cf3df055882bc0b5f95aa9a70b51c4.png

We can also plot latent variable loadings for the species and sites:

o1 <- ordiplot(spmod, ndraws = 0, geom = "text", size = c(3,1),
               errorbar_linewidth = 0.5, errorbar_range = 0.5) + 
  theme(legend.position = "none") 
o2 <- ordiplot(spmod, ndraws = 0, geom = "text", type = "sites", size = c(3,1), 
               errorbar_range = NULL)
o1 + o2 + plot_annotation(tag_levels = "a")
../../_images/9ef43b8ac20d964da8ec5a5729bccd92d81fd905243a88b09bb2ac726b49141e.png

Note that while the model fits the data well, the estimates for the latent variable loadings indicate that the model are trending towards a linear fit with each other indicating the model does not have sufficient information to estimate these latent variables on top of the six environmental predictors. Therefore, we shall fit a reduced model with fewer environmental predictors included.

spmod2 <- update(spmod, newX = spider$x[,c("soil.dry","herb.layer")],
                 cores = 4)
spmod2
Family: neg_binomial 
 With parameters: kappa 
Model type: gllvm with 2 latent variables
  Number of species: 12
  Number of sites: 28
  Number of predictors: 3

Model run on 4 chains with 4000 iterations per chain (2000 warmup).

No parameters with Rhat > 1.01 or Neff/N < 0.05

No warnings are returned, the model has run acceptably with all Rhats < 1.01 and effective sample sizes > 500. Now we can evaluate fit to data, first by looking at some summary statistics and then by looking at fit for each individual species:

p1b <- pp_check(spmod2, plotfun = "ecdf_overlay", ndraws = 20) +
  labs(x = "Abundance") +
  coord_cartesian(xlim = c(0,750))
p2b <- pp_check(spmod2,  plotfun = "ecdf_overlay", 
         summary_stat = function(x) sum(x>0), ndraws = 20) + 
  labs(x = "Richness")+
  coord_cartesian(xlim = c(0,12))
p1b + p2b + plot_annotation(tag_levels = "a")
../../_images/aeb84940f7a57b4653b2b4bd7e2f50b8ad039630e686452e1f2a3a84f8dfefe0.png

We combine (some of) the posterior predictive check plots for the two models for the paper:

p2 + p2b + plot_annotation(tag_levels = "a")
../../_images/c010835f0af91085f5398b3ecf03d8860845be1a06633c297f4f1bb9a384116e.png
multi_pp_check(spmod2, plotfun = "ecdf_overlay")
../../_images/a953e0805065733c6a69c1e82613f0841541e690ba9e09e15c457bcedd21ec67.png

We can also plot the environmental effects for each species:

envplot(spmod2, widths = c(1.2,1))
../../_images/eb5cd6ea46435de4d0e7358618d7c0ba1b5343ecd4c1d7389cef9e7e37b1063e.png

We can also plot latent variable loadings for the species and sites:

o1v2 <- ordiplot(spmod2, ndraws = 0, geom = "text", size = c(3,1),
                 errorbar_linewidth = 0.5, errorbar_range = 0.5) + 
  theme(legend.position = "none")
o2v2 <- ordiplot(spmod2, ndraws = 0, geom = "text", type = "sites", size = c(3,1), 
                 errorbar_range = NULL)
o1v2 + o2v2 + plot_annotation(tag_levels = "a")
../../_images/d14ecd2b5c22ba12370a3d45a2b0c3fa1d8b08e88906edfae337de3e364c57ff.png

Comparing the latent variable plots for the two models shows clearly that the reduced model loses the linear trend in latent variable loadings.

o1 + o2 + o1v2 + o2v2 + plot_annotation(tag_levels = "a")
../../_images/887c7723b7921a9e029384e836e554d2abf7bce1ea2d438111e4a104da795a61.png

References#

Warton et al (2015) So many variables: joint modeling in community ecology. Trends in Ecology & Evolution, 30:766-779. DOI: 10.1016/j.tree.2015.09.007.

Wilkinson et al (2021) Defining and evaluating predictions of joint species distribution models. Methods in Ecology and Evolution, 12:394-404. DOI: 10.1111/2041-210X.13518.

Vehtari, A., Gelman, A., and Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing. 27(5), 1413–1432. DOI: 10.1007/s11222-016-9696-4.