# MixSIAR Script Example (Lake)

#### March 10, 2016

Here we step through the Lake 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_lake.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 lake example script directly with:

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

## Lake Example

The Lake Example data is simulated based on Francis et al. 2011 and we include it as a example of how MixSIAR can include a continuous effect. Here we examine the diet of zooplankton in 21 lakes using:

• 2 biotracers ($$\delta^{13}$$C, $$\delta^{15}$$N)
• 1 continuous effect (Secchi Depth : Mixed Layer Depth)
• Raw source data

MixSIAR fits a continuous covariate as a linear regression in ILR/transform-space. Two terms are fit for the proportion of each source: an intercept and a slope. The plot uses the posterior median estimates of the intercept and slope, and the lines are curved because of the ILR-transform back into p-space. For details, or if you would like to make modifications, see ?plot_continous_var.R.

Fitting a model with a continuous effect is more complex than the categorical fixed/random effects and can be a bit finicky.

### Load MixSIAR package

library(MixSIAR)

### Load mixture data

See ?load_mix_data for details.

The lake consumer data has 1 covariate, which we fit as a continuous effect (cont_effects="Secchi.Mixed"). There are no fixed/random effects (factors=NULL, fac_random=NULL, fac_nested=NULL).

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

iso_names=c("d13C","d15N"),
factors=NULL,
fac_random=NULL,
fac_nested=NULL,
cont_effects="Secchi.Mixed")

### Load source data

See ?load_source_data for details.

We have no fixed/random effects in the model (source_factors=NULL), and we do not have concentration dependence data (conc_dep=FALSE). We have the original “raw” source data, not means and SDs (data_type="raw").

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

source_factors=NULL,
conc_dep=FALSE,
data_type="raw",
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", "lake_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?
# 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).

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)
## [1] 3.053018

### 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)
# default "UNINFORMATIVE" / GENERALIST prior (alpha = 1)
plot_prior(alpha.prior=1,source)

### Write JAGS model file

In the Lake Example we demo the “Residual only” error option. The differences between “Residual * Process”, “Residual only”, and “Process only” are explained in Stock and Semmens (in revision).

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

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:

jags.1 <- run_model(run="normal", mix, source, discr, model_filename,
alpha.prior = 1, resid_err, process_err)

### Analyze diagnostics and output

MixSIAR fits a continuous covariate as a linear regression in ILR/transform-space. Two terms are fit for the proportion of each source: an intercept and a slope. The plot uses the posterior median estimates of the intercept and slope, and the lines are curved because of the ILR-transform back into p-space. For details, or if you would like to make modifications, see ?plot_continous_var. See ?output_JAGS for output options.

The other posterior plots MixSIAR produces for a continuous effect show the estimated diet for the minimum, median, and maximum individuals.

output_JAGS(jags.1, mix, source, output_options)