MixSIAR
MixSIAR copied to clipboard
Published
20 hours ago •
brianstock
brianstock
Producing posterior plots with continuous and two fixed effects
Hello Brian and MixSIAR users, I am estimating diet proportions for an animal using two isotopes (13C and 15N) and 4 diet sources. Ideally, given the geography and biology of the population, I'd like use a model with 2 fixed factors ('forage.zone' and 'lact.stat') and account for individual differences in body condition with a continuous effect 'ifbf.per'. I have two questions: (1) When I run the model with this structure (2 fixed effects and 1 continuous effect), I receive the following warning: "Failed to set trace monitor for ilr.global..." I don't know how to interpret that. (2) I would like to produce a figure like Figures 5 or 6 in the Peer J paper, but while accounting for the fixed and random effects. So, in this case, produce posterior distributions for each combination of 'lact.stat' x 'forage.zone' as a function of 'ifbf.per'. Included in code below is Brian's code for producing the dataframe used to generate Figure 6, but adapted for this case. It will not, however, account for the interacting fixed factors, and I don't know how to access that information if it exists in the jags output. Perhaps I have an incorrect assumption about the model structure or what to expect.
I've never posted a reprex, so here's a first attempt. I will also add that I know this not a well-resolved mixing system, but I would like to learn how to do this anyway:
#----------------------------------------------------------------------
library(MixSIAR)
library(wrapr)
#bring in data
d_mix <- wrapr::build_frame(
"id" , "d13C", "d15N", "forage.zone", "lact.stat" , "ifbf.per" |
"MCH 2022-35", -23.44, 2.669 , "W" , "non-lactating", 6.713 |
"MCH 2022-37", -22.94, 2.833 , "W" , "non-lactating", 18.19 |
"MCH 2022-36", -23.26, 3.431 , "W" , "non-lactating", 17.18 |
"MCH 2022-34", -23.18, 1.529 , "W" , "lactating" , 6.789 |
"MCH 2022-27", -23.2 , 1.078 , "W" , "lactating" , 6.088 |
"MCH 2022-61", -23.31, 2.137 , "E" , "lactating" , 6.088 |
"MCH 2022-62", -22.68, 3.578 , "E" , "non-lactating", 8.757 |
"MCH 2020-02", -23.28, 3.501 , "W" , "lactating" , 6.789 |
"MCH 2012-15", -23.03, 2.179 , "W" , "non-lactating", 14.52 |
"MCH 2021-11", -23.47, 2.712 , "W" , "non-lactating", 11.38 |
"MCH 2021-01", -23.06, 3.179 , "E" , "non-lactating", 12.97 |
"MCH 2021-04", -23.47, 2.083 , "E" , "lactating" , 6.088 |
"MCH 2021-03", -23 , 3.949 , "E" , "lactating" , 7 |
"MCH 2021-16", -22.62, 1.724 , "W" , "lactating" , 7.491 |
"MCH 2022-32", -24.24, 3.053 , "W" , "lactating" , 6.088 |
"MCH 2022-30", -23.59, 2.073 , "W" , "non-lactating", 8.845 |
"MCH 2022-28", -23.52, 1.786 , "W" , "non-lactating", 5.783 |
"MCH 2022-26", -23.62, 1.847 , "W" , "non-lactating", 8.192 |
"MCH 2022-58", -23.27, 2.003 , "E" , "lactating" , 7.69 |
"MCH 2022-63", -23.11, 2.26 , "E" , "non-lactating", 14.77 |
"MCH 2022-41", -22.6 , 1.088 , "W" , "non-lactating", 6.37 |
"MCH 2022-39", -23.41, 0.8646, "W" , "lactating" , 6.088 |
"MCH 2022-46", -22.52, 4.815 , "E" , "non-lactating", 8.148 |
"MCH 2022-49", -22.02, 3.187 , "E" , "non-lactating", 7.81 |
"MCH 2022-57", -23.06, 3.591 , "W" , "non-lactating", 16.2 |
"MCH 2022-29", -23.29, 1.691 , "W" , "lactating" , 6.283 |
"MCH 2022-52", -23.35, 0.9461, "W" , "lactating" , 6.088 |
"MCH 2022-53", -22.71, 2.405 , "W" , "non-lactating", 12.14 |
"MCH 2022-54", -22.99, 3.31 , "W" , "non-lactating", 9.765 |
"MCH 2022-55", -23.65, 1.963 , "W" , "lactating" , 6.088 |
"MCH 2022-56", -23.31, 1.493 , "W" , "lactating" , 6.088 |
"MCH 2022-31", -22.88, 1.24 , "W" , "lactating" , 6.27 |
"MCH 2022-60", -22.74, 4.147 , "E" , "lactating" , 6.088 |
"MCH 2021-08", -23.01, 3.341 , "E" , "non-lactating", 11.08 |
"MCH 2021-10", -22.93, 3.102 , "E" , "non-lactating", 6.894 |
"MCH 2021-07", -23.07, 3.626 , "E" , "non-lactating", 13.89 |
"MCH 2021-09", -22.41, 3.142 , "E" , "non-lactating", 7.739 )
write.csv(d_mix ,'data/isotopes_covariates.csv', row.names = FALSE)
#bring in factors
d_disc <- wrapr::build_frame(
"source.ftagg" , "Meand13C", "SDd13C", "Meand15N", "SDd15N" |
"fung" , 2.61 , 0.61 , 4.23 , 1.05 |
"gram" , 2.61 , 0.61 , 4.23 , 1.05 |
"lich" , 2.61 , 0.61 , 4.23 , 1.05 |
"shrub.moss" , 2.61 , 0.61 , 4.23 , 1.05 )
write.csv(d_disc,'data/disc_factors.csv', row.names = FALSE)
#bring in source data
d_source <- wrapr::build_frame(
"source.ftagg" , "Meand13C", "SDd13C", "Meand15N", "SDd15N", "d13C.ftagg.N", "Concd13C", "Concd15N", "n" |
"fung" , -22.45 , 1.796 , 4.957 , 0.4065 , 4L , 0.4615 , 0.0305 , 10000L |
"gram" , -26.91 , 1.34 , 1.828 , 1.092 , 12L , 0.4492 , 0.01733 , 10000L |
"lich" , -23.97 , 2.025 , -5.002 , 1.756 , 26L , 0.4785 , 0.004692 , 10000L |
"shrub.moss" , -27.55 , 0.905 , -3.655 , 1.759 , 20L , 0.4803 , 0.0125 , 10000L )
write.csv(d_source,'data/sources.csv', row.names = FALSE)
#load mixture into mixsiar
mix = load_mix_data(filename='data/isotopes_covariates.csv'
iso_names=c("d13C","d15N"),
factors=c('forage.zone', 'lact.stat'),
fac_random = c(F, F),
fac_nested = c(F, F),
cont_effects = 'ifbf.per')
#load source data into mixsiar
source=load_source_data(filename='data/disc_factors.csv',
source_factors = NULL,
conc_dep = T,
data_type = "means",
mix)
#load disc data into mixsiar
discr = load_discr_data(filename="data/disc_factors.csv", mix)
#jags model
model_filename <- 'mixmods/mod1.txt'
resid_err <- F
process_err <- T#choosing this because the animal feeds like a filter feeders
write_JAGS_model(model_filename, resid_err, process_err, mix, source)
jags.1 <- run_model(run='normal',mix,source,discr,model_filename,alpha.prior =1,resid_err,process_err)
#view outputs
output_JAGS(jags.1, mix, source, output_options=list(summary_save = TRUE, summary_name ='mixmod_outputs/proportions',
sup_post = T, plot_post_save_pdf = T, plot_post_name ='mixmod_outputs/posteriors',
sup_pairs = T, plot_pairs_save_pdf = F, plot_pairs_name = NULL, sup_xy
= TRUE, plot_xy_save_pdf = F, plot_xy_name = "xy_plot", gelman = TRUE, heidel =
FALSE, geweke = TRUE, diag_save = TRUE, diag_name ='mixmod_outputs/diagnostics',,
indiv_effect =F, plot_post_save_png = F, plot_pairs_save_png = FALSE, plot_xy_save_png =
FALSE, diag_save_ggmcmc = F, return_obj=T))
#if you want to see posteriors for the run
g.post <- output_posteriors(jags.1, mix, source)
#to process posteriors to look at continuous effect
#code adapted from the PeerJ alligator example
#creates a dataframe "df"
{
# Process posteriors to make Figure 6
# Have to
R2jags::attach.jags(jags.1)
n.sources <- source$n.sources
source_names <- source$source_names
calc_eps <- function(f){
n.sources <- length(f)
gam <- rep(1/n.sources,n.sources)
phi <- rep(0,n.sources)
phi[1] <- 1
sqrt(sum((f-gam)^2))/sqrt(sum((phi-gam)^2))
}
ce=1
label <- mix$cont_effects[ce]
cont <- mix$CE[[ce]]
ilr.cont <- get(paste("ilr.cont",ce,sep=""))
n.plot = 200
chain.len = dim(p.global)[1]
Cont1.plot <- seq(from=round(min(cont),1), to=round(max(cont),1), length.out=n.plot)
ilr.plot <- array(NA,dim=c(n.plot, n.sources-1, chain.len))
for(src in 1:n.sources-1){
for(i in 1:n.plot){
ilr.plot[i,src,] <- ilr.global[,src] + ilr.cont[,src]*Cont1.plot[i]
}
}
# Transform every draw from ILR-space to p-space
e <- matrix(rep(0,n.sources*(n.sources-1)),nrow=n.sources,ncol=(n.sources-1))
for(i in 1:(n.sources-1)){
e[,i] <- exp(c(rep(sqrt(1/(i*(i+1))),i),-sqrt(i/(i+1)),rep(0,n.sources-i-1)))
e[,i] <- e[,i]/sum(e[,i])
}
# dummy variables for inverse ILR calculation
cross <- array(data=NA,dim=c(n.plot, chain.len, n.sources, n.sources-1))
tmp <- array(data=NA,dim=c(n.plot, chain.len, n.sources))
p.plot <- array(data=NA,dim=c(n.plot, chain.len, n.sources))
for(i in 1:n.plot){
for(d in 1:chain.len){
for(j in 1:(n.sources-1)){
cross[i,d,,j] <- (e[,j]^ilr.plot[i,j,d])/sum(e[,j]^ilr.plot[i,j,d]);
}
for(src in 1:n.sources){
tmp[i,d,src] <- prod(cross[i,d,src,]);
}
for(src in 1:n.sources){
p.plot[i,d,src] <- tmp[i,d,src]/sum(tmp[i,d,]);
}
}
}
# now take quantiles, after ILR transform of every draw
# get_high <- function(x){return(quantile(x,.975))} # 95% CI
# get_low <- function(x){return(quantile(x,.025))}
get_high <- function(x){return(quantile(x,.95))} # 90% CI
get_low <- function(x){return(quantile(x,.05))}
# get_high <- function(x){return(quantile(x,.75))} # 50% CI
# get_low <- function(x){return(quantile(x,.25))}
p.low <- apply(p.plot, c(1,3), get_low)
p.high <- apply(p.plot, c(1,3), get_high)
p.median <- apply(p.plot, c(1,3), median)
colnames(p.median) <- source_names
Cont1.plot <- Cont1.plot*mix$CE_scale + mix$CE_center # transform Cont1.plot (x-axis) back to the original scale
df <- data.frame(reshape2::melt(p.median)[,2:3], rep(Cont1.plot,n.sources), reshape2::melt(p.low)[,3], reshape2::melt(p.high)[,3])
colnames(df) <- c("source","median","ifbf.per","low","high")}
