Supplement

A BUGS code for correlation model

model{

   for(i in 1:N){
      dat[i,1:4]~dmnorm(mu[species[i],], R[,]) 
   }

   for(j in 1:Nspecies){	
      mu[j,1]~dnorm(gamaw, tauaw) # mu(aw) 
      mu[j,2]~dnorm(gambw, taubw) # mu(bw) 
      mu[j,3]~dnorm(gamaf, tauaf) # mu(af)  
      mu[j,4]~dnorm(gambf, taubf) # mu(bf) 
   }
   # Prior distributions	
   gamaw~dnorm(-4.59,2.11)
   tauaw~dgamma(6.25,5.3708)
   gambw~dnorm(1.11,81.47)
   tau_bw~dgamma(3,1.2)
   gamaf~dnorm(-5,0.25) 
   tauaf~dgamma(3,12)
   gambf~dnorm(1,2.5) 
   taubf~dgamma(3,1.2)

   R[1:4,1:4]~dwish(omega[1:4,1:4],4)	
   covar[1:4,1:4]<-inverse(R[,])
   for(i in 1:4){
      for(j in 1:4){
         corr[i,j]<-covar[i,j]/(sqrt(covar[i,i])*sqrt(covar[j,j]))
      }
   }
}

list(Nspecies=3, N=70,
omega=structure(.Data=c(0.001,0,0,0,0,0.001,0,0,0,0,0.001,0,0,0,0,0.001), 
.Dim=c(4,4)))


B BUGS code for model without correlation structure


model{
	
   for(i in 1:N){
      dat[i,1]~dnorm(mu[species[i],1], r[1]) #aw 
      dat[i,2]~dnorm(mu[species[i],2], r[2]) #bw
      dat[i,3]~dnorm(mu[species[i],3], r[3]) #af
      dat[i,4]~dnorm(mu[species[i],4], r[4]) #bf
   }

   for(j in 1:Nspecies){	
      mu[j,1]~dnorm(gamaw, tauaw) #mu(aw)
      mu[j,2]~dnorm(gambw, taubw) #mu(bw)
      mu[j,3]~dnorm(gamaf, tauaf) #mu(af)
      mu[j,4]~dnorm(gambf, taubf) #mu(bf)
   }

   # Prior distributions
   gamaw~dnorm(-4.59,2.11)
   tauaw~dgamma(6.25,5.3708)
   gambw~dnorm(1.1079,81.47)
   tau_bw~dgamma(3,1.2)
   gamaf~dnorm(-5,0.25) 
   tauaf~dgamma(3,12)
   gambf~dnorm(1,2.5) 
   taubf~dgamma(3,1.2)

   n<-4
   v<-0.001
   a<-n/2
   b<-v/2

   for(j in 1:4){
      r[j]~dgamma(a,b)
      var[j]<-1/r[j]
   }
}

list(Nspecies=3, N=70)