# Measuring covariance among multiple y-variables using a dirlichet model, implemented in JAGS in R

I have data on 3 species abundances as fractional compositions. Within each sample, the sum of all species abundances is 1. Essentially, each observation is a proportional abundance. Here is an example of this type of data using R:

y1 <- rnorm(100, 100, 5)
y2 <- y1*2 + rnorm(100,0,10)
y3 <- rnorm(100,50,5)
spp.y <- data.frame(y1,y2,y3)
tot <- rowSums(spp.y)
spp.y <- spp.y / tot
spp.y <- as.matrix(spp.y)


In this example, the abundances of y1 and y2 are correlated. y1 and y2 are uncorrelated with y3. I would like to specify a model that can capture and quantify the magnitude and sign of this correlation, while accounting for the fact that these observations are constrained by summing to 1. To do this I specified a dirlichet model in JAGS, the multivariate extension of the beta model:

Setup JAGS data as an R object, jags.data:

jags.data <- list(y = spp.y,
N = nrow(spp.y),
N.spp = ncol(spp.y))


Specify dirlichet model as an R object, jags.model:

jags.model = "
model {
#setup priors for each species
for(j in 1:N.spp){
m0[j] ~ dgamma(1.0E-3, 1.0E-3) #intercept prior
}

#implement dirlichet
for(i in 1:N){
for(j in 1:N.spp){
log(a0[i,j]) <- m0[j] ## eventually add linear model here
}
y[i,1:N.spp] ~ ddirch(a0[i,1:N.spp])
}

} #close model loop.
"


Fit the model in R using the runjags package:

library(runjags)
jags.out <- run.jags(jags.model,
data=jags.data,
burnin = 100,
sample = 200,
n.chains=3,
monitor=c('m0'))


This works well. The intercepts capture the within species relative abundances.

colMeans(spp.y)
a.m <- exp(summary(jags.out)[,"Median"])
a.m/sum(a.m)


However, this model assumes the species abundances are independent, which they are not. Furthermore, I want to capture and quantify the correlation between y1 and y2. I believe this is possible by adding covariance among species abundances, but I do not know how to implement this in JAGS. Any advice (and especially code!) would be super appreciated.