Ex 5: Isopod

Here we step through the Isopod Example using the script version of MixSIAR. For a demonstration using the GUI version, see the MixSIAR Manual. For a thorough walkthrough of how to use MixSIAR in a script, see the Wolves Example, which provides more commentary and explanation.

For a clean, runnable .R script, look at mixsiar_script_isopod.R in the example_scripts folder of the MixSIAR package install:

Isopod Example

The Isopod Example is from Galloway et al. 2014 and demonstrates MixSIAR applied to an 8-dimensional fatty acid dataset. Here the mixture data are isopod polyunsaturated fatty acid (PUFA) profiles, with 5 replicates at each of 6 sites in Puget Sound, WA:

8 biotracers (carbon 16.4w3, 18.2w6, 18.3w3, 18.4w3, 20.4w6, 20.5w3, 22.5w3, 22.6w3)
1 random effect (Site)
source data as means and SDs

Here we treat Site as a random effect. This makes sense if we are interested in the overall population and think of Site as a nuisance factor. Fitting Site as a fixed effect would make more sense if we were interested specifically in the diet at each Site, as opposed to the overall population diet and variability between Sites. This differs from the analysis in Galloway et al. 2014.

Fatty acid data greatly increase the number of biotracers beyond the typical 2 stable isotopes, d13C and d15N, which gives the mixing model power to resolve more sources. We caution, however, that using fatty acid data is not a panacea for the “underdetermined” problem (# sources > # biotracers + 1). As the number of sources increases, the “uninformative” prior \((\alpha=1)\) has greater influence, even if there are many more biotracers than sources. See the Cladocera Example prior with 7 sources and 22 biotracers.

Load MixSIAR package

library(MixSIAR)

Load mixture data

See ?load_mix_data for details.

The isopod consumer data has 1 covariate (factors="Site"), which we fit as a random effect (fac_random=TRUE). “Site” is not nested within another factor (fac_nested=FALSE). There are no continuous effects (cont_effects=NULL).

# Replace the system.file call with the path to your file
mix.filename <- system.file("extdata", "isopod_consumer.csv", package = "MixSIAR")

mix <- load_mix_data(filename=mix.filename,
                     iso_names=c("c16.4w3","c18.2w6","c18.3w3","c18.4w3","c20.4w6","c20.5w3","c22.5w3","c22.6w3"),
                     factors="Site",
                     fac_random=TRUE,
                     fac_nested=FALSE,
                     cont_effects=NULL)

Load source data

See ?load_source_data for details.

We do not have any fixed/random/continuous effects or concentration dependence in the source data (source_factors=NULL, conc_dep=FALSE). We only have source means, SD, and sample size—not the original “raw” data (data_type="means").

# Replace the system.file call with the path to your file
source.filename <- system.file("extdata", "isopod_sources.csv", package = "MixSIAR")

source <- load_source_data(filename=source.filename,
                           source_factors=NULL,
                           conc_dep=FALSE,
                           data_type="means",
                           mix)

Load discrimination data

See ?load_discr_data for details.

Note that Galloway et al. 2014 conducted feeding trials to create a “resource library”. In the mixing model, the sources are actually consumers fed exclusively each of the sources. This allowed them to set the discrimination = 0 (see isopod_discrimination.csv).

# Replace the system.file call with the path to your file
discr.filename <- system.file("extdata", "isopod_discrimination.csv", package = "MixSIAR")

discr <- load_discr_data(filename=discr.filename, mix)

Plot data

This is your chance to check:

Are the data loaded correctly?
Is your mixture data in the source polygon?
Are one or more of your sources confounded/hidden?

When there are more than 2 biotracers, MixSIAR currently plots every pairwise combination. Here, that means \({8 \choose 2} = 28\) plots are produced. In the future, MixSIAR will offer non-metric multidimensional scaling (NMDS) plots for these cases.

# Make an isospace plot
plot_data(filename="isospace_plot", plot_save_pdf=TRUE, plot_save_png=FALSE, mix,source,discr)

Plot prior

Define your prior, and then plot using “plot_prior”

RED = your prior
DARK GREY = “uninformative”/generalist (alpha = 1)
LIGHT GREY = “uninformative” Jeffrey’s prior (alpha = 1/n.sources)

# default "UNINFORMATIVE" / GENERALIST prior (alpha = 1)
plot_prior(alpha.prior=1,source)

Write JAGS model file

In the Isopod Example we demo the “Residual only” error option. The differences between “Residual * Process”, “Residual only”, and “Process only” are explained in Stock and Semmens (2016).

# Write the JAGS model file
model_filename <- "MixSIAR_model.txt"
resid_err <- TRUE
process_err <- FALSE
write_JAGS_model(model_filename, resid_err, process_err, mix, source)

Run model

First use run = "test" to check if 1) the data are loaded correctly and 2) the model is specified correctly:

jags.1 <- run_model(run="test", mix, source, discr, model_filename)

After a test run works, increase the MCMC run to a value that may converge:

jags.1 <- run_model(run="normal", mix, source, discr, model_filename)

Analyze diagnostics and output

See ?output_JAGS for details.

output_JAGS(jags.1, mix, source, output_options)

Since we fit Site as a random effect, we can make inference on diet at the overall population level (p.global, posterior plot) as well as at individual sites (Site CP, posterior plot).