# Estimate beta binomial distribution

I have a dataset as follows:

ni-xi 76 48 38 19 20 14 26 29 39 45 38 36 32 34 26 21 23 12 24 18 14 61 35 21 32 20 89 30 34 35

xi    0  0 16 26  1  8 23 19 24 16 16 19  0  0  0  6  8  3  5  9  5  0 34 17 13  0  1 11 13 22


I want to model this dataset as a beta-binomial distribution: $$x_i \sim \mathcal{Binomial}(p, n_i)$$ where $$p \sim \mathcal{Beta}(α, β)$$.

I use VGAM package to do the fitting.

fit=vglm(cbind(ni-xi, xi) ~ 1, betabinomial.ab, trace = TRUE)
coef(fit)

VGLM    linear loop  1 :  loglikelihood = -95.9867
VGLM    linear loop  2 :  loglikelihood = -95.6633
VGLM    linear loop  3 :  loglikelihood = -95.6483
VGLM    linear loop  4 :  loglikelihood = -95.6477
VGLM    linear loop  5 :  loglikelihood = -95.6477
VGLM    linear loop  6 :  loglikelihood = -95.6477


However, I get the results as follows:

(Intercept):1 (Intercept):2
0.9006006    -0.3794863


Can anyone tell me why I get the negative parameters? How can i solve this problem.

• Welcome to our site! Please take a moment to read our faq and consult the markup help, then go back and edit your question to make it more readable. In the meantime, it is quite clear your data won't be fit well: you have a large scatter of values plus a concentration of zeros. Your approach appears doomed. You might perhaps want to ask a more fundamental question, such as what would be an appropriate model. If you do, please tell us where your data come from and what analyses you are really interested in.
– whuber
Mar 15, 2013 at 22:03

I think you're supposed to use Coef instead of coef ... that seems like a bad design on the part of the package authors, but maybe they have a good reason for that.

ni <- scan()
76 48 38 19 20 14 26 29 39 45 38 36 32 34 26 21 23 12 24 18 14 61 35 21 32 20 89 30 34 35

xi <- scan()
0 0 16 26 1 8 23 19 24 16 16 19 0 0 0 6 8 3 5 9 5 0 34 17 13 0 1 11 13 22

library(VGAM)
fit=vglm(cbind(ni, xi) ~ 1, betabinomial.ab, trace = TRUE)
coef(fit) # not what you want
Coef(fit) # probably what you want


When you run this, you get

> coef(fit) # not what you want
(Intercept):1 (Intercept):2
0.9006006    -0.3794863
> Coef(fit) # probably what you want
shape1    shape2
2.4610809 0.6842128


Reference page 132 of http://cran.r-project.org/web/packages/VGAM/VGAM.pdf

• Thanks. Do you think the beta binomial is a good model to fit the above data set? Actually the above data set can also be assumed to come from two binomial discription, the p is close to 0 and 0.5 repectively. Most of the data in this dataset follows the binomial discription which p is close to 0. Is there a better mothod to fit this dataset (estimate the p close to 0)? Mar 17, 2013 at 1:14