Analyzing mixing systems using a new generation of Bayesian tracer mixing models

The ongoing evolution of tracer mixing models has resulted in a confusing array of software tools that differ in terms of data inputs, model assumptions, and associated analytic products. Here we introduce MixSIAR, an inclusive, rich, and flexible Bayesian tracer (e.g., stable isotope) mixing model framework implemented as an open-source R package. Using MixSIAR as a foundation, we provide guidance for the implementation of mixing model analyses. We begin by outlining the practical differences between mixture data error structure formulations and relate these error structures to common mixing model study designs in ecology. Because Bayesian mixing models afford the option to specify informative priors on source proportion contributions, we outline methods for establishing prior distributions and discuss the influence of prior specification on model outputs. We also discuss the options available for source data inputs (raw data versus summary statistics) and provide guidance for combining sources. We then describe a key advantage of MixSIAR over previous mixing model software—the ability to include fixed and random effects as covariates explaining variability in mixture proportions and calculate relative support for multiple models via information criteria. We present a case study of Alligator mississippiensis diet partitioning to demonstrate the power of this approach. Finally, we conclude with a discussion of limitations to mixing model applications. Through MixSIAR, we have consolidated the disparate array of mixing model tools into a single platform, diversified the set of available parameterizations, and provided developers a platform upon which to continue improving mixing model analyses in the future.


Preamble
MixSIAR incorporates many advances to mixing models since previous software packages (IsoSource, MixSIR, SIAR, IsotopeR) were written. While these advances were published elsewhere (Table 1), the novelty of MixSIAR is the integration of these into one framework. MixSIAR accomplishes this via the write_JAGS_model function, which constructs a JAGS model file (MixSIAR_model.txt) given the user's data structure and desired model options. This is not trivial because the model components interact. For example, the choice to include "Individual" as a random effect impacts the choice of error structure on the mixture. The type of source data (raw vs. mean/SD/n) determines whether you can fit the source data with covariance, and if this propagates through to the mixture covariance. As such, we describe the MixSIAR equations in a modular format. For the same reason, we also recommend that users save the JAGS model file and include it as a supplement to any publication, as this file specifies the equations used (there is no one "MixSIAR model").
Throughout, we use bold font to indicate multivariate variables (vectors, matrices), capital letters to indicate the total number of tracers (datapoints, sources, etc.), and lower case letters as indices. For instance, p is a vector of proportions of length K (total number of sources), and p k indicates the proportion of the kth source.

Source data
Unlike MixSIR and SIAR, MixSIAR fits the source data hierarchically (also referred to as "fully Bayesian"). In other words, the model admits that the source means/SDs come from a sample and are not the truth. Therefore, the souce means/SDs used in the mixture likelihood, µ s jk , are allowed to deviate from the source sample means/SDs (and the amount of deviation depends on the source sample sizes; more deviation allowed for lower sample sizes). Note that in systems with lots of mixture data and few source data, the likelihood may be maximized by fitting the source means far from their sample means.
In order of preference (and model complexity), MixSIAR gives users three options to fit source data: 1. If a user has raw source data, MixSIAR includes covariance between tracers (preferred because it is the most complete "fully Bayesian" model). 2. If a user does not have raw source data and only provides summary statistics (mean, SD, and sample size), MixSIAR must assume tracers are independent (no covariance). For large datasets (500+ mixture data points), switching to mean/SD/n can reduce model runtime. This model is also "fully Bayesian" because it also estimates the 'true' source means and variances used in the mixture likelihood. 3. Users can effectively turn off source fitting by using the mean/SD/n option and changing the source sample size to an arbitrarily large number (i.e. set n = 10000). This will fix the source means at their sample means. In poorly resolved mixing systems (many sources, large SD, low n, few tracers), this can help the model converge. The resulting model is not "fully Bayesian", matching previous mixing model software (MixSIR, SIAR).

Raw source data
The data for tracer j, source k, Y s jk , are fit hierarchically as in Hopkins and Ferguson (2012) and Parnell et al. (2013): The fitted source means, µ s jk , are used to construct the mixture likelihood. The priors on the source means are: The source covariance matrices, Σ s k , are constructed as: Option #1 should not be used for compositional tracer (e.g. fatty acid profile, FA) data, because the source data will not conform to the normality assumption. Therefore, we advise users with FA data to use option #2 described below. This is ok because the observed mixture tracer values are a (weighted) sum of random variables, which should be normally distributed according to the Central Limit Theorem (CLT). While FA data are not independent since each tracer must sum to 1, they are likely in the "sufficiently weakly correlated" class such that the CLT still holds. Article S2 demonstrates this via simulations.
Note: The diffuse Normal and Inv-Gamma priors on the source means, µ s jk , and variances, ω s jk 2 , work well for stable isotope and fatty acid tracers, because their order of magnitude is 10 −1 -10 1 . For other tracer types that are of larger orders of magnitude (e.g. element concentrations used in sediment mixing can be 10 3 -10 5 , Nosrati et al. 2014), these priors would not work. Instead of using the data to set the prior (i.e. by setting the prior mean equal to the sample mean), we scale the data so the same prior can be used regardless of data scale.
MixSIAR normalizes (subtract mean, divide by SD) the mixture and source tracer data before running the model. Normalizing the tracer data does not affect the proportion estimates, p k , but does affect users seeking to plot the posterior predictive distribution for their data. For each tracer, we calculate the pooled mean and SD of the mix and source data, then subtract these pooled means and SDs from the mix and source data, and divide by the pooled SD. See lines 226-269 of run_model.R.

Mean / SD / n source data
In the event that a user does not have raw source data and only provides summary statistics (mean, variance, and sample size), we cannot fit the above model with covariance. Instead, we fit the source parameters µ s jk and Σ s k as in Ward et al. (2010): where: m jk = tracer j sample mean for source k (data), s 2 jk = tracer j sample variance for source k (data), n k = source k sample size (data), µ s jk = tracer j mean for source k (parameter), ω s jk 2 = tracer j variance for source k (parameter), Σ s k = source k covariance matrix (calculated from ω s jk terms). Then, µ s jk and Σ s k are used in the mixture likelihood as for raw source data.

Source data by factor
If a user includes a factor as a fixed or random effect on the mixture data, MixSIAR allows the user to also include this factor on the source data. In this case, the source data for each factor level are fit independently and used in the mixture likelihood for mix datapoints in the same factor level. For example, in the wolves dataset (Semmens et al. 2009), the source data from Region 1 is used to fit the mixture data from Region 1, etc.

Fixed / random / continuous effects
In all cases, the overall ("global") source proportions, p k , are drawn from: where p and α are vectors of length K, the number of sources. By default, the "uninformative" prior, α k = 1 for all k, is used. Users can specify their own informative prior using the run_model function. We then transform the global source proportions into ILR-space parameters, β 0 (now a vector of length k − 1), following Egozcue (2003): MixSIAR fits both fixed and random effects as offsets from the overall intercepts, β 0 , in ILR-space.

Fixed effects
The offset for the first level of a factor is set to 0 to avoid identifiability issues (i.e. first level, β 1 k (1), becomes the intercept, β 0 ): Then for the remaining L − 1 factor levels, the offset for the kth source for level l receives the prior: To get the proportion vector for the ith mixture, p i , we add the offsets for the factor level corresponding to mixture i, β 1 (l i ), to the intercept, and then back-transform into p-space using the inverse ILR: A second fixed effect can be added in the same way: With two fixed effects the intercept corresponds to the first level of factor 1 and the first level of factor 2.

Random effects
Random effects are added in much the same way, as offsets from the overall proportions with mean = 0:

Choosing between fixed and random effects
MixSIAR fits both fixed and random effects as offsets from the overall intercepts, β 0 , in ILR-space. We recognize that the terms "fixed" and "random" effects are unclear (Gelman 2005), and in Gelman's constant versus varying terminology, both fixed and random effects in MixSIAR are varying (different for each factor level l). In MixSIAR, for a categorical factor with L levels in a model with K sources:

Fixed effect Random effect
(Effective) parameters added (L − 1)(K − 1) < L(K − 1) + 1 > K Relationship between levels Independent Hierarchical Equations There are two practical distinctions between these models: 1. Number of (effective) parameters added: the fixed effects version generally has more, but it depends on the dataset. The number of parameters will be greater in the random effects version, but the effective number of parameters is often lower because they share information (Gelman et al. 2004). The random effects model adds one parameter for each factor, γ 1 , but the effective number of β 1 k is between L (one for each level) and 1 (mean for all levels). Using the wolves example (Semmens et al. 2009), estimating Pack offsets as fixed effects results in (8 − 1)(3 − 1) = 14 additional parameters, while the random effects model adds somewhere between 3 and 8(3 − 1) + 1 = 17. Thus, if the factor has many levels, it may be better to use the random effects model. If the factor has few levels (< 5), however, it can be difficult to estimate the random effect variance term (γ 2 1 ). If the factor has only 2 levels (e.g. Sex), γ 2 1 cannot be estimated and the factor should be treated as a fixed effect.
2. Independence: the random effects model draws offsets from a shared distribution, which makes sense if the factor levels are related. Since hierarchical structure is common in biological systems, random effects often make sense. If the factor levels are truly independent, then treating the factor as a fixed effect may be best.

Continuous effects
We add continuous effects as linear terms in ILR-space. Before model fitting, we standardize the covariate (subtract mean and divide by standard deviation, see load_mix_data).
To get the source proportions for mixture i, we multiply the linear terms, β 1 , by the value of the covariate for mixture i, x i :

Nested vs. non-nested factors
When two categorial factors are included, the user must tell MixSIAR whether the factors are independent or nested within each other. This affects the calculation of proportions at the factor level. To use the wolves example, each Pack (β 2 ) is nested within a Region (β 1 ). Therefore, to calculate the proportions for pack l: where β 1 (l) is the Region offset where Pack l is found, and β 2 (l) is the Pack offset from the Region mean.
If we change this example such that factor 2 is not nested within factor 1 (e.g. Species that are found in each Region, Region and Species treated independently), then MixSIAR calculates the proportions for the lth level of factor 2 without adding the offset for factor 1:

Mixture mean
MixSIAR assumes mass balance and calculates the mixture mean for each datapoint, µ m ij , as a convex combination of the source proportions, p i , times the fitted source means, µ s jk , adjusted by the mean TDF, λ jk , and concentration of tracer j in source k, q jk : This is the same as previous mixing models (IsoSource, MixSIR, SIAR, IsotopeR, etc.)

Mixture variance
There are three user options for the mixture variance. For motivation, description, and simulation test results, see Stock and Semmens (2016). In brief: 1. Process x Residual error (default) a. with covariance (if raw source data) b. without covariance (if source means/SD/n)

Residual error
Situations where the true variation in source tracer values is not reflected in the mixture data (e.g. integrated sampling), or it does not make sense to think of the consumers (mixtures) sampling individual prey (source) items (e.g. oysters filter feeding, sediment fingerprinting). For large datasets (500+ mixture data points), switching to the residual error model can reduce model runtime.

Process error
Required when there is only one mixture datapoint (not possible to fit mixture variance term), or one mixture datapoint per fixed/random effect level. Including "ID" or "Individual" as a fixed/random effect is an important example where this applies.

Process x Residual error (no covariance)
When there is no information on the covariance of the sources (i.e. user inputs mean/SD/n data), the off-diagonal entries of Σ s k are 0. Then the mixture variance is: E.g. for j=2,

Residual error only
Sometimes the true variation in source tracer values is not reflected in the mixture data (e.g. integrated sampling), or it does not make sense to think of the consumers (mixtures) sampling individual prey (source) items (e.g. oysters filter feeding, sediment fingerprinting). In these cases, the source and mixture data do not follow the standard assumptions about the mixing process (see discussion of "process" vs. "residual" error in Stock and Semmens (2016)). Since there is no information about the source variance, MixSIAR directly fits the mixture variance, Σ, as a residual error term with a Wishart prior: where I is the identity matrix, and j is the number of tracers.

Process error only
When there is only one mixture datapoint (or one mixture datapoint per fixed effect level), it is not possible to fit a mixture variance term. In order to define the likelihood of the one mixture datapoint, MixSIAR assumes that the mixture variance is defined by the proportions and the source variances (i.e. the distribution of mixture tracer data is the mathematical result of adding k independent normal random variables, the sources):

Mixture likelihood
Once the mixture mean and covariances are calculated as above, the likelihood for the data from consumer i and tracer j, Y ij is: Y ij ∼ N ormal µ m ij , Σ i  (2002) Residual error Unexplained variability of mixture data Parnell et al. (2010)