Ex 1: Wolves

Load MixSIAR package

library(MixSIAR)

Load mixture data

Load the mixture data, i.e. your:

• Consumer isotope values (trophic ecology / diet)

• Mixed sediment/water tracer values (sediment/hydrology fingerprinting)

  filename: name of the CSV file with mix/consumer data

  iso_names: column headings of the tracers/isotopes you’d like to use

  factors: vector of random/fixed effect column headings in ‘filename’. NULL if no factors.

  fac_random: vector of TRUE/FALSE, TRUE if factor is random effect, FALSE if fixed effect. NULL if no factors.

  fac_nested: vector of TRUE/FALSE, TRUE if factor is nested within the other. Only applies if 2 factors. NULL otherwise.

  cont_effects: column headings of any continuous effects

The wolves consumer data has 2 covariates: Region and Pack, where Pack is nested within Region (fac_nested=c(FALSE,TRUE)). By “nested”, we mean that all wolves in a given pack are in the same region - each pack is entirely within one region. This is an excellent example of hierarchical structure, fit with 2 random effects (fac_random=c(TRUE,TRUE)).

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

# Load the mixture/consumer data
mix <- load_mix_data(filename=mix.filename, 
                     iso_names=c("d13C","d15N"), 
                     factors=c("Region","Pack"), 
                     fac_random=c(TRUE,TRUE), 
                     fac_nested=c(FALSE,TRUE), 
                     cont_effects=NULL)

Load source data

Load the source data, i.e. your:

• Source isotope values (trophic ecology / diet)

• Sediment/water source tracer values (sediment/hydrology fingerprinting)

  filename: name of the CSV file with source data

  source_factors: column headings of random/fixed effects you have source data by

  conc_dep: TRUE or FALSE, do you have concentration dependence data in the file?

  data_type: “means” or “raw”, is your source data as means+SD, or do you have raw data?

If you look at wolves_sources.csv, you will see that each Region has different isotope values - this is specified with source_factors="Region". We do not have concentration dependence data here, so conc_dep=FALSE. We only have source summary statistics (Mean, SD, and sample size), not the original “raw”" data, so data_type="means". Note that wolves_sources.csv has a column titled “n”" with the sample size of each source estimate. This must be in your source data file when you run your data!

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

# Load the source data
source <- load_source_data(filename=source.filename,
                           source_factors="Region", 
                           conc_dep=FALSE, 
                           data_type="means", 
                           mix)

Load discrimination data

Load the discrimination data, i.e. your:

• Trophic Enrichment Factor (TEF) / fractionation values (trophic ecology/diet)

• xxxxxxxx (sediment/hydrology fingerprinting)

  filename: name of the CSV file with discrimination data

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

# Load the discrimination/TDF data
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?

  filename: name you’d like MixSIAR to save the isospace plot as (extension will be added automatically)

  plot_save_pdf: TRUE or FALSE, should MixSIAR save the plot as a .pdf?

  plot_save_png: TRUE or FALSE, should MixSIAR save the plot as a .png?

You should always look at the isospace plot—this is a good check that the data is loaded correctly, and that the isospace geometry makes sense. If the mixture data are well outside the source polygon, you have a serious violation of mixing model assumptions, and it must be true that either 1) You’re missing a source, or 2) You’re using an incorrect discrimination factor. MixSIAR, like SIAR, can fit a residual error term and thus always find a solution even if it is nonsensical.

Also note that the MixSIAR isospace plot adds the discrimination means AND SDs to the raw source values. This is because model uses the source + discrimination values to fit the mixture data, calculated as \(\sqrt{\sigma^2_{source} + \sigma^2_{discr}}\), under the assumption of independence. Error bars indicate \(\pm\) 1 SD.

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

Calculate convex hull area

If 2 isotopes/tracers, calculate normalized surface area of the convex hull polygon(s) as in Brett (2014).

Note 1: discrimination SD is added to the source SD (see ?calc_area for details)

Note 2: If source data are by factor (as in wolf ex), computes area for each polygon (one for each of 3 regions in wolf ex)

# Calculate the convex hull area, standardized by source variance
calc_area(source=source,mix=mix,discr=discr)

## [1]  0.8750158 13.9726701 13.6440547

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)

Bayesian analyses require priors, and MixSIAR includes a plot_prior function to plot the prior on the mixture (diet) proportions (at the highest hierarchical level, p.global). The prior represents our knowledge about the proportions before we consider the biotracer data. A natural tendency is to want a flat/“uninformative” prior, where all values between 0 and 1 are equally likely. However, because proportions are not independent, there is no truly uninformative prior (e.g. the histograms are not flat). The best we can do with the Dirichlet distribution is set \(\alpha\) = c(1,1,1), which is uninformative on the simplex. In other words, all combinations of the proportions are equally likely. See the section titled “Constructing informative Bayesian priors” in the forthcoming MixSIAR paper.

Because the mean of the “uninformative” prior, \(\alpha\) = c(1,1,1), is \(\frac{1}{n.sources}\), we also call it the generalist prior. This reflects two facts: 1) it is not really uninformative, and 2) for weakly informative data it shifts the posterior towards a generalist diet, \(p_1 = p_2 = p_3 = \frac{1}{3}\). The amount of shift depends on the informativeness of the data (quality, quantity, and polygon geometry).

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

Write JAGS model file

Write the JAGS model file (define model structure). The model will be saved as model_filename (“MixSIAR_model.txt” is default, but you may want to change if you create many different models).

There are 3 error term options available:

1. Residual * Process (resid_err = TRUE, process_err = TRUE)
2. Residual only (resid_err = TRUE, process_err = FALSE)
3. Process only (resid_err = FALSE, process_err = TRUE)

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

Note: If you have only 1 mix datapoint (or 1 mix datapoint per factor level, i.e. region or pack in this example), you have no information about the mixture/consumer variability. In this case, we use the original MixSIR error model (which does not fit a residual error term). This is the same behavior as siarsolo in SIAR.

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

Run model

Choose one of the MCMC run options:

You can also set custom MCMC parameters, e.g:

run <- list(chainLength=200000, burn=150000, thin=50, chains=3, calcDIC=TRUE)

Good idea to use run = "test" first 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)

jags.1 will be an rjags object where you can access the MCMC chains for plotting, aggregating sources a posteriori, etc.

Analyze diagnostics and output

First you can set output options like file names, plot file types, etc. (see ?output_JAGS for details).

output_options <- list(summary_save = TRUE,
                       summary_name = "summary_statistics",
                       sup_post = FALSE,
                       plot_post_save_pdf = TRUE,
                       plot_post_name = "posterior_density",
                       sup_pairs = FALSE,
                       plot_pairs_save_pdf = TRUE,
                       plot_pairs_name = "pairs_plot",
                       sup_xy = TRUE,
                       plot_xy_save_pdf = FALSE,
                       plot_xy_name = "xy_plot",
                       gelman = TRUE,
                       heidel = FALSE,
                       geweke = TRUE,
                       diag_save = TRUE,
                       diag_name = "diagnostics",
                       indiv_effect = FALSE,
                       plot_post_save_png = FALSE,
                       plot_pairs_save_png = FALSE,
                       plot_xy_save_png = FALSE,
                       diag_save_ggmcmc = FALSE)

Then you can call output_JAGS to process diagnostics, summary statistics, and create posterior density plots:

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

For a thorough explanation of the output from output_JAGS, see the Wolves Example section of the MixSIAR Manual. You will also find examples of accessing the MCMC chains for post hoc plotting and analysis there.