# MixSIAR Script Example (Killer Whale)

#### Brian Stock

#### March 10, 2016

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:

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

    ## [1] "/tmp/Rtmpm2SefP/Rinst74a7f7f1e2c/MixSIAR/example_scripts"

You can run the killer whale example script directly with:

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

## Killer Whale Example

The Killer Whale Example demonstrates the difference informative priors make, and illustrates how to construct them:

- 2 biotracers (\(\delta^{13}\)C, \(\delta^{15}\)N)
- Source data as means + SDs
- Informative vs. uninformative priors

### Load MixSIAR package

    library(MixSIAR)

### Load mixture data

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

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

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

### Load source data

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

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

    source <- load_source_data(filename=source.filename,
                               source_factors=NULL,
                               conc_dep=FALSE,
                               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", "killerwhale_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?

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.

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

### Run the model with uninformative 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,alpha.prior = 1, resid_err, process_err)
    # jags.uninf <- run_model(run="very long",mix,source,discr,model_filename,alpha.prior = 1, resid_err, process_err)

    # Process diagnostics, summary stats, and posterior plots
    output_JAGS(jags.uninf, mix, source)

### Run the model with informative prior

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, resid_err, process_err)
    # jags.inf <- run_model(run="very long",mix,source,discr,model_filename,alpha.prior=kw.alpha, resid_err, process_err)

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