model
{

 for (i in 1:N) {
   r[i] ~ dbin(p[i], n[i]);

# alternative link functions (try one at a time):   
  logit(p[i])   <- alpha.star + beta*(x[i]-mean(x[])); 
#  probit(p[i])  <- alpha.star + beta*(x[i]-mean(x[]));      
#  cloglog(p[i]) <- alpha.star + beta*(x[i] - mean(x[]));   
      
   r.hat[i] <- p[i]*n[i];  # fitted values
 }
 alpha.star ~ dnorm(0.0, 1.0E-3);
 beta      ~ dnorm(0.0, 1.0E-3);
 alpha     <- alpha.star - beta*mean(x[]);

}


data
list(x = c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.8610, 1.8839),
     n = c(59, 60, 62, 56, 63, 59, 62, 60),
     r = c(6, 13, 18, 28, 52, 53, 61, 60), N = 8
)


# Here is a version of the model that lets you see how WinBUGS is calculating
# the log likelihood ;

model
{

 #  constants used in computing normalizing constants of log likelihood

 for (i in 1:N) {
	lfactn[i] <- logfact(n[i])
	lfactr[i] <- logfact(r[i])
	lfactnminr[i] <- logfact(n[i] - r[i])
}

 # saturated log-likelihoods ;
 for (i in 1:(N-1)) {
   llike.sat[i] <- lfactn[i] - lfactr[i] - lfactnminr[i] +  
          r[i]*log(r[i]/n[i]) + (n[i]-r[i])*log(1-r[i]/n[i]);
 }
 llike.sat[N] <- lfactn[N] - lfactr[N] - lfactnminr[N] + 0 ;  
                      # assign this one separately to avoid error due to ;
                      # trying to evaluate log(0) ;

 for (i in 1:N) {
   r[i] ~ dbin(p[i], n[i]);

# alternative link functions (try one at a time):   
#  logit(p[i])   <- alpha.star + beta*(x[i]-mean(x[])); 
#  probit(p[i])  <- alpha.star + beta*(x[i]-mean(x[]));      
  cloglog(p[i]) <- alpha.star + beta*(x[i] - mean(x[]));   

# log likelihood for sample i ;
  llike[i]     <- lfactn[i] - lfactr[i] - lfactnminr[i] +  
              r[i]*log(p[i]) + (n[i]-r[i])*log(1-p[i]);  
      
   r.hat[i] <- p[i]*n[i];  # fitted values
 }
 alpha.star ~ dnorm(0.0, 1.0E-3);
 beta      ~ dnorm(0.0, 1.0E-3);
 alpha     <- alpha.star - beta*mean(x[]);

 sumllik <- sum(llike[])       # -2 times this is what WinBUGS calls the deviance
 D <- 2 * (sum(llike.sat[]) - sumllik);    # frequentist deviance
}