brianstock
commented
10 months ago 3000 posterior draws, which you can summarize (mean, median, mode, SD) or plot (eg histogram, density, box plot) as you like. See that when you do these yourself they match the output summary statistics.
On Mon, Jan 8, 2024, 9:17 AM StudentSaskia @.***> wrote:
Hi,
I want to see if the output for my 3 Raja species is significantly different (so does RJM significantly eat more shrimp than RJH for example). However I am having trouble finding the posterior output data. I know it is supposed to be in Jags.1$BUGSoutput$p.global, however this is a matrix of 3000 x 5 in my model. I get that the 5 are the 5 prey species, however where does the 3000 come from?
This is my code: Define the file path as a character string
mix.ray <- "CorrelationLengthSizeClass3.csv" Load the mixture data
mix <- read.csv(mix.ray, header = TRUE) Load the mixture/consumer data
mix <- load_mix_data(filename= mix.ray, iso_names=c("d13C","d15N"), factors=c("Sex","sclass"), fac_random=c(FALSE,FALSE), fac_nested=c(FALSE,FALSE), cont_effects="Length..cm.") Replace the system.file call with the path to your file
source.ray <- "ray_sources3.csv" source <- read.csv(source.ray, header = TRUE) Load the source data
source <- load_source_data(filename=source.ray, source_factors=NULL, conc_dep=FALSE, data_type = "means", mix) load discrimination factors
discr.ray <- "ray_discrimination3.csv" discr <- read.csv(discr.ray,header = TRUE) discr <- load_discr_data(filename="ray_discrimination3.csv", mix)
make isospace plot
plot_data(filename="isospace_plot", plot_save_pdf=TRUE, plot_save_png=FALSE, mix,source,discr)
calculate normalized surface area
if(mix$n.iso==2) calc_area(source=source,mix=mix,discr=discr)
plot informative prior
kw.alpha <- c(58,68, 24, 3, 72) kw.alpha <- kw.alpha*length(kw.alpha)/138 kw.alpha[which(kw.alpha==0)] <- 0.01 plot_prior(alpha.prior=kw.alpha,source=source,plot_save_pdf=TRUE, plot_save_png=FALSE,filename="prior_plot")
write JAGS model
?write_JAGS_model model_filename <- "MixSIAR_model.txt" resid_err <- TRUE process_err <- FALSE write_JAGS_model(model_filename, resid_err, process_err, mix, source)
Run JAGS model
?run_model jags.1 <- run_model(run="long",mix,source,discr,model_filename, alpha.prior = 1,resid_err,process_err)
maximize output
getOption("max.print") options(max.print = 100000) getOption("max.print")
process output
output_JAGS(jags.1,mix,source, output_options = list(summary_save = TRUE, summary_name = "summary_statistics", sup_post = FALSE, plot_post_save_pdf = FALSE, plot_post_name = "posterior_density", sup_pairs = TRUE, plot_pairs_save_pdf = TRUE, plot_pairs_name = "pairs_plot", sup_xy = TRUE, plot_xy_save_pdf = TRUE, plot_xy_name = "xy_plot", gelman = TRUE, heidel = FALSE, geweke = TRUE, diag_save = TRUE, diag_name = "diagnostics", indiv_effect = FALSE, plot_post_save_png = FALSE, plot_pairs_save_png = FALSE, plot_xy_save_png = TRUE, diag_save_ggmcmc = TRUE))
Thanks in advance!
— Reply to this email directly, view it on GitHub https://github.com/brianstock/MixSIAR/issues/365, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABHDA46Y3DGFVHN6GPCFPJ3YNOTRFAVCNFSM6AAAAABBRBWEI2VHI2DSMVQWIX3LMV43ASLTON2WKOZSGA3DSOJUGA3DEMQ . You are receiving this because you are subscribed to this thread.Message ID: @.***>
StudentSaskia
Hi,
I want to see if the output for my 3 Raja species is significantly different (so does RJM significantly eat more shrimp than RJH for example). However I am having trouble finding the posterior output data. I know it is supposed to be in Jags.1$BUGSoutput$p.global, however this is a matrix of 3000 x 5 in my model. I get that the 5 are the 5 prey species, however where does the 3000 come from?
This is my code:
Define the file path as a character string
mix.ray <- "CorrelationLengthSizeClass3.csv"
Load the mixture data
mix <- read.csv(mix.ray, header = TRUE)
Load the mixture/consumer data
mix <- load_mix_data(filename= mix.ray, iso_names=c("d13C","d15N"), factors=c("Sex","sclass"), fac_random=c(FALSE,FALSE), fac_nested=c(FALSE,FALSE), cont_effects="Length..cm.")
Replace the system.file call with the path to your file
source.ray <- "ray_sources3.csv" source <- read.csv(source.ray, header = TRUE)
Load the source data
source <- load_source_data(filename=source.ray, source_factors=NULL, conc_dep=FALSE, data_type = "means", mix)
load discrimination factors
discr.ray <- "ray_discrimination3.csv" discr <- read.csv(discr.ray,header = TRUE) discr <- load_discr_data(filename="ray_discrimination3.csv", mix)
make isospace plot
plot_data(filename="isospace_plot", plot_save_pdf=TRUE, plot_save_png=FALSE, mix,source,discr)
calculate normalized surface area
if(mix$n.iso==2) calc_area(source=source,mix=mix,discr=discr)
plot informative prior
kw.alpha <- c(58,68, 24, 3, 72) kw.alpha <- kw.alpha*length(kw.alpha)/138 kw.alpha[which(kw.alpha==0)] <- 0.01 plot_prior(alpha.prior=kw.alpha,source=source,plot_save_pdf=TRUE, plot_save_png=FALSE,filename="prior_plot")
write JAGS model
?write_JAGS_model model_filename <- "MixSIAR_model.txt" resid_err <- TRUE process_err <- FALSE write_JAGS_model(model_filename, resid_err, process_err, mix, source)
Run JAGS model
?run_model jags.1 <- run_model(run="long",mix,source,discr,model_filename, alpha.prior = 1,resid_err,process_err)
maximize output
getOption("max.print") options(max.print = 100000) getOption("max.print")
process output
output_JAGS(jags.1,mix,source, output_options = list(summary_save = TRUE, summary_name = "summary_statistics", sup_post = FALSE, plot_post_save_pdf = FALSE, plot_post_name = "posterior_density", sup_pairs = TRUE, plot_pairs_save_pdf = TRUE, plot_pairs_name = "pairs_plot", sup_xy = TRUE, plot_xy_save_pdf = TRUE, plot_xy_name = "xy_plot", gelman = TRUE, heidel = FALSE, geweke = TRUE, diag_save = TRUE, diag_name = "diagnostics", indiv_effect = FALSE, plot_post_save_png = FALSE, plot_pairs_save_png = FALSE, plot_xy_save_png = TRUE, diag_save_ggmcmc = TRUE))
Thanks in advance!