var
 X[I,J,K],      # observations
 n[I,J],        # total for each covariate pattern
 E[I,J,K],      # fitted values
 OlogOE[I,J,K], # O log O/E
 G2,            # goodness-of-fit statistic
 mu[I,J,K],     # Poisson means
 phi[I,J,K],    # exp (beta[k] ' x[i,j])
 p[I,J,K],      # fitted probabilities
 lambda[I,J],   # baseline rates in each covariate strata
 alpha[K],      # factor for food = 2,3,4,5
 beta[I,K],     # factor for lakes = 2,3,4, for each food
 b[I,K],        # factor for lakes = 2,3,4, relative to food 1, centred
 gamma[J,K],    # factor for size = 2,  for each food
 g[J,K];        # factor for size = 2,  relative to food 1, centred
model {
# PRIORS
   alpha[1] <- 0;       # zero contrast for baseline food
   for (k in 2:K){
      alpha[k] ~ dnorm(0,0.00001);  # vague priors
   }
   # Loop around lakes:
   for (k in 1:K){
      beta[1,k] <- 0;   # corner-point contrast with first lake 
   } 
   for (i in 2:I) {     
      beta[i,1] <- 0;   # zero contrast for baseline food
      for (k in 2:K) {
         beta[i,k] ~ dnorm(0,0.00001); # vague priors
      }
   }
   # Loop around sizes:
   for (k in 1:K){
      gamma[1,k] <- 0;  # corner-point contrast with first size 
   }
   for (j in 2:J) {     
      gamma[j,1] <- 0;  # zero contrast for baseline food
      for ( k in 2:K){
         gamma[j,k] ~ dnorm(0,0.00001); # vague priors
      }
   }

# LIKELIHOOD

   for (i in 1:I) {        # loop around lakes
      for (j in 1:J) {     # loop around sizes

   # Multinomial response
         X[i,j,] ~ dmulti( p[i,j,] , n[i,j]  );
         for (k in 1:K) {  # loop around foods
            p[i,j,k]        <- phi[i,j,k] / sum(phi[i,j,]);
            log(phi[i,j,k]) <- alpha[k] + beta[i,k]  + gamma[j,k];
         }
      }  
   }

# TRANSFORM OUTPUT TO ENABLE COMPARISON WITH AGRESTI'S RESULTS

   for (k in 1:K) {        # loop around foods
      for (i in 1:I) {     # loop around lakes
         b[i,k] <- beta[i,k] - mean(beta[,k]);   # sum to zero constraint
      }
      for (j in 1:J) {     # loop around sizes
         g[j,k] <- gamma[j,k] - mean(gamma[,k]); # sum to zero constraint
      }
   }

# FITTED VALUES
   for (i in 1:I) {     # loop around lakes
      for (j in 1:J) {     # loop around sizes
         for (k in 1:K) {     # loop around foods
            E[i,j,k] <- p[i,j,k] * n[i,j];
            OlogOE[i,j,k] <- X[i,j,k] * log( X[i,j,k] / E[i,j,k] );
         }
      }  
   }
   G2 <- 2 * sum( OlogOE[,,] );
}