# log-likelihood for a beta-binomial in R

Although I may ask a question that is already solved (but I found none that would explicitly refer to this), I would like to know, if (and I am no statistician) I am right with programming the log-likelihood function of a beta-binomial distribution in R as follows:

sum(lgamma(theta) - lgamma(Pi * theta) - lgamma((1 - Pi) * theta) + lgamma(n + 1) -
lgamma(Y + 1) - lgamma((n - Y) + 1) + lgamma(Y + Pi * theta) +
lgamma(n - Y + (1 - Pi) * theta) - lgamma(n + theta))


Here, theta is used as an overdispersion parameter. n is the number of trials, Y are the actual successes (from the real data), and Pi is the probability of success as derived from my model in JAGS. I (naively) took the density function as stated in Bolker (2008) and tried to take the logarithm of it. After that, I used this to compute an AIC and it gave reasonable numbers, but I am not sure!

P.S. I know there are caveats about using AIC, especially in a bayesian context, but I just want to know if I am right in the way to derive the answer to this, although trying it may be foolish from a philosophical point of view in the first place...

• Just found the aodml-function from package "aods3", which implements a beta-binomial regression from a frequentist point of view. The log-likelihood provided by aodml differs only 0.11 from my version, and the AIC only by 0.09. – Woosah Aug 9 '14 at 11:30