How to extract posterior data from output_posteriors object? about mixsiar HOT 1 OPEN

I am having this same issue. In the output options for my model with one fixed categorical variable (Year) and one continuous variable (Age), I can get plots of proportion by age for each year, but I would like to be able to modify these in ggplot. However, when I extract the code using a similar set of functions as the alligator code, I get the age and proportion plot for all years combined. How would I extract posterior data for each level of my categorical variable?

Here is my code to replicate extracting the posteriors for the continuous variable with all levels of my categorical variable combined:

R2jags::attach.jags(jags.2)
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]
}
}

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

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,]);
}
}
}

get_high <- function(x){return(quantile(x,.95))} # 90% CI
get_low <- function(x){return(quantile(x,.05))}
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")

from mixsiar.