1
$\begingroup$

There was a recent publication discussing the advantages of the zero-inflated dirichlet for microbiome count data which is compositional (you are modeling a matrix of species relative abundance data as a dependent variable, where the relative abundances within a sample sum to 1), and has lots of zeros (many taxa are not observed in over half of samples).

I am sold on this in principle. Authors provide an R package to fit these models. However, the group I work with is most interested in using Bayesian statistics for a suite of reasons. As a result, I need to translate this into the JAGS/BUGS language.

Fitting a standard dirichlet is easy enough and I copy data generation and fitting code below using JAGS and the runjags package in R. To do this I first have to add a 1 to all observations (the dirichlet can't handle hard zeros as is) before I normalize by rowSums to get proportions. Any advice on how to zero-inflate this model would be extremely welcome.

#clear environment, load packages.
rm(list=ls())
library(runjags)

#data generation.-----
#3 species abundances
set.seed(1234)
N <- 200
N.spp <- 3
mu <- rnorm(N.spp, 300, sd = 10)

#sequence depth per sample (essentially sampling 'effort').
depth <- 2000

#get a zero inflation probability.
z.p <- runif(N.spp,0.2,0.3)

#Generate matrix of abundances. Round because these need to be counts.
spp.list <- list()
for(i in 1:N){
  spp.list[[i]] <- round(rnorm(N.spp,mu,sd=10) * rbinom(N.spp,1,z.p))
}
spp.dat <- data.frame(do.call(rbind,spp.list))
other <- depth - rowSums(spp.dat)
spp.dat <- data.frame(other, spp.dat)
colnames(spp.dat) <- c('other','spp.1','spp.2','spp.3')

#normalize.
spp.dat <- spp.dat + 1
spp.dat <- spp.dat / rowSums(spp.dat)

#specify JAGS dirichlet model without zero inflation.----
jags.model = "
model {
  alpha ~ dnorm(0, 1.0E-3) 
  mu[1] <- 0
  for (j in 2:N.spp) {
    mu[j] ~ dnorm(0, 1.0E-3)
  }

  #implement dirlichet
  for (i in 1:N){
    for(j in 1:N.spp){
      log(a0[i,j]) <- alpha + mu[j]
    }
    y[i,1:N.spp] ~ ddirch(a0[i,1:N.spp]) 
  }
}
"

#setup jags data object, run model.----
#jags data list.
jd <- list(y=as.matrix(spp.dat), N.spp = ncol(spp.dat), N=N)

#fit jags model with run.jags function from runjags.
jags.fit <- run.jags(model = jags.model,
                     data = jd,
                     n.chains = 3,
                     adapt = 200,
                     burnin = 500,
                     sample = 500,
                     monitor = c('alpha','mu'))
summary(jags.fit)

Once this is run it is pretty easy to show the model over-estimates the abundances of the zero inflated species.

out <- summary(jags.fit)
alpha <- out[1,4]
mu <- out[2:nrow(out),4]
mu <- mu + alpha
fitted <- boot::inv.logit(mu)
observed <- colMeans(spp.dat)
data.frame(fitted,observed)
         fitted   observed
mu[1] 0.8570854 0.89392964
mu[2] 0.1876097 0.03061128
mu[3] 0.1830003 0.02834331
mu[4] 0.2061653 0.04711577
$\endgroup$
  • $\begingroup$ There seems to be a statistical question embedded here, but otherwise requests for coding advice are likely to be off-topic here. Please emphasise what you are asking that doesn't hinge on how to code it. $\endgroup$ – Nick Cox Oct 2 '18 at 15:04
  • $\begingroup$ @NickCox i have removed the word "code". This is an applied statistics questions however, and implementation is central to answering it. It requires understanding both the statistical problem and how Bayesian models are fit using MCMC sampling software, making this more than just a coding question. Not sure how to get around this, many applied statistics questions are implementation questions. $\endgroup$ – colin Oct 2 '18 at 15:19
  • $\begingroup$ Indeed, and many such are asked and answered on Stack Overflow. As statistical expertise is required to answer this, this one seems to fall this side of the fence. Thanks for the edit. $\endgroup$ – Nick Cox Oct 2 '18 at 15:40

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.