Suppose we have an outcome variable $y_{ji}$ which is a count of behaviors performed by group $j$ in round $i$, for $j = 1,...,n$ and $i = 1,...,8$. The outcome $y_{ji}$ counts are non-independent within each group $j$. Possible values for $y$ are bounded between 0 and 12, and because the behaviors appear to follow an almost "all or nothing" kind of pattern, the outcome is a u-shaped distribution:

enter image description here

My question is: how should we set up the beta-binomial model when the data has a multilevel / random-effects component?

I was planning on modelling the outcome data $y_{ji}$ with a Beta-binomial($m$, $\alpha$, $\beta$) distribution, where $y_{ji}$ are binomial observations with a number of trials $m = 12$, and where the probabilities $\pi_1,...,\pi_n$ follow a Beta($\alpha$, $\beta$) distribution. This way the U-shaped outcome distribution could be modeled with $0 < \alpha < 1$ and $0 < \beta < 1$.

In a linear random-effects model I would model a random intercept (or more) varying by group, and this would account for the non-independent error terms within each group's 8 rounds of behavior. Is it possible to do that with a beta-binomial model?

Here is what I've tried so far:


A hierarchical / random-effects GLMM package in R which says it can model beta-binomial by modelling the fixed-effects portion with binomial(link=logit), and the "random-effects" portion with Beta(link=logit).

question: If I understand correctly, the random-effects portion is modelling the dispersion parameter (they call it $\alpha$...). Is this "random-effect" the same as the random effect in a HLM / MLM / linear random-effects model, in that it would account for the within-group non-independence?


Another package that has an explicit betabin function that models the beta-binomial mean parameter $\mu$ as a function of any number of explanatory variables. The dispersion parameter $\phi$ can only be modeled with one explanatory variable.

question: is $\phi$ being treated as the "random effect"? If so, I have the same question as above.


Dr. Bolker is still working on implementing the beta-binomial, so I can't use this.


Perhaps I should just model it myself?

Thank you all for any help (pointers to papers or books would be appreciated).

  • $\begingroup$ I'm given to understand that Stan can generate posterior samples for any log-posterior you can code up in R. (I have yet to use it myself.) $\endgroup$
    – Cyan
    Commented Jun 20, 2013 at 16:29
  • $\begingroup$ I don't doubt that it can handle it (so could Jags, I believe), but I guess I'm really wondering what it means that the dispersion parameter is a 'random effect,' and whether that is the same kind of random effect that we're used to from Linear Mixed Models. Any ideas? $\endgroup$ Commented Jun 22, 2013 at 3:33
  • $\begingroup$ Stan is explicitly built for the purpose of computing hierarchical models: Hamiltonian Monte Carlo is generally more efficient in sampling from the hierarchical posterior than typical MCMC routines. $\endgroup$
    – Sycorax
    Commented Jan 5, 2014 at 15:18
  • $\begingroup$ As for the Bayesian way, you can use an R package "Rgbp" that produces independent posterior samples of all the model parameters via an acceptance-rejection sampling. $\endgroup$
    – Tak
    Commented Nov 20, 2016 at 20:04

1 Answer 1


In linear random effects models, the additional source of variability due to the random effect results in an additive increase in the total variance. In the beta-binomial model the additional variability is accounted for by a multiplicative overdispersion factor $\phi$. The random effect here is modeled implicitly.

\begin{align} & y_{ji}|\pi_{j} \sim \textit{Ber}(\pi_{j}) , \quad \pi_j \sim Beta(\mu,\rho) \nonumber \\ & E(\pi_{j}) = \mu \quad \textbf{Var}(\pi_{j}) = \rho\mu(1-\mu) \end{align}

$\mu$ and $\rho$ are the mean and intra-class correlation and can be reparametrized in terms of $\alpha$, $\beta$ of the Beta distribution.

Together the mean and variance of the marginal proportions $y_{j.}$ look like

\begin{align} E(y_{j.}) &= m\pi_{j} & \\ \text{Var}(y_{j .}) &= m\mu(1-\mu) + m(m-1)\rho\mu(1-\mu) & \\ & = m\mu(1-\mu)\phi & \nonumber \end{align}

The effect of assuming that $\pi_j$ is beta distributed is an overdispersion $\phi = (1+(m-1)\rho)$.

If your goal is inference on the mean proportions, method of moments estimates for $\mu$ and $\rho$ are quite effective, see [1]. These are particularly good estimates, if the number of observations per group are balanced.


[1] Kleinman, Joel C. "Proportions with extraneous variance: single and independent samples." Journal of the American Statistical Association 68.341 (1973): 46-54.

  • $\begingroup$ Shouldn't the parameters of beta distribution be α & β instead μ& ρ here ? $\endgroup$
    – sree22
    Commented Nov 9, 2017 at 4:33
  • $\begingroup$ This is a known reparametrization of beta-bernoulli distributions in terms of mean and overdispersion parameters $\endgroup$
    – user3303
    Commented Nov 19, 2017 at 23:08

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.