Merge branch 'master' of https://github.com/andrewcparnell/simms_course

.gitignore

@@ -44,3 +44,5 @@ Src/*.o
 *.aux
 *.log
 *.synctex.gz
+MixSIAR_model.txt
+isospace_plot_1_2.pdf

ap_notes/complex_simms/complex_SIMMs.Rmd

@@ -6,8 +6,17 @@ output:
   beamer_presentation:
     includes:
       in_header: ../header.tex
+editor_options: 
+  chunk_output_type: console
+classoption: "aspectratio=169"
 ---
 
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(dev = 'pdf', fig.height = 5)
+par(mar=c(3,3,2,1), mgp=c(2,.7,0), tck=-.01,las=1)
+```
+
+
 ## Learning outcomes
 
 - Understand how covariates might be included in a SIMM
@@ -99,7 +108,7 @@ model {
 
 ## The Geese data
 
-\includegraphics[width=\textwidth]{Geese2IsoSpace.png}
+\includegraphics[width=0.8\textwidth]{Geese2IsoSpace.png}
 
 ## A Fourier basis
 
@@ -155,6 +164,7 @@ output=coda.samples(model=model,variable.names=c('p','beta','mu_f','sigma_f'),n.
 gelman.diag(output,multivariate = FALSE)
 ```
 
+\small 
 ```{r,eval=FALSE}
 modelstring ='
 model {
@@ -177,7 +187,7 @@ Full script in `run_geese_harmonic.R` file
 
 ## Output
 
-\includegraphics[width=\textwidth]{Geese2Harmonic.png}
+\includegraphics[width=0.8\textwidth]{Geese2Harmonic.png}
 
 ## Comparing models
 

extra_code/mixsiar_examples/cladocera.R

@@ -0,0 +1,75 @@
+## ---- eval=FALSE--------------------------------------------------------------
+#  library(MixSIAR)
+#  mixsiar.dir <- find.package("MixSIAR")
+#  paste0(mixsiar.dir,"/example_scripts")
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  source(paste0(mixsiar.dir,"/example_scripts/mixsiar_script_cladocera.R"))
+
+## -----------------------------------------------------------------------------
+library(MixSIAR)
+library(tidyr) # For pivoting later in new output_jags
+library(ggplot2)
+library(GGally)
+
+# My new output function
+source("extra_code/output_jags_new.R")
+
+## -----------------------------------------------------------------------------
+# Replace the system.file call with the path to your file
+mix.filename <- system.file("extdata", "cladocera_consumer.csv", package = "MixSIAR")
+
+mix <- load_mix_data(filename=mix.filename,
+                     iso_names=c("c14.0","c16.0","c16.1w9","c16.1w7","c16.2w4",
+                                 "c16.3w3","c16.4w3","c17.0","c18.0","c18.1w9",
+                                 "c18.1w7","c18.2w6","c18.3w6","c18.3w3","c18.4w3",
+                                 "c18.5w3","c20.0","c22.0","c20.4w6","c20.5w3",
+                                 "c22.6w3","BrFA"),
+                     factors="id",
+                     fac_random=FALSE,
+                     fac_nested=FALSE,
+                     cont_effects=NULL)
+
+## -----------------------------------------------------------------------------
+# Replace the system.file call with the path to your file
+source.filename <- system.file("extdata", "cladocera_sources.csv", package = "MixSIAR")
+
+source <- load_source_data(filename=source.filename,
+                           source_factors=NULL,
+                           conc_dep=FALSE,
+                           data_type="means",
+                           mix)
+
+## -----------------------------------------------------------------------------
+# Replace the system.file call with the path to your file
+discr.filename <- system.file("extdata", "cladocera_discrimination.csv", package = "MixSIAR")
+
+discr <- load_discr_data(filename=discr.filename, mix)
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  # default "UNINFORMATIVE" / GENERALIST prior (alpha = 1)
+#  plot_prior(alpha.prior=1,source)
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  # Write the JAGS model file
+ model_filename <- "MixSIAR_model.txt"
+ resid_err <- FALSE
+ process_err <- TRUE
+ write_JAGS_model(model_filename, resid_err, process_err, mix, source)
+
+## ---- eval=FALSE--------------------------------------------------------------
+ jags.1 <- run_model(run="test", mix, source, discr, model_filename)
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  jags.1 <- run_model(run="normal", mix, source, discr, model_filename)
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  output_JAGS(jags.1, mix, source, output_options)
+ output_JAGS(jags.1, mix = mix, source = source,
+             c('summary_diagnostics',
+               'summary_statistics',
+               'summary_quantiles',
+               'plot_global',
+               'plot_global_matrix',
+               'plot_factors'))
+

extra_code/mixsiar_examples/geese.R

@@ -0,0 +1,84 @@
+## ---- eval=FALSE--------------------------------------------------------------
+#  library(MixSIAR)
+#  mixsiar.dir <- find.package("MixSIAR")
+#  paste0(mixsiar.dir,"/example_scripts")
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  source(paste0(mixsiar.dir,"/example_scripts/mixsiar_script_geese.R"))
+
+## -----------------------------------------------------------------------------
+library(MixSIAR)
+library(tidyr) # For pivoting later in new output_jags
+library(ggplot2)
+library(GGally)
+
+# My new output function
+source("extra_code/output_jags_new.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)
+
+## -----------------------------------------------------------------------------
+# 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)
+
+## -----------------------------------------------------------------------------
+# 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)
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  # Make an isospace plot
+#  plot_data(filename="isospace_plot", plot_save_pdf=TRUE, plot_save_png=FALSE, mix,source,discr)
+
+## -----------------------------------------------------------------------------
+# Calculate the convex hull area, standardized by source variance
+calc_area(source=source,mix=mix,discr=discr)
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  # default "UNINFORMATIVE" / GENERALIST prior (alpha = 1)
+#  plot_prior(alpha.prior=1,source)
+
+## ---- 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)
+
+## ---- eval=FALSE--------------------------------------------------------------
+ jags.1 <- run_model(run="test", mix, source, discr, model_filename,
+                     alpha.prior = 1, resid_err, process_err)
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  jags.1 <- run_model(run="short", mix, source, discr, model_filename,
+#                      alpha.prior = 1, resid_err, process_err)
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  output_JAGS(jags.1, mix, source, output_options)
+my_output <- output_JAGS(jags.1, mix = mix, source = source,
+             c('summary_diagnostics',
+               'summary_statistics',
+               'summary_quantiles',
+               'plot_global',
+               'plot_global_matrix',
+               'plot_factors'))
+
+# To get the full list of parameters use:
+ rownames(my_output$summary_statistics)
+

extra_code/mixsiar_examples/isopod.R

@@ -0,0 +1,73 @@
+## ---- eval=FALSE--------------------------------------------------------------
+#  library(MixSIAR)
+#  mixsiar.dir <- find.package("MixSIAR")
+#  paste0(mixsiar.dir,"/example_scripts")
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  source(paste0(mixsiar.dir,"/example_scripts/mixsiar_script_isopod.R"))
+
+## -----------------------------------------------------------------------------
+library(MixSIAR)
+library(tidyr) # For pivoting later in new output_jags
+library(ggplot2)
+library(GGally)
+
+# My new output function
+source("extra_code/output_jags_new.R")
+
+## -----------------------------------------------------------------------------
+# Replace the system.file call with the path to your file
+mix.filename <- system.file("extdata", "isopod_consumer.csv", package = "MixSIAR")
+
+mix <- load_mix_data(filename=mix.filename,
+                     iso_names=c("c16.4w3","c18.2w6","c18.3w3","c18.4w3","c20.4w6","c20.5w3","c22.5w3","c22.6w3"),
+                     factors="Site",
+                     fac_random=TRUE,
+                     fac_nested=FALSE,
+                     cont_effects=NULL)
+
+## -----------------------------------------------------------------------------
+# Replace the system.file call with the path to your file
+source.filename <- system.file("extdata", "isopod_sources.csv", package = "MixSIAR")
+
+source <- load_source_data(filename=source.filename,
+                           source_factors=NULL,
+                           conc_dep=FALSE,
+                           data_type="means",
+                           mix)
+
+## -----------------------------------------------------------------------------
+# Replace the system.file call with the path to your file
+discr.filename <- system.file("extdata", "isopod_discrimination.csv", package = "MixSIAR")
+
+discr <- load_discr_data(filename=discr.filename, mix)
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  # Make an isospace plot
+#  plot_data(filename="isospace_plot", plot_save_pdf=TRUE, plot_save_png=FALSE, mix,source,discr)
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  # default "UNINFORMATIVE" / GENERALIST prior (alpha = 1)
+#  plot_prior(alpha.prior=1,source)
+
+## ---- 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)
+
+## ---- eval=FALSE--------------------------------------------------------------
+ jags.1 <- run_model(run="test", mix, source, discr, model_filename)
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  jags.1 <- run_model(run="normal", mix, source, discr, model_filename)
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  output_JAGS(jags.1, mix, source, output_options)
+ output_JAGS(jags.1, mix = mix, source = source,
+             c('summary_diagnostics',
+               'summary_statistics',
+               'summary_quantiles',
+               'plot_global',
+               'plot_factors'))

extra_code/mixsiar_examples/killerwhale.R

@@ -0,0 +1,107 @@
+## ---- eval=FALSE--------------------------------------------------------------
+#  library(MixSIAR)
+#  mixsiar.dir <- find.package("MixSIAR")
+#  paste0(mixsiar.dir,"/example_scripts")
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  source(paste0(mixsiar.dir,"/example_scripts/mixsiar_script_killerwhale.R"))
+
+## -----------------------------------------------------------------------------
+library(MixSIAR)
+library(tidyr) # For pivoting later in new output_jags
+library(ggplot2)
+library(GGally)
+
+# My new output function
+source("extra_code/output_jags_new.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)
+
+## -----------------------------------------------------------------------------
+# 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)
+
+## -----------------------------------------------------------------------------
+# 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)
+
+## ---- eval=FALSE--------------------------------------------------------------
+#  # Make an isospace plot
+#  plot_data(filename="isospace_plot", plot_save_pdf=TRUE, plot_save_png=FALSE, mix,source,discr)
+
+## ---- 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)
+ output_JAGS(jags.uninf, mix = mix, source = source,
+             c('summary_diagnostics',
+               'summary_statistics',
+               'summary_quantiles',
+               'plot_global_matrix',
+               'plot_global'))
+
+
+## ---- 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)
+output_JAGS(jags.inf, mix = mix, source = source,
+            c('summary_diagnostics',
+              'summary_statistics',
+              'summary_quantiles',
+              'plot_global_matrix',
+              'plot_global'))