The Winner's Curse: Why it Happens and What to Do About It

Introduction

Imagine we’ve run an A/B test to figure out which subject line is the best at making people open emails. This test involves ten different subject lines, each of which is sent to 1,000 people. The results are shown below.

Of course, if we actually observed results in such a clear pattern we would be deeply suspicious. But you can imagine we’ve simply labeled the variants in ascending order of observed performance; the perfect spacing is irrelevant to the point of this post; all that matters is there is some spread in observed performance.

When we plug these results into my A/B testing calculator, we see that most of these comparisons are statistically significant. For example, although there is no stat sig difference between variants I and J, or between H and J, there is between G and J, and between F and J, and so forth. Certainly J is one of the best variants, if not the best. Similarly, A is clearly one of the worst variants, if not the worst.

If we were to declare variant J “the winner”, and use it for all emails going forward, what kind of open rate might we expect? The maximum likelihood estimate is simply the observed open rate: 19%. Is that the best prediction of future performance? Well, I can virtually guarantee the performance will be worse than 19%, a phenomenon known as the Winner’s Curse. In this post, we discuss why the Maximum Likelihood Estimate actually does not work very well as a point estimate when we have several experiment groups.

Regression to the Mean

The Winner’s Curse is also known as “regression to the mean”, a phenomenon originally studied by Sir Francis Galton in the 1800s in the context of human stature. Galton observed that tall parents tended to have shorter children—not short, but shorter than the parents. And short parents tended to have taller children—not tall, but taller than the parents. He coined the term “regression to the mean” to describe this phenomenon.

And it’s a good thing this phenomenon exists! Imagine it were the opposite, that tall parents had taller children, and short parents had shorter children. The human race would bifurcate in a few generations, which clearly does not happen. Regression to the mean is a stabilizing force in evolution.

But why does it happen? It’s because height is a combination of genes and random chance, and these things balance out. But it’s easier to understand in a more familiar context: test scores.

Suppose a large group of students take a test. Maybe it’s an intro to statistics course with 1,000 students who are taking the midterm exam. Some of the students do very well, some of the students do pretty poorly, and many students are somewhere in the middle. We might decide the students who do very well are promising candidates for data science careers, and so we monitor their performance on the final. While some of the students still do well, most of the students do worse on the final than they had on the midterm. Meanwhile, many of the students who had done poorly on the midterm, do better on the final. The students we originally thought promising do worse, and the seemingly hopeless students do better: regression to the mean.

To understand this phenomenon, we need only consider why the students who did well on the midterm did so well. Certainly most of them have a pretty good understanding of the material, but many of them also just got lucky: they guessed correctly on the multiple choice questions, and coincidentally the problems asked were very similar to the problems they had studied. The students who did the best had a strong understanding of the material and got lucky. On the final, these students still had a strong understanding of the material, but were not, on average, as lucky as before. Most people don’t get lucky twice in a row. So the overall performance appeared to decrease.

Similarly, the students who did the worst on the midterm struggled with the material, but were also unlucky. They may have studied just as much as the other students, but happened to study problems that didn’t appear in the midterm. But it is just as unlikely to be unlucky twice in a row as it is to be lucky, so on average the performance improves.

We can isolate this phenomenon as: “if you find a group of lucky people, they will tend to be less lucky the next time”. And of course the statement works swapping lucky with unlucky. In contrast, if you start with just a group of people, with no preference towards the lucky or unlucky, we wouldn’t expect to see any particular shift in luck over time. It’s the initial selection of lucky people that creates this phenomenon, and whenever we grab the best (or worst) of a bunch, where “best” is at least partially influenced by luck, we are setting ourselves up for regression to the mean.

This shows up all over, and it goes by many names: regression to the mean, winner’s curse, sophomore slump—it’s the same phenomenon. Imagine there is some promising NBA player in his first year. Everyone thinks he’s the next Lebron James. But then the next season he doesn’t do as well. He still does well, but not as well as people expected, so everyone is disappointed. A musician comes out with an amazing album and becomes really popular, but then their second album is a flop. You find a new restaurant you love, but then the food isn’t as good the next time you go. In all of these cases, we’re selecting someone or something as being especially good, and in so doing we are condemning ourselves to disappointment.

And this phenomenon is exacerbated the more things are involved. How many new musicians release albums per year? For them to stand out, it isn’t enough for them to be excellent musicians, they also need to have the right sound at the right time. That’s luck! And the more musicians there are, the bigger the difference between the luckiest and unluckiest. The more monkeys typing randomly, the more likely one is to type Shakespeare. Mathematically, if the “luck” exhibited by person $i$ is $u_i \sim F$, then the luck associated with the luckiest person, $\xi_n = \max\{ u_1, \ldots, u_n \}$ will tend to increase with $n$.

The larger $n,$ the larger a role luck plays in the top (or bottom) performers. In an A/B test, the more experiment groups we have, the more likely we are to observe the Winner’s Curse.

Empirical Bayes and Hierarchical Models

Going back to our original example, it’s clear that while subject line J is probably pretty good, we should expect some kind of regression to the mean. When we communicate to the business what kind of performance they can expect, we should temper expectations. But how much?

The Empirical Bayes methodology is a simple and appealing approach. We postulate that the “true” open rates (which we do not observe) come from some distribution, such as a Beta distribution, with unknown shape parameters: $\pi^{(k)} \sim \mathrm{Beta}(\alpha, \beta)$. Then the observed opens in each group has a binomial distribution, conditional on the open rate: $y^{(k)} \; | \; \pi^{(k)} \sim \mathrm{Binom}(m^{(k)}, \pi^{(k)})$, where $m^{(k)}$ is the number of emails sent using subject line $i$, and $y^{(k)}$ is the corresponding number of opens. The observed open rate is $y^{(k)} / m^{(k)}$, which is probably close to but distinct from the true open rate $\pi^{(k)}$.

The Empirical Bayes method starts by estimating $\alpha$ and $\beta$ from the data, and then using the posterior distribution $\pi^{(k)} \; | \; y^{(k)}$ to predict the open rate going forward. It can be shown that the open rate so reported will be a weighted average of the observed performance for that subject line and the overall average across all subject lines. As a result, we will tend to shrink both the worst performances and the best performances towards the average, exactly agreeing with the regression to the mean we anticipate.

The “Empirical” part is that we are estimating the parameters of the prior distribution, $\alpha$ and $\beta$ directly from the data. Bayesian folks don’t like this terminology because it makes it sound like other Bayesian approaches aren’t empirical. But I’d rather not get hung up on semantics. It’s a neat approach!

We can do even better. The subject lines might fall into categories, in which case it makes sense to shrink estimates based only on other subject lines in the same category, or giving more weight to those subject lines. Or more generally if we have metadata on the different subject lines, model the prior distribution as depending on these features. In this approach, we can even predict how untested subject lines will perform! See Bayesian Data Analysis by Gelman et al for details.

But What About Maximum Likelihood Estimates?

The narrative above is sensible enough, but it was shocking in the 1950s when Charles Stein first published an argument like the above. At the time, the statistics community was enamored with approaches based on maximum likelihood estimates (MLEs). In our subject line example, the MLEs are simply the observed open rates. MLEs are often simple and intuitive, and have many useful properties like consistency and asymptotic normality. But Stein constructed an estimator, in the context of linear regression, that was uniformly better than MLE in typical situations. From a modern perspective, we would consider this an early form of regularization, and we can (loosely) summarize Stein’s conclusion as: regularization almost always improves point estimates.

Since then, regularization (and Bayesian approaches like Empirical Bayes and Hierarchical Modeling that accomplish the same result) has been embraced in statistics (e.g. ridge regression, the lasso) and especially in Machine Learning, with techniques such as drop out and early stopping in neural networks.

Okay, so regularization almost always improves point estimates. What about hypothesis tests? What about confidence intervals? I think here the answer is more complicated. In the context of A/B tests where the sample size is much larger than the number of parameters in our models, I would hesitate to change my approach for assessing statistical significance or reporting confidence intervals.

With regards to hypothesis testing, we need to ask if a regularization-based test statistic has better power than those we would normally use. Under the typical assumptions of “classical” statistics, e.g. the number of parameters is fixed as the sample size increases, then we ask: do we have simple null and alternative hypotheses? If so, then the likelihood ratio test is most powerful. If the alternative hypothesis is compound, then the score test is most powerful against local alternatives. If the null hypothesis is compound, then the maximum likelihood ratio statistic is asymptotically sufficient and thus most powerful, and is also asymptotically equivalent to the score and Wald tests. These are classic results in statistical inference, and regularization plays no role here. But in a high dimensional setting where we are looking for a few promising candidates in a sea of null effects, a lasso-based approach will likely work very well. That’s usually not the situation in A/B testing.

With regards to confidence intervals, it’s a similar story. A good confidence interval is characterized by having the correct coverage above all else. Confidence intervals based on inverting hypothesis tests automatically have the correct coverage (assuming the hypothesis test does). So I always make my confidence interval based on my hypothesis test, not based on my point estimate. (I certainly do not compute a confidence interval as the point estimate plus or minus two standard errors.)

I do worry that under certain circumstances we might end up adjusting the point estimate so that it is nowhere near the center of the confidence interval. It conceivably could even be outside of the interval! I feel like if this happened in practice, it’s an indication we have data issues, or that our model assumptions (e.g. that all open rates are drawn from the same Beta distribution) are wrong. I just can’t imagine this happening under normal circumstances.

As we discussed above, this regression to the mean phenomenon is exacerbated the more experiment groups we have. I probably wouldn’t bother with two or three experiment groups, but I would definitely consider it for five or more. It certainly does no harm to calculate both sets of point estimates and see how different they are!

Estimating the Prior Parameters

Added February 20, 2022.

I haven’t really been able to find a detailed explanation of how to estimate the prior parameters $\alpha$ and $\beta$ from the data. I’ve seen some hand-wavy explanations. I’ve seen instances where the Binomial distribution is approximated as Gaussian. I’ll sketch out here what I know so far.

Recall: $$ \begin{align} y^{(k)} \; | \; \pi^{(k)} &\sim \mathrm{Binom}\left(m^{(k)}, \pi^{(k)}\right) \\\ \pi^{(k)} &\sim \mathrm{Beta}(\alpha, \beta), \quad k=1, \ldots, K \end{align} $$

As derived in Theory of Point Estimation by Erich Leo Lehmann and George Casella, section 4.6, the marginal distribution of $\alpha$ and $\beta$ is: $$ L(\alpha, \beta) = \Pi_{k=1}^K {m^{(k)} \choose y^{(k)}} \frac{\Gamma(\alpha + \beta) \Gamma(\alpha + y^{(k)}) \Gamma(\beta + m^{(k)} - y^{(k)})}{\Gamma(\alpha) \Gamma(\beta) \Gamma(\alpha + \beta + m^{(k)})}. $$ Taking the logarithm and omitting terms which do not depend on $\alpha$ or $\beta$, the log-likelihood is: $$ \begin{align} \ell(\alpha, \beta) &= \sum_{k=1}^K \left[ \log \Gamma(\alpha + \beta) + \log \Gamma(\alpha + y^{(k)}) + \log \Gamma(\beta + m^{(k)} - y^{(k)}) \right. \\\ &\phantom{=} \hspace{20pt} \left. - \log \Gamma(\alpha) - \log \Gamma(\beta) - \log \Gamma(\alpha + \beta + m^{(k)}) \right] \\\ &= K \left[ \log \Gamma(\alpha + \beta) - \log \Gamma(\alpha) - \log \Gamma(\beta) \right] \\\ &\phantom{=} + \sum_{k=1}^K \left[ \log \Gamma(\alpha + y^{(k)}) + \log \Gamma(\beta + m^{(k)} - y^{(k)} ) \right. \\\ &\phantom{=} \hspace{30pt} \left. - \log \Gamma(\alpha + \beta + m^{(k)}) \right]. \end{align} $$

Now, according to Convex Optimization by Stephen Boyd and Lieven Vandenberghe, section 3.5.1, $\Gamma(x)$ is log-convex for $x \geq 1$. That doesn’t seem to help us much, because we have differences of convex functions that we wish to maximize. This point seems to be completely ignored both by Lehmann and Casella, and in Empirical Bayes Methods by J.S. Maritz and T. Lwin.

Using the relationship between Beta and Gamma functions, we can re-write the log-likelihood as $$ \ell(\alpha, \beta) = -K \log B(\alpha, \beta) + \sum_{k=1}^K \log \left( B(\alpha + y^{(k)}, \beta + m^{(k)} - y^{(k)}) \right). $$ Although this isn’t in Convex Optimization, a little Googling suggests that the Beta function is also log-convex, so the first term of this is concave, but the other terms are convex. So this doesn’t help much either.

Given a good initial guess of $\alpha$ and $\beta$, perhaps a local maximum would work well. That’s what the texts above do, without even acknowledging the issue. It is a little frustrating, though, when the function you’re trying to maximize is not concave.

We may be able to use a convex-concave procedure to improve on a simple local search. We write the problem as: $$ \begin{array}{cl} \text{minimize} & f(\alpha, \beta) - g(\alpha, \beta), \quad (1) \end{array} $$ where $$ \begin{align} f(\alpha, \beta) &= K \log B(\alpha, \beta), \\\ g(\alpha, \beta) &= \sum_{k=1}^K \log \left( B(\alpha + y^{(k)}, \beta + m^{(k)} - y^{(k)}) \right) \end{align} $$ are both convex. We form a linear approximation to $g$ and use that to approximately solve $(1)$. Given estimates $\alpha_0$ and $\beta_0$, we approximate $$ g(\alpha, \beta) \approx \hat{g}(\alpha, \beta) = g(\alpha_0, \beta_0) + \nabla g(\alpha_0, \beta_0)^T \left[ \begin{array}{c} \alpha - \alpha_0 \\\ \beta - \beta_0 \end{array} \right], $$ and then solve $$ \begin{array}{cl} \text{minimize} & f(\alpha, \beta) - \hat{g}(\alpha, \beta). \end{array} $$ The solution to this problem is used to create a new approximation to $g$, and we continue this process until convergence.

This method requires a good initial guess of $\alpha$ and $\beta$ to work well. We could use a Gaussian approximation to the Binomial distribution and fit an Empirical Bayes model. Based on the estimated (Gaussian) prior, we find a Beta distribution that is a good fit, either matching the mean and variance, or by sampling data from the Gaussian distribution and finding a Maximum Likelihood Beta distribution.

Is this better than local search? I’d like to simulate some data and calculate the Bayes error rate associated with (a) the true $\alpha$ and $\beta$ (which is the best the error rate could be), (b) the Gaussian approximation, (c) local search, and (d) the convex-concave procedure outlined above. It may be the case that although MLE for the Beta-Binomial distribution is not a convex optimization problem, and is therefore subject to non-global local optima, the local search still works just as well.

Global Optimization of the Prior Parameters

Added October 14, 2026.

Above, I worried that the marginal log-likelihood is not concave, so local search might stall at a local maximum. I proposed a simulation comparing (a) the true $\alpha$ and $\beta$, (b) a Gaussian approximation, (c) local search, and (d) the convex-concave procedure. Here is that simulation, with a fifth contender: (e) a method that finds the global maximum and proves it has done so. The punchline is that local maxima exist but are rare, and when local search lands on one, the resulting estimates are no worse. Lehmann and Casella were right not to worry. What does matter is making sure local search converges.

The non-concavity lives in one dimension

A little algebra reveals more structure than the difference-of-convex decomposition above. Write $x^{\overline{n}} = x (x + 1) \cdots (x + n - 1) = \Gamma(x + n) / \Gamma(x)$ for the rising factorial. Then the log-likelihood becomes $$ \ell(\alpha, \beta) = \sum_{k=1}^K \left[ \log \alpha^{\overline{y^{(k)}}} + \log \beta^{\overline{m^{(k)} - y^{(k)}}} - \log (\alpha + \beta)^{\overline{m^{(k)}}} \right]. $$

Each term $\log x^{\overline{n}} = \sum_{i=0}^{n-1} \log(x + i)$ is concave in $x$. The negative log-likelihood therefore splits as $$ -\ell(\alpha, \beta) = A(\alpha) + B(\beta) + C(\alpha + \beta), $$

where $A(\alpha) = -\sum_k \log \alpha^{\overline{y^{(k)}}}$ and $B(\beta) = -\sum_k \log \beta^{\overline{m^{(k)} - y^{(k)}}}$ are convex, and $C(s) = \sum_k \log s^{\overline{m^{(k)}}}$ is concave. Without $C$, minimizing $-\ell$ would be a convex problem. And $C$ depends on $\alpha$ and $\beta$ only through their sum, $s = \alpha + \beta$. The non-convexity lives in a single dimension.

Bounding the optimum with weak duality

Write the troublesome part as a supremum over a dummy variable, then swap the infimum and the supremum. Treat $s$ as a separate variable, linked to $\alpha$ and $\beta$ by the constraint $s = \alpha + \beta$. Since $\sup_\nu \nu (\alpha + \beta - s)$ is zero when the constraint holds and $+\infty$ otherwise, minimizing $-\ell$ is the same as calculating $$ p^\star = \inf_{\alpha, \beta, s} \sup_\nu \left[ A(\alpha) + B(\beta) + C(s) + \nu (\alpha + \beta - s) \right]. $$

Swapping the infimum and supremum can only decrease the result (weak duality), so $p^\star \geq \sup_\nu g(\nu)$, where $$ g(\nu) = \inf_\alpha \left[ A(\alpha) + \nu \alpha \right] + \inf_\beta \left[ B(\beta) + \nu \beta \right] + \inf_s \left[ C(s) - \nu s \right]. $$

The dual function splits into three one-dimensional problems. The first two are convex and easy. The third is $-\infty$ for every $\nu$, because $C(s) \to -\infty$ as $s \to 0$. Over the full range of $s$, the bound is useless.

It becomes useful when we restrict $s$ to an interval $[s_\mathrm{lo}, s_\mathrm{hi}]$. A concave function minus a linear function attains its minimum over an interval at an endpoint, so the third term becomes the smaller of two numbers. In fact, the resulting bound equals the optimal value of a convex problem: the original problem with $C$ replaced by its chord. (The chord is the line segment connecting $C(s_\mathrm{lo})$ and $C(s_\mathrm{hi})$.) The chord lies below the concave function $C$, so the bound never exceeds the true minimum over the interval. It is also tight when the interval is narrow: the gap is at most

$$ \frac{(s_\mathrm{hi} - s_\mathrm{lo})^2}{8} \sum_{k=1}^K \sum_{l=0}^{m^{(k)} - 1} \frac{1}{(s_\mathrm{lo} + l)^2}, $$

the standard error bound for linear interpolation, since $|C^{\prime\prime}(s)| = \sum_k \sum_l (s + l)^{-2}$.

Branch and bound

Narrow intervals give tight bounds, which suggests branch and bound (see Boyd’s EE364b notes). Divide the range of $s$ into intervals. For each interval, calculate the lower bound above, and calculate an upper bound by evaluating $-\ell$ at any point in the interval. An interval whose lower bound exceeds the best upper bound found so far cannot contain the global minimum, so we discard it. Split the interval with the smallest lower bound in two, and repeat. Once every remaining interval has a lower bound within some tolerance of the best upper bound, we have found the global optimum, along with a certificate proving it.

Two practical modifications make this work well. First, complete pooling (shrinking every group all the way to the grand mean) corresponds to $\alpha + \beta \to \infty$, a point at infinity. Yet complete pooling is the maximum likelihood estimate for 28% of the simulated data sets below. So we work instead with the prior mean $\mu = \alpha / (\alpha + \beta)$ and $\theta = 1 / (\alpha + \beta)$, where complete pooling corresponds to $\theta = 0$. The same structure survives: the negative log-likelihood becomes

$$ F(\mu, \theta) + \sum_{k=1}^K \sum_{l=0}^{m^{(k)} - 1} \log(1 + l \theta), $$

where

$$ F(\mu, \theta) = -\sum_{k=1}^K \left[ \sum_{i=0}^{y^{(k)} - 1} \log(\mu + i \theta) + \sum_{j=0}^{m^{(k)} - y^{(k)} - 1} \log(1 - \mu + j \theta) \right] $$

is jointly convex (a sum of negative logarithms of affine functions), and the second term is concave in $\theta$ alone. We branch on $\theta \in [0, 1000]$.

Second, we avoid solving a two-dimensional convex problem for every interval. Instead, we minimize $F$ over $\mu$ for fixed $\theta$ (a one-dimensional root-finding problem), giving the profile $h(\theta) = \min_\mu F(\mu, \theta)$. Partial minimization preserves convexity, so $h$ is convex, and its tangent lines at the interval endpoints are valid lower bounds. The envelope theorem gives the derivative $h^\prime(\theta)$ at no extra cost. Each interval then requires only two function evaluations.

In the 24,000 simulated data sets below, branch and bound certified the global optimum (to within $10^{-6}$ in log-likelihood) every time. A brute-force search over a fine grid never found anything better. It took a median of 0.9 seconds per data set in Python (maximum 11 seconds), compared to 2 milliseconds for local search.

Simulation design

We simulated 24 scenarios, 1,000 data sets each, crossing:

• Number of groups: $K = 3$ or $K = 15$.
• Traffic: exactly 2,000 people per group; approximately equal (uniform between 1,800 and 2,200); or unequal (log-uniform between 100 and 20,000).
• Baseline conversion rate: 2% or 50%.
• True rates: either drawn from a Beta prior with $\alpha + \beta = 2000$, or all equal to the baseline except one randomly chosen winner. In the first case the model is correct, and the prior standard deviation equals the binomial standard error for a group of 2,000 people. In the second case the winner converts at 3% instead of 2%, or 54% instead of 50% (3 to 4 standard errors at 2,000 people), and the Beta prior is misspecified.

We compared the following estimates of $\alpha$ and $\beta$, each plugged into the posterior mean $(\alpha + y^{(k)}) / (\alpha + \beta + m^{(k)})$:

• (a) True prior, available only when the true rates come from a Beta distribution.
• (b) Method of moments (MoM), the Gaussian approximation. Estimate the prior mean by the pooled conversion rate, and the prior variance by the variance of the observed rates minus the average binomial variance. Then choose the Beta distribution with that mean and variance. If the variance estimate is negative, pool completely.
• (c) Local search (LS): L-BFGS-B with scipy’s default tolerances over $(\log \alpha, \log \beta)$, started either at the MoM estimate or at $\alpha = \beta = 1$. We also tried the MoM start with much tighter tolerances.
• (d) Convex-concave procedure (CCP), as sketched above, started at the MoM estimate.
• (e) Certified global optimum, from branch and bound.

We also report no pooling (each group’s observed rate) as a reference. We measure accuracy by the mean squared error (MSE) of the estimated conversion rates, averaged over groups and data sets.

How often does local search find the global optimum?

We say a method recovered the global optimum if its log-likelihood is within 0.001 of the certified optimum. Otherwise, one of two things happened. Either it converged to a different local maximum (separated from the global maximum by a valley in the profile likelihood), or it stopped short. A method stops short when it is still climbing toward the global maximum but declares convergence too early. Table 1 summarizes how often each failure occurred.

Table 1. Percent of data sets. Each row pools 4,000 data sets (both baseline rates and both kinds of true rates). The last four columns give the percent of data sets in which each method stopped short.

Two things stand out. First, multiple local maxima arose only when traffic differed sharply between groups, and even then in under 4% of data sets. With equal or approximately equal traffic, we never saw one. (I suspect that with exactly equal group sizes the maximum is always unique, but I have only proven this for tiny groups.) Every method landed on the wrong local maximum in about the same data sets: 1.2% with three groups and unequal traffic, and at most 0.2% otherwise.

Second, local search stopped short far more often than it landed in the wrong place. Starting from MoM with default tolerances, it stopped short in 18% to 54% of data sets. Near complete pooling the profile likelihood is often nearly flat, and complete pooling sits at infinity in the $(\log \alpha, \log \beta)$ parametrization. The gradient becomes tiny long before the optimizer arrives. Tighter tolerances fix most of this.

Does the global optimum give better estimates?

Table 2 measures each method against the frequentist default: taking the observed conversion rates at face value. Four observations stand out.

Table 2. Impact to MSE relative to the observed conversion rates, in percent. Negative values mean a reduction in MSE (better than the observed rates). Each entry averages the impact across the two baseline rates. Monte Carlo 95% intervals for individual entries range from ±0.5 points (15 groups) to ±4 points (3 groups). Every column in a row uses the same data sets, so differences between columns are far more precise: within ±0.8 points, except the MoM column (up to ±1.9).

Empirical Bayes helps more with more groups. With the true rates drawn from a Beta prior and similar traffic, the global optimum cut MSE by 22% with three groups and by 40% with fifteen. With unequal traffic, the reductions were 65% and 82%. More groups mean a more precise estimate of the prior, and so better-calibrated shrinkage. Unequal traffic magnifies the benefit because the small, noisy groups have the most to gain from borrowing strength.

A Beta prior works best when the rates in fact come from a Beta distribution. Empirical Bayes rests on a modeling assumption, and the “winner” scenarios violate it. With three groups, similar traffic, and a single winner, every method increased MSE relative to the observed rates, by 14% (MoM) to 31% (local search from $\alpha = \beta = 1$). With fifteen groups the overall MSE still fell by 45%, but the average hides where the errors went. The global optimum cut MSE for the fourteen groups that did not win by 75%, and roughly quadrupled it for the winner itself. (With three groups, the winner’s MSE rose by about 80%.) Shrinkage cures the winner’s curse for the groups that did not win, and overcorrects for the one that did. Since the winner is usually the group we intend to ship, that is not a comfortable trade.

All methods for estimating the prior lead to roughly the same MSE. Within each row, the global optimum, the CCP, and the three flavors of local search agree to within about 3 points. MoM is the exception. It matched or beat maximum likelihood when traffic was similar, but trailed by 5 to 7 points in three of the four unequal-traffic rows, since it gives small, noisy groups as much weight as large ones. Multiple local maxima, as we saw in Table 1, are rare in these scenarios, and when they do occur they hardly matter.

When local search lands on a different local maximum, the global optimum does not help. Across the 61 data sets where local search converged to the wrong local maximum, the global optimum gave the more accurate estimates only 41% to 45% of the time (depending on the method). The average squared errors were within 4% of each other. In 80% of these data sets, the global maximum was complete pooling, while the competing local maximum shrank each group only partway. The two solutions had nearly equal likelihoods (a median difference of 0.1 in log-likelihood). The data could not tell them apart, and neither was reliably closer to the truth.

Here is one such data set, from the “winner” scenario with three groups and unequal traffic.

The global maximum is complete pooling: every group gets the pooled rate of 50.3%. Local search from the MoM estimate (and the CCP) found a different local maximum at $\alpha + \beta \approx 794$, with a log-likelihood 0.9 lower.

Stopping short cost more than local maxima did. Starting from MoM with default tolerances, local search gave up 1.7 to 2.8 points relative to the global optimum for 15 groups with unequal traffic (78.9% versus 81.7%, for example). Equivalently, its MSE was 11% to 15% higher. Tightening the tolerances closed the gap to 0.1 points. Starting from $\alpha = \beta = 1$ gave up 2.1 to 3.0 points for 15 groups with similar traffic. The CCP was within 1 point of the global optimum in every row, and within 0.1 points in all but the three-group, unequal-traffic rows.

Where we know the true prior, using it gives the largest reduction in MSE. With three groups and similar traffic, the true prior cut MSE by 50%, compared to 22% for the estimated prior; with unequal traffic, 86% compared to 65%. With fifteen groups the estimated prior came tantalizingly close: 40% to 41% versus 50% with similar traffic, and 82% versus 85% with unequal traffic. The price of estimating the prior dwarfs the price of imperfect optimization.

Conclusions

Back to the question I asked in 2022: is global optimization better than local search? In this simulation, no. (We tested one prior strength, with signal and noise comparable at 2,000 people per group; other regimes may differ.) Local search works just as well, provided it converges. Some recommendations:

• Use a parametrization in which complete pooling is a finite point, such as $(\mu, \theta)$ above, or at least compare your answer to the completely pooled fit.
• Tighten the convergence tolerances, or check the gradient at the solution. Stopping short cost far more accuracy than local maxima did.
• If you want a guarantee, branch and bound certifies the global optimum in about a second. The inf-sup trick applies to any problem whose non-convexity enters through a small number of variables.
• Spend your worry on the prior, not the optimizer. With few groups, an estimated prior is far worse than a known one. And when the Beta assumption is wrong, as with a single true winner, shrinkage can hurt the group we care about most.