# Trouble specifying a hierarchical dirichlet model in JAGS

I have a sampling design where samples (cores) are taken within plots. Those plots are then nested within sites. There are multiple sites. I would like to get a hierarchical site-level estimate of observation y, observed at the core level. Furthermore, my y observation is a multivariate relative abundance. It is the relative abundance of 3 species, and the total of those 3 relative abundances always sums to 1. Here is R code to simulate data, dat like this:

#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

#simulate from dirichlet distribution.
site.list <- list()
for(i in 1:n.site){
plot.list <- list()
for(j in 1:n.plot){
cores <- DirichletReg::rdirichlet(n.core,y[i,])
cores <- data.frame(cbind(rep(LETTERS[i],n.core),rep(LETTERS[j],n.core),LETTERS[1:n.core],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)
dat[,4:6] <- sapply(dat[4:6],as.character)
dat[,4:6] <- sapply(dat[4:6],as.numeric)
y <- dat[,4:6]

#core_plot, plot_site factors.
dat$plotID <- paste0(dat$siteID,'_',dat$plotID) core_site = dat$siteID
core_plot = dat$plotID plot_site = dat[seq(1, nrow(dat), 3), ]$siteID

Ignoring plots, I can get site level means using the following JAGS model:

#JAGS model for site only case.----
jags.model1 = "
model {
#get site level means.
for(i in 1:N.core){
for(j in 1:N.spp){
log(core.hat[i,j]) <- site_mu[core_site[i],j]
}
y[i,1:N.spp] ~ ddirch(core.hat[i,1:N.spp])
}
#prior
for(i in 1:N.site){
for(j in 1:N.spp){
site_mu[i,j] ~ dnorm(0,1E-3)
}
}
}" #end jags model.

#JAGS data for site only case.----
jd1 <- list(y=as.matrix(y), N.site = n.site, N.core = nrow(y), N.spp = ncol(y),
core_site = as.factor(core_site))

#4. Run JAGS model.----
test <- run.jags(model = jags.model1,
data = jd1,
n.chains = 3,
monitor = c('site_mu'),
burnin = 1000,
sample = 1000)

This model fits well, as evidenced by the following plot of predicted vs. observed:

out <- summary(test)
spp.sum <- matrix(out[,4], nrow = n.site, ncol = ncol(y))
spp.sum <- boot::inv.logit(spp.sum) / rowSums(boot::inv.logit(spp.sum))
plot(as.vector(truth) ~ as.vector(spp.sum));mod<-lm(as.vector(truth) ~ as.vector(spp.sum));abline(mod)
mtext(paste0('R2=',round(summary(mod)$r.squared,3)), side = 3, adj = 0.05, line = -2) This begins to fail when I add an intermediate plot-level hierarchy. Here is the JAGS model with a plot-level hierarchy as well as code to specify a JAGS data object and fit the model: # JAGS model for hierarchical plot-site case.---- jags.model2 = " model { #get plot level means. for(i in 1:N.core){ for(j in 1:N.spp){ log(core.hat[i,j]) <- plot_mu[core_plot[i],j] } y[i,1:N.spp] ~ ddirch(core.hat[i,1:N.spp]) } #get site level means. for(i in 1:N.plot){ for(j in 1:N.spp){ log(plot.hat[i,j]) <- site_mu[plot_site[i],j] } plot_mu[i,1:N.spp] ~ ddirch(plot.hat[i,1:N.spp]) } #prior for(i in 1:N.site){ for(j in 1:N.spp){ site_mu[i,j] ~ dnorm(0,1E-3) } } }" #end jags model. #JAGS data for hierarchical plot-site case.---- jd2 <- list(y=as.matrix(y), N.site = n.site, N.plot = length(plot_site), N.core = nrow(y), N.spp = ncol(y), core_plot = as.factor(core_plot), plot_site=as.factor(plot_site)) #Run hierarchical plot-site model.---- test <- run.jags(model = jags.model2, data = jd2, n.chains = 3, monitor = c('site_mu'), adapt = 200, burnin = 1000, sample = 1000) This model does not converge, and does not capture the site-level means. The correlation between predicted vs. observed is actually negative: out <- summary(test) spp.sum <- matrix(out[,4], nrow = n.site, ncol = ncol(y)) spp.sum <- boot::inv.logit(spp.sum) / rowSums(boot::inv.logit(spp.sum)) spp.sum;truth plot(as.vector(truth) ~ as.vector(spp.sum));mod<-lm(as.vector(truth) ~ as.vector(spp.sum));abline(mod) mtext(paste0('R2=',round(summary(mod)$r.squared,3)), side = 3, adj = 0.05, line = -2)

My question is: Where am I going wrong in the hierarchical model? How can I capture thee true site-level means while accounting for the hierarchical sampling design?

It seems that the way in which you are thinking about the Dirichlet distribution is missing something crucial, which has affected both your data simulation and the model itself. Remember that the response of a Dirichlet distribution is a vector of probabilities (summing to 1), but the parameter vector is NOT a vector of probabilities but strictly positive numbers that may go over 1. In fact, the variance is inversely related to the sum of the parameter vector - if the parameters are all less than 1 then there is huge variance in the response. Think of a simplified example using the Beta distribution (of which the Dirichlet is a generalisation) - the following three parameterisations have equal mean but very different variance:

curve(dbeta(x,500,500), col='green')

In effect, your data simulation is forcing the response variable to do something like the red line - the variance is locked at a (very high) quantity because you are using a parameterisation where the vector must sum to 1. You need to add parameters to the data generation process to control the variance associated with each hierarchical partition:

#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:
plot_alpha <- 10
core_alpha <- 20

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

# Examine between-core quantiles of y by site:
for(i in 1:n.site){
print(apply(y[dat$siteID==levels(dat$siteID)[i],],2,quantile))
}

#core_plot, plot_site factors.
dat$plotID <- factor(paste0(dat$siteID,'_',dat$plotID)) core_site = dat$siteID
core_plot = dat$plotID # There seems to be an error in this code: # plot_site = dat[seq(1, nrow(dat), 3), ]$siteID
# Corrected:
plot_site <- unique(dat[,c('siteID','plotID')])$siteID #JAGS data for hierarchical plot-site case.---- jd2 <- list(y=as.matrix(y), N.site = n.site, N.plot = length(plot_site), N.core = nrow(y), N.spp = ncol(y), core_plot = as.factor(core_plot), plot_site=as.factor(plot_site)) This allows the similarity of cores-within-plots to be adjusted (core_alpha) relative to the similarity of plots-within-sites (plot_alpha), which you can now see reflected in the simulated data. The model can then be re-written to reflect the data generation: jags.model2 = " model { # Observations (single set per core): for(i in 1:N.core){ y[i,1:N.spp] ~ ddirch(plot_mu[core_plot[i],1:N.spp] * core_alpha) } # Plot means: for(i in 1:N.plot){ plot_mu[i,1:N.spp] ~ ddirch(site_mu[plot_site[i],1:N.spp] * plot_alpha) } # Site means: for(i in 1:N.site){ site_mu[i,1:N.spp] ~ ddirch(priors[1:N.spp]) } # Priors: for(s in 1:N.spp){ priors[s] <- 1 } plot_alpha <- 10 #~ dgamma(0.01, 0.01) core_alpha <- 20 #~ dgamma(0.01, 0.01) }" #end jags model. library('runjags') test <- run.jags(model = jags.model2, data = jd2, n.chains = 2, monitor = c('site_mu', 'core_alpha', 'plot_alpha'), adapt = 1000, burnin = 1000, sample = 1000) test out <- summary(test, var='site_mu') spp.sum <- matrix(out[,4], nrow = n.site, ncol = ncol(y)) spp.sum;truth plot(as.vector(truth) ~ as.vector(spp.sum));mod<-lm(as.vector(truth) ~ as.vector(spp.sum));abline(mod) mtext(paste0('R2=',round(summary(mod)$r.squared,3)), side = 3, adj = 0.05, line = -2)

Now the estimated site_mu values seem to follow the simulation parameters reasonably well, although you would obviously need to check for convergence and effective sample size (1000 iterations is not enough).

Notice however that I have cheated, and fixed the alpha values in the model to the values that are known from the simulation parameters. This is because the model is very sensitive to the priors for these parameters (some examples commented out), and tends to shrink the variance associated with sites-within-plots to the point where the model estimates effectively no variation at this level and consequently gives unrealistically small confidence intervals for site_mu. In reality you would not be able to fix these parameters at known values, so this model may be difficult to implement unless you have a much larger dataset in real life and/or you have good priors for these variance partitions. At the very least you would have to compare your results with simpler models, and investigate the complexity of the model that it is possible to fit to your data. In the end ignoring the site-within-plot hierarchy may be the best option.

For the record it is also possible to implement this hierarchical model as a GLM - but I tried this and the model actually performs worse:

jags.model3 <- "
model {

for(i in 1:N.core){
# Observations (single set per core):
y[i,1:N.spp] ~ ddirch(plot_mu[i,1:N.spp] * core_alpha)

# A GLM-type formulation (random effect of plot and fixed effect of site):
for(s in 1:N.spp){
log(plot_mu[i,s]) <- intercept[s] + site_effect[siteID[i],s] + plot_effect[plotID[i],s]
}
}

# Random effect of plot (tau could differ between species if desired):
for(s in 1:N.spp){
for(p in 1:N.plot){
plot_effect[p,s] ~ dnorm(0, plot_tau)
}
}

# Prior for site effects (zero for the first site):
for(s in 1:N.spp){
site_effect[1,s] <- 0
for(t in 2:N.site){
site_effect[t,s] ~ dnorm(0, 10^-6)
}

# Needed only to monitor site_mu:
for(t in 1:N.site){
log(theta[t,s]) <- intercept[s] + site_effect[t,s]
site_mu[t,s] <- theta[t,s] / sum(theta[t,1:N.spp])
}
}

# Prior for intercepts (zero for the first species):
intercept[1] <- 0
for(s in 2:N.spp){
intercept[s] ~ dnorm(0, 10^-6)
}

# Prior for core alpha:
core_alpha ~ dgamma(0.01, 0.01)
# Prior for plot tau:
plot_tau ~ dgamma(0.01, 0.01)

}" #end jags model.

jd3 <- list(y=as.matrix(y), N.core=nrow(dat), N.site = length(levels(dat$siteID)), N.plot = length(levels(dat$plotID)), N.spp = ncol(y),
siteID = dat$siteID, plotID = dat$plotID)

test <- run.jags(model = jags.model3,
data = jd3,
n.chains = 2,
monitor = c('site_mu', 'plot_tau', 'core_alpha'),
modules = 'glm',