A Worked Gamma-Poisson Allocation Example

From Event Counts to Allocation Probabilities

When a trial outcome is a count of events observed over follow-up, the count alone does not describe the evidence about a treatment’s event rate. Four events during 10 person-years and four events during 100 person-years imply different rates. A Gamma-Poisson model incorporates both the observed event count and its exposure, then gives a posterior distribution for each treatment’s rate. Thompson sampling uses those distributions to calculate the probability that each treatment has the best rate.

Throughout this tutorial, events are adverse outcomes, so a smaller event rate is preferred. The calculations would reverse the selection direction if larger rates represented better outcomes. All examples are hypothetical.

This model assumes that, within each treatment, events follow a Poisson process with a constant rate over the observed exposure; exposure is known; and the prior distributions across treatments are independent. It also assumes that the event counts are recorded and classified correctly. These conditions can be questionable if rates change materially over time, outcomes are clustered, or the count variance is much larger than the Poisson model permits. The basic update is useful only if the model is a reasonable description of the data.

The Gamma-Poisson Update

For a single treatment, the likelihood of an observed count \(y\), given exposure \(E\) and rate \(\lambda\), is proportional to \(\lambda^y e^{-E\lambda}\). Multiplying by the Gamma prior density gives a posterior kernel proportional to \(\lambda^{a+y-1}e^{-(b+E)\lambda}\). Therefore, the posterior is again Gamma:

The posterior mean is \((a+y)/(b+E)\) events per person-year. The prior shape \(a\) is updated by the event count, while the prior rate parameter \(b\) is updated by exposure. Exposure does not increase the shape: it contributes to the posterior rate parameter.

Best-Treatment Probabilities

For \(K\) treatments with independent posterior rates, define \(\pi_j\) as the posterior probability that treatment \(j\) has the smallest event rate. With continuous Gamma distributions, exact ties have probability zero. If each allocation is made by drawing one rate from each treatment’s posterior and choosing the treatment with the smallest draw, then the probability of allocating to treatment \(j\) is \(\pi_j\).

When all posterior shapes are 2, the calculation for three treatments has a convenient closed form. Write their posterior rate parameters as \(r_1,r_2,r_3\), so the posterior density for treatment \(j\) is \(r_j^2x e^{-r_jx}\), and its survival probability at \(x\) is \(e^{-r_jx}(1+r_jx)\). Conditioning on the value drawn for treatment \(j\), multiply its density by the probabilities that each other rate exceeds that value, then integrate:

The last equality follows from \(\int_0^\infty x^m e^{-Tx}\,dx=m!/T^{m+1}\) for \(T>0\), applied for \(m=1,2,3\). This result is specific to three independent Gamma posteriors with shape 2; it is not a general formula for arbitrary shapes or numbers of treatments.

A Two-Treatment Calculation With Unequal Shapes

The shape-2 formula is convenient, but posterior shapes need not match. For two independent posteriors, the best-treatment probability can always be expressed by integrating one posterior density against the other posterior’s survival function. In a simple unequal-shape case, that calculation is still exact.

Computation, Diagnostics, and Limits

For general Gamma shapes, or for many treatments, posterior best probabilities can be estimated by simulation. On each iteration, independently draw one rate from each treatment’s posterior and record which draw is smallest. The estimated probability for treatment \(j\) is the fraction of iterations on which it has the smallest draw. A sufficiently large simulation makes Monte Carlo error small, but it does not correct a mistaken model, prior, exposure, or event count.

Several checks are especially useful:

• Check the parameterization. In the shape-rate convention, posterior parameters are \(a+y\) and \(b+E\). A shape-scale implementation uses a different second parameter and can silently produce incorrect draws.
• Check the exposure units. If exposure is entered in days but the prior rate parameter is in person-years, convert the exposure before updating.
• Check that probabilities sum to one. Exact probabilities should sum to one; simulated estimates should be close, with small deviations attributable to Monte Carlo error.
• Check the selection direction. For adverse events, select the smallest sampled rate. Selecting the largest reverses the intended clinical objective.
• Check the Poisson assumption. Strong time variation, clustering, or excess variability can make a constant-rate Poisson model inadequate. Posterior precision can then be overstated.

Posterior best probability differs from both a plug-in greedy rule and predictive response probability. A greedy rule selects according to point summaries, such as the smallest posterior mean, and does not preserve the uncertainty-based randomization of Thompson sampling. A predictive response probability concerns the chance of an outcome for a future participant; it need not equal the probability that a treatment has the best underlying rate. The posterior probability here answers a specific question: which treatment has the smallest event rate, given the model and observed data?

These calculations rely on independent treatment posteriors and a specified clinical direction in which lower rates are better. They do not resolve whether the prior is appropriate or whether a rate alone captures clinical value. When posterior probabilities are computed numerically, compare estimates across larger simulation runs and, where possible, verify a small case with an exact calculation such as the examples above. A stable numerical answer is reassuring about computation, not proof that the model assumptions hold.