# Understanding the multilevel / random-effects beta-binomial regression model

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:

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:

hglm

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?

aod

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.

bayesian

Perhaps I should just model it myself?

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

• 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.)
– Cyan
Jun 20, 2013 at 16:29
• 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? Jun 22, 2013 at 3:33
• 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.
– Sycorax
Jan 5, 2014 at 15:18
• 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.
– Tak
Nov 20, 2016 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.

References

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

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