jsdmstan: An R package for fitting joint species distribution models in Stan

jsdmstan: An R package for fitting joint species distribution models in Stan | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results jsdmstan: An R package for fitting joint species distribution models in Stan View ORCID Profile Fiona Seaton doi: https://doi.org/10.1101/2025.11.10.687559 Fiona Seaton 1 UK Centre for Ecology & Hydrology , Lancaster, UK Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Fiona Seaton For correspondence: fseaton{at}ceh.ac.uk Abstract Full Text Info/History Metrics Data/Code Preview PDF Abstract Joint species distribution models (JSDMs) have become an increasingly utilised tool for modelling and predicting change within communities of species across environmental gradients. Here we present a new R package for the fitting of JSDMs using the Bayesian probabilistic programming language Stan. The jsdmstan package can model species responses to environmental covariates and also species interactions through either full specification of the species covariance matrix or through a latent variable formulation. It also provides tools for simulating joint species distribution data according to those models. It supports specification of prior distributions of all parameters in the model, as well as access to the full suite of Stan diagnostics. The ability of this package to fit these models is demonstrated upon two real data sets, one on tree species in a survey of broadleaved woodlands and another on dune spider populations. The models are able to successfully recover population characteristics such as richness and show ecologically interesting results regarding residual species correlations and responses to environmental factors. The jsdmstan package provides a user-friendly interface for fitting JSDMs, with tools available to better understand both the consequences of the model assumptions and how well the model is performing. Introduction Joint species distribution models (JSDMs) involve modelling the distribution of multiple species simultaneously, and have been growing in popularity amongst ecologists ( Warton et al., 2015 ). The main difference between JSDMs and standard species distribution models is that the covariance between species is modelled at the same time as the response of each species to the environmental covariates ( Pollock et al., 2014 ). This allows for information to be borrowed across species, such that the responses of rare species to the environment may be related to the responses of common species and therefore be more biologically realistic. However, the efficacy of this will depend upon the similarity of the niches of the different species ( Erickson & Smith, 2023 ). Modelling the species variance-covariance matrix can be done in a variety of ways and some JSDM methods, i.e. those incorporating a latent variable approach, also allow for model based ordination ( Hui et al., 2015 ). JSDMs can, however, be high dimensional models, with each model potentially incorporating hundreds or thousands of parameters – generally constrained such that the models are identifiable. This complexity can lead to difficulties in understanding if the model fitting process has occurred successfully, and whether the model fit is truly representing the data generating process. Within this paper we present the jsdmstan package, which incorporates both JSDM data simulation and model fitting using Stan to enable greater understanding of the constraints and features of JSDMs through access to a wider suite of diagnostics and exploratory tools. Stan is a probabilistic programming language used to specify statistical models, and can be used to fit these statistical models using a variety of Bayesian algorithms including Hamiltonian Monte Carlo ( Carpenter et al., 2017 ). The default algorithm is a dynamic Hamiltonian Monte Carlo algorithm, incorporating a series of improvements upon the original No-U-Turn Sampler ( Betancourt, 2017 ). This provides several advantages over Gibbs samplers, adopted within other Bayesian software packages such as BUGS and JAGS, including greater efficiency and robustness in exploring the entire parameter space even in very high dimensional models. Stan also provides diagnostics that show when the entire parameter space has not been explored, and that when the model structure or fitting procedure needs to be revisited. This is particularly important in the case of JSDMs due to their high complexity and dimensionality, which can lead to poor model fits and pathological behaviour within the fitting algorithms. Therefore, the use of Stan to fit these models provides robust information for ecological inference. Underlying modelling approach Within jsdmstan each JSDM is fit such that the site ( i ) by species ( j ) matrix ( M ij ) is a function of a species specific intercept ( β i0 ), an environmental effect represented by species specific coefficients multiplied against the environmental covariate matrix and the species covariance ( u ij ). Optionally, an additional site-specific intercept ( α i ) can be added to adjust for total abundance or richness. This is represented in Equation 1 : Where g (·)is the link function, is the transpose of the environmental covariate vector x i , and for each taxon jβ 0j is an intercept and β j is a vector of regression coefficients related to measured site varying predictors. There are two ways of representing the species covariance that can be used within jsdmstan. The first is that the whole species covariance matrix can be modelled, referred to as a Multivariate Generalised Linear Mixed Model (MGLMM) throughout the jsdmstan package. Here the species covariance matrix is treated as a random draw from a multivariate normal distribution with mean 0 and covariance matrix ∑ , shown in equation 2 . 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. Alternatively, therefore, we also include Generalised Linear Latent Variable Models (GLLVMs) in which the species covariance, as included in equation 1 ( u ij ), is specified as a linear function of a set of latent variables ( z i ), equation 3 . The latent variables are assumed to draw from standard normal distributions, i.e. mean 0 and standard deviation of 1. The number of latent variables is specified by the user and is generally much lower than the total number of species. In order to make GLLVMs identifiable the upper triangle of the species loadings to latent variables ( λ i ) must be fixed to zero, i.e. the first species must only load onto the first latent variable, the second species onto the first two latent variables and so on. Also, as these latent variables have the same ecological meaning whether they are positive or negative, as the sampling can lead to sign switching, which is fixed post-sampling by jsdmstan so that the signs of the latent variables do not switch within the results. Within jsdmstan, both the GLLVM and MGLMM methods model the response of species to environmental covariates either as each response independently drawing from a prior distribution (the default), or constrained by a covariance matrix between the environmental covariates. This second option 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 β i (a j by k matrix, assuming k environmental covariates) being a random draw from a multivariate normal with mean 0 and covariance matrix (of dimension k by k ) Ф . Package description The jsdmstan R package provides a variety of functions for both the simulation of species community data and the fitting of JSDMs to species community data. It interfaces with the rstan package for model fitting and the bayesplot package for plotting the model outputs ( Gabry & Mahr, 2022 ; Stan Development Team, 2025 ). It can also interface with the loo package for leave-one-out cross validation of the model ( Vehtari et al., 2020 ). The key functions are jsdm_sim_data which simulates data according to a specified JSDM structure and stan_jsdm which takes the data input and fits the appropriate JSDM in Stan to estimate model parameters. Both of these functions work with both a MGLMM approach and a GLLVM approach, and allow for specification of as many sites, species and environmental covariates as the user requires. The community of species can be analysed using a variety of different response families, currently including the Gaussian, Poisson, negative binomial, Bernoulli, binomial, zero-inflated Poisson and zero-inflated negative binomial families. There are a suite of functions that can be used to examine the results of each model, including summary calls and a variety of plotting functions. Workflow Simulation of data . Data can be simulated using jsdm_sim_data , or with the aliases mglmm_sim_data and gllvm_sim_data which simulate data with the appropriate method. The number of species, sites, and environmental covariates can be specified by the user, in addition to directly supplying the function with environmental data to simulate the species data from. The prior distribution of the parameters used to simulate the data can also be changed through a call to jsdm_prior() within the simulation function. Fitting the models . Models are fit with the stan_jsdm function, or alternatively the stan_mglmm and stan_gllvm aliases of this function. The Stan code used to fit each model can be accessed through supplying the relevant arguments to jsdm_stancode() , and the prior distributions used within the model can be modified through a call to jsdm_prior() within the stan_jsdm function . Summarising the model fits . Textual summaries. The summary and print functions return the estimated mean and quantiles of each parameter as well as the R-hat statistic and effective sample size of the bulk (Bulk-ESS) and tail (Tail-ESS) of each parameter estimate. Graphical summaries. The default plot behaviour is to call mcmc_combo from the bayesplot package and return a traceplot and density plot for a selection of parameters. Graphical MCMC checks are provided by the mcmc_plot function, which interfaces to the mcmc_ family of functions within bayesplot. The envplot function plots the estimated environmental effects upon each species. The ordiplot function is specific to GLLVM models and plots the species or sites scores against the latent variables. The corrplot function is specific to MGLMM models and plots the modelled species correlations. Predicting from the model . 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 the corresponding probability distribution to produce realisations on the original scale. Graphical posterior predictive checks . These are supported through the pp_check function, which interfaces to the ppc_ family of functions within bayesplot. This family of functions provides a graphical way to check your posterior against the data used within the model to evaluate model fit. 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 summary statistic can be changed, as can whether it is calculated for every species or every site. The complementary multi_pp_check function plots a graphical posterior predictive check for each species. Examples To demonstrate the potential usage of the jsdmstan package, we will now describe its application to two real ecological datasets. The first dataset is of records of tree genera occurrence from a survey of broadleaved woodlands where multiple plots were measured in each site. This allows demonstration of the usage of jsdmstan for binomial data where the number of trials is greater than one, and we take a MGLMM approach of estimating the entire covariance matrix. The second case study uses the widely analysed dune spider dataset, which has been used to demonstrate the usage of GLLVM models across many fitting methodologies and packages (e.g. Hui et al., 2015 ; Popovic et al., 2022 ; van der Veen et al., 2021 ). Here, we apply a GLLVM approach to this dataset for comparability with other methods, and discuss how jsdmstan deals with situations where there is a low amount of data for the complexity of the model. Broadleaved woodland data Here we use data from a 1971 survey of broadleaved woodlands across Great Britain, which surveyed up to 16 plots per site for vegetation and soil characteristics. We model all 45 tree genera in response to soil pH, soil organic matter and three climate variables averaged over the 30 years previous to survey (annual rainfall, minimum winter temperature and maximum summer temperature). Full details on this dataset are available in Wood et al. (2015) , and this data is available as part of the jsdmstan package under the name “bunce71”. The model was a MGLMM model, with all species assumed to have a binomial distribution with the number of trials equal to the number of plots per site and assuming no covariance structure within the environmental effect parameterisation (that is, the β j s were taken to be independent of each other). All variables were scaled to have a mean of zero and standard deviation of one before model fitting. The model was fit with default jsdmstan parameters and took approximately thirty minutes over four cores on a windows laptop. All R-hat values were under 1.01 and the number of effective samples for every parameter in both the bulk and the tail of the distribution was above 1000. Graphical posterior checks indicate a good fit to the data, with both the number of genera per site and the number of sites per genus being recovered well by the model ( Figure 1 ). It can be seen that most sites have between 10-15 genera present, while most genera appear in less than half the sites. The more common genera have more clearly resolved environmental responses and species covariance compared to the less common genera ( Figure 2 ). The model can clearly find biological responses such as the preference of ash ( Fraxinus excelsior – the only species present in the Fraxinus genus) for high pH soils and the presence of pines upon high carbon content soils ( Figure 2a ). The residual species correlations can be interpreted in terms of the successional gradient of the woodlands. For example, birch (Betula) is usually present early in succession while ash is present mid-to-late succession and so these species are generally negatively correlated with each other and instead found more often with genera typical of their usual successional stage ( Figure 2b ), e.g. birches are negatively associated with elms (Ulmus) and maples (Acer) while ash is positively associated with these ( Atkinson, 1992 ; Thomas, 2016 ; Thomas et al., 2018 ). We cannot accurately disentangle from these results whether these are true biotic associations or represent shared response to environmental drivers, but it is likely that these represent shared response to time since disturbance or other metrics relevant to successional stage that we did not include in our model. Download figure Open in new tab Figure 1. Model recovery of the number of genera per site (a) and the number of sites per genus (b). The dark blue line is the data and the light blue lines are a selection of twenty random draws from the model estimates. Close correspondence between the data and the model draws demonstrates the good fit of the model. Download figure Open in new tab Figure 2. The estimated effects of soil pH, and SOM upon each genus (a) and the estimated correlations between genera for Betula and Fraxinus (b). Genera are listed in order of frequency of appearance, with the genera that appear in the most sites at the top of the figure. The estimates are shown as the median (circle), 50% interval (thick blue line) and the 90% interval (thin blue line). Dune spider data Here we use data from a pitfall trap survey of spiders in Dutch dunes to demonstrate fitting a GLLVM. This data consists of 12 spider species measured over 28 dune sites with 6 environmental covariates measured at the same sites and is available within the mvabund R package ( van der Aart & Smeenk-Enserink, 1975 ; Wang et al., 2012 ). The model was run with two latent variables (from equation 3 , z ={ z 1 z 2 }), a negative binomial response and default jsdmstan parameters. This took approximately 2 minutes on a windows laptop over four cores, with all parameters having R-hat less than 1.01 and effective sample sizes of the bulk and tail being over 1000. Initial attempts to run this model with all six environmental covariates showed that the loadings onto the latent variables for both the species and the sites trended towards a diagonal with high uncertainty around individual estimates ( Figure 3a-b ). This indicates that the model was unable to fully resolve the latent variables, likely due to there not being enough information in the data to separate out the latent variable effects from the environmental covariate effects. Therefore, the model was run again with two environmental covariates (soil.dry and herb.layer) and two latent variables. Again, this model had all parameters with R-hat less than 1.01 and effective sample sizes of the bulk and tail being over 1000. This model resolved the latent variable effects, with no diagonal loadings apparent ( Figure 3c-d ). It is worth noting that both models successfully recovered the distribution of abundance and richness across the sites ( Figure 4 ), indicating that fit to data alone is insufficient to detect all potential model fit issues. Download figure Open in new tab Figure 3. The species (a, c) and site (b, d) mean loadings onto the two latent variables for the model with six environmental covariates (a, b) and the model with two environmental covariates (c, d). For species only the 50% range of the loadings are also shown as lines upon the graph. Download figure Open in new tab Figure 4. Model recovery of species richness for the model with all environmental covariates (a) and the model with two environmental covariates (b). The dark blue line is the data and the light blue lines are draws from the model estimate. Ecological interpretation of the second model suggests that the separation out of pardmont, arctperi, alopfabr, alopacce within this second model could be due to a shared positive response to moss. As moss is no longer included as an environmental predictor, the latent variables are free to take up the explanatory power of the moss gradient, which is consistent with the understanding of latent variables generally being measures of missing environmental variables ( Vallé et al., 2024 ). Balancing the ability to include multiple environmental predictors and more general latent factors in low data contexts is a challenge, and the most appropriate choice will depend on the intended use of any particular model. Conclusions The jsdmstan package offers a new way to fit joint species distribution models, through using the Stan probabilistic programming language with advanced model diagnostics and graphical checks incorporated. It supports multiple ways of specifying the species interactions, including both specification of the full covariance matrix and specification of the covariance matrix through a latent variable model. Future directions for this package include incorporating spatio-temporal and other random effect information into the model, as well as working on speed improvements to model fitting and comparison. We expect that the jsdmstan package will provide ecologists with the tools to fit a variety of joint species distribution models; therefore achieving greater understanding of both how well these models are describing the data and how environmental effects and species associations are affecting species communities. Data and code availability The Bunce woodland data is available from the Environmental Information Data Centre under ( Kirby et al., 2013c , 2013a , 2013b ). The climatic data is available to download from the NERC EDS Centre for Environmental Data Analysis (Met Office et al., 2022 ). The Bunce woodland data used in this paper are available within the jsdmstan package presented here, and the spider data is available within the mvabund R package ( Wang et al., 2012 ). The jsdmstan R package is available at https://github.com/NERC-CEH/jsdmstan ( Seaton, 2025 ), and the R code for the case studies included in the Supplementary Material below and at https://github.com/NERC-CEH/ds-toolbox-notebook-jsdmstan . Supplementary Material: Code for Case Studies library(jsdmstan) library(dplyr) library(ggplot2) library(patchwork) set.seed(3757892) ## ----read in bunce71 data data(“bunce71”) ## ----scale b71 data bunce71_env <-mutate(bunce71$env, across(pH:WinterMinTemp, scale)) bunce71_abund 0), decreasing = TRUE))] ## ----bw model fit, cache = TRUE bwmod <-stan_mglmm(∼Rainfall + SummerMaxTemp + WinterMinTemp + pH + SOM, data = bunce71_env, Y = bunce71_abund, Ntrials = bunce71_env$Nplots, family = “binomial”, cores = 4) bwmod ## ----bw pp_check pp_check(bwmod) p1 0), ndraws = 20) + labs(x = “#Genera per site”) p2 0), ndraws = 20) + labs(x = “#Sites per genus”) p1+p2 + plot_annotation(tag_levels = “a”) ## ----bw multi_pp_check multi_pp_check(bwmod, plotfun = “ecdf_overlay”, species = 1:9) multi_pp_check(bwmod, plotfun = “ecdf_overlay”, species = 10:18) multi_pp_check(bwmod, plotfun = “ecdf_overlay”, species = 19:27) multi_pp_check(bwmod, plotfun = “ecdf_overlay”, species = 28:36) multi_pp_check(bwmod, plotfun = “ecdf_overlay”, species = 37:45) ## ----bw envplot envplot(bwmod, widths = c(1.5,1,1)) ## ----bw corrs corrplot(bwmod, species = c(“Betula”,”Fraxinus”,”Quercus”,”Corylus”)) 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)”)) ## ----get spider data data(“spider”, package = “mvabund”) spider$x <-scale(spider$x) ## ----run spider mod spmod <-stan_gllvm(X = spider$x, Y = spider$abund, D = 2, family = “neg”, cores = 4) spmod ## ----sp pp_check p1 <-pp_check(spmod, plotfun = “ecdf_overlay”, ndraws = 20) + labs(x = “Abundance”) + coord_cartesian(xlim = c(0,750)) p2 0), ndraws = 20) + labs(x = “Richness”)+ coord_cartesian(xlim = c(0,12)) p1 + p2 + plot_annotation(tag_levels = “a”) ## ----sp multi_pp_check multi_pp_check(spmod, plotfun = “ecdf_overlay”) ## ----sp envplot envplot(spmod, widths = c(1.2,1,1)) ## ----sp ordiplot 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”) ## ----update sp mod with fewer covariates spmod2 <-update(spmod, newX = spider$x[,c(“soil.dry”,”herb.layer”)], cores = 4) spmod2 ## ----sp2 pp_check p1b <-pp_check(spmod2, plotfun = “ecdf_overlay”, ndraws = 20) + labs(x = “Abundance”) + coord_cartesian(xlim = c(0,750)) p2b 0), ndraws = 20) + labs(x = “Richness”)+ coord_cartesian(xlim = c(0,12)) p1b + p2b + plot_annotation(tag_levels = “a”) p2 + p2b + plot_annotation(tag_levels = “a”) ## ----sp2 multi_pp_check multi_pp_check(spmod2, plotfun = “ecdf_overlay”) ## ----sp2 envplot envplot(spmod2, widths = c(1.2,1)) ## ----sp2 ordiplot 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) ## ----sp all ordiplots o1 + o2 + o1v2 + o2v2 + plot_annotation(tag_levels = “a”) Acknowledgements This work would not have been possible without the Stan forums community, particularly Jim Downie who developed Stan code for a GLLVM. Thanks are also due to Adam Kimberley and Simon Smart who offered advice for understanding the woodland data and model results, and to Pete Henrys for his feedback on the modelling approach. This research was supported by NERC, through the UKCEH National Capability for UK Challenges Programme NE/Y006208/1. Funder Information Declared Natural Environment Research Council , NE/Y006208/1 Footnotes https://doi.org/10.5281/zenodo.17396119 References ↵ Atkinson , M. D. ( 1992 ). Betula Pendula Roth (B. Verrucosa Ehrh.) and B. Pubescens Ehrh . Journal of Ecology , 80 ( 4 ), 837 – 870 . doi: 10.2307/2260870 OpenUrl CrossRef Web of Science ↵ Betancourt , M. ( 2017 ). A Conceptual Introduction to Hamiltonian Monte Carlo . http://arxiv.org/abs/1701.02434 ↵ Carpenter , B. , Gelman , A. , Hoffman , M. D. , Lee , D. , Goodrich , B. , Betancourt , M. , Brubaker , M. A. , Guo , J. , Li , P. , & Riddell , A. ( 2017 ). Stan: A probabilistic programming language . Journal of Statistical Software , 76 ( 1 ), 1 – 32 . doi: 10.18637/jss.v076.i01 OpenUrl CrossRef PubMed ↵ Erickson , K. D. , & Smith , A. B. ( 2023 ). Modeling the rarest of the rare: A comparison between multii species distribution models, ensembles of small models, and singleispecies models at extremely low sample sizes . Ecography , 2023 ( 6 ), 861 – 875 . doi: 10.1111/ecog.06500 OpenUrl CrossRef ↵ Gabry , J. , & Mahr , T. ( 2022 ). bayesplot: Plotting for Bayesian Models . https://mc-stan.org/bayesplot/ ↵ Hui , F. K. C. , Taskinen , S. , Pledger , S. , Foster , S. D. , & Warton , D. I. ( 2015 ). Model-based approaches to unconstrained ordination . Methods in Ecology and Evolution , 6 ( 4 ), 399 – 411 . doi: 10.1111/2041-210X.12236 OpenUrl CrossRef ↵ Kirby , K. J. , Smart , S. M. , Black , H. I. J. , Bunce , R. G. H. , Corney , P. M. , Smithers , R. J. , & Shaw , M. W. ( 2013a ). Woodlands survey site information 1971-2001 [Dataset]. NERC Environmental Information Data Centre . doi: 10.5285/d6409d40-58fe-4fa7-b7c8-71a105b965b4 OpenUrl CrossRef ↵ Kirby , K. J. , Smart , S. M. , Black , H. I. J. , Bunce , R. G. H. , Corney , P. M. , Smithers , R. J. , & Shaw , M. W. ( 2013b ). Woodlands survey soil data 1971-2001 [Dataset]. NERC Environmental Information Data Centre . doi: 10.5285/fb1e474d-456b-42a9-9a10-a02c35af10d2 OpenUrl CrossRef ↵ Kirby , K. J. , Smart , S. M. , Black , H. I. J. , Bunce , R. G. H. , Corney , P. M. , Smithers , R. J. , & Shaw , M. W. ( 2013c ). Woodlands survey tree diameter data 1971-2001 [Dataset]. NERC Environmental Information Data Centre . doi: 10.5285/4d93f9ac-68e3-49cf-8a41-4d02a7ead81a OpenUrl CrossRef ↵ Met Office , Hollis , D. , McCarthy , M. , Kendon , M. , & Legg , T. ( 2022 ). HadUK-Grid Gridded Climate Observations on a 1km grid over the UK, v1.1.0.0 (1836-2021) [Dataset] . NERC EDS Centre for Environmental Data Analysis . doi: 10.5285/bbca3267dc7d4219af484976734c9527 OpenUrl CrossRef ↵ Pollock , L. J. , Tingley , R. , Morris , W. K. , Golding , N. , O’Hara , R. B. , Parris , K. M. , Vesk , P. A. , & Mccarthy , M. A. ( 2014 ). Understanding co-occurrence by modelling species simultaneously with a Joint Species Distribution Model (JSDM) . Methods in Ecology and Evolution , 5 ( 5 ), 397 – 406 . doi: 10.1111/2041-210X.12180 OpenUrl CrossRef ↵ Popovic , G. C. , Hui , F. K. C. , & Warton , D. I. ( 2022 ). Fast modelibased ordination with copulas . Methods in Ecology and Evolution , 13 ( 1 ), 194 – 202 . doi: 10.1111/2041-210X.13733 OpenUrl CrossRef ↵ Seaton , F.M. ( 2025 ). NERC-CEH/jsdmstan: v0.4.3 (v0.4.3) . Zenodo . doi: 10.5281/zenodo.17396119 OpenUrl CrossRef ↵ Stan Development Team . ( 2025 ). RStan: The R interface to Stan . https://mc-stan.org/ ↵ Thomas , P. A. ( 2016 ). Biological Flora of the British Isles: Fraxinus excelsior . Journal of Ecology , 104 ( 4 ), 1158 – 1209 . doi: 10.1111/1365-2745.12566 OpenUrl CrossRef ↵ Thomas , P. A. , Stone , D. , & La Porta , N. ( 2018 ). Biological Flora of the British Isles: Ulmus glabra . Journal of Ecology , 106 ( 4 ), 1724 – 1766 . doi: 10.1111/1365-2745.12994 OpenUrl CrossRef ↵ Vallé , C. , Poggiato , G. , Thuiller , W. , Jiguet , F. , Princé , K. , & Le Viol , I. ( 2024 ). Species associations in joint species distribution models: From missing variables to conditional predictions . Journal of Biogeography , 51 ( 2 ), 311 – 324 . doi: 10.1111/jbi.14752 OpenUrl CrossRef ↵ van der Aart , P. J. M. , & Smeenk-Enserink , N. ( 1975 ). Correlations between distributions of hunting spiders (Lycos-idae, Ctenidae) and environmental characteristics in a dune area . Netherlands Journal of Zoology , 25 , 1 – 45 . OpenUrl ↵ van der Veen , B. , Hui , F. K. C. , Hovstad , K. A. , Solbu , E. B. , & O’Hara , R. B. ( 2021 ). Model-based ordination for species with unequal niche widths . Methods in Ecology and Evolution , 12 ( 7 ), 1288 – 1300 . doi: 10.1111/2041-210X.13595 OpenUrl CrossRef ↵ Vehtari , A. , Gabry , J. , Magnusson , M. , Yao , Y. , Bürkner , P.-C. , Paananen , T. , & Gelman , A. ( 2020 ). loo: Efficient leave-one-out cross-validation and WAIC for Bayesian models . https://mcstan.org/loo ↵ Wang , Y. , Naumann , U. , Wright , S. T. , & Warton , D. I. ( 2012 ). mvabund—An R package for model-based analysis of multivariate abundance data . Methods in Ecology and Evolution , 3 ( 3 ), 471 – 474 . doi: 10.1111/j.2041-210X.2012.00190.x OpenUrl CrossRef ↵ Warton , D. I. , Blanchet , F. G. , O’Hara , R. B. , Ovaskainen , O. , Taskinen , S. , Walker , S. C. , & Hui , F. K. C. ( 2015 ). So many variables: Joint modeling in community ecology . Trends in Ecology & Evolution , 30 ( 12 ), 766 – 779 . doi: 10.1016/j.tree.2015.09.007 OpenUrl CrossRef PubMed ↵ Wood , C. M. , Smart , S. M. , & Bunce , R. G. H. ( 2015 ). Woodland Survey of Great Britain 1971–2001 . Earth System Science Data , 7 ( 2 ), 203 – 214 . doi: 10.5194/essd-7-203-2015 OpenUrl CrossRef View the discussion thread. Back to top Previous Next Posted November 11, 2025. Download PDF Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following jsdmstan: An R package for fitting joint species distribution models in Stan Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share jsdmstan: An R package for fitting joint species distribution models in Stan Fiona Seaton bioRxiv 2025.11.10.687559; doi: https://doi.org/10.1101/2025.11.10.687559 Share This Article: Copy Citation Tools jsdmstan: An R package for fitting joint species distribution models in Stan Fiona Seaton bioRxiv 2025.11.10.687559; doi: https://doi.org/10.1101/2025.11.10.687559 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Ecology Subject Areas All Articles Animal Behavior and Cognition (7635) Biochemistry (17690) Bioengineering (13892) Bioinformatics (41936) Biophysics (21451) Cancer Biology (18588) Cell Biology (25499) Clinical Trials (138) Developmental Biology (13378) Ecology (19899) Epidemiology (2067) Evolutionary Biology (24320) Genetics (15609) Genomics (22506) Immunology (17736) Microbiology (40394) Molecular Biology (17181) Neuroscience (88603) Paleontology (666) Pathology (2832) Pharmacology and Toxicology (4824) Physiology (7641) Plant Biology (15152) Scientific Communication and Education (2045) Synthetic Biology (4294) Systems Biology (9825) Zoology (2271)