var 
  pattern[Npattern,T], # response pattern e.g. (1,1) or (2,4) etc. 
  Ncum[Npattern,T],    # cumulative total
  response[N,T],       # response for patient i in period t
  p[N,T,(Ncut+1)],     # prob of response j for patient i in period t 
  Q[N,T,Ncut],         # cumulative prob of response worse than j
                       # for patient i in period t
  group[N],            # treatment group (1=AB; 2=BA)
  mu[G,T],             # logistic mean for group g & period t
  treat[2,T], beta,    # treatment effect
  period[2,T], pi,     # period effect
  carry[2,T], kappa,   # carryover effect
  a[Ncut],             # cut points for latent response variable
  b[N],                # subject random effect
  tau,                 # precision of subject effects
  sigma, log.sigma;
data {
   # Construct individual response data from contingency table
   for (i in 1:Ncum[1,1]) { 
      group[i] <- 1; for (t in 1:T) { response[i,t] <- pattern[1,t] }
   }
   for (i in (Ncum[1,1]+1):Ncum[1,2]) { 
      group[i] <- 2; for (t in 1:T) { response[i,t] <- pattern[1,t] }
   }

   for (k in 2:Npattern) {
      for(i in (Ncum[k-1,2]+1):Ncum[k,1]) {
         group[i] <- 1; for (t in 1:T) { response[i,t] <- pattern[k,t] }
      }
      for(i in (Ncum[k,1]+1):Ncum[k,2]) {
         group[i] <- 2; for (t in 1:T) { response[i,t] <- pattern[k,t] }
      }
   }
}
model {
  for (i in 1:N) {
     for (t in 1:T) {
        response[i,t] ~ dordered.logit(-mu[group[i],t] + b[i], a[1:Ncut])
     }
     # Random effects
     b[i] ~ dnorm(0.0, tau);
  }

#
# Fixed effects
#
  for (g in 1:G) {
     for(t in 1:T) { 
        # logistic mean for group i in period t
        mu[g,t] <- beta*treat[g,t]/2 + pi*period[g,t]/2 + kappa*carry[g,t]; 
     }
  }                                                             
  beta ~ dnorm(0, 1.0E-06);
  pi ~ dnorm(0, 1.0E-06);
  kappa ~ dnorm(0, 1.0E-06);

# ordered cut points for underlying continuous latent variable  
  for(i in 1:3) {
     a[i] ~ dnorm(0, 1.0E-6);
  }
 
  tau ~ dgamma(1.0E-3, 1.0E-3);
  sigma <- sqrt(1/tau);
  log.sigma <- log(sigma);

}