4
$\begingroup$

I understand how to set up a binomial probability distribution. I'm trying to extend my understanding to the beta-binomial. On Wikipedia, there is a beta-binomial example given at https://en.wikipedia.org/wiki/Beta-binomial_distribution. I can follow the computations to get m1 and m2 and from there (6.23 and 42.31 respectively, n=12), hat-alpha and hat-beta (using the method of moments -- 34.1350 and 31.6085 respectively). I can reproduce in the example table the values for the row "Predicted (Binomial p = 0.519215)", but I can't figure out how to reproduce the values for the row "Predicted (Beta-Binomial)". How were those values computed?

Here's the data from the table (the number of male children among the first 12 children of family size 13 in 6115 families taken from hospital records in 19th century Saxony, so x ranges from 0-12):

Number of males:

0 1 2 3 4 5 6 7 8 9 10 11 12

Observed:

3 24 104 286 670 1033 1343 1112 829 478 181 45 7

Beta:

2.3 22.6 104.8 310.9 655.7 1036.2 1257.9 1182.1 853.6 461.9 177.9 43.8 5.2
$\endgroup$

1 Answer 1

2
$\begingroup$

Your example data is contained in the R package vcd, so I will use it from there. First, the betabinomial distribution has two parameters which we can estimate by maximum likelihood. Let us do that in R:

data(Saxony, package="vcd")
    
library(bbmle)
# negative loglikelihood function:
    
bb_nloglik0 <- function(x) {
        function(alpha, beta) {
             -sum(extraDistr::dbbinom(seq(from=0, to=length(x)-1, 
               by=1), length(x)-1, alpha, beta, log=TRUE)*x)
            }
    }
    
bb_nloglik <- bb_nloglik0(Saxony)  

Then finding the maximum likelihood estimates (by minimizing the negative log likelihood defined above)

mod <- bbmle::mle2(bb_nloglik, list(alpha=50, beta=50))
mod.prof <-bbmle::profile(mod)
coef(mod)
       alpha     beta 
    34.10599 31.58112 

For good measure we can also find confidence intervals by likelihood profiling, even if you didn't ask:

confint(mod.prof)
             2.5 %   97.5 %
    alpha 27.38399 44.53472
    beta  25.35719 41.23744

And, finally, we can calculate the expected number of boys under the estimated beta-binomial model, by calculating the point probabilities and multiplying by N:

    round(rbind(Obs=Saxony,  Exp=extraDistr::dbbinom(0:12, 12, 
                 coef(mod)[1], coef(mod)[2])*N), 1)
          0    1     2     3     4      5    6      7     8     9    10   11  12
    Obs 3.0 24.0 104.0 286.0 670.0 1033.0 1343 1112.0 829.0 478.0 181.0 45.0 7.0
    Exp 2.3 22.6 104.8 310.9 655.7 1036.2 1258 1182.2 853.6 461.9 177.9 43.8 5.2
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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