MixSIAR icon indicating copy to clipboard operation
MixSIAR copied to clipboard

posterior extraction of group level means from combine_sources

Open klato1 opened this issue 3 years ago • 3 comments

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?

klato1 avatar Aug 17 '21 22:08 klato1

Please post/email your code and saved model objects, otherwise it's hard to help.

brianstock avatar Aug 17 '21 23:08 brianstock

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)]

klato1 avatar Aug 17 '21 23:08 klato1

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)

klato1 avatar Aug 17 '21 23:08 klato1