Access to posterior output

[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!

[brianstock]

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 <#365>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABHDA46Y3DGFVHN6GPCFPJ3YNOTRFAVCNFSM6AAAAABBRBWEI2VHI2DSMVQWIX3LMV43ASLTON2WKOZSGA3DSOJUGA3DEMQ> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[StudentSaskia]

Thanks for the fast response! Is there a way to match the 3000 draws to each prey and consumer? Or are the only statistics to work with the summary?