3
$\begingroup$

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!

$\endgroup$
1
$\begingroup$

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)
$\endgroup$

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.