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

Here we step through the Geese 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_geese.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 geese example script directly with:
```{r, eval=FALSE}
source(paste0(mixsiar.dir,"/example_scripts/mixsiar_script_geese.R"))
```

## Geese Example

The "Geese Example" uses data from [Inger et al. (2006)](https://doi.org/10.1111/j.1365-2656.2006.01142.x) 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)](https://doi.org/10.1002/env.2221):

+ 2 biotracers ($\delta^{13}$C, $\delta^{15}$N)
+ 1 **fixed effect** (Group)
+ Source data as means and SDs
+ **Concentration dependence**
+ Residual only error

### Load MixSIAR package

```{r}
library(MixSIAR)
```

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

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

mix <- load_mix_data(filename=mix.filename,
					 iso_names=c("d13C","d15N"),
					 factors="Group",
					 fac_random=FALSE,
					 fac_nested=FALSE,
					 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"`).

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

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

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

### Calculate convex hull area

Calculate normalized surface area of the convex hull polygon(s) as in [Brett (2014)](https://www.int-res.com/articles/suppl/m514p001_supp.pdf).

**Note 1:** discrimination SD is added to the source SD (see ?calc_area for details)

```{r}
# Calculate the convex hull area, standardized by source variance
calc_area(source=source,mix=mix,discr=discr)
```

### Plot prior

Define your prior, and then plot using "plot_prior"

+ RED = your prior
+ DARK GREY = "uninformative"/generalist (alpha = 1)
+ LIGHT GREY = "uninformative" Jeffrey's prior (alpha = 1/n.sources)

```{r, eval=FALSE}
# default "UNINFORMATIVE" / GENERALIST prior (alpha = 1)
plot_prior(alpha.prior=1,source)
```

### 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)](https://doi.org/10.1002/ecy.1517).

```{r, eval=FALSE}
# Write the JAGS model file
model_filename <- "MixSIAR_model.txt"
resid_err <- TRUE
process_err <- FALSE
write_JAGS_model(model_filename, resid_err, process_err, mix, source)
```

### 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:
```{r, eval=FALSE}
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:
```{r, eval=FALSE}
jags.1 <- run_model(run="short", mix, source, discr, model_filename,
                    alpha.prior = 1, resid_err, process_err)
```

### Analyze diagnostics and output

See ?output_JAGS for details.

```{r, eval=FALSE}
output_JAGS(jags.1, mix, source, output_options)
```

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