MixSIAR
MixSIAR copied to clipboard
posterior extraction of group level means from combine_sources
Hi Brian,
I ran a model with species as a fixed effect (2 species/groups) and individual as a random effect. I am really only interested right now in extracting the fixed effect posterior estimates, however when I 'combine_sources', I don't quite understand what the posterior columns are giving me. There are estimates for every species, individual, and source (before combining), and columns 'p.global[1]' and 'p.global[2]', but how can I extract the posterior values for the combined sources for each species?
Please post/email your code and saved model objects, otherwise it's hard to help.
Sure, thanks Brian
mix<-load_mix_data(mix.filename, iso_names=c("d13C","d15N", "d34S"), #these are the stable isotope names in the mixture file factors= c("Species", "TagID"), fac_random=c(FALSE, TRUE), #set whether factors are random or fixed fac_nested=c(FALSE,FALSE),# nested factors (i.e. sex between locations) cont_effects = NULL) # continuous covariates
#load source data source <- load_source_data(filename=source.filename, source_factors=NULL,# if sources match factor in mixture file then call by mix$file conc_dep=TRUE, data_type="means", # means, SD, or sample size mix) #output from load_mix_data
#load discrimination factors discr <- load_discr_data(filename=discr.filename, mix)
model_filename <- "34s_model.txt" # Name of the JAGS model file resid_err <- TRUE process_err <- TRUE write_JAGS_model(model_filename, resid_err, process_err, mix, source) #write the model
species.34S <- run_model(run="very long", mix, source, discr, model_filename)
combined <- MixSIAR::combine_sources(jags.1 = species.34S, mix = mix, source = source, alpha.prior=1, groups=list(marine = c("bivalve", "fish_crust"), urban=c("wheat","fastfoodmeat")))
apply(combined$post, 2, median) df <- as.data.frame(combined$post)]
The column names include: 'deviance', 'fac2.sig', 'ilr.fac1[1,1]', 'ilr.fac1[2,1]', 'ilr.fac1[1,2]'.......'p.global[1]', 'p.global[2]', 'p.both[1,1,1]', 'p.both[1,2,1]'....... I gather that the 'p.both' is the output for each individual, but I am really looking for the group (fixed effect) means, which is the first number in the 'p.both' column title (ie. group = 1 or 2)