cosimmr: an R package for fast fitting of Stable Isotope Mixing Models with covariates

1. 1.

   The study of animal diets and the proportional contribution that different foods make to their diets is an important task in ecology. Stable Isotope Mixing Models (SIMMs) are an important tool for studying an animal’s diet and understanding how the animal interacts with its environment.

2. 2.

   We present cosimmr, a new R package designed to include covariates when estimating diet proportions in SIMMs, with simple functions to produce plots and summary statistics. The inclusion of covariates allows for users to perform a more in-depth analysis of their system and to gain new insights into the diets of the organisms being studied.

3. 3.

   A common problem with the previous generation of SIMMs is that they are very slow to produce a posterior distribution of dietary estimates, especially for more complex model structures, such as when covariates are included. The widely-used Markov chain Monte Carlo (MCMC) algorithm used by many traditional SIMMs often requires a very large number of iterations to reach convergence. In contrast, cosimmr uses Fixed Form Variational Bayes (FFVB), which we demonstrate gives up to an order of magnitude speed improvement with no discernible loss of accuracy.

4. 4.

   We provide a full mathematical description of the model, which includes corrections for trophic discrimination and concentration dependence, and evaluate its performance against the state of the art MixSIAR model. Whilst MCMC is guaranteed to converge to the posterior distribution in the long term, FFVB converges to an approximation of the posterior distribution, which may lead to sub-optimal performance. However we show that the package produces equivalent results in a fraction of the time for all the examples on which we test. The package is designed to be user-friendly and is based on the existing simmr framework.

Keywords: bayesian statistics; mixing models; stable isotopes; trophic ecology; variational bayes;

Stable Isotope Mixing Models (SIMMs) are commonly used in ecology to study the proportional contribution that different foods make to an animal’s diet (Phillips, 2012). This information can be important as it allows scientists to look at diet, which resources are important for different species (McDonald et al., 2020), and consequently niche overlap and competition (Teixeira et al., 2021; Aksu et al., 2023), as well as being useful in looking at trophic position and energy flow in an ecosystem (Manlick and Newsome, 2022). These models have been extensively used by ecologists over the past 20 years with recent papers revealing the foraging behaviour in Dugongs (Thibault et al., 2024), overlap in trophic niche between native and non-native species of carp (Aksu et al., 2023), and assessment of nursery areas used by the scalloped hammerhead shark (Paez-Rosas et al., 2024). SIMMs have been shown to produce results comparable to direct observation (Swan et al., 2020). The approach relies on the fact that the stable isotopes of several elements, but most usefully those of hydrogen($\delta^{2}H$), carbon ($\delta^{13}C$), nitrogen ($\delta^{15}N$) and sulphur ($\delta^{34}S$), are incorporated into animal tissues from the diet in a predicable manner (Inger and Bearhop, 2008). Thus, if the isotope ratios of potential dietary items are known then animal diets can be reconstructed from the stable isotope ratios from proteinaceous tissues using SIMMs. cosimmr is a new R package (R Core Team, 2021) developed to allow for the fast running of SIMMs, especially but not limited to those that include covariates. It has been designed to be easy to use for non-expert R users, with S3 classes used throughout. SIMMs are widely used and cited in other fields, such as geology (Munoz et al., 2019) and pollution studies (Zaryab et al., 2022), amongst many others. In other scientific areas, SIMMs and similar models can be referred to in the literature as ‘source apportionment models’ (Hopke, 1991), ‘end member analysis’ (Hooper et al., 1990), or ‘mass balance analysis’ (Miller et al., 1972). Further discussion of these models can be found in Govan et al. (2023).

The basic mathematical equation for a statistical SIMM is:

Here $y$ is the mixture (consumer tissues) value (for example, the $\delta^{13}C$ or $\delta^{15}N$ values for the species we wish to study), $p_{k}$ are the proportions contributed by each source (dietary item) $k$ (of $K$ total sources), $s_{k}$ is the source tracer value for source $k$, and $\epsilon$ is a residual term. The parameters $p_{k}$ are usually the main focus of scientific interest. These models are commonly expanded in diet analyses to allow for processes that cause the mixture and source tracer values to differ besides the source proportions, such as Trophic Discrimination Factors (TDFs; Inger and Bearhop, 2008) and concentration dependence (Phillips and Koch, 2002). Other expansions include process error on the dietary proportions (Moore and Semmens, 2008) as well as hierarchical source fitting (Ward et al., 2010). The models are further made richer by incorporating random effects on the source values (Semmens et al., 2009). Here, whilst we include many of these extensions, our focus is on the inclusion of covariate dependence on the proportions $p_{k}$. The restriction that these must sum to unity (i.e. a simplex) makes their estimation more complex, and specialist link functions are required to map their values on to covariates. We provide a more detailed explanation of the mathematical model behind cosimmr in Section 3.

Modern SIMMs are fitted using the standard tools of Bayesian inference. Most commonly this involves using Markov chain Monte Carlo (MCMC) to obtain samples from the posterior. For complex models with covariates this can be extremely slow, with models requiring millions of posterior samples and taking several days to converge, if they converge at all. By contrast in cosimmr we use Fixed Form Variational Bayes (FFVB), specifically Gaussian Variational Bayes with Cholesky Decomposed Variance (Titsias and Lázaro-Gredilla, 2014; Tan and Nott, 2018). FFVB is an optimisation-based algorithm which works by first defining the form of the posterior distribution (here a multivariate normal distribution), and then minimising the Kullback-Leibler divergence between this VB approximation and the true posterior distribution. More details on the assumed distributions used in cosimmr can be found in Section 4. The main advantage of using FFVB, and thus of cosimmr, is that it works by optimisation rather than sampling, and therefore can be much faster to produce a posterior distribution. However, convergence issues can still occur (Yao et al., 2018). Our approach gives cosimmr a key advantage over other packages, which tend to use JAGS (Plummer, 2003) to implement the MCMC algorithm. Additionally, cosimmr runs using C++ code via Rcpp (Eddelbuettel and François, 2011) which allows for an additional speed boost.

SIMMs have been a popular method for studying animal diets, with thousands of citations across different papers, and many packages have been developed in order to make this easier. Some of the most popular packages are listed below. A summary is also provided in Table 1. A more detailed comparison of cosimmr, simmr, and MixSIAR is provided in Table 2.

• •

  Isosource (Phillips and Gregg, 2003) was one of the earliest developed packages for running SIMMs. It worked by simulating possible values of each proportion to produce many potential combinations of proportions. Importantly, it lacked an explicit statistical basis.

• •

  MixSIR (Moore and Semmens, 2008) adopted a Bayesian framework and allowed for the inclusion of Trophic Discrimination Factors (TDFs). MixSIR utilised Importance Sampling, generating many samples of possible proportion combinations and calculating importance weights to produce the final posterior sample. MixSIR also allowed for the inclusion of prior information and allowed uncertainty to be incorporated into SIMMs.

• •

  SIAR (Parnell et al., 2010) was developed as an R package. It utilised Markov chain Monte Carlo (MCMC) sampling. SIAR also included a residual component $\epsilon$ in the model. SIAR is no longer updated and simmr was developed to replace it.

• •

  simmr (Govan et al., 2023) is an R package that follows a Bayesian framework, and provides the option of using either JAGS (Just Another Gibbs Sampler; Plummer, 2003) or Fixed Form Variational Bayes (FFVB; Tran et al., 2021) for running the models. simmr allows for the inclusion of concentration dependence and Trophic Discrimination Factors but does not allow for covariates on the dietary proportions.

• •

  FRUITS (Fernandes et al., 2014) allowed for the inclusion of concentration dependence and prior information in a Bayesian framework. FRUITS is encoded in Visual Basic and runs via the FRUITS computer programme. This package simplifies the incorporation of prior information.

• •

  IsotopeR (Hopkins III and Ferguson, 2012) adopted a Bayesian framework and allowed for inclusion of concentration dependence and TDFs, and also includes covariance. It uses MCMC for running the models.

• •

  MixSIAR (Stock et al., 2018) is an R package that fits models in JAGS. MixSIAR allows for the inclusion of covariates as fixed, random, or continuous effects. It fits the source means hierarchically, either using raw data or sample statistics (means, variances, and sample sizes).

TDFs (Trophic Discrimination Factors) account for the fact that consumers may differentially lose ‘light’ versions of isotopes with respect to ‘heavy’ versions during the process of assimilation (Inger and Bearhop, 2008). Thus, while TDFs are important in ecological applications, they have less relevance to geological or pollution-based applications as these same processes do not occur (although similar processes such as isotopic fractionation may need to be accounted for). TDFs can be calculated in the lab, or calculated mathematically, such as in Greer et al. (2015). Alternatively SIDER (Healy et al., 2018) is an R package that allows for estimates of TDFs based on phylogenetic relatedness of species. Estimates can also be obtained from the scientific literature.

Concentration dependence accounts for the fact that different food sources can contribute proportionally different amounts of each isotope (Phillips and Koch, 2002). The standard two-isotope model assumes all food sources contribute both isotopes equally. However, there are often occasions where a food source can be rich in one isotope and poor in another, thus not contributing equally to both. Instead, concentration dependence assumes a source’s contribution to each isotope is proportional to the mass of the food source times the elemental concentration of the isotope within the food source. Inclusion of concentration depedence facilitates conversion between consideration of either the total mass of food sources assimilated and the mass of specific elements within them.

The inclusion of covariates in the SIMM allows users to avoid pseudo-replication (Hurlbert, 1984), because if a covariate is important to the diet proportions, its exclusion violates the assumption that all mixtures are independently and identically distributed. Including covariates in SIMMs allows users to determine the potentially causal relationships between covariates and diet proportions. Although the model returns these coefficients, they are only available in a transformed space (via the link function) and not directly interpretable. We have designed cosimmr to produce interpretable output in ‘coefficient space’ where users can determine the direction of the relationship and evaluate uncertainty, and also in ‘proportion space’ ($p$-space) which allows users to see the effect of the covariate directly on the dietary proportions. These tools are defined via a predict function that to allow the user to predict dietary proportions based on combinations of covariates that may not necessarily be present in the original dataset. We follow the ‘tidyverse’ (Wickham et al., 2019) style guide, with the Snake case naming convention and S3 classes used throughout. Being able to evaluate the model at new values of the covariates allows for a more detailed picture to be seen. Uncertainty intervals, in the form of Bayesian credible intervals, are provided for all estimated quantities. Advanced users have access to the full posterior distributions as created by the FFVB algorithm.

The data required by cosimmr can be illustrated by Figure 1. This example is discussed more fully in Section 7.3. This figure shows an ‘iso-space’ plot generated by cosimmr. It shows each individual plotted in iso-space, with the axes representing $\delta^{13}C$ and $\delta^{15}N$. The diet sources (Marine and Freshwater in this case) are also plotted on the graph. It is important that all individuals lie within the mixing polygon created by the source means. If they do not, it indicates that the mixing system does not follow the model assumptions. Possible reasons for observations lying outside the mixing polygon include: issues with data collection; inaccurate TDFs; or missing food sources, amongst others. However, note that the mixing polygon vertices are sample means subject to sampling error, so there is some uncertainty in their exact position. Hierarchical Bayesian models allow for the source means to deviate from the source sample means by maximizing the likelihood of the source and mixture tracer data together (Ward et al., 2010; Hopkins III and Ferguson, 2012; Stock et al., 2018).

The fundamental SIMM we fit can be written as:

The source and TDF random effects add an additional burden of $2KJ$ parameters into the model which can cause a significant computational slowdown. Moore and Semmens (2008) proposed proposed marginalising across these parameters to produce a more complex, but more computationally tractable likelihood:

The remaining prior distribution is that of the $p_{ik}$ terms which must retain the constraint that $\sum_{k}p_{ik}=1$, but also allow for the terms to be dependent on the covariates $\mathbf{x}_{i}$. We use a Centralised Log-Ratio (CLR; Aitchison (1986)) link function so that:

This prior has the advantage that the resulting terms $f_{ik}$ are unconstrained and made to depend directly on $\mathbf{x}_{i}$. We model this dependence linearly, but extensions that capture more nuanced dependence seem like a fruitful avenue for future research. We thus set: $f_{ik}=\mathbf{x}_{i}^{T}\mathbf{\beta}_{k}$. In other words, we can write the proportion for individual $i$ consuming food $k$ as $p_{ik}=CLR(\mathbf{x}_{i}^{T}\mathbf{\beta}_{k})$. where the parameters $\mathbf{\beta}_{k}$ model the dependence of the covariates across source $k$. We require a further prior distribution on $\mathbf{\beta}_{k}$ to ensure identifiability. By default we thus set $\beta_{kl}\sim N(0,1)$. In certain circumstances where prior information is available on the $\beta_{kl}$ values we may use an informative prior of the form $\beta_{kl}\sim N(\mu_{\beta,kl},\sigma_{\beta,kl}^{2})$. The prior distribution for $\beta$ can be changed by the user in cosimmr via the cosimmr_ffvb function.

Fixed Form Variational Bayes (FFVB) is an optimisation-based algorithm that aims to approximate the posterior distribution of a Bayesian model in a pre-defined form (Salimans and Knowles, 2013). It aims to finds the parameters of the ‘closest’ probability distribution to that of the true posterior. Unlike traditional MCMC sampling methods, the greedy nature of the optimisation can usually find this approximate posterior far faster. In cosimmr we use a sub-type of FFVB known as Gaussian Variational Bayes with Cholesky decomposed covariance (Titsias and Lázaro-Gredilla, 2014; Tan and Nott, 2018). To avoid becoming diverted in mathematical detail, we defer a full description of our approach to Appendix A. However here we provide an intuitive guide to how the fitting process works.

Our model assumes that the joint posterior distribution of all the parameters is multivariate normal. The parameters for our model are $\mathbf{\beta}$, and $\mathbf{\sigma^{2}}$, representing the regression parameters across sources and covariates, and the residual variances across tracers. Recall that $p$ is a deterministic function of $\beta$ so not included in the algorithm. Since the variances are all restricted to be positive, we model these on the log scale. We thus write $\mathbf{\theta}=\{\mathbf{\beta},\mathbf{\log(\sigma^{2}})\}=\{\beta_{11},\ldots,\beta_{KL},\log(\sigma^{2}_{1}),\ldots,\log(\sigma^{2}_{J})\}$ as the set of parameters for which we want a posterior.

For FFVB we need to define the form of the posterior distribution. We use:

where $\mathbf{\mu}_{\theta}$ and $\mathbf{\Sigma}_{\mathbf{\theta}}$ are the mean and covariance matrix of the approximated posterior distribution. To avoid the positive semi-definite constraints on $\mathbf{\Sigma}_{\theta}$ we model the Cholesky decomposition of this matrix so that $\mathbf{\Sigma}_{\mathbf{\theta}}=\mathbf{LL}^{T}$. Together these terms are vectorised and written as: $\mathbf{\lambda}=\{\mbox{vec}(\mathbf{\mu}_{\mathbf{\theta}}),\mbox{vec}(\mathbf{L)}\}$. The goal of the algorithm is to provide a set of optimal values for $\mathbf{\lambda}$ which captures this posterior distribution. The main steps of the algorithm are as follows:

• •

  Get starting values and use these to get an estimate of the difference between the posterior and variational posterior.

• •

  Calculate the gradient of this difference and use this to update $\mathbf{\lambda}$. The gradient of the posterior is calculated using automatic differentation and the gradient of the variational posterior is calculated manually.

• •

  Use new values of $\mathbf{\lambda}$ in place of starting values and repeat until stopping condition is met.

The algorithm requires several user-set hyper-parameters for fitting the model. The main hyper-parameters are: the patience $P$, which determines when the algorithm stops; and $S$ which provides the number of parameter samples used at each stage of the algorithm. For other algorithm parameters, we have used the default values from Tran et al. (2021). These parameters include: the fixed (beta_1 and beta_2) and adaptive (eps_0) learning rates; the size of the window to use when calculating stopping conditions (t_W); and the threshold for exploring the learning space before a the learning rate is decreased (tau).The stopping conditions work by calculating a moving average of the lower bound (the difference between the log of the posterior and the log of the variational posterior). When the moving average does not improve after $P$ iterations then the algorithm stops. All hyper-parameters can be changed by the user if they wish when running the FFVB algorithm through cosimmr, though we have provided reasonable defaults which should work in most circumstances

The use of FFVB in covariate-dependent SIMMs is novel, and confers an advantage over MCMC due to the increase in speed. This method is flexible and can be extended in the future, to include hierarchical source fitting, raw data, as well as random effects, features which are currently available in other SIMM software. We see up to an order of magnitude speed increase when comparing FFVB to MCMC, with comparable results produced.

The cosimmr package is available via CRAN. The first step to use it is to install and load the package:

R> install.packages("cosimmr")
R> library(cosimmr)

The user must provide mixture data ($y$), source means ($\mu_{s,kj}$), and source standard deviations ($\sigma_{s,kj}$) to run cosimmr. TDFs, Concentration Dependence, and any covariates are not necessary to run the model but they may need to be included in order for the model to be ecologically valid. The cosimmr_load function can be used to read the data into R. This ensures the data is loaded in the correct format. For illustration purposes we will use an artificially generated dataset, though see later sections for real case studies:

R> y = matrix(c(5, 5.1, 4.7, 3.6, 3.2, 0, -1, -2, -3, -7,
+               3.1, 5.6, 3.6, 4.7, 1.3, 1, -4, -3, -7, -9),
+               ncol = 2)
R> colnames(y) = c("iso1", "iso2")
R> mu_s = matrix(c(-10, 0, 10, -10, 10, 0), ncol = 2, nrow = 3)
R> sigma_s = matrix(c(1, 1, 1, 1, 1, 1), ncol = 2, nrow = 3)
R> s_names = c("A", "B", "C")
R> x = c(1.6, 1.7, 2.1, 2.5, 1.1, 3.7, 4.5, 6.8, 7.1, 7.7)

This dataset contains measurements of two isotopic ratios, ‘iso1’ and ‘iso2’, as well as three food sources named ‘A’, ‘B’, and ‘C’. There is one continuous covariate, named ‘x’. These data can then be loaded into cosimmr using cosimmr_load to create an object of class cosimmr_in:

R> cosimmr_in_1 = cosimmr_load(
+                              formula = y ~ x,
+                              source_names = s_names,
+                              source_means = mu_s,
+                              source_sds = sigma_s
)

As discussed in the introduction, it is recommended that these data are plotted on an ‘iso-space’ plot before modelling. The iso-space plot shows tracer values of the mixtures as well as sources, with each axis representing a tracer. These tracers are often isotope ratios. It is important that the mixture data lies within the polygon formed when the sources are joined with straight lines (this polygon is referred to as the mixing polygon). If the mixtures do not lie within this polygon it indicates that there is an issue - potential reasons are that TDFs are inappropriate or a food source has been omitted. The polygon vertices are subject to sampling error and this could be another potential source of error. The iso-space plot can be generated by cosimmr by running the following code:

R> plot(cosimmr_in_1, col_by_cov = TRUE, cov_name = "x")

The resulting plot can be seen in Figure 2. It shows the mixtures lie within the mixing polygon. These data can then be run through the cosimmr_ffvb function to produce an output of class cosimmr_out:

R> cosimmr_out_1 = cosimmr_ffvb(cosimmr_in_1)

cosimmr has built-in functions to produce summaries of the model run. Both graphical and textual summaries can be produced, as shown below.

R> summary(cosimmr_out_1, type = "statistics")

Summary for Observation 1

         mean    sd
P(A)    0.093 0.026
P(B)    0.424 0.052
P(C)    0.482 0.045
sd_iso1 1.449 1.133
sd_iso2 1.790 1.437

In this example, we have specified that we wish to produce ‘statistics’. As we have not specified an observation then the function defaults to returning statistics for observation 1. Any/multiple individuals can be selected and ‘statistics’, ‘quantiles’ or ‘correlations’ can be produced for each individual. The ‘statistics’ summary produces a table of the means and standard deviations for the estimates of the proportion of each food eaten by the individual. An estimate of the marginal residual error of the isotope ratios is produced. In this summary we can see that individual 1’s diet is mostly composed of foods B and C, with A contributing very little to their diet. This matches the observations we can make from the iso-space plot, where individual 1 lies at (5, 3.1), equidistant from food sources B and C and quite far away from food source A. The ‘quantiles’ summary produces the 2.5%, 25%, 50%, 75% and 97.5% quantiles for the same values as provided by the ‘statistics’ summary. The ‘correlations’ option produces the correlation values between each source and the residual error of the isotope ratios.

Graphical plots can be produced in cosimmr. For example, below we create a proportion plot for observation 1. We can see the plot in Figure 3. This shows the range of the proportion estimates for each food source for individual 1. We also create a ‘covariates plot’ which shows the change in the proportion of foods consumed as the covariate changes. This plot can be seen in Figure 4.

R> plot(cosimmr_out_1, type = c("prop_histogram", "covariates_plot"),
+                      obs = 1, cov_name = "x")

Another important function within cosimmr is the ability to predict values based on covariate values, using the predict function, as illustrated below:

R> x_pred = data.frame(x = c(3, 5))
R> pred_out = predict(cosimmr_out_1, x_pred)
R> summary(pred_out, type = "statistics", obs = c(1,2))

Summary for Observation 1

      mean    sd
P(A) 0.181 0.033
P(B) 0.345 0.041
P(C) 0.474 0.039

Summary for Observation 2

      mean    sd
P(A) 0.393 0.043
P(B) 0.221 0.054
P(C) 0.386 0.054

Here we can select values of the covariate for which we have no observations with those values, add these to the predict function along with out cosimmr_out object, and return a new object which can then be used in the summary and plot functions as with a cosimmr_out object. In this example we return dietary predictions for observations with covariate values 3 and 5, and can see what their estimated diet would look like if there were individuals with these predicted values.

We can use the prior_viz function to visualise how the posterior has changed from the prior. This plot overlays the posterior distribution and the prior distribution to see how it has changed or if the posterior has not changed much from the prior. This plot can be seen in Figure 5. We can see that our posterior estimates have changed from the prior estimate. One figure is shown here but it is recommended that users create plots for multiple individuals when running a model.

For convenience a summary of the main cosimmr functions is presented in Table 4.

In this section we use simulated data to verify that cosimmr returns valid estimates when a run is performed. We simulate data from the model using a variety of different data set sizes, varying $N$, $J$, $K$ and $L$, and changing the main parameters in the model. We evaluate the performance of the model by looking at how often the posterior distribution obtained by cosimmr matches the true values. The code for performing the runs in this section can be found at \urlhttps://github.com/emmagovan/cosimmr_paper/.

The values selected for several different simulations are presented in Table 5. We run low ($N=50$, $J=2$, $K=3$, $L=2$), medium ($N=200$, $J=3$, $K=4$, $L=5$) and high ($N=500$, $J=4$, $K=5$, $L=10$) versions to capture a range of scenarios, where N = no. of individuals, J = no. of tracers, K = no. of food sources, and L = no. of covariates. For each of these we simulate data using the default prior distribution of $\beta_{kl}\sim N(0,1)$ and with $\sigma\sim Ga(1,1)$. We set $\mu_{s}\sim U(-10,10)$, and $\sigma_{s}\sim U(0,2)$. TDFs and Concentration dependence are ignored for this example.

After running each model, we produce posterior uncertainty intervals at the 50% level and calculate the proportion of posterior samples inside these values. As an example, the posterior predictive plot for tracer 1 in the Medium run is shown in Figure 6. Posterior predictive plots for other simulations and tracers are found in the Appendix B (Figures 20, 21 and 22).

The results show that cosimmr produces accurate estimates for the posterior values, even with more complex data and increased numbers of tracers and covariates.

We show an example posterior distribution of $\mathbf{\beta}$ in Figure 7 for the ‘Low’ run. From Figure 7 we can see that cosimmr is performing well and producing $\mathbf{\beta}$ values similar to the original values used to generate data for this example, as illustrated by the red lines visible on the plots. This holds true for the Medium and High examples. Plots are omitted to avoid repetition.

We now perform a direct comparison of cosimmr and MixSIAR to evaluate both the accuracy of the FFVB posterior and check the computational gains. We provide three case studies: the first using the Geese data of Inger et al. (2006), for which we include a single categorical covariate (Group number). The second uses the Isopod data of Galloway et al. (2014) which contains 8 tracers and a single covariate (Site). The third uses the Alligator data of Nifong et al. (2015), for which we provide a detailed model comparison across 8 different potential covariate panels. In each case we compare the posterior distributions of the parameters, the posterior predictive performance, and the computational speed differences. All of the data for our model fits is available in the cosimmr package, available at \urlhttps://github.com/emmagovan/cosimmr and on CRAN, and the code for running the models is available at \urlhttps://github.com/emmagovan/cosimmr_paper.

Fixed Form Variational Bayes is a novel technique within the field of Stable Isotope Mixing Models. It is an optimization-based algorithm which contrasts with the sampling-based approaches traditionally used in SIMMs, such as MCMC. FFVB works by estimating a variational posterior which approximates the true posterior. Through the examples presented in this paper we have demonstrated that it performs as well as MCMC in terms of results produced, while also offering a significant speed improvement of up to one order of magnitude.

The use of FFVB over MCMC allows users to run more complex models in a shorter time. Alternatively it may allow users to compare more models across differing covariates with a view to finding one that matches the data best. We believe that this speed advantage (without loss of accuracy) is an important development for SIMMs. It is important to remember that the FFVB method only ever produces an approximation of the posterior so model checking through, e.g. posterior predictive distributions, is especially important.

We have introduced the package cosimmr to implement these new methods. The inclusion of covariates allows users to avoid violating the assumption of IID data. The package contains functions that allow users to make predictions for combinations of covariates not found or recorded during data collection, allowing for a deeper understanding of the system being studied. We have shown that cosimmr is demonstrably faster than previous packages (due to FFVB), while returning comparable results. We have designed it to be user-friendly for non-expert users, with built-in summary, plotting and predict functions. It is important that users are aware plots are generally created for an individual observation, with the option to specify the individual(s) for which to create plots built into the package functions. Other functions help users to see which of their covariates are having an impact on the diet of the animal being studied, and how well the model fits the data.

Future work on cosimmr could include allowing for random effects, hierarchical modelling, or source fitting - these are all options that are currently available in MixSIAR but could also be implemented using FFVB in order to speed up model fitting. The challenge here is ensuring the optimization still converges satisfactorily despite the additional parameters. We plan to implement these options in a future evolution of the package.

The results in this paper were obtained using R 4.4.1. R is available from the Comprehensive R Archive Network (CRAN) at \urlhttps://CRAN.R-project.org/. C++ compiler gcc version 13.1.0

This publication has emanated from research conducted with the financial support of Science Foundation Ireland under Grant number SFI/12/RC/2289_P2_Parnell to ACP. For the purpose of Open Access, the authors have applied a CC BY public copyright licence to any Author Accepted Manuscript version arising from this submission.

Emma Govan and Andrew C. Parnell conceived the ideas, designed methodology and led the writing of the manuscript; All authors contributed critically to the drafts and gave final approval for publication.