Dirichlet / Gamma model...
model {
for (i in 1:npupil) {
Goals[i] ~ dcat(p[School[i],])
}
for (j in 1:nschool) {
for (k in 1:3) {
p[j,k] <- q[j,k]/sum(q[j,])
q[j,k] ~ dgamma(a[k], 1)
}
}
for (k in 1:3) {
a[k] ~ dgamma(1, 0.001)
p.pop[k] <- a[k]/sum(a[]) # population mean of p[,k]
}
dummy <- Gender[1] # stop WinBUGS complaining about unused variable
}


Data:
Click arrow for data - same data used for all models in this example

Inits:
list(a=c(1,1,1))

Poor convergence. Click arrow for trace plots of a[], which are highly cross-correlated.

No difference between schools...
model {
for (i in 1:npupil) {
Goals[i] ~ dcat(p[i,])
for (k in 1:3) {
p[i,k] <- q[i,k]/sum(q[i,])
log(q[i,k]) <- a[i,k]
}
a[i,1] <- b[1] + b.boy*Gender[i]
a[i,2] <- b[2]
a[i,3] <- 0
}
b[1] ~ dnorm(0, 0.0001)
b[2] ~ dnorm(0, 0.0001)
b.boy ~ dnorm(0, 0.0001)
or.boy <- exp(b.boy)

qboy[1] <- exp(b[1] + b.boy); qboy[2] <- exp(b[2]); qboy[3] <- 1
qgirl[1] <- exp(b[1]); qgirl[2] <- exp(b[2]); qgirl[3] <- 1

# Probabilities of preferring 1) sports 2) popularity 3) grades for boys and girls separately
for (k in 1:3) {
   p.boy[k] <- qboy[k]/sum(qboy[])
   p.girl[k] <- qgirl[k]/sum(qgirl[])
}
dummy <- School[1] + nschool
}

Inits:
list(b.boy=0, b=c(0,0))


   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   b.boy   0.9855   0.2476   0.005047   0.5125   0.9823   1.473   4003   50000
   or.boy   2.763   0.6999   0.01422   1.669   2.671   4.364   4003   50000
   
Dbar = post.mean of -2logL; Dhat = -2LogL at post.mean of stochastic nodes
   Dbar   Dhat   pD   DIC   
Goals   957.489   954.473   3.016   960.506   
total   957.489   954.473   3.016   960.506   
   
Independent school effects...
model {
for (i in 1:npupil) {
Goals[i] ~ dcat(p[i,])
for (k in 1:3) {
p[i,k] <- q[i,k]/sum(q[i,])
log(q[i,k]) <- a[i,k]
}
a[i,1] <- b[School[i], 1] + b.boy*Gender[i]
a[i,2] <- b[School[i], 2]
a[i,3] <- 0
}
for (j in 1:nschool) {
b[j,1] ~ dnorm(0, 0.0001)
b[j,2] ~ dnorm(0, 0.0001)
b[j,3] <- 0
for (k in 1:3) {    
egirl[j,k] <- exp(b[j,k])
}
eboy[j,1] <- exp(b[j,1] + b.boy); for (k in 2:3) {eboy[j,k] <- exp(b[j,k])}
for (k in 1:3) {
p.girl[j,k] <- egirl[j,k] / sum(egirl[j,])
p.boy[j,k] <- eboy[j,k] / sum(eboy[j,])
}
}
b.boy ~ dnorm(0, 0.0001)
or.boy <- exp(b.boy)
}

Inits:
list(b.boy=0, b=structure(.Data=c(0, 0, NA, 0, 0, NA, 0, 0, NA, 0, 0, NA, 0, 0, NA, 0, 0, NA, 0, 0, NA, 0, 0, NA, 0, 0, NA),.Dim=c(9,3)))

   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   b.boy   1.088   0.2639   0.005174   0.5781   1.084   1.61   4002   50000
   or.boy   3.074   0.8356   0.01662   1.783   2.956   5.005   4002   50000

   Dbar   Dhat   pD   DIC   
Goals   936.567   917.409   19.158   955.725   
total   936.567   917.409   19.158   955.725   

Hierarchical school effects...
model {
for (i in 1:npupil) {
Goals[i] ~ dcat(p[i,])
for (k in 1:3) {
p[i,k] <- q[i,k]/sum(q[i,])
log(q[i,k]) <- a[i,k]
}
a[i,1] <- b[School[i], 1] + b.boy*Gender[i]
a[i,2] <- b[School[i], 2]
a[i,3] <- 0
}
for (j in 1:nschool) {
b[j,1] ~ dnorm(mu[1], tau[1])
mub2[j] <- mu[2] + s[2]/s[1]*cor*(b[j,1] - mu[1])
b[j,2] ~ dnorm(mub2[j], taub2)
b[j,3] <- 0
for (k in 1:3) {    
egirl[j,k] <- exp(b[j,k])
}
eboy[j,1] <- exp(b[j,1] + b.boy); for (k in 2:3) { eboy[j,k] <- exp(b[j,k]) }
for (k in 1:3) {
p.girl[j,k] <- egirl[j,k] / sum(egirl[j,])
p.boy[j,k] <- eboy[j,k] / sum(eboy[j,])
}
}
vb2 <- (1 - cor*cor)*v[2]
tau[1] <- 1/v[1]
taub2 <- 1/vb2
for (k in 1:2) {
mu[k] ~ dnorm(0, 0.0001)
v[k] <- s[k]*s[k]
s[k] ~ dunif(0, 100)
}
cor ~ dunif(0, 1)
b.boy ~ dnorm(0, 0.0001)
or.boy <- exp(b.boy)
}

Inits:
list(mu=c(0,0), s=c(1,1), cor=0.5)

   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   b.boy   1.053   0.2549   0.003429   0.5596   1.052   1.56   4002   200000
   or.boy   2.961   0.7715   0.01031   1.75   2.862   4.758   4002   200000

   Dbar   Dhat   pD   DIC   
Goals   937.938   924.479   13.459   951.397   
total   937.938   924.479   13.459   951.397

Monitor p.girl and p.boy to obtain the data for the graph in Figure 10.5, which was produced in R.