---
title: "Ex 4: Killer Whale"
output: html_vignette
vignette: >
  %\VignetteIndexEntry{Ex 4: Killer Whale}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc} 
---

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](https://github.com/brianstock/MixSIAR/blob/master/inst/mixsiar_manual_small.pdf). For a thorough walkthrough of how to use MixSIAR in a script, see the [Wolves Example](https://brianstock.github.io/MixSIAR/articles/wolves_ex.html), 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:
```{r, eval=FALSE}
library(MixSIAR)
mixsiar.dir <- find.package("MixSIAR")
paste0(mixsiar.dir,"/example_scripts")
```

You can run the killer whale example script directly with:
```{r, eval=FALSE}
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

```{r}
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`).

```{r}
# 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"`).

```{r}
# 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.

```{r}
# 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.**

```{r, eval=FALSE}
# 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

```{r, eval=FALSE}
# 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)
```

### 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](https://doi.org/10.1111/j.1461-0248.2008.01163.x). 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.

```{r, eval=FALSE}
# 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](https://github.com/brianstock/MixSIAR/blob/master/Manual/kw_prior_plot.pdf). Compare the model results with the [uninformative prior](https://github.com/brianstock/MixSIAR/blob/master/Manual/kw_uninf_posterior_density_diet_p_global.pdf) with those from the [informative prior](https://github.com/brianstock/MixSIAR/blob/master/Manual/kw_inf_posterior_density_diet_p_global.pdf).