beta-binomial distribution with R

Question

I am studying an experiment of the kind: Let $n_{ij}$ be the number of fetuses, $X_{ij}$ the number of responses i.e. the number of fetuses with a malformation in the jth litter of the ith dose level for j=1,...,25 and i=1,...,5 . Then, $p_{ij}$ is the probability of response of in the jth litter of the ith dose level and hence we have: $$ P(X_{ij} = x_{ij}|p_{ij}) \sim Bin(n_{ij},p_{ij}) $$ But the probability of response $p_{ij}$ follow a beta distribution hence $$ P(p_{ij})=B^{-1}(\alpha_i , \beta_i )p_{ij}^{x_{ij}}(1-p_{ij})^{n_{ij}-x_{ij}} $$ and hence, at the end, $X_{ij}$ follow a beta-binomial distribution.

My problem is that I have to generate the number of responses $X_{ij}$ but I'm having some troubles.

The data that I have are all the $n_{ij}$ and $p_1, p_2, p_3, p_4, p_5$ (and this probabilities follow a logistic model) i.e. the probability response for each dose group, hence I don't have $p_ij$.

What should I do? I think that I should first generate $p_{ij}$ using the fact that they follow a beta distribution. But in which way should I do? How to estimate the parameter $\alpha_i$ and $\beta_i$?

Maybe someone has some ideas.. Thank you in advance!

Answer 1

In a Bayesian approach, you would not estimate $\alpha_i$ and $\beta_i$, but these would be prior beliefs supplied by you.

I am not so sure what $i$ and $j$ in your problem exactly represent, but here is some basline code for the Beta-Binomial distribution that hopefully gets you started. It is similar in spirit to a Gibbs sampler:

rm(list=ls())
library(VGAM)

R <- 10000
T <- 20000
x <- matrix(NA,R+T,2)

n <- 10
alpha <- 7
beta <- 19

x[1,] <- c(1,.5)

for(i in 2:(R+T)) {
  x[i,2] <- rbeta(1,alpha,beta)
  x[i,1] <- rbinom(1,n,x[i,2])
}
x <- x[(R+1):(R+T),]

plot(table(x[,1])/T, col="sienna4",type="p",pch=20)
betabinomialdensity <- function(x,n,alpha,beta) {choose(n, x)*beta(alpha+x,beta+n-x)/beta(alpha,beta)}
points(0:n,betabinomialdensity(0:n,n,alpha,beta),type="o", pch=22, lty=2, col="red")     

AvgDraws = mean(x[,1])
TheoreticalExpectedValue = n*alpha/(alpha+beta)