0
$\begingroup$

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: enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.