EconBase
← Back to paper

Empirical Bayes Selection for Value Maximization

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

62,714 characters · 10 sections · 59 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Empirical Bayes Selection for Value Maximization

\authornote{Both authors contributed equally to this research.} \email{[email removed]} \orcid{0000-0002-3040-9662} \authornotemark[1] \email{[email removed]} \orcid{0000-0002-3911-182X}

abstractWe study the problem of selecting the best $m$ units from a set of $n$ as $m / n \to \alpha \in (0, 1)$, where noisy, heteroskedastic measurements of the units' true values are available and the decision-maker wishes to maximize the aggregate true value of the units selected. Given a parametric prior distribution, the empirical Bayes decision rule incurs $\mathcal{O}_p(n^{-1})$ regret relative to the Bayesian oracle that knows the true prior. More generally, if the error in the estimated prior is of order $\mathcal{O}_p(r_n)$, regret is $\mathcal{O}_p(r_n^2)$. In this sense selection of the best units is fundamentally easier than estimation of their values. We show this regret bound is sharp in the parametric case, by giving an example in which it is attained. Using priors calibrated from a dataset of over four thousand internet experiments, we confirm that empirical Bayes methods perform well in detecting the best treatments with only a modest number of experiments.
CCSXML<ccs2012> <concept> <concept_id>10002950.10003648.10003662.10003667</concept_id> <concept_desc>Mathematics of computing Density estimation</concept_desc> <concept_significance>300</concept_significance> </concept> <concept> <concept_id>10002944.10011123.10011131</concept_id> <concept_desc>General and reference Experimentation</concept_desc> <concept_significance>500</concept_significance> </concept> <concept> <concept_id>10002951.10003227.10003241</concept_id> <concept_desc>Information systems Decision support systems</concept_desc> <concept_significance>100</concept_significance> </concept> </ccs2012>

\ccsdesc[300]{Mathematics of computing Density estimation} \ccsdesc[500]{General and reference Experimentation} \ccsdesc[100]{Information systems Decision support systems}

\received{3 February 2025} \received[revised]{23 May 2025} \received[accepted]{23 May 2025}

\ifdefempty{https://doi.org/10.5281/zenodo.15538137}{ \begingroup\raggedrightKDD Availability Link:\\ The source code of the simulations that generate the plots has been made publicly available at \url{https://doi.org/10.5281/zenodo.15538137}. \endgroup }

Introduction

In many important scientific and economic applications, decision-makers are presented with data on the performance of $n$ units, from which they must select a strict subset for further investigation or treatment. Examples include identifying the best teachers, hospitals, or athletes (brown2008season,dimick2010ranking,chetty2014measuringII); genes associated with particular outcomes (efron2002empirical); drug candidates (yu2020bayes); or in the application of this paper, internet experiments. Each unit is associated with an unobserved true value, which is measured with heteroskedastic noise. The constraint that only $m < n$ units can be selected arises naturally when the decision-maker has limited resources to devote to the chosen units, and must restrict attention to the most promising candidates.

A desirable feature of a selection procedure is that the aggregate value of its selections is close to the maximum attainable value. Understanding how different selection procedures perform in this respect enables decision-makers to assess the quality of their decisions in the preceding applications. The empirical Bayes approach to this question involves estimating the unknown, prior distribution from which the true values are drawn, and selecting units with the highest estimated posterior means. We show that if the prior distribution is known to lie within some parametric class, empirical Bayes incurs regret of order $\mathcal{O}_p(n^{-1})$\footnote{$\mathcal{O}_p(\cdot)$ is the stochastic $\mathcal{O}$ notation commonly used in statistics, as defined in van2000asymptotic. We write $X_n = \mathcal{O}_p(a_n)$ if $X_n / a_n$ is bounded in probability.} relative to the oracle Bayes decision rule in which the prior distribution is known.\footnote{We refer to the oracle Bayes decision rule rather than simply the Bayes decision rule throughout, to emphasize that the prior is unknown to the decision-maker.} This is faster than the usual $n^{-1/2}$ parametric rate of convergence from the central limit theorem. In this sense selection is fundamentally easier than estimation: picking a set of units with low regret is easier than pinning down the precise values of those units. This generalizes directly to the nonparametric case: regret converges to zero at the square of the rate that estimation error in the prior converges to zero.

The basic intuition for this result follows. First, mistakes, whether of inclusion or exclusion, are only likely to happen for those units whose true values are sufficiently close to a critical threshold. Units comfortably (or below) above that threshold will be correctly selected (or omitted) with high probability (i.e.\ with probability converging to $1$, abbreviated as w.h.p.). Second, even those mistakes cannot be too costly, as units incorrectly included or excluded are likely to be marginal—almost good enough to be selected, or almost bad enough to be omitted. The regret, which is the product of two terms corresponding to these factors, will therefore be second-order small. We show that our $\mathcal{O}_p(n^{-1})$ bound is sharp in the parametric case, by constructing an example in which regret is at least $C n^{-1}$ with non-vanishing probability for some positive constant $C$.

We illustrate this result with simulations based on internet experimentation data, where units correspond to experiments, and values to treatment effects. Heteroskedasticity arises because experiments vary in sample size. Technology companies may wish to identify a subset of best- or worst-performing experiments for further investigation, in the former case, as candidates to launch to production, or in the latter case, as candidates to stop early. When the follow-up investigation incurs some cost, it may only be feasible to select a strict subset of experiments for more analysis. In this application, as in others, the cost of mistakes depends on their magnitude—that is, on the difference between the aggregate value of the units selected and the units that should have been selected. We simulate true effects from a scale mixture of mean-zero Gaussians calibrated on this dataset, and evaluate the regret of the empirical Bayes approach for selecting the top $10\%$ of experiments. Consistent with our theoretical results, we find that regret is $\mathcal{O}_p(n^{-1})$. By comparison, identifying the set of the top $10\%$ experiments with all misclassifications being equally penalized regardless of their magnitude, or estimating the treatment effects of the selected experiments, or estimating the prior distribution itself, are all structurally harder problems, each of which only exhibits convergence at the usual parametric rate.

Related Work

Our work builds on several large and active strands of the statistics and econometrics literature. Foundational work introducing and developing the empirical Bayes approach to statistics includes robbins1956empirical,kiefer1956consistency,robbins1964empirical,efron1973stein. Applications of the selection problem have proliferated, as the problem of discerning between units which perform well or poorly on the basis of noisy, heteroskedastic measurements describes many real-world settings of interest. Previous work has studied identifying the best teachers (kane2008does,jacob2008can, harris2014skills,chetty2014measuringI,gilraine2020new), the best medical facilities (thomas1994empirical,goldstein1996league,dimick2010ranking,hull2018estimating), the best baseball players (efron1977stein,brown2008season); differentially expressed genes (efron2002empirical,smyth2004linear); promising drug candidates (yu2020bayes); geographic areas associated with the greatest intergenerational mobility (bergman2019creating) or mortality (marshall1991mapping), and employers exhibiting the most evidence of discrimination (kline2021reasonable). Internet experiments are particularly well-suited to empirical Bayes methods (deng2015objective,goldberg2017decision,azevedo2019empirical,coey2019improving,azevedo2020b,guo2020empirical) as datasets are often large enough for accurate estimation of flexibly-specified priors, and the experiment-level sampling error is typically close to normally distributed. For these applications, the aggregate value of the selected units will often be an important component of the decision-maker's utility function. Our results provide theoretical and empirical support for selection based on such methods.

The literature on post-selection inference, including dahiya1974estimation,cohen1989two,gupta2002multiple,fithian2014optimal,hung2019rank,andrews2019inference,guo2021inference, also studies selection problems, but differs from the present work in that its chief focus is estimating the values, differences or ranks of the selected units, rather than analyzing the regret associated with the selection. dahiya1974estimation,cohen1989two provide estimates for the value of a selection unit. gupta2002multiple,fithian2014optimal,hung2019rank,andrews2019inference,guo2021inference largely aim at frequentist inferences. While the notion of regret we consider averages over draws from the distribution of units' true values, an alternative line of inquiry beyond the scope of this paper would be to characterize admissible and minimax decision rules for the frequentist analog of the regret we define, considering the units' values as fixed constants.

gu2020invidious,mogstad2024inference both study similar selection problems to the one we analyze. gu2020invidious take an empirical Bayes approach to selecting the best units while controlling the marginal false discovery rate; mogstad2024inference assert frequentist control over the familywise error rate, which amounts to a zero-one loss based on the correctness of the ranks. Both consider loss functions different from ours. In their frameworks, mistakenly selecting or omitting any unit incurs a discrete cost, whereas in ours the cost of mistakenly selecting or omitting a marginal unit near the selection threshold is small. We view these as complementary perspectives. While in some decision problems mistakes may be undesirable per se, the aggregate performance of the selected units is typically still of interest. For teacher evaluations, for example, policy-makers may rightly be concerned with guarantees over the number of teachers who are incorrectly fired (mogstad2022comment), but may also wish to understand how well their selection procedure is performing from the students' perspective, in terms of aggregate teacher “value-added”. In other contexts, as in internet experimentation or drug discovery, the aggregate value of the selection is the primary concern, and it is harder to justify caring about the number of mistakes per se.

Closely related to our paper is chen2022empirical, which notes the importance of empirical Bayes top-$m$ selection to various social science applications, and derives regret bounds for the problem. Those rate results are more favorable than the ones we present in cases where we can recover the posterior mean but not the prior fast. However in other cases, e.g.\ the parametric case, chen2022empirical bound regret by a term converging slower than $n^{-1/2}$, while we prove $n^{-1}$ convergence and show that rate cannot in general be improved upon. We summarize this comparison in (ref). Furthermore, while the nonparametric rate of convergence of estimated priors to the truth is generally only logarithmic even for optimal procedures carroll1988optimal,fan1991optimal, we may often observe a faster rate of convergence (e.g.\ see (ref)) in practice, which our result translates to a tighter bound on the regret.

table*[table* omitted — 913 chars of source]

The bound in chen2022empirical goes through the mean squared error of the posterior means, which would be more pessimistic if the posterior mean of irrelevant items (e.g. those almost always or never selected) are hard to estimate. Meanwhile, our method relies on controlling the number of mistakes by estimating the distribution, which is likely more challenging than minimizing the mean squared error, leading to potential pessimism in a different way.

Our selection problem may remind readers of the multi-armed bandit literature that studies the problem of identifying the top $m$ arms with a certain probability, e.g.\ shang2020bai,chen2017top. However to target the probability of correct selection is to consider discontinuous loss functions similar to ones in gu2020invidious,mogstad2024inference. We also note that our selection problem is non-sequential which leads to new challenges, as poor choices of parameter in the prior cannot be overcome with additional samples in long run.

Finally, our work is related to the compound decision framework introduced in robbins1951asymptotically, in which a simple decision is made for each unit and the overall loss is the sum of the loss from each individual decision. Convergence and rate results are available for empirical Bayes as applied to compound decision problems, e.g.\ hannan1965rate,van1977empirical,zhang1997empirical,gupta2005empirical,polyanskiy2021sharp, but as weinstein2021permutation observes this framework is rather restrictive and does not encompass the value maximization problem studied here of selecting the best $m$ of $n$ units. weinstein2021permutation generalizes further to a class of simultaneous decision problems that are permutation invariant, which encompasses our problem of selecting $m$ units. However the optimal frequentist solution requires knowledge of the empirical c.d.f.\ (e.c.d.f.), or equivalently the order statistics, of the true effects $\mu_i$. Instead of studying the performance under pathological choices of $\mu_i$ as would be required in a minimax analysis, we take a more Bayesian approach to enable an analysis of regret without knowledge of the order statistics of the $\mu_i$'s.

Our Contribution

Our main contribution relative to this existing literature is to provide the first sharp regret bounds for parametric empirical Bayes selection, and to show how these same ideas extend to control regret in the nonparametric case. Our empirical work on internet experimentation complements this by verifying that regret is quantitatively modest in practice, given a reasonable number experiments of moderate precision. Together, our theoretical and empirical results suggest optimism for empirical Bayes approaches to selection when the decision-maker is primarily concerned with maximizing the aggregate value of the selected units, as opposed to correctly classifying the top units or estimating their values.

The Top-\texorpdfstring{$m$}{m} Selection Problem

Setup

There are $n$ units, each of which is associated with a unobserved true value $\mu_i \in \mathbb{R}$ and an observed noise standard deviation $\sigma_i > 0$.\footnote{We assume $\sigma_i$ to be known, in line with past applications of empirical Bayes methods, e.g.\ weinstein2018group,guo2020empirical,deng2021post.} The $\mu_i$ and $\sigma_i$ are distributed independently from each other and independently across experiments.\footnote{Independence of $\mu_i$ and $\sigma_i$ is a common maintained assumption in empirical Bayes methods, but may be unrealistic in some applications. chen2022empirical treats this topic in detail.} Their unknown marginal distributions are denoted $G_0$ and $H_0$, \[ (\mu_i, \sigma_i) \sim G_0 \times H_0. \] We consider nondegenerate prior distributions $G$ belong to a potentially nonparametric family $\mathcal{M}$, that forms a metric space with $1$-Wasserstein distance $W_1(\cdot, \cdot)$. We assume the family is not misspecified, i.e. the family includes the truth $G_0$. For each unit $i$, the decision-maker observes a measurement $X_i \in \mathbb{R}$, which is distributed as \[ X_i \mid \mu_i, \sigma_i \sim \mathcal{N}(\mu_i, \sigma_i^2). \] The decision-maker must choose $m$ units for some $m < n$. Their average utility given the index set of choices $J \subset \{1, 2, \ldots, n\}$ is $U(J) = \frac{1}{n}\sum_{i = 1}^n \mathds{1}(i \in J) \mu_i$. Let

equation[equation omitted — 221 chars of source]

denote the posterior mean of $\mu_i$ given $X_i,\sigma_i$, assuming that the prior distribution of $\mu_i$ is $G$, where $\phi(\cdot)$ is the probability density function (p.d.f.) of a standard Gaussian. The true posterior mean is $f_{G_0, \sigma_i}(X_i)$. An estimator $\widehat{G} = \widehat{G}(X_1, \ldots, X_n, \allowbreak \sigma_1, \ldots, \sigma_n)$ of $G_0$ is available, where $\widehat{G}$ converges to $G_0$ at some rate $r_n$ in $1$-Wasserstein distance. It is used as the empirical Bayes prior, and in constructing posterior mean estimates, $f_{\widehat{G}, \sigma_i}(X_i)$. For simplicity we denote $f_{G_0, \sigma_i}(X_i)$ and $f_{\widehat{G}, \sigma_i}(X_i)$, the oracle and empirical Bayes posterior means for unit $i$, as $\theta_i$ and $\widehat\theta_i$ respectively. We can view $\theta_i$ as being drawn i.i.d.\ from a distribution, which we denote $P$.

Given the observed data, an oracle Bayesian decision-maker maximizes expected utility (where the expectation is with respect to the posterior distribution over the unknown true values) by selecting the $m$ units with the highest values of $\theta_i$, breaking ties randomly. The empirical Bayes decision-maker mimics this rule, by selecting the $m$ units with the highest values of $\widehat\theta_i$, breaking ties randomly. Letting $J_{\mathrm{EB}}$ and $J_{\mathrm{Bayes}}$ be the empirical Bayes and oracle Bayes choice sets, the regret from empirical relative to oracle Bayes is

align[align omitted — 522 chars of source]

where $\mathds{1}(\cdot)$ is the indicator function.\footnote{We can also consider the loss $U(J) - U(J_{\mathrm{EB}})$ for some other choice of benchmark $J$. However, other natural choices of $J$ may require oracle knowledge of the order statistics of the $\mu_i$'s, e.g.\ when $J$ is optimal among the class of permutation invariant choice sets (weinstein2021permutation).} We aim to characterize how quickly $\mathcal{R}$ converges to zero as $n \to \infty$ and $m / n \to \alpha \in (0, 1)$.\footnote{For simplicity of exposition, we only consider fixed $m$. We expect our main results to extend to the case where $m$ is allowed to be mildly data-driven, with $m/n \to \alpha$ in probability, as when selecting units with positive posterior means.} The proof that $\mathcal{R} = \mathcal{O}_p(r_n^2)$ proceeds by bounding the regret $\mathcal{R}$ by the product of two terms: the proportion of mistakes, and the maximum possible magnitude of the loss caused by a mistake. We show that each of these terms are of the same order as the estimation error in $\widehat{G}$, and consequently regret must be second-order small, i.e.\ $\mathcal{O}_p(r_n) \cdot \mathcal{O}_p(r_n) = \mathcal{O}_p(r_n^2)$.

Note that in the homoskedastic case when all variances are equal, $\sigma_1 = \cdots = \sigma_n$, the posterior mean of $\mu_i$ is monotone in $X_i$ for any choice of prior (efron2011tweedie,koenker2014convex). Hence the oracle Bayes selection rule amounts to selecting the top-$m$ observations ordered by $X_i$, and this selection problem is trivial.

Establishing a Convergence Bound

To establish the convergence bound, we enlist (ref).

assumption$W_1(\widehat{G}, G_0) = \mathcal{O}_p(r_n)$ for some sequence $(r_n)_{n \in \mathbb{N}}$ with $r_n \ge n^{-1/2}$ for all $n$.
assumptionThe support of the distribution $H_0$ of $\sigma_i$ is compact and bounded away from $0$.

(ref) will be satisfied under mild conditions by the maximum likelihood estimator when $\mathcal{M}$ is parameterized by some finite-dimensional parameter $\eta$, with $r_n = n^{-1/2}$ (keener2010theoretical).\footnote{The convergence rate holds for estimating $\eta$, but this translates to common families such as finite mixture families parameterized only by their weights and canonical exponential families with finite variance. For finite mixture families, note that the $1$-Wasserstein metric can be bounded by the product of the $L_1$-norm of the weight parameters and the maximum $1$-Wasserstein metric between any two mixture components. For exponential families, see (ref) in (ref).} This includes the commonly used “normal-normal” model, in which the prior is $\mathcal{N}(m_g, v_g)$ for unknown $m_g, v_g$ which are estimated by maximum likelihood. In the nonparametric case we will generally obtain slower rates of convergence, as allowed for by this assumption. For example, from the existing literature on convergence rates for deconvolution problems:

itemize• if the prior takes on some finite but unknown number of values, chen1995optimal shows that the best possible convergence rate for estimating the prior in the $1$-Wasserstein metric\footnote{The $1$-Wasserstein metric is a natural choice here. We need a statistical distance that reflects the metric on the observation space, as our regret is tied to that metric as well. We also do not require the distributions to have the same support, or more precisely, be absolutely continuous with respect to $G_0$. Other common statistical distances such as total variation distance or Kullback--Leibler divergence do not meet these two desiderata.} is $n^{-1/4}$; • if the prior has a density function with $k$ bounded derivatives, carroll1988optimal shows that the fastest rate of convergence of any estimator of the prior is $(\log n)^{-k/2}$ in the $L_1$-norm of the p.d.f. Furthermore, if we assume the support of $G$ is bounded, the same rate of convergence applies to the $1$-Wasserstein metric.

(ref) states that there are non-trivial upper and lower bounds on the precision with which the true values are measured, as would be the case in experiments with sample sizes bounded below and above.

Under (ref), our main result that $\mathcal{R} = \mathcal{O}_p(r_n^2)$ follows. We give a brief overview of the proof strategy behind this theorem, before establishing supporting lemmas and giving the proof itself. Regret arises because our estimated posterior means, $\widehat\theta_i$, are different from their oracle Bayes counterparts, $\theta_i$, and consequently the top units ranked by the former may differ from the top units ranked by the latter. The difference between $\widehat\theta_i$ and $\theta_i$ is the error relative to oracle Bayes shrinkage for observation $i$. We show that regret can be bounded above by the product of the maximum magnitude of this shrinkage error from mistakes (whether of inclusion or exclusion) and the proportion of such mistakes. It suffices to show that each of these terms is $\mathcal{O}_p(r_n)$. Using the facts that the posterior mean function $f_{G,\sigma}(X_i)$ is sufficiently well-behaved around the $G_0$ when the observations $X_i$ belong to a compact set ((ref)), and that the $X_i$'s associated with all mistakes lie within a compact set w.h.p.\ ((ref)), we can show that the maximum magnitude of the shrinkage error from mistakes is bounded above by a constant times $W_1(\widehat{G}, G_0)$, and hence is $\mathcal{O}_p(r_n)$ by (ref). Next we argue that for any neighborhood around $P^{-1}(1 - \frac{m}{n})$ shrinking slower than $r_n$, the true values associated with mistakes will lie within that neighborhood w.h.p. This allows us to control the proportion of mistakes, and conclude that they are also $\mathcal{O}_p(r_n)$.

The following lemma is a key preliminary result, establishing continuity of both the posterior mean function $f_{G, \sigma}(X)$ and its inverse, and will be used to establish (ref). The existence of the inverse follows immediately from a classic result by efron2011tweedie.

lemmaUnder (ref), the posterior mean function $f_{G, \sigma}(X)$ and its inverse $f^{-1}_{G, \sigma}(X)$ are both continuous in $(G, \sigma, X) \in \mathcal{M} \times \supp(H_0) \times \mathbb{R}$.
proofWe first prove that $f_{G, \sigma}(X)$ is continuous in $(G, \sigma, X)$. From (ref), \begin{align} f_{G, \sigma}(X) &= \frac{\int \mu \phi\left(\frac{X - \mu}{\sigma}\right) \,dG}{\int \phi\left(\frac{X - \mu}{\sigma}\right) \,dG} \coloneqq \frac{h_1(G, \sigma, X)}{h_0(G, \sigma, X)}, \end{align} where $\phi(\cdot)$ is the p.d.f.\ of a standard Gaussian. As $h_0 > 0$, it suffices to show that $h_0$ and $h_1$ are continuous. Suppose we have a sequence $(G_k, \sigma_k, X_k) \to (G^*, \sigma^*, X^*)$ as $k \to \infty$. For $h_1$, we wish to show that \[ \int \mu \phi\left(\frac{X_k - \mu}{\sigma_k}\right) \,dG_k \to \int \mu \phi\left(\frac{X^* - \mu}{\sigma^*}\right) \,dG^*. \] Note that the function sequence $\mu \phi\left(\frac{X_k - \mu}{\sigma_k}\right)$ converges uniformly to $\mu \phi\left(\frac{X^* - \mu}{\sigma^*}\right).$\footnote{For any $\varepsilon > 0$, there exists a compact interval $C_\varepsilon$ such that $\mu \phi\left(\frac{X_k - \mu}{\sigma_k}\right) < \varepsilon$ on $C_\varepsilon^c$. The function sequence itself is equicontinuous and converges pointwise, so it also converges uniformly within $C_\varepsilon$. Hence for any $\varepsilon > 0$ there is sufficiently large $k$ such that $\mu \phi\left(\frac{X_k - \mu}{\sigma_k}\right)$ is within $\varepsilon$ of $\mu \phi\left(\frac{X^* - \mu}{\sigma^*}\right)$ pointwise.} Hence for any $\varepsilon > 0$, when $k$ is sufficiently large, we have \begin{align} & \left| \int \mu \phi\left(\frac{X_k - \mu}{\sigma_k}\right) \,dG_k - \int \mu \phi\left(\frac{X^* - \mu}{\sigma^*}\right) \,dG_k \right| \nonumber \\ \le& \sup_{\mu \in \mathbb{R}} \left| \mu \phi\left(\frac{X_k - \mu}{\sigma_k}\right) - \int \mu \phi\left(\frac{X^* - \mu}{\sigma^*}\right) \right| \nonumber \\ <& \varepsilon / 2. \end{align} Note also that $\mu \phi\left(\frac{X^* - \mu}{\sigma^*}\right)$ is Lipschitz. Therefore by Kantorovich--Rubinstein duality, when $k$ is sufficiently large, $W_1(G_k, G^*)$ is sufficiently small and \begin{equation} \left|\int \mu \phi\left(\frac{X^* - \mu}{\sigma^*}\right) \,dG_k - \int \mu \phi\left(\frac{X^* - \mu}{\sigma^*}\right) \,dG^*\right| < \varepsilon / 2. \end{equation} Summing up (ref) and (ref) yields the convergence of the numerator. The proof for $h_0$ is almost identical, as $\phi\left(\frac{X_k - \mu}{\sigma_k}\right)$ converges uniformly to $\phi\left(\frac{X^* - \mu}{\sigma^*}\right)$ and $\phi\left(\frac{X^* - \mu}{\sigma^*}\right)$ is Lipschitz. The continuity of $f_{G, \sigma}^{-1}(X)$ follows from the continuity of $f_{G, \sigma}(X)$ by (ref).

The next lemma states that posterior mean function $f_{G, \sigma}(X)$ is locally Lipschitz around $G_0$, uniformly in $(\sigma, X) \in \supp(H_0) \times W$, for any compact $W$. This will be used in (ref) to bound the shrinkage error by a constant times the estimation error in the prior parameter, $\widehat{G}$.

lemmaSuppose (ref) holds. Then for any compact $W \subset \mathbb{R}$, there exist positive constants $K$, $\delta$ such that for all $(G, \sigma, X) \in \mathcal{M} \times \supp(H_0) \times W$, we have $|f_{G, \sigma}(X) - f_{G_0, \sigma}(X)| \le K W_1(G, G_0)$ whenever $W_1(G, G_0) < \delta$.
proofUsing the same definition for $h_0$ and $h_1$ as in (ref), we have \begin{align*} & |f_{G, \sigma}(X) - f_{G_0, \sigma}(X)| \\ =& \left| \frac{h_1(G, \sigma, X)}{h_0(G, \sigma, X)} - \frac{h_1(G_0, \sigma, X)}{h_0(G_0, \sigma, X)} \right| \\ \le& \frac{\left| h_1(G, \sigma, X) - h_1(G_0, \sigma, X) \right|}{h_0(G, \sigma, X)} \\ & \quad + \frac{\left| h_1(G_0, \sigma, X) \right| \cdot \left| h_0(G_0, \sigma, X) - h_0(G, \sigma, X) \right|}{h_0(G, \sigma, X) h_0(G_0, \sigma, X)}. \end{align*} It remains to show the following claims for when $G$ is in a sufficiently small neighborhood of $G_0$: \begin{itemize} • $h_0$ is bounded away from $0$: Since $h_0$ is continuous from the proof of (ref), $h_0(G_0, \sigma, X)$ is strictly positive and $\supp(H_0) \times W$ is compact, for sufficiently small $\delta$, $h_0(G, \sigma, X)$ is bounded away from $0$ whenever $W_1(G, G_0) < \delta$. • $h_1$ is Lipschitz in $G$ with a Lipschitz constant that does not depend on $\sigma$ or $X$: The integrand in $h_1$ is $\mu \phi\left(\frac{X - \mu}{\sigma}\right)$, a Lipschitz function in $\mu$. Since this Lipschitz constant is a continuous function in $\sigma, X$ and $\supp(H_0) \times W$ is compact, $\mu \phi\left(\frac{X - \mu}{\sigma}\right)$ is uniformly Lipschitz. The function $h_1$ is Lipschitz in $G$ again by Kantorovich--Rubinstein duality. • $h_1(G_0, \sigma, X)$ is bounded: From the proof of (ref), $h_1$ and thus $h_1(G_0, \cdot, \cdot)$ are bounded. So $h_1(G_0, \sigma, X)$ is bounded since $\supp(H_0) \times W$ is compact. \qedhere \end{itemize}

We will apply (ref) on a specific $W$ which contains all of the observations corresponding to mistakes made by empirical Bayes selection w.h.p. This is the subject of the following lemma. We use $\triangle$ to denote symmetric difference, so $J_{\mathrm{Bayes}} \triangle J_{\mathrm{EB}}$ is the index set of all mistakes.

lemmaIf (ref) hold, there exists a compact set $W$ such that $X_i \in W$ for all $i \in J_{\mathrm{Bayes}} \triangle J_{\mathrm{EB}}$ w.h.p.
proofOracle Bayes selection essentially thresholds on $\theta^*$, the $m$-th largest order statistic of the $\theta_i$'s. For any $c > 0$, this threshold lands in $(P^{-1}(1 - \frac{m}{n}) - c, P^{-1}(1 - \frac{m}{n}) + c)$ w.h.p.\ and hence $(P^{-1}(1 - \alpha)- c, P^{-1}(1 - \alpha) + c)$ w.h.p., by van2000asymptotic. Let $V$ be a open ball of $G_0$ in $1$-Wasserstein whose radius is fixed but to be determined later. By (ref), $\widehat{G}$ lies in $V$ w.h.p. For any $i \notin J_{\mathrm{EB}}$, under the high probability event $A_n$ defined as $A_n = \{\theta^* \in (P^{-1}(1 - \alpha) - c, P^{-1}(1 - \alpha) + c)\} \cap \{\widehat{G} \in V\}$, \begin{align} \widehat\theta_i & \le \min_{i' \in J_{\mathrm{EB}}} \widehat\theta_{i'} \nonumber \\ & = \min_{i' \in J_{\mathrm{EB}}} f_{\widehat{G}, \sigma_{i'}} \circ f^{-1}_{G_0, \sigma_{i'}}(\theta_{i'}) \nonumber \\ & \le \min_{i' \in J_{\mathrm{EB}}} \max_{\sigma' \in \supp(H_0)} f_{\widehat{G}, \sigma'} \circ f^{-1}_{G_0, \sigma'}(\theta_{i'}) \\ & = \max_{\sigma' \in \supp(H_0)} f_{\widehat{G}, \sigma'} \circ f^{-1}_{G_0, \sigma'}\left(\min_{i' \in J_{\mathrm{EB}}}\theta_{i'}\right) \nonumber \\ & \le \max_{\sigma' \in \supp(H_0)} f_{\widehat{G}, \sigma'} \circ f^{-1}_{G_0, \sigma'}\left(\min_{i' \in J_{\mathrm{Bayes}}}\theta_i\right) \\ & \le \max_{\sigma' \in \supp(H_0)} f_{\widehat{G}, \sigma'} \circ f^{-1}_{G_0, \sigma'}(P^{-1}(1 - \alpha) + c) \\ X_i & \le \max_{\sigma' \in \supp(H_0)} f^{-1}_{\widehat{G}, \sigma_i} \circ f_{\widehat{G}, \sigma'} \circ f^{-1}_{G_0, \sigma'}(P^{-1}(1 - \alpha) + c) \\ & \le \max_{\sigma', \sigma\in \supp(H_0)} f^{-1}_{\widehat{G}, \sigma”} \circ f_{\widehat{G}, \sigma'} \circ f^{-1}_{G_0, \sigma'}(P^{-1}(1 - \alpha) + c) \end{align} By (ref), the function $\max_{\sigma' \in \supp(H_0)} f_{\widehat{G}, \sigma'} \circ f^{-1}_{\widehat{G}, \sigma'}$ is a strictly increasing function, so we can move the minimum inside in (ref). Also applying $f^{-1}_{\widehat{G}, \sigma_i}$ to both sides of (ref) yields (ref). Note that the maximand in (ref) is a composition of functions in $\widehat{G}, \sigma', \sigma''$ that are continuous by (ref), hence also a continuous function itself. By maximum theorem and the compactness of $\supp(H_0)$, (ref) is continuous in $\widehat{G}$ and locally bounded. In other words, for a sufficiently small open ball in $1$-Wasserstein, $V$, centered at $G_0$, (ref) is bounded by some constant. On the other hand, for any $i \in J_{\mathrm{Bayes}}$, under the event $A_n$, \begin{align*} X_i &= f^{-1}_{G_0, \sigma_i}(\theta_i) \\ &\ge f^{-1}_{G_0, \sigma_i}(P^{-1}(1 - \alpha) - c) \\ &\ge \min_{\sigma' \in \supp(H_0)} f^{-1}_{G_0, \sigma'}(P^{-1}(1 - \alpha) - c). \end{align*} Together, there exists a constant bounded interval that contains all $i$ in $J_{\mathrm{Bayes}} \setminus J_{\mathrm{EB}}$ under the event $A_n$. Likewise, there is also a constant bounded interval contains all $i$ in $J_{\mathrm{EB}} \setminus J_{\mathrm{Bayes}}$ under the event $A_n$. Taking $W$ to be the union of these two intervals completes the proof.

With these preliminaries we can prove our main result.

restatable{theorem}{rnsquared} If (ref) hold, then $\mathcal{R} = \mathcal{O}_p(r_n^2)$, the square of the rate of convergence for estimating the prior.
proofWe first decompose an upper bound for $\mathcal{R}$ into two components. \begin{align} \mathcal{R} &\le \frac{1}{n}\sum_{i = 1}^n (\mathds{1}(i \in J_{\mathrm{Bayes}}) - \mathds{1}(i \in J_{\mathrm{EB}})) \theta_i \nonumber \\ & \qquad - \frac{1}{n}\sum_{i = 1}^n (\mathds{1}(i \in J_{\mathrm{Bayes}}) - \mathds{1}(i \in J_{\mathrm{EB}})) \widehat\theta_i \\ &= \frac{1}{n}\sum_{i = 1}^n (\mathds{1}(i \in J_{\mathrm{Bayes}}) - \mathds{1}(i \in J_{\mathrm{EB}})) (\theta_i - \widehat\theta_i) \nonumber \\ &\le \frac{1}{n} \left(\#(J_{\mathrm{Bayes}} \setminus J_{\mathrm{EB}}) + \#(J_{\mathrm{EB}} \setminus J_{\mathrm{Bayes}})\right) \nonumber \\ & \qquad \cdot \max_{i \in J_{\mathrm{Bayes}} \triangle J_{\mathrm{EB}}} |\theta_i - \widehat\theta_i| \\ &= 2 \cdot \underbrace{\frac{1}{n} \#(J_{\mathrm{Bayes}} \setminus J_{\mathrm{EB}})}_{proportion of mistakes} \cdot \underbrace{\max_{i \in J_{\mathrm{Bayes}} \triangle J_{\mathrm{EB}}} |\theta_i - \widehat\theta_i|}_{max magnitude of shrinkage error} , \end{align} where (ref) follows from the fact that $J_{\mathrm{EB}}$ is the set of indices of the $m$ largest $\widehat\theta_i$'s. In (ref), since $\#J_{\mathrm{Bayes}} = \#J_{\mathrm{EB}}$, we have $\#(J_{\mathrm{Bayes}} \setminus J_{\mathrm{EB}}) = \#(J_{\mathrm{EB}} \setminus J_{\mathrm{Bayes}})$. From (ref), it suffices to bound the proportion of mistakes and the maximum magnitude of shrinkage error. We start by bounding the latter term. We denote the event that the observations associated with all mistakes belong in the set $W$ from (ref) and $\widehat{G}$ belongs to a neighborhood $V$ around $G_0$ by $A_n = \{X_i \in W \text { for all } i \in J_{\mathrm{Bayes}} \triangle J_{\mathrm{EB}}\} \cap \{\widehat{G} \in V\}$. By (ref), $\mathbb{P}(A_n) \to 1$, and under this high probability event, we have \begin{align} \max_{i \in J_{\mathrm{Bayes}} \triangle J_{\mathrm{EB}}} |\theta_i - \widehat\theta_i| &= \max_{i \in J_{\mathrm{Bayes}} \triangle J_{\mathrm{EB}}} |f_{G_0, \sigma_i}(X_i) - f_{\widehat{G}, \sigma_i}(X_i)| \nonumber \\ &\le \max_{\substack{X \in W \\ \sigma \in \supp(H_0)}} |f_{G_0, \sigma}(X) - f_{\widehat{G}, \sigma}(X)| \nonumber \\ &\le K W_1(G_0, \widehat{G}). \end{align} (ref) follows from (ref), and is $\mathcal{O}_p(r_n)$ by (ref). Consequently \begin{equation} \max_{i \in J_{\mathrm{Bayes}} \triangle J_{\mathrm{EB}}} |\theta_i - \widehat\theta_i| = \mathcal{O}_p(r_n). \end{equation} Next we bound the proportion of mistakes. Let $\theta^*$ denote the $m$-th largest order statistic of the $\theta_i$'s, and $\theta^{**} = P^{-1}(1 - m/n)$ the $(1 - m/n)^\textrm{th}$ quantile of $P$. For any nondecreasing sequence $(b_n)_{n \in \mathbb{N}}$ with $\lim_{n \to \infty} b_n = \infty$, define the event $B_n$ as \[ B_n = \left\{\max_{i \in J_{\mathrm{Bayes}} \setminus J_{\mathrm{EB}}} \left|\theta_i - \theta^{**}\right| \le b_n r_n \right\}. \] Then \begin{align*} & \frac{1}{n} \#(J_{\mathrm{Bayes}} \setminus J_{\mathrm{EB}}) \\ \le& \mathds{1}(B_n^c) + \mathds{1}(B_n) \frac{1}{n}\#\left\{i: \left|\theta_i - \theta^{**}\right| \le b_n r_n\right\}. \end{align*} We first argue that $\mathbb{P}(B_n^c) \rightarrow 0$ and subsequently that $\frac{1}{n}\#\{i: |\theta_i - P^{-1}(1 - \frac{m}{n})| \le b_n r_n\} = \mathcal{O}_p(b_n r_n)$, giving $\frac{1}{n} \#(J_{\mathrm{Bayes}} \setminus J_{\mathrm{EB}}) = \mathcal{O}_p(b_n r_n)$. Since $(b_n)_{n\in\mathbb{N}}$ was an arbitrary nondecreasing sequence converging to infinity, by (ref) in (ref) this implies $\frac{1}{n} \#(J_{\mathrm{Bayes}} \setminus J_{\mathrm{EB}}) = \mathcal{O}_p(r_n)$. For each $i$ in $J_{\mathrm{Bayes}} \setminus J_{\mathrm{EB}}$, empirical Bayes selection must have excluded it because some other shrinkage estimate was larger, i.e.\ $\widehat\theta_i \le \widehat\theta_{i'}$ for some $i' \in J_{\mathrm{EB}} \setminus J_{\mathrm{Bayes}}$. Hence for each $i \in J_{\mathrm{Bayes}} \setminus J_{\mathrm{EB}}$, there is a $i' \in J_{\mathrm{EB}} \setminus J_{\mathrm{Bayes}}$ such that $\theta^* \le \theta_i \le \theta_{i'} + 2 \max_{i \in J_{\mathrm{Bayes}} \triangle J_{\mathrm{EB}}} |\theta_i - \widehat\theta_i| \le \theta^* + 2 \max_{i \in J_{\mathrm{Bayes}} \triangle J_{\mathrm{EB}}} |\theta_i - \widehat\theta_i|$ and so \begin{equation} \max_{i \in J_{\mathrm{Bayes}} \setminus J_{\mathrm{EB}}} | \theta_i - \theta^*| \le 2 \max_{i \in J_{\mathrm{Bayes}} \triangle J_{\mathrm{EB}}} |\theta_i - \widehat\theta_i|. \end{equation} From the triangle inequality and union bound, we have \begin{align*} & \mathbb{P}(B_n^c) \le \mathbb{P}\left(r_n^{-1}\max_{i \in J_{\mathrm{Bayes}} \setminus J_{\mathrm{EB}}} | \theta_i - \theta^*| > b_n/2 \right) + \\ & \quad \mathbb{P}\left(r_n^{-1} \left|\theta^* - \theta^{**}\right| > b_n/2\right). \end{align*} By (ref) and (ref), $r_n^{-1}\max_{i \in J_{\mathrm{Bayes}} \setminus J_{\mathrm{EB}}} | \theta_i - \theta^*| = \mathcal{O}_p(1)$. By (ref), given standard results on the convergence of sample quantiles (van2000asymptotic), we have $r_n^{-1} | \theta^* - P^{-1}(1 - \frac{m}{n})| = \mathcal{O}_p(r_n^{-1} n^{-1/2})$. As $b_n \rightarrow \infty$, $\mathbb{P}(B_n^c) \to 0$. The probability of $\theta_i$ that falls in $(P^{-1}(1 - \frac{m}{n}) - b_n r_n, P^{-1}(1 - \frac{m}{n}) + b_n r_n)$ is no greater than $ P\left(\theta^{**} + b_n r_n\right) - P\left(\theta^{**} - b_n r_n\right) = \mathcal{O}(b_n r_n) $ by the continuous differentiability of $P$ from (ref). So by Chebyshev's inequality, the proportion $\frac{1}{n}\#\{i: |\theta_i - P^{-1}(1 - \frac{m}{n})| < b_n r_n\}$ is $\mathcal{O}_p(b_n r_n)$, and by arbitrariness of $b_n$, it must also be $\mathcal{O}_p(r_n)$. Because the first part and second parts of (ref) are both $\mathcal{O}_p(r_n)$, it follows that the regret $\mathcal{R}$ is $\mathcal{O}_p(r_n^2)$.

The two main estimation approaches for empirical Bayes are {\em $f$-modeling}, in which a model is specified for the observed outcomes, and {\em $g$-modeling}, in which a model is specified for the unobserved prior (efron2014two). This theorem is consistent with either estimation approach. In the $f$-modelling case, if the estimated distribution for outcomes is consistent with some prior distribution for true effects, i.e.\ falls in the class characterized by guo2020empirical, we can think of $\widehat{G}$ as the prior implicitly specified by deconvolving the estimated observation distribution. For $g$-modelling, we can interpret $\widehat{G}$ directly as the model specified for the unobserved prior.

The bound is also sharp when $r_n = n^{-1/2}$, as shown by our example in (ref).

Sharpness of the Convergence Bound in the Parametric Case

We provide an example where the regret satisfies $\mathcal{R} \ge C n^{-1}$ with non-vanishing probability for some positive constant $C$. Let the location family $G(\eta) = \mathcal{N}(\eta, 1)$ be the model for the prior, where the scalar location parameter $\eta$ is estimated by maximum likelihood. In our example, the truth is $\eta_0 = 0$. We assume the standard deviation of the noise term is drawn i.i.d.\ from \[ \sigma_i =

cases1 & with probability $1/2$, and \\ 2 & with probability $1/2$;

\] and we will select $m = \lfloor \alpha n \rfloor$ units.

The maximum likelihood estimator $\widehat\eta$ converges to $\eta_0$ at rate $n^{-1/2}$. The oracle Bayes shrunken estimate of the posterior mean is $\theta_i = \frac{1}{\sigma_i^2 + 1} X_i$ and the empirical Bayes estimate is $\widehat\theta_i = \frac{1}{\sigma_i^2 + 1} X_i + \frac{\sigma_i^2}{\sigma_i^2 + 1} \widehat\eta$. In particular, the magnitude of $\widehat\theta_i - \theta_i = \frac{\sigma_i^2}{\sigma_i^2 + 1} \widehat\eta$ increases with $\sigma_i$. In our setting $\theta_i$ and $\sigma_i$ are measurable with respect to the Lebesgue measure and the counting measure, respectively. The density of $(\theta_i, \sigma_i)$ with respect to the product measure is then $\frac{1}{2} \cdot \sqrt{2} \phi(\sqrt{2}\theta_i)$ for $\sigma_i = 1$ and $\frac{1}{2} \cdot \sqrt{5}\phi(\sqrt{5}\theta_i)$ for $\sigma_i = 2$. If we condition on $\theta^*$ and $\theta_i < \theta^*$, then $\theta_i$ are i.i.d. In fact $(\theta_i, \sigma_i) \mid \theta^*, \theta_i < \theta^*$ is i.i.d.\ with density

equation[equation omitted — 387 chars of source]

with respect to the product measure, where $\Phi(\cdot)$ is the c.d.f.\ of a standard Gaussian. Likewise, $(\theta_i, \sigma_i) \mid \theta^*, \theta_i > \theta^*$ is i.i.d.\ with density

equation[equation omitted — 404 chars of source]

Consider a compact interval $[\underline{a}, \bar{a}]$ that contains $P^{-1}(1 - \alpha)$ in its interior. Since $\theta^*$ converges to $P^{-1}(1 - \alpha)$ at rate $n^{-1/2}$, the event $A_n$ where the interval $(\theta^* - c n^{-1/2}, \theta^* + c n^{-1/2})$ is a subset of $[\underline{a}, \bar{a}]$ happens w.h.p.\ for any positive constant $c > 0$. For any such $\theta^*$, the density in (ref) over $(\theta^* - cn^{-1/2}, \theta^*)$ and density in (ref) over $(\theta^*, \theta^* + cn^{-1/2})$ are in some strictly positive bounded interval $[\underline{b}, \bar{b}]$ that does not depend on the value of $\theta^*$.

There are three sets of units of interest:

itemize$K_n = \{i: \theta_i \in (\theta^*, \theta^* + c n^{-1/2}) \text{ and } \sigma_i = 1\}$, • $L_n = \{i: \theta_i \in (\theta^* - d n^{-1/2}, \theta^*) \text{ and } \sigma_i = 2\}$, and • $M_n = \{i: \theta_i \in (\theta^* - d n^{-1/2}, \theta^* - \frac{1}{2} d n^{-1/2}) \text { and } \sigma_i = 2\}$,

where $c, d$ are positive constants to be chosen. There are $\lfloor \alpha n \rfloor - 1$ realizations of $\theta_i$ greater than $\theta^*$ and $n - \lfloor \alpha n \rfloor$ realizations of $\theta_i$ smaller than $\theta^*$. Conditional on $\theta^*$, the cardinalities are consequently binomially distributed with

align*[align* omitted — 522 chars of source]

Marginalizing over the event $A_n$ gives the same observation but removes the dependence on $\theta^*$. Hence for some constants $c_K, c_L, c_M > 0$, we have w.h.p.\ $\#K_n \ge c_K n^{1/2}$, $\#L_n \ge c_L n^{1/2}$, and $\#M_n \ge c_M n^{1/2}$. Furthermore $d > 0$ can be chosen sufficiently small such that $\# K_n > \# L_n$ w.h.p.

The rest of the argument focuses on the event where $\widehat\eta > \frac{10}{3} (c + d) n^{-1/2}$, which occurs with non-vanishing probability. Under this event, since $\widehat\eta > 0$, empirical Bayes selection will only mistakenly select units with $\sigma_i = 2$ in place of other units with $\sigma_i = 1$. In particular, for any $i$ in $K_n$ and $i'$ in $L_n$, we have $\theta_i > \theta_{i'}$ but \[ \widehat\theta_i < \theta^* + cn^{-1/2} + \frac{1}{2}\widehat\eta < \theta^* - dn^{-1/2} + \frac{4}{5}\widehat\eta < \widehat\theta_{i'}. \] So w.h.p.\ at least $\min(\# K_n, \# L_n) = \#L_n$ mistakes were made. In fact since $L_n$ consists of units immediately smaller than $\theta^*$ and the relative ordering of all units with $\sigma_i = 2$ does not change, all of $\# L_n$ will be mistakenly selected. This incurs a regret of at least \[ \frac{1}{n} \sum_{i \in L_n} (\theta^* - \theta_i) \ge \frac{1}{n} \sum_{i \in M_n} (\theta^* - \theta_i) \ge \frac{1}{n} \# M_n \cdot \frac{1}{2} d n^{-1/2} \ge \frac{1}{2} c_M d n^{-1}, \] with high probability.

Top-\texorpdfstring{$m$}{m} Selection in Simulation

We illustrate (ref) with a realistic simulation, based on the Upworthy dataset of internet experiments conducted between 2013 and 2015.\footnote{From the publicly accessible Upworthy Research Archive (matias2021upworthy) which is downloadable at \url{https://osf.io/jd64p/}.} The dataset contains a list of experiments, along with effect sizes and standard errors. For the prior $G_0$, we fit a normal scale mixture with fixed components, parameterized only by the weights. The data and modelling details are described in (ref), and the notebook to reproduce the simulations and figures is available as an artifact\footnote{Also available on \url{https://github.com/facebookresearch/eb-selection}}.

We simulate a variety of signal-to-ratio regimes, and choices of family for the prior. In increasing order of flexibility, these are:

inlinelist• the family of normal priors, • the family of scale mixtures of normals, and • the family of all distributions.

The priors for these cases are estimated using the {\tt ebnm} R package (willwerscheid2021ebnm). In particular, the normal scale mixture is estimated using adaptive shrinkage as in stephens2016fdr, and the fully nonparametric case is estimated by nonparametric maximum likelihood estimator (NPMLE) (kiefer1956consistency). Henceforth we refer to these three estimators as EB-NN, EB-NSM, and EB-NPMLE. This enables comparison of the performance of empirical Bayes methods under misspecification (when the restrictive EB-NN estimator is used), under a parsimonious and well-specified model (EB-NSM), and under a highly flexibly and well-specified model (EB-NPMLE). For the distribution $H_0$, we use the empirical distribution of standard errors in the dataset.

Top-$m$ selection here corresponds to selecting a subset of experiments, given a constraint on the subset size. We pick $m = \lfloor 0.1 n \rfloor$ and vary $n$, the number of simulated experiments, showing the distribution of regret for each choice of $n$. For each $n$ we run $\num{1000}$ iterations of the selection simulation. In each iteration, we

enumerate• independently draw $n$ true treatment effects $\mu_i \sim G_0$ and noise standard deviations $\sigma_i \sim H_0$; • generate the $n$ observations $X_i$, where $X_i \mid \mu_i, \sigma_i \sim \mathcal{N}(\mu_i, \sigma_i^2)$; • fit three models for the prior distribution of treatment effects $\mu_i$: EB-NN, EB-NSM, and EB-NPMLE; • compute the choice sets $J_{\mathrm{Bayes}}$, $J_{\mathrm{EB-NN}}$, $J_{\mathrm{EB-NSM}}$, $J_{\mathrm{EB-NPMLE}}$, and $J_{\mathrm{UN}}$ corresponding to the oracle Bayes posterior mean estimators, the three empirical Bayes posterior mean estimators, and the unshrunk $X_i$; • compute the regret relative to oracle Bayes selections, $\mathcal{R}_{M} = \frac{1}{n} \sum_{i = 1}^n (\mathds{1}(i \in J_{\mathrm{Bayes}}) - \mathds{1}(i \in J_{\mathrm{M}}))\theta_i$ for $M = \text{EB-NN},\allowbreak \text{EB-NSM},\allowbreak \text{EB-NPMLE},\allowbreak \text{UN}$.

To assess performance in lower signal-to-noise regimes, we repeat this exercise with varying levels of sampling error. We use standard errors 1, 2 and 4 times greater than the baseline standard errors, corresponding to signal-to-noise ratios of roughly $1.3$, $0.7$ and $0.3$.

The normal scale mixture is a parametric model once the number of components and the scale parameters are fixed. As {\tt ebnm} fits by maximizing the likelihood, the remaining parameters—the weights—converge at $\mathcal{O}_p(n^{-1/2})$. Hence by (ref), $\mathcal{R}_\text{EB-NSM}$ is $\mathcal{O}_p(n^{-1})$. We have no such guarantees for $\mathcal{R}_\text{EB-NN}$ or $\mathcal{R}_\text{EB-UN}$, corresponding to the misspecified normal prior and the “naive” choice which selects the units with the largest $X_i$'s. EB-NPMLE is highly flexible and not misspecified, although its guaranteed convergence rate is very slow soloff2024multivariate.

(ref) shows regret as a function of the number of experiments, for each selection method and each value of the noise multiplier. As $n$ increases, the mean and $99$th percentile across simulations of $\mathcal{R}_\text{EB-NSM}$ and $\mathcal{R}_\text{EB-NPMLE}$ both exhibit declines consistent with $n^{-1}$ convergence, although the regret associated with the latter is larger, suggesting the NPMLE model incurs a cost from its greater flexibility. With just $\num{1000}$ experiments, the regret of the EB-NSM approach can be as low as $10^{-4}$ times the standard error of the noise. The normal prior performs better than the unshrunk selection procedure, but neither has regret approaching zero. These patterns are consistent across different noise levels, although the regret are lower with less noise, as the oracle prior and estimated prior are closer.

figure[figure omitted — 581 chars of source]

We compute other quantities of interest from (ref), such as the proportion of mistakes in (ref) and the maximum magnitude of shrinkage error in (ref), as well as the $1$-Wasserstein distance between the true prior and the estimated prior in (ref). As expected, we see that the proportion of mistakes, their magnitude, and the 1-Wasserstein distance between the true and estimated prior in the correctly specified EB-NSM model all converge to zero at $n^{-1/2}$. The misspecified EB-NN model and the unshrunk procedure perform poorly in comparison, with the proportion and magnitude of mistakes not converging to zero, or even increasing, with the number of experiments. The most flexible model, EB-NPMLE, performs worse along every dimension than the more parsimonious EB-NSM, although the proportion and magnitude of its mistakes both converge to zero.

figure[figure omitted — 578 chars of source]
figure[figure omitted — 590 chars of source]
figure[figure omitted — 618 chars of source]

Estimated standard error

The simulations above assume the known standard error to be known, which is reasonable for large-scale online experiments where each experiments have million of units. We complement the simulations above to demonstrate how the noise in the estimated standard error will affect the performance of empirical Bayes methods, showing the regret as the number of units increase and the estimation for standard error improves in (ref).

figure[figure omitted — 435 chars of source]

Conclusion

Our results show that empirical Bayes methods perform well in maximizing the aggregate value of the selected units, in the sense that the regret they incur converges to zero faster than the estimation error in the values themselves. This stands in contrast to prior work emphasizing the difficulty of accurately selecting the best units when the decision-maker incurs a discrete loss from each misclassification (e.g.\ lockwood2002uncertainty,lin2006loss,gu2020invidious). This underscores that rather than selection being an inherently difficult problem, it depends on whether misclassification errors should be weighted by their severity in the utility function. Finally, we note that many extensions and variations on this setting are yet to be fully explored, including characterizing the performance of decision rules for the frequentist analog of the Bayesian regret we study, treating the true values of units as non-stochastic;\footnote{Specifically, studying the utility $U(J_{\mathrm{EB}})$ under permutation invariance as outlined in weinstein2021permutation.} improving performance by incorporating unit-specific covariates into the analysis; and extending to an empirical Bayes knapsack problem where the selected units incur heterogeneous costs.

As discussed in (ref), the frequentist optimal solution requires unavailable oracle knowledge of the order statistics of $\mu_i$'s. This implies that the empirical Bayes solution is not optimal in a frequentist sense, with mainly two gaps:

inlinelist• the optimal solution in weinstein2021permutation is the Bayesian solution with a uniform prior on the permutations of $\mu_i$'s, while empirical Bayes uses $\widehat{G}$ instead; • weinstein2021permutation focused on the loss for a specific set of $\mu_i$'s, while our analysis averages this over $G_0$.

We suspect these gaps are small. For the first gap, weinstein2021permutation conjectures that the Bayesian solution using the e.c.d.f.\ of $\mu_i$'s is asymptotically optimally. We believe this can be reasonably recovered as $\widehat{G}$ when the class of priors $\mathcal{M}$ is sufficiently large. For the second gap, weinstein2021permutation showed that minimizing the loss is equivalently to minimizing the loss averaged over a uniform permutation of $\mu_i$'s. Asymptotically this should be close to the loss averaged over $G_0$, our regret $\mathcal{R}$. Putting this together, both the solution and loss function are similar between the frequentist and the empirical Bayes settings, hinting at some loose frequentist optimality of the empirical Bayes approach.

acksWe thank Eytan Bakshy, Kevin Chen, Matt Goldman, Daniel Jiang, Jelena Markovic, Sepehr Akhavan Masouleh, Jingang Miao, Adam Obeng, Alex Peysakovich, Okke Schrijvers, Yanyi Song, Daniel Ting and Mark Tygert for their comments and suggestions.