Ex 2: Geese

Here we step through the Geese 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_geese.R in the example_scripts folder of the MixSIAR package install:

library(MixSIAR)
mixsiar.dir <- find.package("MixSIAR")
paste0(mixsiar.dir,"/example_scripts")

You can run the geese example script directly with:

source(paste0(mixsiar.dir,"/example_scripts/mixsiar_script_geese.R"))

Geese Example

The “Geese Example” uses data from Inger et al. (2006) of 251 wintering geese feeding on terrestrial grasses, Zostera spp., Enteromorpha spp., and Ulva lactuca. This is the same data included as a demo in SIAR and in Parnell et al. (2013):

2 biotracers (δ¹³C, δ¹⁵N)
1 fixed effect (Group)
Source data as means and SDs
Concentration dependence
Residual only error

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:

run ==	Chain Length	Burn-in	Thin	# Chains
“test”	1,000	500	1	3
“very short”	10,000	5,000	5	3
“short”	50,000	25,000	25	3
“normal”	100,000	50,000	50	3
“long”	300,000	200,000	100	3
“very long”	1,000,000	500,000	500	3
“extreme”	3,000,000	1,500,000	500	3

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.

- Geese Example