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:
You can run the geese example script directly with:
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):
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 (
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” (
See ?load_discr_data for details.
This is your chance to check:
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)
##  20.27097
Define your prior, and then plot using “plot_prior”
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).
Choose one of the MCMC run options:
|run ==||Chain Length||Burn-in||Thin||# Chains|
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:
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.