Why Dependence Matters for Thompson Sampling
Thompson sampling for treatment allocation uses uncertainty about which treatment has the greatest response. When treatment effects share information through a common control, a shared model parameter, or hierarchical shrinkage, their posterior uncertainties are generally dependent. The probability that one treatment is best then depends on the joint posterior distribution, not just on each treatment’s posterior mean and variance considered separately.
For example, two treatment effects may tend to rise and fall together because both depend on an uncertain common baseline. Their individual posterior distributions may be quite wide while their difference is comparatively precise. Conversely, negative dependence can make their difference more variable. Sampling each treatment independently would discard this structure and can change the probability of selecting each arm.
In a continuous model, exact ties ordinarily have probability zero. In a discrete model or a model with point masses, ties may have positive probability; the allocation rule must then specify how tied maxima are handled. Here the focus is on continuous treatment posteriors. Once the \(q_k\) values have been computed, ordinary Thompson sampling assigns the next participant to arm \(k\) with probability \(q_k\), conditional on the current trial data and the model.
Computing Best-Arm Probabilities From Joint Draws
Let \(\boldsymbol{\theta}^{(m)}=(\theta_1^{(m)},\ldots,\theta_K^{(m)})\), for \(m=1,\ldots,M\), be draws from the joint posterior given \(D\). Each draw must represent a coherent set of treatment effects from the same posterior draw, including shared parameters and latent effects. For each draw, identify the arm with the largest value. The fraction of draws in which arm \(k\) is largest estimates \(q_k\):
This is a Monte Carlo estimate of a posterior probability, not a proportion of trial participants who responded. The treatment-effect draws are simulated from the posterior; they are not additional observed patients. When the draws are independent, the indicator in this expression is a Bernoulli random variable with probability \(q_k\), so the approximate Monte Carlo standard error is
If the draws come from a correlated Markov chain, this formula using the raw number of draws can understate uncertainty. The sequence of indicators for whether arm \(k\) is best is itself correlated; an effective sample size for that indicator, or a batch-based Monte Carlo error estimate, is more appropriate.
Worked Example: Two Correlated Normal Treatment Effects
Suppose a hypothetical trial compares two treatments, with higher response effect preferred. Conditional on the observed data, let their effects have the joint representation \(\theta_1=0.10+U+E_1\) and \(\theta_2=U+E_2\), where \(U\sim N(0,0.04)\), \(E_1\sim N(0,0.01)\), and \(E_2\sim N(0,0.01)\) are mutually independent. The effects are measured in response-rate units. Find the posterior probability that treatment 1 is best.
The shared component \(U\) makes the treatment effects positively correlated. Their difference is \(\theta_1-\theta_2=0.10+E_1-E_2\), because the common component cancels. Thus the difference has mean \(0.10\) and variance
Equivalently, each treatment effect has variance \(0.01+0.04=0.05\), and their covariance is \(0.04\), giving \(0.05+0.05-2(0.04)=0.02\). Therefore,
The posterior best-arm probabilities are about \(0.7602\) for treatment 1 and \(1-0.7602=0.2398\) for treatment 2. If one incorrectly sampled from the two marginal posteriors independently, the variance of the difference would be \(0.05+0.05=0.10\), giving \(\Phi(0.10/\sqrt{0.10})=\Phi(0.3162)\approx0.6241\). That is a different probability because independent sampling removes the posterior covariance.
Hierarchical Posteriors and the Joint-Draw Procedure
A hierarchical model often gives treatment effects a shared distribution, such as \(\theta_j\mid\mu,\tau\sim N(\mu,\tau^2)\), with a posterior for the common mean \(\mu\) and possibly for the between-treatment scale \(\tau\). Conditional on a shared parameter, arm-specific effects may be independent, but after integrating over uncertainty in shared parameters they need not be independent. A posterior draw must therefore include a draw of the shared parameters and the corresponding arm effects. Pairing the arm effects from the same draw preserves the dependence induced by the hierarchy.
The procedure requires a proper joint posterior under the specified model and data, a way to obtain representative draws from it, and a clearly defined direction of preference. If larger values are not better, the parameter must first be expressed on a scale where the stated maximization is appropriate. The model must also reflect the trial’s dependence structure; preserving the dependence in an incorrect model does not make the resulting probabilities reliable.
For each Monte Carlo iteration, draw any common parameters or latent quantities from their posterior distribution.
Generate the complete vector of arm effects conditional on that same shared state, following the fitted model.
Identify the arm with the largest sampled effect and add one to its best-arm count.
Divide each arm’s count by the number of draws; then assess Monte Carlo error, accounting for serial dependence if applicable.
Worked Example: Counting Maxima Across Three Arms
Consider a hypothetical three-arm trial. A hierarchical posterior sampler has produced the following 12 joint draws of treatment effects, in response-rate units. The columns in each row belong together: they were generated in the same iteration, including the shared posterior quantities.
| Draw | Arm A | Arm B | Arm C | Largest effect |
|---|---|---|---|---|
| 1 | 1.2 | 0.7 | 0.1 | A |
| 2 | 0.3 | 0.8 | 0.5 | B |
| 3 | 0.6 | 0.4 | 0.9 | C |
| 4 | 0.9 | 0.2 | 0.4 | A |
| 5 | 0.4 | 1.1 | 0.7 | B |
| 6 | 0.5 | 0.6 | 0.2 | B |
| 7 | 1.0 | 0.9 | 0.8 | A |
| 8 | 0.2 | 0.3 | 0.6 | C |
| 9 | 0.8 | 0.7 | 0.1 | A |
| 10 | 0.1 | 0.5 | 0.4 | B |
| 11 | 0.7 | 0.6 | 0.8 | C |
| 12 | 0.9 | 0.4 | 0.5 | A |
Arm A is largest in draws 1, 4, 7, 9, and 12, giving 5 maxima. Arm B is largest in draws 2, 5, 6, and 10, giving 4; Arm C is largest in draws 3, 8, and 11, giving 3. Therefore the estimated probabilities are
They sum to 1 because each draw has exactly one largest effect. If these were the only available draws, the estimates would be very imprecise; the example illustrates the calculation, not a recommended simulation size.
A Larger Simulation and Its Monte Carlo Error
In practice, the joint draw may come from a hierarchical model fitted to trial summaries or participant-level outcomes. To make the calculation concrete, suppose a hypothetical three-arm trial has observed response-effect summaries \(y=(0.35,0.10,-0.10)\), each with standard error \(0.20\). A model assigns the arm effects a common hierarchical distribution and samples its shared parameters and arm effects jointly. Assume for this illustration that the sampler has converged and has yielded 20,000 effectively independent joint draws. The specific posterior probabilities below are hypothetical simulation output, not values inferred from these summaries alone.
Worked Example: Estimating Probabilities From 20,000 Joint Draws
Suppose Arm A is the maximum in 11,440 draws, Arm B in 6,020, and Arm C in 2,540. The counts sum to 20,000. Dividing by that total gives
For Arm A, the independent-draw Monte Carlo standard error is \(\sqrt{0.5720(1-0.5720)/20000}\approx0.00350\). For Arm B it is \(\sqrt{0.3010(1-0.3010)/20000}\approx0.00324\), and for Arm C it is \(\sqrt{0.1270(1-0.1270)/20000}\approx0.00235\). These quantify simulation error in estimating the posterior probabilities, conditional on the model and data. They do not quantify uncertainty about whether the model is appropriate, nor do they turn a posterior probability into a frequentist confidence statement.
Under ordinary Thompson sampling, the next participant would be assigned to A, B, or C with probabilities \(0.5720\), \(0.3010\), and \(0.1270\), respectively, conditional on the current data. The largest probability does not imply that Arm A is certain to be best; it means that it was best in about 57.2% of the joint posterior draws.
Why the Joint Draw Is the Essential Unit
The justification for counting maxima is direct. Under the joint posterior, the event that arm \(k\) is best is represented by an indicator that equals 1 when \(\theta_k\) is the largest and 0 otherwise. Its posterior expectation is exactly the event probability:
A Monte Carlo average of these indicators therefore estimates \(q_k\). This argument depends on drawing the vector from the joint posterior. If one instead shuffles each arm’s draws separately, the resulting rows combine values that may never occur together under the posterior. The computed maxima then refer to an artificial independent distribution, not the fitted model.
A shared additive component is one special case: it induces positive correlation in arm effects but cancels from every pairwise difference, so it does not change which arm is largest in a draw. Correlation is not automatically consequential in the same direction for every model. The relevant question is how the dependence changes the joint distribution of contrasts and, ultimately, the ordering of the effects.
Diagnostics, Alternatives, and Limitations
| Approach | What it computes | When it is appropriate or problematic |
|---|---|---|
| Joint posterior draws | Best-arm frequencies using complete vectors from the fitted joint posterior | Appropriate for correlated or hierarchical posteriors when draws represent the specified joint model |
| Independent arm draws | Best-arm frequencies after drawing each arm separately from its marginal posterior | Appropriate only when the joint posterior factorizes across arms; otherwise it imposes a different dependence structure |
| Plug-in greedy choice | Selects the arm with the largest posterior mean | Does not calculate posterior best-arm probabilities and does not perform Thompson sampling |
Several checks help detect problems. Verify that every iteration contains one value for every arm from the same posterior state. Confirm that the maximum is selected on the intended response scale and that no arm labels are shifted or reordered. Check that the estimated probabilities sum to 1, up to rounding, when there is one winner per draw. Increase the number of draws and compare estimates; large changes indicate inadequate Monte Carlo precision. For Markov chain draws, inspect convergence and account for autocorrelation in the best-arm indicators, not only in the individual treatment effects.
The method also has substantive limitations. Best-arm probabilities depend on the likelihood, prior, hierarchical structure, and observed data. A model that assumes exchangeability among treatments may shrink effects toward a common distribution in ways that are questionable when treatments are not genuinely comparable. Further, the probability of having the largest modeled response parameter is not by itself a statement about clinical importance or predictive response probability. The introductory material on best-arm probability versus predictive response probability distinguishes those quantities; the calculation here concerns the former.
Common Mistakes
- Multiplying or sampling marginal posteriors as if they were independent. This is justified only if the joint posterior actually factorizes. Otherwise it discards the dependence that the joint draws must preserve.
- Choosing the arm with the largest posterior mean and calling that Thompson sampling. A largest-mean rule is a plug-in choice. Thompson sampling randomizes according to the posterior probability of being best.
- Treating a small Monte Carlo error as evidence that the model is correct. Monte Carlo precision measures simulation noise under the fitted posterior; it does not diagnose prior sensitivity or model misspecification.
- Pairing arm values from different iterations. A row of treatment effects must be a joint draw. Rearranging columns independently changes the dependence and can change the estimated allocation probabilities.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- For two arms, a joint posterior gives \(\theta_1-\theta_2\sim N(0.06,0.04)\), where the second argument is the variance. What is the posterior probability that Arm 1 is best?
- In 5,000 independent joint draws, Arm A is largest 2,500 times and Arm B is largest 1,500 times. What is the estimated best-arm probability for the remaining arm?
- Why can independently shuffling the columns of a set of hierarchical posterior draws give incorrect best-arm probabilities?
- A simulation reports \(\widehat q_A=0.40\) based on 10,000 independent draws. Give its approximate Monte Carlo standard error and state what it measures.
Solutions. (1) The standard deviation of the difference is \(\sqrt{0.04}=0.20\), so \(P(\theta_1>\theta_2\mid D)=\Phi(0.06/0.20)=\Phi(0.30)\approx0.6179\). (2) The remaining arm wins \(5000-2500-1500=1000\) draws, so its estimated probability is \(1000/5000=0.2000\). (3) Shuffling breaks the within-draw dependence induced by shared parameters; the newly paired values may not represent draws from the joint posterior, so their maxima estimate probabilities under an artificial distribution. (4) The approximate standard error is \(\sqrt{0.40(0.60)/10000}=\sqrt{0.000024}\approx0.00490\). It measures Monte Carlo error in estimating the posterior best-arm probability under independent posterior draws, not uncertainty about model adequacy or a frequentist sampling error for the trial estimate.