Deriving the Beta-Binomial Predictive Distribution

Why Integrate Over the Response Rate?

In a single-arm trial with a binary endpoint, the interim data inform an unknown response rate. That uncertainty matters when describing outcomes among patients who have not yet been observed: the response rate is not known, so a predictive distribution should account for its posterior uncertainty. The beta-binomial distribution does this by averaging a future binomial distribution over the interim posterior for the response rate.

This tutorial derives the distribution of the number of responders in a fixed group of future patients, conditional on the interim observations. It does not make a terminal decision; the distribution is the quantity from which later decision summaries can be calculated.

Definition: Suppose the response probability $p$ has a Beta$(a,b)$ prior, where $a>0$ and $b>0$. After observing $r$ responses among $n$ patients, let $Y$ be the number of responses among a further $m$ patients. Under the shared-response-rate model below, the conditional distribution of $Y$ given the interim data is beta-binomial with updated parameters $a+r$, $b+n-r$, and $m$.

Model and Interim Posterior

Let each patient’s endpoint be a binary response indicator. The model assumes that, conditional on one common response probability $p$, patient outcomes are independent and have response probability $p$. The prior is $p\sim\mathrm{Beta}(a,b)$, with density proportional to $p^{a-1}(1-p)^{b-1}$ for $0<p<1$. The interim sample contains $n$ patients, of whom $r$ responded.

The likelihood for the interim data, apart from a factor that does not depend on $p$, is $p^r(1-p)^{n-r}$. Multiplication by the prior gives a posterior density proportional to $p^{a+r-1}(1-p)^{b+n-r-1}$. Thus the interim posterior is $p\mid r,n\sim\mathrm{Beta}(a+r,b+n-r)$. Write $\alpha=a+r$, $\beta=b+n-r$, and $S=\alpha+\beta$ for the updated parameters and their sum.

For the future group, assume the number $m$ of patients is fixed, their binary outcomes are observed in the same way, and they share the same $p$ as the interim patients. Conditional on $p$, $Y\sim\mathrm{Binomial}(m,p)$. The predictive distribution averages this conditional binomial distribution over the posterior uncertainty in $p$.

Deriving the Beta-Binomial Probability

For $y=0,1,\ldots,m$, the probability of exactly $y$ future responses, conditional on the interim data, is the posterior average of the binomial probability. The beta function is $B(u,v)=\int_0^1 p^{u-1}(1-p)^{v-1}\,dp$, for $u,v>0$. A Beta$(\alpha,\beta)$ density is $p^{\alpha-1}(1-p)^{\beta-1}/B(\alpha,\beta)$.

\begin{aligned} P(Y=y\mid r,n) &=\int_0^1 P(Y=y\mid p)\,\pi(p\mid r,n)\,dp\\ &=\int_0^1 {m\choose y}p^y(1-p)^{m-y} \frac{p^{\alpha-1}(1-p)^{\beta-1}}{B(\alpha,\beta)}\,dp\\ &={m\choose y}\frac{1}{B(\alpha,\beta)} \int_0^1 p^{\alpha+y-1}(1-p)^{\beta+m-y-1}\,dp\\ &={m\choose y}\frac{B(\alpha+y,\beta+m-y)}{B(\alpha,\beta)}. \end{aligned}

The integral in the third line is the beta function with arguments $\alpha+y$ and $\beta+m-y$, giving the final expression. Equivalently, using $B(u,v)=\Gamma(u)\Gamma(v)/\Gamma(u+v)$, the probability can be written ${m\choose y}\Gamma(\alpha+y)\Gamma(\beta+m-y)\Gamma(S)/[\Gamma(\alpha)\Gamma(\beta)\Gamma(S+m)]$. All probabilities are nonnegative. They sum to one because, for each $p$, the binomial probabilities sum to one; integrating that sum against the posterior density also gives one.

Predictive probability mass function: For $y=0,\ldots,m$, $$ P(Y=y\mid r,n)={m\choose y}\frac{B(a+r+y,\ b+n-r+m-y)}{B(a+r,\ b+n-r)}. $$ The support is the set of integer counts from zero through $m$. The distribution is conditional on the interim data and the specified prior and model.

Worked Probability Calculations

Worked Example: Two Future Patients

Suppose a hypothetical trial uses a Beta$(2,3)$ prior. At the interim analysis, 4 of 10 patients have responded. Find the predictive probabilities for the number of responses among the next $m=2$ patients.

The posterior is Beta$(2+4,3+10-4)=\mathrm{Beta}(6,9)$. Thus $\alpha=6$, $\beta=9$, and $S=15$. For a beta-binomial distribution with $m=2$, the rising-factorial form of the probability is ${2\choose y}(\alpha)_y(\beta)_{2-y}/(S)_2$, where $(q)_k=q(q+1)\cdots(q+k-1)$ and $(q)_0=1$. The denominator is $(15)_2=15(16)=240$.

For zero responses, $P(Y=0\mid r,n)=9(10)/240=90/240=0.3750$. For one response, $P(Y=1\mid r,n)=2(6)(9)/240=108/240=0.4500$. For two responses, $P(Y=2\mid r,n)=6(7)/240=42/240=0.1750$. These add to $0.3750+0.4500+0.1750=1.000$, as required. The expected number of responses is $2(6/15)=0.8000$ patients.

Worked Example: The Full Distribution for Four Future Patients

Consider another hypothetical interim result: a Beta$(2,2)$ prior and 3 responses among 8 patients. Calculate the complete predictive distribution for $m=4$ additional patients.

The interim posterior is Beta$(5,7)$, so $\alpha=5$, $\beta=7$, and $S=12$. For integer parameters, use $P(Y=y\mid r,n)={4\choose y}(5)_y(7)_{4-y}/(12)_4$. The common denominator is $(12)_4=12(13)(14)(15)=32{,}760$.

For $y=0$, the numerator is $1(1)(7\cdot8\cdot9\cdot10)=5{,}040$, so the probability is $5{,}040/32{,}760=0.1538$. For $y=1$, it is $4(5)(7\cdot8\cdot9)=10{,}080$, giving $0.3077$. For $y=2$, it is $6(5\cdot6)(7\cdot8)=10{,}080$, also $0.3077$. For $y=3$, it is $4(5\cdot6\cdot7)(7)=5{,}880$, giving $0.1795$. For $y=4$, it is $5\cdot6\cdot7\cdot8=1{,}680$, giving $0.05128$. The numerators sum to $32{,}760$, and the rounded probabilities sum to approximately 1.000.

The predictive mean is $m\alpha/S=4(5/12)=1.667$ responses among these four patients. The distribution assigns probability to each possible count, not just to the mean; the mean need not itself be a possible observed count.

Predictive Mean, Variance, and Posterior Dependence

The predictive mean follows by conditioning on $p$: the conditional mean is $mp$, so its posterior average is $m\alpha/S$. The variance has two components. Conditional on $p$, the binomial variance is $mp(1-p)$; across possible posterior values of $p$, the conditional mean $mp$ also varies.

\begin{aligned} E(Y\mid r,n)&=E\{E(Y\mid p,r,n)\mid r,n\} =mE(p\mid r,n)=\frac{m\alpha}{S},\\ \operatorname{Var}(Y\mid r,n) &=E\{\operatorname{Var}(Y\mid p,r,n)\mid r,n\} +\operatorname{Var}\{E(Y\mid p,r,n)\mid r,n\}\\ &=mE\{p(1-p)\mid r,n\}+m^2\operatorname{Var}(p\mid r,n)\\ &=\frac{m\alpha\beta(S+m)}{S^2(S+1)}. \end{aligned}

For the last equality, a Beta$(\alpha,\beta)$ variable has mean $\alpha/S$ and variance $\alpha\beta/[S^2(S+1)]$. Also $E[p(1-p)]=\alpha\beta/[S(S+1)]$. Substituting these into the variance decomposition gives $m\alpha\beta/[S(S+1)]+m^2\alpha\beta/[S^2(S+1)]$, which factors to the displayed result. This derivation makes clear why parameter uncertainty contributes to predictive variation.

Worked Example: Comparing Predictive and Plug-In Variance

Use the Beta$(5,7)$ posterior from the preceding example and $m=4$. Its mean response rate is $5/12$. Compare the beta-binomial predictive variance with the variance from a binomial model that fixes $p$ at this posterior mean.

The beta-binomial variance is $4(5)(7)(12+4)/[12^2(12+1)]=2{,}240/1{,}872=1.197$ responses squared. The plug-in binomial variance is $4(5/12)(7/12)=140/144=0.9722$ responses squared. The beta-binomial predictive distribution has the larger variance because it includes uncertainty about $p$, while the plug-in calculation treats $p=5/12$ as known.

This does not mean that the posterior mean is an unreasonable summary. It means that using it as if it were the true, fixed response rate discards uncertainty relevant to the distribution of future counts.

Integrating over a shared uncertain $p$ also induces dependence among future outcomes after $p$ is marginalized out. For two distinct future patients with binary response indicators $X_1$ and $X_2$, conditional independence gives $E(X_1X_2\mid r,n)=E(p^2\mid r,n)$. Therefore $\operatorname{Cov}(X_1,X_2\mid r,n)=\operatorname{Var}(p\mid r,n)>0$. The future indicators are independent given $p$, but generally not independent given only the interim data.

Worked Example: More Predictive Probability at the Extremes

Suppose the interim posterior is Beta$(3,3)$, and consider $m=3$ future patients. The posterior mean response rate is $1/2$. Compare the probability of no future responses under the beta-binomial model with the calculation that fixes $p$ at that mean.

Integrating over the posterior gives $P(Y=0\mid r,n)=B(3,6)/B(3,3)=(3\cdot4\cdot5)/(6\cdot7\cdot8)=60/336=0.1786$. The fixed-rate binomial calculation is $(1-1/2)^3=0.1250$. The beta-binomial probability is larger because some posterior values of $p$ are below its mean, making a run of zero responses more plausible. A corresponding increase in the probability of high counts is possible as well; the full distribution, rather than one tail alone, describes this spread.

Assumptions, Diagnostics, and Alternatives

The derivation requires a fixed number $m$ of future outcomes, a proper Beta prior with positive parameters, and a common response probability for interim and future patients. Conditional on that probability, outcomes must follow the assumed independent Bernoulli model. The beta-binomial calculation is exact under these assumptions; it does not correct a misspecified response model.

Before relying on the distribution, check that the endpoint definition is consistent across patients and that the interim count $r$ and denominator $n$ are correct. Confirm that the future cohort size is defined as a count of patients whose outcomes will be observed, and that the trial population, treatment, follow-up, and response ascertainment are plausibly stable. If patients differ systematically in prognosis, outcomes are clustered, response rates change over time, or an outcome is missing in a way related to response, a single common $p$ and simple binomial likelihood may not represent the data-generating process. A richer model may be needed.

A useful alternative is the plug-in binomial distribution $\mathrm{Binomial}(m,\bar p)$, where $\bar p=\alpha/S$ is the posterior mean. It is simpler and has the same predictive mean, but ignores posterior uncertainty and generally understates predictive variance. The beta-binomial integrates that uncertainty while retaining a common-rate model. Other alternatives, such as models with patient-level predictors or multiple response rates, can represent heterogeneity but require additional assumptions and parameters.

A common mistake is to say, “The posterior mean is $0.42$, so future outcomes are binomial with $p=0.42$.” That substitutes one estimate for an unknown parameter and removes uncertainty. The correction is to average the conditional binomial probabilities over the entire posterior, using the beta-binomial formula. Conversely, a beta-binomial model is not automatically more realistic: if the shared-$p$ assumption is unsuitable, integrating a beta posterior does not fix that failure.

Key takeaway: Given $r$ responses among $n$ patients, a Beta$(a,b)$ prior yields a Beta$(a+r,b+n-r)$ posterior. Integrating a future binomial likelihood over this posterior produces a beta-binomial distribution that reflects both future sampling variation and uncertainty about the response rate.

Check Your Understanding

Answer each question, then compare with the brief solutions.

With a Beta$(1,2)$ prior and 2 responses among 5 patients, what is the interim posterior?
Under a Beta$(4,6)$ posterior, what is the probability of zero responses among two future patients?
Why does integrating over $p$ generally produce greater predictive variance than fixing $p$ at its posterior mean?
For a beta-binomial model, what is the support of the future responder count $Y$ when $m=5$?
Name one assumption whose failure could make the shared-rate beta-binomial model unsuitable.

Solutions. (1) The posterior is Beta$(1+2,2+5-2)=\mathrm{Beta}(3,5)$. (2) The probability is $B(4,8)/B(4,6)=(6\cdot7)/(10\cdot11)=42/110=0.3818$. (3) The integrated variance includes variation in the conditional mean $mp$ across posterior values of $p$, in addition to binomial variation at a fixed $p$. (4) $Y$ can take the integer values $0,1,2,3,4,5$. (5) Examples include materially changing response rates over time, dependence between patients not represented by the shared-$p$ model, or outcome ascertainment that differs across patients.