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 Jun 20 '13 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$ – Christopher Poile Jun 22 '13 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 Jan 5 '14 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 Nov 20 '16 at 20:04

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 Nov 9 '17 at 4:33
  • $\begingroup$ This is a known reparametrization of beta-bernoulli distributions in terms of mean and overdispersion parameters $\endgroup$ – user3303 Nov 19 '17 at 23:08

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.