0
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – whuber Mar 15 '13 at 22:03
1
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ 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)? $\endgroup$ – user22062 Mar 17 '13 at 1:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.