# Excellent fit, zero convergence hierarchical dirichlet model in JAGS

I am fitting a hierarchical dirichlet model to some data in JAGS. My samples (referred to as cores in the code) are observations of the relative abundance of 3 species. The rows of the matrix always sum to 1. Samples (cores) are nested within plots, and plots are nested within sites. I am trying to hierarchically estimate plot and site means based on core observations.

I have built a dirichlet model to hierarchically fit these means. The model estimates means very well, however chains do not converge. This is a problem as the whole reason I am fitting this model in JAGS is to get a robust estimate of uncertainty. Any ideas on how to improve model convergence here would be greatly appreciated. I have tried running this with 1 million burnin iterations, which just confirmed it will never converge. Data generation depends on the DirichletReg R package and fitting the model depends on the runjags R package.

Generate pseudo data:

#1. generate pseudo data.----
#3 species relative abundances across 4 sites.
spp1 <- c(0.2,0.8,0.5,0.1)
spp2 <- c(0.3,0.1,0.4,0.1)
spp3 <- 1 - (spp1 + spp2)
y <- data.frame(spp1,spp2,spp3)
y <- as.matrix(y)
site_sd <- 0.02
truth <- y

#get some numbevr of pltos and cores within plots.
n.site <- nrow(y)
n.plot <- 15
n.core <- 10

## Variability partitions:
site_alpha <- n.plot
plot_alpha <- n.core
core_alpha <- 1000

## Modified simulation to reflect variability partitions:
site.list <- vector('list', length=n.site)
for(i in 1:n.site){
plot_mu <- DirichletReg::rdirichlet(1, y[i,]*plot_alpha)[1,]
plot.list <- vector('list', length=n.plot)
for(j in 1:n.plot){
cores <- DirichletReg::rdirichlet(n.core, plot_mu*core_alpha)
# Re-written to avoid conversion of numeric to character/factor:
cores <- cbind(data.frame(rep(LETTERS[i],n.core),rep(LETTERS[j],n.core),LETTERS[1:n.core], stringsAsFactors=TRUE), cores)
colnames(cores) <- c('siteID','plotID','coreID','spp1','spp2','spp3')
plot.list[[j]] <- cores
}
site.list[[i]] <- do.call(rbind,plot.list)
}
dat <- do.call(rbind,site.list)
y <- dat[,4:6]

#core_plot, plot_site factors.
dat$$plotID <- factor(paste0(dat$$siteID,'_',dat$$plotID)) core_plot = dat$$plotID
plot_site <- unique(dat[,c('siteID','plotID')])$siteID  Specify JAGS model as an object in R: #jags model.---- jags.model1 = " model { #plot means. for(i in 1:N.core){ for(j in 1:N.spp){ log(core.hat[i,j]) <- plot_mu[core_plot[i],j] + plot.intercept } y[i,1:N.spp] ~ ddirch(core.hat[i,1:N.spp]) } #link plot and site scales. for(i in 1:N.plot){ plot_mu[i,1:N.spp] <- log(plot_l1[i,1:N.spp]) } #site means. for(i in 1:N.plot){ for(j in 1:N.spp){ log(plot.hat[i,j]) <- site_mu[plot_site[i],j] + site.intercept } plot_l1[i,1:N.spp] ~ ddirch(plot.hat[i,1:N.spp]) } #priors. plot.intercept ~ dnorm(0,1E-3) site.intercept ~ dnorm(0,1E-3) for(i in 1:N.site){ site_mu[i,1] <- 0 for(j in 2:N.spp){ site_mu[i,j] ~ dnorm(0,1E-3) } } }" #end jags model.  Specify JAGS data object, fit model using the runjags package. #JAGS data for hierarchical plot-means model.---- jd1 <- list(y=as.matrix(y), N.plot = length(plot_site), N.site = n.site, N.core = nrow(y), N.spp = ncol(y), core_plot=as.factor(core_plot), plot_site=as.factor(plot_site)) #run the hierarchical jags model.---- jmod <- run.jags(model = jags.model1, data = jd1, n.chains = 3, monitor = c('plot.intercept','site.intercept','plot_mu','site_mu'), adapt = 1000, burnin = 4000, sample = 3000)  Summarize model output and plot fit: #Summarize model output.---- #check the site parameters. plot.int <- summary(jmod, vars = 'plot.intercept') site.int <- summary(jmod, vars = 'site.intercept') plot.mu <- summary(jmod, vars = 'plot_mu') site.mu <- summary(jmod, vars = 'site_mu') #plot(mod, vars = 'site_mu') #Get fitted plot and site means. plot.sum <- matrix(plot.mu[,4], nrow = n.site*n.plot, ncol = ncol(y)) + plot.int[,4] plot.sum <- exp(plot.sum) / rowSums(exp(plot.sum)) site.sum <- matrix(site.mu[,4], nrow = n.site , ncol = ncol(y)) + site.int[,4] site.sum <- exp(site.sum) / rowSums(exp(site.sum)) #Get true plot and site means. plot.truth <- aggregate( . ~ plotID, data = dat, FUN = mean) site.truth <- aggregate( . ~ siteID, data = dat, FUN = mean) plot.truth <- as.matrix(plot.truth[,4:6]) site.truth <- as.matrix(site.truth[,4:6]) #plot plot and site fits.---- par(mfrow=c(1,2)) #plot plot-level fit. plot(as.vector(plot.truth) ~ as.vector(plot.sum));mod<-lm(as.vector(plot.truth) ~ as.vector(plot.sum));abline(mod, lwd =2) mtext(paste0('R2=',round(summary(mod)$r.squared,2)), side = 3, line = -1.5, adj =0.05)
abline(0,1, col='purple',lty = 2)

#plot site level fit.
plot(as.vector(site.truth) ~ as.vector(site.sum));mod<-lm(as.vector(site.truth) ~ as.vector(site.sum));abline(mod, lwd =2)
mtext(paste0('R2=',round(summary(mod)\$r.squared,2)), side = 3, line = -1.5, adj =0.05)
abline(0,1, col='purple',lty = 2)


Fit: