Ex 2: Geese

Load MixSIAR package

library(MixSIAR)

Load mixture data

See ?load_mix_data for details.

The geese consumer data has 1 covariate (factors="Group"), which we fit as a fixed effect (fac_random=FALSE). We choose to treat Group as a fixed effect instead of a random effect here because we are interested in the diet of each group separately and NOT in the overall diet. “Group” 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", "geese_consumer.csv", package = "MixSIAR")

mix <- load_mix_data(filename=mix.filename,
                     iso_names=c("d13C","d15N"),
                     factors="Group",
                     fac_random=FALSE,
                     fac_nested=FALSE,
                     cont_effects=NULL)

Load source data

See ?load_source_data for details.

If you look at geese_sources.csv, you will see that our geese source data are not by Group (source_factors=NULL), but we DO have concentration dependence data (conc_dep=TRUE). We only have source means, SD, and sample size—not the original “raw” (data_type="means").

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

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

Load discrimination data

See ?load_discr_data for details.

# Replace the system.file call with the path to your file
discr.filename <- system.file("extdata", "geese_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?

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

Calculate convex hull area

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)

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

## [1] 20.27097

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 Geese 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

Choose one of the MCMC run options:

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, 
                    alpha.prior = 1, resid_err, process_err)

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

jags.1 <- run_model(run="short", mix, source, discr, model_filename,
                    alpha.prior = 1, resid_err, process_err)

Analyze diagnostics and output

See ?output_JAGS for details.

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

Note that there is no global/overall estimated diet—this is because we fit Group as a fixed effect instead of a random effect.