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:

mixsiar.dir <- find.package("MixSIAR")

You can run the geese example script directly with:


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):

Load MixSIAR package

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).

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").

Load discrimination data

See ?load_discr_data for details.

Plot data

This is your chance to check:

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)

## [1] 20.27097

Plot prior

Define your prior, and then plot using “plot_prior”

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).

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:

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

Analyze diagnostics and output

See ?output_JAGS for details.

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