Function to fit a beta-binomial exchangeable model

fit.model.R

fit.model <- function(data){
   # data is a list with two components 
   # - y: binomial counts
   # - n: binomial sample sizes
    require(LearnBayes)
    mode <- laplace(betabinexch, c(1, 1), 
                    cbind(data$y, data$n))$mode
    eta <- exp(mode[1]) / (1 + exp(mode[1]))
    K <- exp(mode[2])
    list(eta=eta, K=K, 
         d=data.frame(data, est=(data$y + K * eta) / (data$n + K)))
  }