Correlated binomials in JAGS

[Arnaud]

Dear JAGS users,

I am searching a way to produce correlated random numbers sampled from several binomial distributions in JAGS.
In R, I use normal copula to correlate samples from several binomials, however, I do not know how to do that in JAGS.

Any help would be appreciated,

[Martyn Plummer]

Something like this:

Model file "copula.bug"

model {

   ## Set up prior for bivariate normal.

   # Mean vector
   zero[1] <- 0
   zero[2] <- 0

   # Variance-covariance matrix
   Sigma.X[1,1] <- 1
   Sigma.X[1,2] <- rho
   Sigma.X[2,1] <- rho
   Sigma.X[2,2] <- 1

   # Precision matrix
   Tau.X <- inverse(Sigma.X)

   ## Correlated bivariate normal
   X[1:2] ~ dmnorm(zero, Tau.X)

   ## Correlated uniforms
   p[1] <- pnorm(X[1], 0, 1)
   p[2] <- pnorm(X[2], 0, 1)

   ## Set up breaks for transformation to binomial
   for (i in 1:N) {
      breaks[i] <- pbin(i-1, prob, N)
   }

   ## Transform to binomial
   Y[1] ~ dinterval(p[1], breaks)
   Y[2] ~ dinterval(p[2], breaks)
}

R script:

library(rjags)
m <- jags.model("copula.bug", data=list(rho = 0.5, N=5, prob=0.3))
s <- coda.samples(m, "Y", n.iter=10000)[[1]]

# Marginal distributions are correct
table(s[,1])
table(s[,2])
table(rbinom(10000, prob=0.3, size=5))

# Binomials are correlated
cor(s[,1], s[,2])

[Arnaud]

Thanks Martyn, I will have to take time to study your answer !!
Cheers

[Arnaud]

Thanks a lot, it's exactly what I need.
Much appreciated !

[Arnaud]

Hi Martyn,

Long time since I tried to use your suggestion.
Here is a slightly modified version that fit my needs

model(

for (year in 1:5){
   # Mean vector
   for (i in 1:length(prob[,year])){
     zero[year,i] <- 0
   }

   # Variance-covariance matrix
   for (i in 1:length(prob[,year])){
     for(j in 1:length(prob[,year])){
       Sigma.X[year,i,j] <- ifelse(i == j, 1, rho)
     }
   }

   ## Correlated bivariate normal
   X[year, 1:length(prob[,year])] ~ dmnorm(zero[year,], inverse(Sigma.X[year,,]))  #inverse(Sigma.X[year,,]) is the precision matrix

   ## Correlated uniforms
   for (i in 1:length(prob[,year])){
     p[year,i] <- pnorm(X[year,i], 0, 1)
   }

   ## Set up breaks for transformation to binomial
   for (j in 1:length(prob[,year])){
     for (i in 1:max(size[,year])) {
        breaks[year,j,i] <- ifelse(i > size[j,year], 999, pbin(i-1, prob[j,year], size[j,year]))
     }

     ## Transform to binomial
     n[year,j] ~ dinterval(p[year,j], breaks[year,j,1:size[j,year]])
     probSample[year,j] <- n[year,j] / size[j,year]
   }
  }
#  ## fitting
# for (year in 1:5){
#  obsProb[year] ~ dnorm(probSample[year,3], 1000)
#  }
 }
)

The R script:

size <- matrix(rep(50,15), nrow=3)
prob <- matrix(rep(c(0.2,0.5,0.8), each=5), nrow=3, byrow=T)
obsProb <- rep(0.8,5)

library(rjags)
m <- jags.model("C:/temp/modelCOPULA.bug", data=list(rho = 0.7, size=size, prob=prob, obsProb=obsProb))
s <- coda.samples(m, c("probSample"), n.iter=10000)[[1]]

cor(s[,substring(colnames(s),1,13) == "probSample[2,"])  #test correlation for year 2

par(mfrow=c(2,2))
hist(s[,substring(colnames(s),1,14) == "probSample[2,1"], breaks=20, xlim=c(0,1)); hist(s[,substring(colnames(s),1,14) == "probSample[2,2"], breaks=20, xlim=c(0,1)); hist(s[,substring(colnames(s),1,14) == "probSample[2,3"], breaks=20, xlim=c(0,1))

For each of my 5 years, the model gives me three proportions sampled from correlated binomials (each of them having a defined mean). The correlation between the binomial is given.

In the current version, the "fitting" portion of the code is not active. In this case, no problem, the correlation among the proportion sampled for the binomial correspond to the one given as rho (always just a little lower).

Now, here is my problem. In activating the "fitting" portion of the code, I include observed proportions for each year and I make a link with the sampled one (including a normal error). In that case, the correlation no longer exists (just run the code removing the comments in front of the "fitting" part and look at the correlation at year 2, for example).

How can I do to keep the correlation in that case ?
In fact, I thought that the fitting would have updated the information in the three binomials, but keeping the correlation.

Any idea ??

Arnaud

[Arnaud]

Dears JAGS users,

Sorry to insist on this subject, but I hope finding a way to solve my problem.
Thanks to Martyn Plummer, I have a solution to produce correlated samples from binomials distributions (see above).
However, when fitting samples from one of the distributions to some "observed values", the correlation between binomials samples disappear.
In fact, it seems that the distribution of the binomial for which I provide "observed" values is updated (perfect), but samples from the other binomial distributions are no longer correlated.

Is their a way to maintain this correlation even when fitting to "observed" data?

Thanks

Arnaud