Here we step through the Killer Whale 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_killerwhale.R in the example_scripts folder of the MixSIAR package install:
You can run the killer whale example script directly with:
The Killer Whale Example demonstrates the difference informative priors make, and illustrates how to construct them:
See ?load_mix_data for details.
The killer whale consumer data has no covariates—we are estimating the diet of one population where all individuals are considered identical (factors=NULL, fac_random=NULL, fac_nested=NULL, cont_effects=NULL).
See ?load_source_data for details.
We do not have any fixed/random/continuous effects or concentration dependence in the source data (source_factors=NULL, conc_dep=FALSE). We only have source means, SD, and sample size—not the original “raw” data (data_type="means").
See ?load_discr_data for details.
This is your chance to check:
Notice that the high degree of correlation places the sources on a line, which makes it difficult for the model to determine if killer whales are eating Chinook vs. Coho, or Sockeye vs. Steelhead vs. Chum. We can give the model additional information to resolve this problem by using an informative prior.
# Plot uninformative prior
plot_prior(alpha.prior=1, source, filename = "prior_plot_kw_uninf")
# Define model structure and write JAGS model file
model_filename <- "MixSIAR_model_kw_uninf.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 the JAGS model ("very long" took ~5 min)
jags.uninf <- run_model(run="test",mix,source,discr,model_filename)
# jags.uninf <- run_model(run="very long",mix,source,discr,model_filename)
# Process diagnostics, summary stats, and posterior plots
output_JAGS(jags.uninf, mix, source)From Semmens et al. (in prep):
One of the benefits to conducting mixture models in a Bayesian framework is that information from other data sources can be included via informative prior distributions Moore and Semmens 2008. Suppose we have collected 14 fecal samples from the killer whale population we are studying, and we find 10 chinook, 1 chum, 0 coho, 0 sockeye, and 3 steelhead. The sum of the Dirichlet hyperparameters roughly correspond to prior sample size, so one approach would be to construct a prior with \(\alpha\) = (10, 1, 0, 0, 3). A downside of this prior is that a sample size of 14 represents a very informative prior, with much of the parameter space given very little weight. Keeping the relative contributions the same, the hyperparameters can be rescaled to have the same mean, but different variance. An alternative is to scale the prior to have a weight of 5, equal to the same weight as the “uninformative” prior \(\alpha\) = (1, 1, 1, 1, 1). This prior could be represented as \(\alpha=\left(\frac{10*5}{14},\frac{1*5}{14},\frac{0*5}{14},\frac{0*5}{14},\frac{3*5}{14}\right)=\left(3.57,0.36,0.01,0.01,1.07\right)\).
You cannot have \(\alpha\) parameters = 0, so while we had no coho and sockeye in our fecal samples, we set their α parameters = 0.01. Remember that the sources are sorted alphabetically.
# Our 14 fecal samples were 10, 1, 0, 0, 3
kw.alpha <- c(10,1,0,0,3)
# Generate alpha hyperparameters scaling sum(alpha)=n.sources
kw.alpha <- kw.alpha*length(kw.alpha)/sum(kw.alpha)
# the Dirichlet hyperparameters for the alpha.prior cannot be 0 (but can set = .01)
kw.alpha[which(kw.alpha==0)] <- 0.01
# Plot your informative prior
plot_prior(alpha.prior=kw.alpha,
source=source,
plot_save_pdf=TRUE,
plot_save_png=FALSE,
filename="prior_plot_kw_inf")
# Define model structure and write JAGS model file
model_filename <- "MixSIAR_model_kw_inf.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 the JAGS model ("very long" took ~5 min)
jags.inf <- run_model(run="test",mix,source,discr,model_filename,alpha.prior=kw.alpha)
# jags.inf <- run_model(run="very long",mix,source,discr,model_filename,alpha.prior=kw.alpha)
# Process diagnostics, summary stats, and posterior plots
output_JAGS(jags.inf, mix, source)The plot_prior function shows your prior next to the uninformative prior. Compare the model results with the uninformative prior with those from the informative prior.