Likelihood Inference for Beta-Binomial Dispersion

Why Dispersion Inference Needs Special Care

A beta-binomial model allows successes within a group to be positively associated, so its success count can vary more than a binomial count with the same mean. Inference about that extra variation is not quite an ordinary one-parameter likelihood problem: the no-overdispersion model occurs at the edge of the allowed parameter space. This tutorial develops likelihood-based intervals and tests while making that boundary explicit.

Suppose there are $n$ independent groups. Group $i$ has $m_i$ trials and an observed success count $Y_i=y_i$. Write $p_i$ for its mean success probability and $\rho$ for the common intraclass correlation, or dispersion parameter. The mean probabilities may be specified directly or through a regression model; in the latter case, its regression coefficients are nuisance parameters for inference about $\rho$.

Definition: For $0<p_i<1$ and $0<\rho<1$, the beta-binomial model for group $i$ is $Y_i\sim\operatorname{BetaBinomial}(m_i,\alpha_i,\beta_i)$, where $\alpha_i=p_i(1-\rho)/\rho$ and $\beta_i=(1-p_i)(1-\rho)/\rho$. Its probability mass function is $$ f_i(y_i;p_i,\rho)= \binom{m_i}{y_i} \frac{B(y_i+\alpha_i,\;m_i-y_i+\beta_i)} {B(\alpha_i,\beta_i)},\qquad y_i=0,\ldots,m_i. $$ At $\rho=0$, define the model by its binomial limit: $Y_i\sim\operatorname{Binomial}(m_i,p_i)$. The parameter space for dispersion is $0\leq\rho<1$.

The independent-group likelihood is $L=\prod_{i=1}^n f_i(y_i;p_i,\rho)$, with log-likelihood $\ell=\sum_i\log f_i$. If $p_i=p_i(\gamma)$ is determined by regression coefficients $\gamma$, inference about $\rho$ must account for estimating $\gamma$, rather than pretending the probabilities are known. A group with $m_i=1$ has a Bernoulli marginal distribution that does not depend on $\rho$; such a group contributes no direct information about dispersion.

Building and Profiling the Likelihood

For computation, the beta-binomial mass can also be written with the rising factorial $(a)_k=a(a+1)\cdots(a+k-1)$, where $(a)_0=1$:

f_i(y_i;p_i,\rho)= \binom{m_i}{y_i} \frac{(\alpha_i)_{y_i}(\beta_i)_{m_i-y_i}} {(\alpha_i+\beta_i)_{m_i}}.

This form makes clear that each observed group contributes its own trial count $m_i$ and probability $p_i$. In general, one maximizes numerically over $\rho$ and the nuisance parameters. For each fixed $\rho$, let $\hat\gamma_\rho$ maximize the log-likelihood over $\gamma$. The profile log-likelihood for dispersion is $\ell_p(\rho)=\ell(\hat\gamma_\rho,\rho)$. If the $p_i$ are known, there is no nuisance maximization and $\ell_p(\rho)=\ell(\rho)$.

Worked Example: Likelihood Contributions With Unequal Trial Counts

Consider two independent groups with known mean probability $p=0.4$, common $\rho=0.2$, trial counts $m_1=2$ and $m_2=3$, and observed counts $y_1=y_2=1$. Calculate their joint likelihood.

The beta parameters are $\alpha=0.4(0.8)/0.2=1.6$ and $\beta=0.6(0.8)/0.2=2.4$, so $\alpha+\beta=4$. For $m=2,y=1$, the rising-factorial formula gives $2(1.6)(2.4)/(4\cdot5)=0.384$.

For $m=3,y=1$, it gives $3(1.6)(2.4)(3.4)/(4\cdot5\cdot6)=39.168/120=0.3264$. Independence therefore gives joint likelihood $0.384(0.3264)=0.1253376$. If $\rho=0$ instead, the two binomial probabilities are $2(0.4)(0.6)=0.48$ and $3(0.4)(0.6)^2=0.432$, with joint likelihood $0.20736$. These are likelihood values at specified parameter values, not a comparison of fitted models: a likelihood-ratio comparison requires maximizing each model as appropriate.

The Binomial Limit

The binomial model is not a separate, unrelated choice of distribution: it is the limit of the beta-binomial model as $\rho$ decreases to zero. This limit is important because it identifies the null hypothesis of no extra dispersion and shows why that null lies at the boundary.

Result: For fixed $m_i$, $p_i\in(0,1)$, and $y_i\in\{0,\ldots,m_i\}$, the beta-binomial probability mass converges as $\rho\downarrow0$ to the binomial mass: $$ \lim_{\rho\downarrow0} f_i(y_i;p_i,\rho) =\binom{m_i}{y_i}p_i^{y_i}(1-p_i)^{m_i-y_i}. $$

Proof. Put $s=\alpha_i+\beta_i=(1-\rho)/\rho$. Then $\alpha_i=p_i s$ and $\beta_i=(1-p_i)s$. By the rising-factorial expression, the ratio of the beta-binomial mass to the binomial coefficient is $(p_i s)_{y_i}((1-p_i)s)_{m_i-y_i}/(s)_{m_i}$. For any fixed nonnegative integer $k$, $(cs)_k/s^k=\prod_{j=0}^{k-1}(c+j/s)\to c^k$ as $s\to\infty$. Applying this finite-product limit to the numerator and denominator, and using $m_i=y_i+(m_i-y_i)$, the ratio tends to $p_i^{y_i}(1-p_i)^{m_i-y_i}$. Since $s\to\infty$ as $\rho\downarrow0$, the stated limit follows. $\square$

Likelihood-Ratio Intervals

A likelihood-based confidence set compares the profile likelihood at a candidate value $\rho_0$ with its maximum. Define the likelihood-ratio statistic $D(\rho_0)=2\{\ell_p(\hat\rho)-\ell_p(\rho_0)\}$, where $\hat\rho$ maximizes the profile likelihood over $0\leq\rho<1$. When $\rho_0$ is an interior value and regularity conditions hold, $D(\rho_0)$ has an asymptotic $\chi^2_1$ distribution under $\rho=\rho_0$. Inverting that test gives the usual approximate 95% profile-likelihood confidence set:

\left\{\rho_0\in(0,1):D(\rho_0)\leq \chi^2_{1,0.95}\right\}, \qquad \chi^2_{1,0.95}\approx3.841.

Here regularity means, in particular, that the candidate value is interior, the model is identifiable, the log-likelihood is sufficiently smooth, and the information for the parameters being estimated is nonsingular. For a confidence set that also considers $\rho_0=0$, the boundary null needs its own calibration; it should not automatically be treated as an interior candidate. The profile likelihood itself remains well-defined at zero by the binomial limit.

Worked Example: A Profile-Likelihood Confidence Set

Suppose $p=0.5$ is known and there are ten independent groups, each with $m_i=2$. Six groups have $y_i=0$ or $2$, and four have $y_i=1$. For a group, the probability of either particular extreme outcome is $(1+\rho)/4$, while the probability of $y_i=1$ is $(1-\rho)/2$. Ignoring factors constant in $\rho$, the log-likelihood is $\ell(\rho)=6\log(1+\rho)+4\log(1-\rho)+C$.

Differentiating gives $\ell'(\rho)=6/(1+\rho)-4/(1-\rho)$. Setting this to zero yields $6(1-\rho)=4(1+\rho)$, hence $\hat\rho=0.2$. The log-likelihood difference from its maximum is $\ell(\rho)-\ell(0.2)=6\log\{(1+\rho)/1.2\}+4\log\{(1-\rho)/0.8\}$.

For positive candidate values, the approximate 95% profile-likelihood set uses $2\{\ell(0.2)-\ell(\rho)\}\leq3.841$. At the boundary, $\rho=0$, this statistic is $-2[6\log(1/1.2)+4\log(1/0.8)]\approx0.403$, so zero is included. Solving $2\{\ell(0.2)-\ell(\rho)\}=3.841$ for the upper endpoint gives approximately $0.709$. Thus the set runs from the boundary value zero to about $0.709$; the lower endpoint is zero because the parameter cannot be negative and the boundary value is not rejected. The interval is wide relative to the point estimate, reflecting the limited information in ten groups.

Testing No Overdispersion at the Boundary

The hypothesis of no overdispersion is $H_0:\rho=0$ against $H_1:\rho>0$. Since the parameter space only allows nonnegative values, zero is a boundary point. Under suitable regularity conditions, the likelihood-ratio statistic for this test does not have the usual $\chi^2_1$ null distribution.

Theorem: Suppose independent groups provide an increasing amount of information, the null nuisance parameters are identifiable and interior, the efficient information for $\rho$ is positive, and the likelihood has the usual local quadratic expansion at $\rho=0$. Suppose also that zero is the only active boundary and the alternative allows $\rho\geq0$. Then under $H_0:\rho=0$, the likelihood-ratio statistic converges in distribution to $$ \frac12\chi^2_0+\frac12\chi^2_1, $$ where $\chi^2_0$ denotes a point mass at zero. The mixture is also called the chi-bar-square distribution for this one-boundary-parameter case.

Proof. After accounting for the regular nuisance parameters, the local log-likelihood for a displacement $h/\sqrt{n}$ in dispersion has the quadratic form $hZ_n-\tfrac12h^2 I+o_p(1)$, where $I>0$ is the efficient information and $Z_n$ converges in distribution to $N(0,I)$ under the null. The alternative restricts $h\geq0$. Maximizing the quadratic over that half-line gives $\hat h=\max(0,Z_n/I)$; under the null the constrained maximum is zero when $Z_n\leq0$. Twice the gain over the null is therefore $\{\max(0,Z_n)\}^2/I+o_p(1)$. Since $Z_n/\sqrt I$ converges to a standard normal variable $Z$, the limit is $\{\max(0,Z)\}^2$. A standard normal variable is positive with probability $1/2$, giving a zero statistic in that half of the cases; conditional on being positive, its square has the $\chi^2_1$ distribution. This proves the stated mixture limit. $\square$

For a 5% test using this mixture, the critical value $c$ satisfies $0.5+0.5P(\chi^2_1\leq c)=0.95$, so $c$ is the 90th percentile of $\chi^2_1$, approximately $2.706$. Equivalently, the asymptotic boundary-null $p$-value for a positive observed statistic $D$ is one half of the usual $\chi^2_1$ upper-tail probability. This calibration assumes the regularity conditions in the theorem; for example, if dispersion is not identifiable under the null, the mixture result need not apply.

Worked Example: Testing the Binomial Model

Let $p=0.5$ be known and again take groups of size two. Suppose ten groups yield eight extreme outcomes ($y=0$ or $2$) and two middle outcomes ($y=1$). Up to a constant, the log-likelihood is $\ell(\rho)=8\log(1+\rho)+2\log(1-\rho)+C$.

The score is $8/(1+\rho)-2/(1-\rho)$. Setting it to zero gives $8(1-\rho)=2(1+\rho)$, so $\hat\rho=0.6$. The likelihood-ratio statistic for $H_0:\rho=0$ is $D=2\{\ell(0.6)-\ell(0)\}=2\{8\log(1.6)+2\log(0.4)\}\approx3.855$.

An ordinary $\chi^2_1$ calibration gives a $p$-value of about $0.050$. The boundary calibration instead gives approximately half that upper-tail probability, about $0.025$, and compares $3.855$ with the 5% critical value $2.706$. It therefore rejects the binomial model at the 5% level under the stated asymptotic boundary conditions. The difference is a consequence of the one-sided parameter space, not a change to the fitted likelihood.

Practical Interpretation and Common Pitfall

Profile likelihood is useful because it respects the parameter space and incorporates the estimation of nuisance parameters. In practice, fit the model at a grid of fixed dispersion values, re-estimate the mean-model parameters at each value, and inspect the resulting profile. The likelihood can be flat near zero or asymmetric, so a symmetric interval based on a standard error may be a poor summary. A likelihood-based set also makes clear when the data are compatible with the binomial limit.

Common mistake: using the ordinary $\chi^2_1$ test at $\rho=0$. The tempting reasoning is that the likelihood-ratio statistic compares models differing by one parameter, so its reference distribution must be $\chi^2_1$. That rule relies on an interior null value. Here the alternative can move away from zero in only one direction; under the null, the fitted displacement is sometimes exactly zero. The resulting point mass changes the reference distribution to the half-$\chi^2_0$, half-$\chi^2_1$ mixture under the stated conditions. Use the ordinary $\chi^2_1$ calibration for regular interior candidate values, and the boundary calibration for the no-overdispersion test.

Key takeaway: Estimate dispersion by profiling the beta-binomial likelihood, but distinguish interior inference from inference at $\rho=0$. The binomial limit is a boundary model, and its likelihood-ratio test generally needs the mixture reference distribution rather than an ordinary one-degree-of-freedom chi-square distribution.

Check Your Understanding

Use the definitions and likelihood arguments above to answer these questions.

Why does a group with one trial provide no information about the beta-binomial dispersion parameter?
State the beta-binomial parameter values $\alpha_i$ and $\beta_i$ in terms of $p_i$ and $\rho$.
What is the asymptotic reference distribution for the likelihood-ratio test of $\rho=0$, and why does it differ from $\chi^2_1$?
When forming a profile-likelihood interval, what must be re-estimated at each fixed value of $\rho$ if the mean probabilities come from a regression model?
Under the boundary mixture, what 5% critical value is used for testing no overdispersion, and how is it obtained?