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.
105,787 characters · 20 sections · 111 citation commands
Empirical Bayes When Estimation Precision Predicts Parameters
{\onehalfspacing
} \onehalfspacing
Applied economists often use empirical Bayes methods to shrink noisy parameter estimates, in hopes of accounting for the imprecision in the estimates and improving subsequent decisions. Many such settings\footnote{Empirical Bayes methods are applicable whenever many parameters for heterogeneous populations are estimated in tandem. These settings include value-added modeling angrist2017leveraging,mountjoy2021returns,chandra2016health,doyle2017evaluating,hull2018estimating,einav2022producing,abaluck2021mortality, place-based effects chyn2021neighborhoods,finkelstein2021place,chetty2018opportunity,chetty2018impacts,diamond2021standard,baum2019microgeography,aloni2023one, discrimination kline2022systemic,kline2023discrimination,rambachan2021identifying,egan2022harry,arnold2022measuring,montiel2021empirical, meta-analysis azevedo2020b,meager2022aggregating,andrews2019identification,elliott2022detecting,wernerfelt2022estimating,dellavigna2022rcts,abadie2023estimating, and correlated random effects in panel data chamberlain1984panel,arellano2009robust,bonhomme2020much,bonhomme2015grouped,liu2020forecasting,giacomini2023robust,bonhommedenis. } can be described by a heteroskedastic Gaussian sequence model with known variances. That is, researchers obtain statistical estimates $Y_i$ and accompanying standard errors $\sigma_i$ for parameters $\theta_i$ associated with units $i=1,\ldots, n$. Motivated by the central limit theorem, we model $Y_i$ as unbiased Gaussian signals on $\theta_i$ with known variances $\sigma_i^2$: \[Y_i \mid \theta_i, \sigma_i \sim \Norm(\theta_i, \sigma_i^2) \quad i=1,\ldots, n. \addtocounter{equation}{1}\tag{\theequation} \label{eq:gaussian_heteroskedastic_location}\] Loosely speaking, empirical Bayes methods improve decisions---e.g., estimating $\theta_i$ or identifying units with high $\theta_i$---by pooling strength across the many estimates $(Y_i, \sigma_i)_{i=1}^n$ and accounting for differing levels of noise $\sigma_i$.
Commonly used empirical Bayes methods often assume precision independence---that the known standard errors $\sigma_i$ do not predict the underlying parameters $\theta_i$ (i.e., $\sigma_i \indep \theta_i$). However, precision independence is economically questionable and empirically rejected in many contexts. Inappropriately imposing it can harm empirical Bayes decisions, possibly even making them underperform decisions without shrinkage. Motivated by these concerns, this paper introduces and analyzes empirical Bayes methods that allow for precision dependence.
To be concrete, our empirical application bergman2019creating uses empirical Bayes methods to shrink raw economic mobility estimates $(Y_i, \sigma_i)$ of low-income children, curated by chetty2018opportunity. Here, $\theta_i$ represents true unobserved economic mobility of low-income children from Census tract $i$. In this context, precision independence assumes that the standard errors of these estimates do not predict true economic mobility. However, more upwardly mobile Census tracts tend to have noisier estimates, in part because they contain fewer low-income households. Consequently, the standard errors $\sigma_i$ and true mobility $\theta_i$ are positively correlated.
In this context, imposing precision independence can be costly for decision-making. bergman2019creating select high-mobility Census tracts by choosing those with high empirical Bayes posterior means (i.e., shrinkage estimates). Under precision independence, empirical Bayes methods shrink all estimates to their unconditional mean (i.e., $\E[\theta_i]$) and shrink noisier estimates more aggressively. If $\theta_i$ and $\sigma_i$ are positively correlated, such shrinkage tends to systematically underestimate true mobility of high-$\sigma_i$ tracts. This can harm subsequent selection decisions, if we wish to target high-mobility---hence disproportionately high-$\sigma_i$---tracts.\footnote{For a few measures of economic mobility where precision independence is severely violated, we find that screening on conventional estimates selects less economically mobile tracts, on average, than screening on the unshrunk estimates. Fortunately, for the measure of economic mobility (mean income rank pooling over all demographic groups whose parents are at the 25\th percentile of household income) used in bergman2019creating, the violation of precision independence is sufficiently mild, so that screening on these empirical Bayes shrinkage estimates still outperforms screening on the raw estimates. } In contrast, screening on shrinkage estimates computed by our methods selects substantially more mobile tracts.
To introduce empirical Bayes methods, let us return to the Gaussian model (ref). Under this setup, empirical Bayes methods are rationalized as approximations of unknown optimal decisions. Assume that $(\theta_i, \sigma_i)$ are drawn randomly from some distribution. Then the optimal, infeasible decisions take the form of Bayes decision rules for an oracle Bayesian, whose prior is the unknown distribution of $ (\theta_i, \sigma_i)$. Empirical Bayes methods emulate these {oracle decisions} by estimating the {oracle's prior} from the data. For instance, shrinkage estimation, discussed so far, corresponds to using the estimated posterior means of $\theta_i$ given $(Y_i, \sigma_i)$ as a decision rule for predicting $\theta_1,\ldots,\theta_n$. Under this backdrop, precision independence further simplifies the problem of estimating the oracle's prior, but introduces poor performance when it fails to hold.
This paper has two contributions. First, we propose a flexible but tractable framework for modeling precision dependence that nests various proposals in the literature. Our methods are then natural estimation strategies under this framework. (ref) models $\theta_i \mid \sigma_i$ as a conditional location-scale family,\footnote {A location-scale family with shape $G$, indexed by location $m$ and scale $s$, is a set of distributions with cumulative distribution functions (CDFs) $F_ {m,s} (t) = G \pr{\frac{t-m}{s}}$ as $m$ and $s$ vary. For instance, the family $\Norm(m,s^2)$ is location-scale with shape $G(t) = \Phi(t)$, for $\Phi$ the standard Gaussian CDF.} controlled by $\sigma_i$-dependent {location} hyperparameter $m_0 (\sigma) = \E [\theta \mid \sigma]$ and scale hyperparameter $s_0^2(\sigma) = \var(\theta \mid \sigma)$. Under this assumption, different values of $\sigma_i$ translate, compress, or dilate the distribution $\theta_i \mid \sigma_i$, but the underlying {shape} $G_0$ of this distribution is constant over $\sigma_i$. This model subsumes precision independence as the special case where the location and scale parameters are constant functions of $\sigma_i$.
This model naturally gives rise to a family of conditional location-scale empirical Bayes methods---which we call close---by estimating the hyperparameters $(m_0 (\sigma), s_0(\sigma),G_0)$. The close framework also makes estimating these objects highly tractable. The location and scale hyperparameters $m_0 (\cdot), s_0(\cdot)$ can be written as conditional moments of $Y \mid \sigma$, reducing their estimation to learning conditional expectation functions. Subsequently, given $ (m_0(\cdot), s_0(\cdot))$, it is possible to normalize the data $(Y_i, \sigma_i)$ so as to remove precision dependence. After normalization, one could then apply conventional empirical Bayes methods to estimate the remaining hyperparameter $G_0$.
The close framework unifies and generalizes several proposals in the literature kline2023discrimination,weinstein2018group,george2017mortality,ignatiadis2019covariate. These proposals can be viewed as specific modeling and estimation choices for $ (m_0, s_0, G_0)$. Various subsets of these proposals emphasize a nonparametric perspective for modeling and estimating various components of $(m_0, s_0, G_0)$; thus, a natural way to generalize is to adopt a nonparametric perspective for all of them. In particular, we advocate for using nonparametric regression to estimate $(m_0(\cdot), s_0 (\cdot))$ and for using nonparametric maximum likelihood (npmle) to estimate $G_0$ kiefer1956consistency,jiang2009general,koenker2014convex. We refer to this variant as close-npmle. We view close-\textsc{npmle} as a flexible, minimalist, and computationally efficient default, in the absence of substantive knowledge that motivates further restrictions on $ (m_0, s_0, G_0)$.
The second contribution of the paper is a theoretical analysis of close-npmle in (ref). Our main result ((ref)) establishes that, under the close assumptions, close-npmle emulates the oracle Bayesian as well as possible in terms of squared error loss. Specifically, we establish upper and lower bounds for the squared error Bayes regret for \textsc{close}-\textsc{npmle}. These upper and lower bounds match up to logarithmic factors in the number of observations, indicating that \textsc{close}-\textsc{npmle} attains a regret rate that is approximately minimax optimal. These results extend existing regret guarantees for \textsc{npmle}-based empirical Bayes to account for precision dependence soloff2021multivariate,jiang2020general,jiang2009general,saha2020nonparametric. The key technical difficulty is accounting for estimation error in $m_0$ and $s_0$, which feed into \textsc{npmle} estimation.
We enrich our main result in two additional ways. First, to assess robustness of close-npmle to the close assumption, we study a population version of close-npmle under misspecification of the location-scale model. (ref) finds that its worst-case risk---under arbitrarily different shapes of $\theta_i \mid \sigma_i$ as a function of $\sigma_i$---is within a bounded multiple of the risk of a minimax procedure. Second, we also extend our guarantee for squared error regret to the Bayes regret for two ranking-related decision problems, including the problem of selecting high-mobility tracts in bergman2019creating. (ref) shows that the Bayes regret in squared error dominates the Bayes regret for these other decision problems. Coupled with (ref), this implies that close-\textsc{npmle} has good performance for these ranking-related problems as well.
To illustrate our method, (ref) applies close to two empirical exercises chetty2018opportunity,bergman2019creating. The first exercise is a simulation calibrated to the Opportunity Atlas, the dataset published by chetty2018opportunity. For all 15 measures of economic mobility that we consider, close-npmle improves over all alternative methods and captures over 90% of possible mean-squared error (MSE) gains relative to no shrinkage, whereas conventional empirical Bayes methods capture only 70% on average and as little as 50% for some.
The second exercise evaluates the out-of-sample performance of various procedures for selecting high-mobility Census tracts bergman2019creating, using an out-of-sample validation procedure based on the coupled bootstrap that we introduce oliveira2021unbiased. bergman2019creating use empirical Bayes procedures to select high-mobility Census tracts in Seattle. In an exercise that mimics theirs, we find that close-npmle selects more economically mobile tracts than conventional methods. Conventional methods, on the other hand, frequently select less mobile tracts than screening based on the noisy estimates directly. The improvements of close-npmle over the standard method are on median 2.6 times the value of basic empirical Bayes---that is, the improvements the standard method delivers over screening on the raw estimates $Y_i$ directly. Therefore, for this application, if one finds using the standard empirical Bayes method a worthwhile methodological investment, then the additional gain of using close is likewise meaningful.
We observe estimates $Y_i$ and their standard errors $\sigma_i$ for parameters $\theta_i$, over populations $i \in \br{1,\ldots,n}$. We maintain two assumptions that are standard in the empirical Bayes literature gilraine2020new,jiang2020general,soloff2021multivariate,gu2023invidious,guwalters_eb,eb_hole.
First, we assume throughout that the estimates are conditionally Gaussian with known variances equal to $\sigma_i^2$ and are independent across $i$ (ref). The Gaussian model (ref) is heuristically motivated by a central limit theorem applied to the underlying micro-data. This assumption is not without loss: We ignore the fact that the central limit theorem is only an approximation and treat the Normality as exact. \Copy{seremark}{As a concrete example eb_hole, suppose $\theta_i = \E_ {Q_i} [Y_ {ij}]$ is the population mean of some variable $Y_{ij} \sim Q_i$ drawn from population $Q_i$. A natural estimator $Y_i$ of $\theta_i$ is the sample mean of $Y_{i1},\ldots, Y_{in_i}$. A natural estimate for the variance of $Y_i$ is $\sigma_i^2 = n_i^{-2} \sum_ {j=1}^ {n_i} (Y_ {ij}-Y_i)^2$. By standard arguments, as $n_i \to \infty$, $\smash{\sigma_i^{-1}(Y_i - \theta_i) \dto \Norm(0,1)}$. This heuristically motivates (ref) by replacing \smash{“$\dto$”} with “$\sim$.”\footnote{Note too that $Y_i - \theta_i = O_P(n_i^{-1/2})$ and $\sigma_{i} - n_i^{-1/2}\var_ {Q_i} (Y_{ij}) = O_P(n_i^{-1})$, and so the estimation error in $\sigma_i$ is negligible compared to the estimation error in $Y_i$, thereby heuristically justifying treating the estimated standard error $\sigma_i$ as the true variance of $Y_i$.}}
Second, we assume that $ (\theta_i, \sigma_i)$ are random and sampled i.i.d. from some distribution. Since empirical Bayes methods estimate the distribution of $ (\theta_i,\sigma_i)$, it is natural to think of $(\theta_i, \sigma_i)$ as random. For minor technical reasons, throughout, we condition on $\sigma_ {1:n} = (\sigma_1,\ldots, \sigma_n)$ and treat them as fixed. Thus, we think of $\theta_i$ as drawn independently but not necessarily identically: \[ \theta_{i} \mid \sigma_{i} \overset{\mathrm{i.n.i.d.}}{\sim} G_{(i)}. \addtocounter{equation}{1}\tag{\theequation} \label{eq:eb_sampling} \] Let $P_0 \equiv (G_ { (1)},\ldots, G_ {(n)})$ denote the conditional distribution $\theta_ {1:n} \mid \sigma_{1:n}$.
\Copy{covariates}{Throughout, we focus on a setting without additional covariates $X_i$, returning to accommodating for covariates in the empirical application ((ref)). Our methods generalize immediately to settings with covariates $X_i$---as long as $Y_i \mid X_i, \theta_i, \sigma_i \sim \Norm (\theta_i, \sigma_i^2)$---by treating $X_i$ symmetrically as $\sigma_i$. We focus on $\sigma_i$ since it is always present in heteroskedastic empirical Bayes settings, and it enters the likelihood of $Y_i$ unlike other covariates.} Likewise, for simplicity, we focus on a setting where $(Y_i, \theta_i, \sigma_i)$ are independently distributed: We briefly discuss dependence across $i$ in (ref).
Under these assumptions, empirical Bayes methods are desirable for decision-making: They approximate optimal but infeasible decision rules. To see this, consider a decision problem with loss function $L(\bm{\delta}, \theta_{1:n})$, which evaluates an action $\bm{\delta}$ at a vector of parameters $\theta_{1:n}$. The optimal decision---in terms of expected loss $\E_{P_0}[L(\cdot, \theta_{1:n}) \mid \sigma_{1:n}]$ over $(Y_i,\theta_i) \mid \sigma_i$---chooses actions that minimize the posterior expected loss under prior $P_0$: \[ \bm{\delta}^\star(Y_{1:n}, \sigma_{1:n}; P_0) \in \argmin_{\bm{\delta}} \E_{P_0}[L(\bm{\delta}, \theta_{1:n}) \mid Y_{1:n}, \sigma_{1:n}]. \addtocounter{equation}{1}\tag{\theequation} \label{eq:oracle_bayes} \] For this reason, we refer to $\bm{\delta}^\star$ as the oracle Bayes decision rule, and think of it as the Bayes decision rule for an oracle whose prior is $P_0$. $\bm{\delta}^\star$ is infeasible since we do not know $P_0$. To remedy, empirical Bayes methods seek to approximate the oracle Bayes rule $\bm{\delta}^\star$ . Naturally, one recipe is to plug an estimate $\hat P$ for $P_0$ into (ref):\footnote{To emphasize the distinction between the true expectation with respect to the data-generating process (ref) and a posterior mean taken with respect to some possibly estimated measure $\hat P$, we shall use $\E$ to refer to the former and $\mathbf{E}$ to refer to the latter. Subscripts typically make the distinction clear as well. } \[ \bm{\delta}_{\mathrm{EB}} (Y_{1:n}, \sigma_{1:n}; \hat P) \in \argmin_{\bm{\delta}} \mathbf{E}_{\hat P}[L (\bm{\delta}, \theta_{1:n}) \mid Y_{1:n}, \sigma_{1:n}]. \addtocounter{equation}{1}\tag{\theequation} \label{eq:empirical_bayes_rule} \] For the decision problem where $L(\bm{\delta}, \theta_{1:n}) = \frac{1}{n}\sum_ {i=1}^n (\delta_i - \theta_i)^2$ is mean-squared error, (ref) generates empirical Bayes posterior means $\mathbf{E}_{\hat P} [\theta_i \mid Y_i, \sigma_i]$, often referred to as shrinkage estimates james1992estimation,efron1973stein.
To simplify the estimation of $P_0$, popular empirical Bayes methods often assume precision independence: $\theta_i \indep \sigma_i$, or, equivalently, $G_{(1)} = \cdots = G_{(n)}$ in (ref) and equal to some distribution $G_{(0)}$. For instance, the standard parametric empirical Bayes method models $G_{(i)}$ as i.i.d. Gaussian, $G_{ (0)} \sim \Norm(m_0, s_0^2)$ morris1983parametric. State-of-the-art empirical Bayes methods relax the parametric assumptions on $G_{(0)}$ and estimate $G_{(0)}$ with nonparametric maximum likelihood, or npmle jiang2020general,gilraine2020new,soloff2021multivariate. Henceforth, we refer to these methods as independent-gauss and independent-npmle, respectively. The “independent” here emphasizes precision independence.
Despite its convenience, precision independence may be economically implausible; imposing it may cause empirical Bayes methods to underperform. We illustrate this with an application to the Opportunity Atlas chetty2018opportunity. There, one published measure of economic mobility $\theta_i$ of tract $i$ defines it as the probability that a Black individual becomes relatively high-income (i.e., having family income in the top 20 percentiles nationally) after growing up relatively poor in tract $i$ (i.e., with parents at the 25\th percentile nationally).
Intuitively, Census tracts with more low-income Black households should have {more precise} estimates of $\theta_i$, simply because there is a larger sample size to estimate $\theta_i$. However, it is likely that these tracts are also on average poorer and are thus less economically mobile. Thus, these Census tracts should have smaller $\sigma_i$ but also lower $\theta_i$, meaning that $(\sigma_i, \theta_i)$ are positively correlated.
As this economic intuition predicts, precision independence is readily rejected for this measure of economic mobility. (ref) plots the estimates $Y_i$ against their standard errors, overlaying an estimate of the conditional mean function $m_0(\sigma_i) \equiv \E[\theta_i \mid \sigma_i] = \E[Y_i \mid \sigma_i]$. If $\theta_i$ were independent of $\sigma_i$, then the true conditional mean function $m_0(\sigma_i)$ should be constant. (ref) shows the contrary---tracts with more imprecisely estimated $Y_i$ indeed tend to have higher $\theta_i$.
What happens if we apply empirical Bayes methods that assume precision independence here? (ref) overlays empirical Bayes posterior means on the scatterplot. In the top left panel, independent-gauss shrinks $Y_i$ towards a common estimated mean $\hat m_0$, depicted as the black line. When $\sigma_i$ and $\theta_i$ are positively correlated, estimated posterior means under independent-gauss {systematically undershoot} $\theta_i$ for tracts with imprecise estimates. Similarly, the top right panel of (ref) shows that independent-npmle suffers from the same undershooting. In contrast, the bottom panel of (ref) previews our preferred procedure, close-npmle, which shrinks towards the conditional mean $\E[\theta_i \mid \sigma_i]$, thus avoiding the undershooting.
\Copy{msevsranking}{Nonetheless, posterior means from independent-gauss or independent-npmle may still be better predictors, on average, for $\theta_i$ in mean-squared error than the noisy $Y_i$ james1992estimation,efron1975data. However, the undershooting for large $\sigma_i$ is particularly problematic if one hopes instead to select high-mobility Census tracts based on these posterior means, as do bergman2019creating.} On average, high-mobility tracts are exactly those with high $\sigma_i$. Underestimating mobility for these tracts thus leads to suboptimal selections that may even underperform screening directly based on $Y_i$ mehta2019measuring.
To see this, (ref) zooms into a subregion of (ref) and highlights two Census tracts, one in Englewood, NJ, and one in Richmond, CA---referring to them by tracts $A$ and $B$, respectively. Demographically, tract $A$ is 77% nonwhite according to the 2010 Census, and tract $B$ is 57% nonwhite, contributing to different $\sigma_i$'s. Tract $A$ has a lower raw estimate $Y_i$ than tract $B$ ($Y_A < Y_B$); and tracts with similar $\sigma_i$ to tract $A$, on average, also have lower estimates than those similar to tract $B$ (i.e., $\hat m (\sigma_A) < \hat m (\sigma_B)$). Either gap between the two tracts is substantial.\footnote{Both $Y_B-Y_A$ and $\hat m(\sigma_B) - \hat m(\sigma_A)$ are about five percentage points. For reference, an estimate of the unconditional standard deviation of $\theta_i$ is 3.7 percentage points.} These observations are compelling evidence in favor of $\theta_B > \theta_ {A}$: If one would like to select a Census tract to recommend, then, between $A$ and $B$, one is probably better off recommending tract $B$.
However, independent-gauss shrinks both to an estimate of the unconditional mean, which results in a higher posterior mean estimate for tract $A$. In doing so---fooled by an excessively low shrinkage target for tract $B$---independent-gauss recommends tract $A$ over $B$ instead. In contrast, our preferred method (close-npmle) computes posterior means that preserve the more plausible ordering of the two tracts. We do so by modeling the conditional distribution of $\theta_i \mid \sigma_i$ more flexibly, which we turn to now.
We propose the following conditional location-scale model as a relaxation: For a distribution $G_0$ normalized to have zero mean and unit variance, $\theta_i$ has the following representation
(ref) states that the conditional distribution $\theta \mid \sigma$ depends on $\sigma$ via $m_0(\sigma)$ and $s_0(\sigma)$. The function $m_0(\cdot)$ translates the location of the distribution and the function $s_0(\cdot)$ controls the scaling. The underlying shape of the distribution is governed by $\tau_i \sim G_0$ and is restricted by (ref) to be invariant across different $\sigma_i$ values. By the normalization of $G_0$, we can think of $m_0 (\cdot)$ as the conditional mean of $\theta_i \mid \sigma_i$ and $s_0^2(\cdot)$ as the conditional variance.
Applying the empirical Bayes recipe (ref) amounts to estimating the unknown hyperparameters $(\eta_0, G_0)$. Estimating $\eta_0 = (m_0(\cdot), s_0(\cdot))$ is straightforward, as $\eta_0$ can be written as conditional moments of $Y_i \mid \sigma_i$: \[m_0 (\sigma) = \E[\theta \mid \sigma] = \E[Y \mid \sigma] \quad\text{ and } \quad s_0^2 (\sigma) = \var (\theta \mid \sigma) = \var(Y\mid \sigma) - \sigma^2. \addtocounter{equation}{1}\tag{\theequation} \label{eq:m_s_def} \] Estimating $\eta_0$ thus reduces to estimating conditional expectation functions.
Estimating $G_0$ is more complicated. We do so by normalizing away the precision dependence. Consider transforming $(Y_i, \sigma_i)$ into $(Z_i, \nu_i)$, defined by $Z_i \equiv \frac{Y_i - m_0(\sigma_i)}{s_0(\sigma_i)}$ and $\nu_i \equiv \frac{\sigma_i}{s_0(\sigma_i)}$. Note that (ref) implies that \[ Z_i \mid \tau_i, \nu_i^2 \sim \Norm(\tau_i, \nu_i^2) \quad \tau_i \mid \sigma_i, \nu_i \iid G_0. \addtocounter{equation}{1}\tag{\theequation} \label{eq:location_scale_tau_form} \] (ref) makes clear that, first, the transformed triplet $(Z_i, \tau_i, \nu_i)$ obeys an analogue of the Gaussian model (ref), where $Z_i$ is a noisy Gaussian signal on $\tau_i$ with variance $\nu_i^2$. Second, precision independence holds in (ref), since $\tau_i \mid \nu_i \iid G_0$.
This observation motivates the following strategy: First, estimate $m_0$ and $s_0$ with $\hat m(\cdot)$ and $\hat s(\cdot)$ so as to transform $ (Y_i, \sigma_i)$ into $ (\hat Z_i, \hat \nu_i)$: \[ \hat Z_i = \frac{Y_i - \hat m(\sigma_i)}{\hat s(\sigma_i)} \quad \text {and } \quad \hat\nu_i = \frac{\sigma_i} {\hat s(\sigma_i)}. \addtocounter{equation}{1}\tag{\theequation} \label{eq:transformed_data}\] Second, apply empirical Bayes methods that assume precision independence on $(\hat Z_i, \hat\nu_i)$ to estimate $G_0$. This leads to a family of empirical Bayes strategies that we refer to as conditional location-scale empirical Bayes, or close:
This framework produces a family of empirical Bayes strategies, since (ref) can take different forms that practitioners can plug and play. \Copy{covariatesrec}{When there are additional covariates $X_i$ (independent of the noise $ \frac{Y_i - \theta_i}{\sigma_i}$), researchers can choose instead to model $m_0(\sigma_i, X_i)$ and $s_0(\sigma_i, X_i)$ that include these covariates, and estimate $G_0$ after normalizing by $m_0(\sigma_i, X_i)$ and $s_0(\sigma_i, X_i)$.}
This paper focuses on a particular implementation which we call close-npmle. It uses nonparametric regression for (ref) and npmle for (ref). We recommend this method as a flexible default and primarily analyze it in (ref). We conclude this section with several self-contained discussions on implementations of these two steps, the rationale for (ref) and close-npmle, and other miscellaneous issues.
For (ref), one can exploit (ref) by plugging in estimates of conditional expectation functions. For $\hat \E[\cdot \mid \sigma]$ an estimator of conditional means, we may let $\hat m (\sigma) = \hat\E [Y \mid \sigma]$ and $\hat s^2(\sigma) = \hat\E[(Y - \hat m(\sigma))^2 \mid \sigma] - \sigma^2$. The estimator $\hat\E[\cdot \mid \sigma]$ itself may be nonparametric or based on judiciously chosen parametric models eb_hole. \Copy {supportcomment}{The estimation of $\eta_0$ should also impose known support restrictions on $\eta_0$. For instance, the conditional variance estimate $\hat s_0$ should be nonnegative (see (ref)), and the conditional mean estimate should be within the support of $\theta_i$.} Our subsequent theoretical results simply assume that the estimators for $m_0(\cdot), s_0(\cdot)$ are well-behaved and are uniformly accurate.
For (ref), one could again model $G_0$ nonparametrically or parametrically. As a flexible, performant, and minimalist default in the absence of stronger views on the shape $G_0$, we focus on using npmle to estimate $G_0$ koenker2019comment. Formally, the npmle $\hat G_n$ maximizes the log-likelihood of $\hat Z_i$, whose marginal distribution is the convolution $G_0 \star \Norm(0, \hat\nu_i^2)$: For $\varphi(\cdot)$ the Gaussian probability density function and $\mathcal P(\R)$ the set of all distributions supported on $\R$, we maximize \[ \hat G_n \in \argmax_{G \in \mathcal{P}(\R)} \frac{1}{n} \sum_{i=1}^n \log \int_{-\infty}^\infty \varphi\pr{\frac{\hat Z_i - \tau}{\hat \nu_i}} \frac{1}{\hat \nu_i} \, G(d\tau). \addtocounter{equation}{1}\tag{\theequation} \label{eq:npmle} \] In practice, we approximate $\mathcal P(\R)$ with finitely-supported distributions on a grid in order to compute (ref) koenker2014convex.\footnote{koenker2017rebayes provide an efficient software implementation for (ref), which we use throughout.
In terms of grid choice, theoretically, the only downside of a finer grid is computational burden. Ideally, adjacent grid points should have a sufficiently small and economically insignificant gap between them. In our empirical exercises, since the distribution $G_0$ of $\tau_i$ have zero mean and unit variance, we find that a fine grid within $[-6, 6]$ (e.g., 400 equally spaced grid points), with a coarse grid on $[\min_i \hat Z_i, \max_i \hat Z_i] \setminus [-6, 6]$ (e.g., 100 equally spaced grid points), performs well. Our subsequent theory accommodates an approximate maximizer of the likelihood, and thus accommodates the discretization ((ref)). }
\Copy{closegauss}{On the other hand, a default parametric model for $G_0$ is to simply assume that $G_0 \sim \Norm(0,1)$, which we refer to as close-gauss. This approach amounts to using independent-gauss on the transformed estimates $(Z_i, \nu_i)$, with knowledge that the prior $G_0$ has zero mean and unit variance.} Under this model, the oracle Bayes posterior means are: \[ \delta_{\text{\textsc{close}-\textsc{gauss}}}^*(Y_i, \sigma_i) = \frac{\sigma_i^2}{s_0^2(\sigma_i) + \sigma_i^2} m_0 (\sigma_i) + \frac{s_0^2(\sigma_i)}{s_0^2(\sigma_i) + \sigma_i^2} Y_i. \addtocounter{equation}{1}\tag{\theequation} \label{eq:gaussian_cond_b} \] Despite being rationalized under the assumption $\theta_i \mid\sigma_i \sim \Norm(m_0(\sigma_i), s_0^2(\sigma_i))$, this oracle (ref) enjoys strong robustness properties\footnote{(ref) shows that oracle versions of close-npmle satisfy analogous but weaker robustness properties when the location-scale model fails. } even without the location-scale model (ref) and the assumption that $G_0 \sim \Norm(0,1)$. First, (ref) is the optimal linear-in-$Y$ decision rule for estimating $\theta_i$ in squared error weinstein2018group; second, (ref) is minimax in the sense that it minimizes the worst-case mean squared error over choices of $G_{(1)}, \ldots, G_{(n)}$ among all decision rules (see (ref) for formal statements, respectively). This method performs almost as well as \textsc{close}-\textsc{npmle} in our empirical exercises.
We argue that the location-scale assumption provides a unifying framework for a number of existing methods, and close-npmle is a natural generalization of these methods within this framework. We also briefly speculate how to generalize beyond close-npmle.
Several existing methods can be thought of as implementations of close by making different choices in (ref) and (ref). (ref) summarizes how these methods fit into the close framework. Among these methods, some choose nonparametric models for (ref) and some choose nonparametric models for (ref). For instance, weinstein2018group propose close-gauss, with a partition-based nonparametric estimator for $m_0, s_0^2$. kline2023discrimination consider a scale family $\theta_i = s_0(\sigma_i; \beta) \tau_i$ for some $\tau_i \mid \sigma_i \iid G_0$; they model $s_0(\sigma_i; \beta)$ parametrically, but model $G_0$ flexibly using a log-spline sieve efron2016empirical. george2017mortality propose a fully Bayesian model whose components feature parametric choices for $m_0, s_0$ with $G_0 \sim \Norm(0,1)$.
While the right modeling approach likely depends on the particular empirical context, various subsets of these proposals emphasize being flexible in at least one of the two steps. Thus, absent substantive knowledge that motivates more restrictive assumptions, a natural default that unifies these approaches is to be flexible in both steps. Among nonparametric methods, close-npmle may be particularly attractive due to its minimalism: The npmle is free of tuning parameters koenker2019comment, and tuning parameter choices for nonparametric regression are relatively well-understood calonico2019nprobust,armstrong2018optimal. That said, at a high level, when precision dependence is an issue, any approach that models and estimates $m_0, s_0, G_0$ well is likely to perform well.
\Copy{transforms}{While close-npmle naturally generalizes the existing methods in (ref), one might consider methods that do not impose (ref) and are even more flexible. These methods are potentially more theoretically and computationally cumbersome: For instance, we can show that these flexible methods can no longer transform $Y_i$ into some $Z_i = h(Y_i, \sigma_i)$ so as to exploit precision independence on the transformed model $Z_i \mid \tau (\theta_i, \sigma_i), \sigma_i$.\footnote{This is because transforms that preserve linear exponential family structure are necessarily affine. Exponential family structure is important for empirical Bayes because Tweedie's formula holds efron2011tweedie,efron2022exponential. For an affine transform, the only way for $Z_i = a(\sigma_i) + b(\sigma_i) Y_i$ to satisfy precision independence is if (ref) holds. See (ref) for a precise statement.} In this sense, these methods must depart substantially from those that impose precision independence.}
A natural approach is to estimate npmle locally around $\sigma$ values, and we consider these approaches important venues of future work. One might consider discretizing observed $\sigma_i$ values into bins and apply independent-npmle within each bin.\footnote{Our Monte Carlo exercise in (ref) uses a similar approach to construct a Monte Carlo data-generating process. Thus, the oracle performance in the Monte Carlo is the best-case scenario for the performance of this procedure. There, we find close-npmle performs well relative to the oracle and thus to this procedure ((ref)).} A smoother---but more computationally intensive---alternative is to estimate the posterior at some given $\sigma$ by considering only observations with $\sigma_i \in [\sigma-h, \sigma+h]$ and again use independent-npmle for these observations. For these methods, the number of bins and bandwidth $h$ are tuning parameters. While we anticipate ad hoc choices of tuning parameters to perform well, a proper theoretical analysis likely needs to link tuning choices to smoothness in the conditional distribution $\sigma \mapsto f_{Y\mid \sigma}(\cdot \mid \sigma)$ with respect to certain distributional distances. The corresponding regularity conditions thus seem more complex than smoothness conditions for conditional expectations required by close-\textsc{npmle}.
As a review, we observe $(Y_i, \sigma_i)_{i=1}^n$, where $(\theta_i,\sigma_i)$ satisfies (ref) and $(Y_i, \theta_i, \sigma_i)$ obeys (ref). The procedure close-npmle transforms the data $(Y_i, \sigma_i)$ into $(\hat Z_i, \hat \nu_i)$, with estimated conditional moments $\hat\eta = (\hat m, \hat s)$ for $\eta_0 = (m_0, s_0)$ in (ref). It then estimates $G_0$ via npmle (ref) on $(\hat Z_i, \hat \nu_i)_{i=1}^n$. This section introduces a few statistical guarantees on the performance of close-npmle in terms of regret. To unify presentation, we first review decision theory primitives and introduce regret.
Let $\bm{\delta}(Y_{1:n}, \sigma_{1:n})$ be a decision rule mapping the data $(Y_{1:n}, \sigma_{1:n})$ to actions. Recall that $L(\bm{\delta}, \theta_{1:n})$ denotes a loss function mapping actions and parameters to a scalar. Let $R_{\mathrm{B}} (\bm{\delta}; P_0) = \E_{P_0}[L(\bm{\delta}, \theta_{1:n}) \mid \sigma_{1:n}]$ be the Bayes risk of $\bm{\delta}$ under $P_0$. The oracle Bayes decision rule $\bm{\delta}^\star$ (ref) is optimal in the sense that it minimizes $R_{\mathrm{B}}$. Thus, a natural performance measure for the empirical Bayesian (ref) is the gap between the Bayes risks of $\bm{\delta}_ {\mathrm{EB}}$ and $\bm{\delta}^\star$. We refer to this quantity as Bayes regret:
where the right-hand side integrates over the randomness in $\theta_{1:n}, Y_{1:n}$, and, by extension, $\hat P$. If an empirical Bayes method achieves low Bayes regret, then it successfully imitates the decisions of the oracle Bayesian, and its decisions are thus approximately optimal. Our results show that Bayes regret for close-npmle vanishes quickly as a function of $n$.
In what follows, we use the symbol $C$ to denote a generic positive and finite constant which does not depend on $n$. We use the symbol $C_{x}$ to denote a generic positive and finite constant that depends only on $x$, some parameter(s) that describe the problem. Occurrences of the same symbol $C, C_x$ may not refer to the same constants. Since all expectation or probability statements are with respect to the conditional distribution $P_0$ of $\theta_{1:n} \mid \sigma_{1:n}$, going forward, we treat $\sigma_ {1:n}$ as fixed and simply write $\E[\cdot], \P(\cdot)$ to denote the expectation and probability over $\theta_{1:n} \mid \sigma_{1:n} \sim P_0$; we may omit the subscript $P_0$ or the conditioning on $\sigma_{1:n}$.
Our main result concerns the canonical statistical problem of estimating the parameters $\theta_{1:n}$ under MSE.
\Copy{mseregretdef}{ For (ref), define $\mathrm{MSERegret}_n$ as the excess loss of the empirical Bayes posterior means relative to that of the oracle Bayes posterior means:
where $\theta_i^*$ are the oracle posterior means and $\hat\theta_ {i,G,\eta}$ are the posterior means under a prior parametrized by $(G, \eta)$. } The corresponding Bayes regret (ref) for close-npmle in this decision problem is then the $P_0$-expectation of $\mathrm{MSERegret}_n$: \[ \mathrm{BayesRegret}_n = \E\bk{ \mathrm{MSERegret}_n(\hat G_n, \hat\eta)} = \E_{P_0}\bk{ \frac{1}{n} \sum_{i=1}^n (\theta_i^* - \hat\theta_{i, \hat G_n, \hat\eta})^2 \addtocounter{equation}{1}\tag{\theequation} \label{eq:mse_regret_and_mse} }. \] Equation (ref) additionally notes that expected $\mathrm{MSERegret}_n$ is equal to the expected mean-squared difference between the empirical Bayesian posterior means $\hat\theta_{i, \hat G_n, \hat\eta}$ and their oracle counterparts $\theta_i^*$. Our subsequent results ((ref)) state upper and lower bounds for $\mathrm{BayesRegret}_n$, over a class of data generating processes $\mathcal P_0\ni P_0$. We now introduce and discuss the assumptions on $\mathcal P_0$.
We first assume that $\hat G_n$ is an approximate maximizer of the log-likelihood on the transformed data $(\hat Z_i, \hat\nu_i)$ satisfying some support restrictions. This is not restrictive, as the actual maximizers of the log-likelihood function satisfy it (Proposition 4, soloff2021multivariate). This assumption also accommodates for the fact that the npmle is approximated by a discrete distribution on a grid.
We now state further assumptions on $\mathcal P_0$ beyond (ref). First, we assume that $G_0$ is sufficiently thin-tailed such that its moments grow slowly.\footnote{An equivalent statement to (ref) is that there exists $a_1, a_2 > 0$ such that $\P_{G_0}(|\tau| > t) \le a_1\exp\pr{-a_2 t^\alpha}$ for all $t > 0$. Note that when $\alpha = 2$, $G_0$ is subgaussian, and when $\alpha = 1$, $G_0$ is subexponential vershynin2018high. (ref) is slightly stronger than requiring that all moments exist for $G_0$, and weaker than requiring $G_0$ to have a moment-generating function. Similar tail assumptions feature in the theoretical literature on empirical Bayes soloff2021multivariate,jiang2009general,jiang2020general. } \Copy{alpha} {The thickness of its tail is parametrized by $\alpha \in (0,2]$, which subsequently affects the log factors in (ref).}
Next, (ref) imposes that members of $\mathcal P_0$ have various variance parameters uniformly bounded away from zero and $\infty$. This is a standard assumption in the literature, maintained likewise by jiang2020general and soloff2021multivariate.
Lastly, we require that $m_0(\cdot)$ and $s_0(\cdot)$ satisfy some smoothness restrictions. We also require that $\hat m(\cdot)$ and $\hat s(\cdot)$ satisfy some corresponding regularity conditions. Let $C_{A_1}^p([\sigma_\ell,\sigma_u])$ denote the H\"older class of order $p \ge 1$ with maximal H\"older norm $A_1 > 0$ supported on $ [\sigma_\ell,\sigma_u]$ vaart1996weak.
(ref) is a H\"older smoothness assumption on the conditional moments $m_0$ and $s_0$, which is a standard regularity condition for nonparametric regression. Moreover, it is also a high-level assumption on the quality of the estimation procedure for $(\hat m, \hat s)$. It expects that $\hat m$ and $\hat s$ are accurate in $\norm{\cdot}_\infty$, belong to a class with manageable metric entropy, and obey the bounds for $s_{0}$.\footnote{ (ref)(2) is slightly stronger than an estimation rate requirement $\norm{\hat\eta - \eta_0}_\infty = O_P\pr {n^{-p/(2p+1)}(\log n)^ {\beta_0}}$, in the sense that the probability of large deviations are additionally controlled. Local polynomial smoothing estimators can attain the desired estimation rate of $n^{-p/(2p+1)}(\log n)^{\beta_0}$ in $\norm{\cdot}_\infty$ tsybakov2008introduction,stone1980optimal. Since the data is assumed to be thin-tailed in (ref), such estimators also attain the stronger requirement in (ref)(2).
For (ref)(3), if the estimators $\hat m$ and $\hat s$ are $p$-H\"older smooth almost surely, we can simply take $\mathcal V = C_{A_1'}^p([\sigma_\ell,\sigma_u])$ for some potentially different $A_1'$. This can be achieved in practice by, say, projecting estimated parameters $\tilde \eta$ to $C_{A_1}( [\sigma_\ell, \sigma_u])$ in $\norm {\cdot}_\infty$.
Finally, (ref)(4) also expects the conditional moment estimates $\hat\eta$ to respect the boundedness constraints for $s_0$. This is mainly so that our results are easier to state.
We show in (ref) that a local linear regression estimator (with $\hat s$ suitably truncated) satisfies weaker conditions than (ref)(2)--(4) that are nonetheless sufficient for the conclusion of (ref). }
(ref) specify a class of distributions $\mathcal P_0$ and estimators $\hat\eta = (\hat m(\cdot), \hat s(\cdot))$ regulated by a set of hyperparameters $\mathcal{H} = (\sigma_\ell, \sigma_u, s_\ell, s_u, A_0, A_1, \alpha, \beta_0, p).$ Our subsequent theoretical results are uniform over $\mathcal P_0$ for a fixed $\H$.
Our main result is a non-asymptotic upper bound for (ref): The MSE regret of close-npmle converges to zero no slower than $n^{-\frac{2p}{2p+1}}(\log n)^{C}$.
Second, we give a corresponding minimax lower bound on the regret, which shows that (ref) cannot be improved by more than logarithmic factors.
(ref) continues a recent statistics literature on empirical Bayes methods via npmle, by characterizing the effect of an estimated first-step parameter $\hat\eta$. Our theory hews closely to---and extends---the results in jiang2020general and soloff2021multivariate, which themselves extend earlier results in the homoskedastic setting jiang2009general,saha2020nonparametric. In particular, soloff2021multivariate show that the MSE regret rate is of the form $C (\log n)^ {\beta} \frac{1}{n}$ under precision independence and assumptions similar to ours. In this context, we show that first-step estimation error degrades this regret rate gracefully, and we link the corresponding regret rate to the smoothness of $\eta_0$. The proof of (ref) is deferred to the Online Appendix, but its main ideas are outlined in (ref).
(ref) shows that the rate (ref) is optimal up to logarithmic factors. These logarithmic factors partly reflect inefficiencies in the proof of (ref), but in any case the gap is not large. We prove (ref) by showing that any good posterior mean estimate $\hat\theta_i$ implies a good estimate $\hat m(\sigma_i)$ for $m_0$ for some particular choice of $G_0, \sigma_ {1:n}, s_0^2(\cdot)$. Minimax lower bounds for estimation of $m_0$ tsybakov2008introduction then imply lower bounds for estimation of the oracle posterior means $\theta_i^*$ ignatiadis2019covariate.
We additionally note that these regret upper bounds readily extend to the case where covariates are present and the location-scale assumption (ref) is specified with respect to the additional covariates $X_i$: \[\theta_i \mid \sigma_i, X_i \sim G_0\pr{\frac{\cdot - m_0(X_i, \sigma_i)}{s_0(X_i, \sigma_i)}}, \addtocounter{equation}{1}\tag{\theequation} \label{eq:location-scale-covariates} \] under smoothness assumptions on $(m_0, s_0, \hat m, \hat s)$ analogous to (ref). The resulting convergence rate would reflect the dimensionality of the covariates, and the term $n^ {- \frac{2p}{2p+1}}$ would be replaced with $n^ {- \frac{2p}{2p+1+d}}$, where $d$ is the dimension of $X$.
Taken together, (ref) are statistical optimality guarantees for close-npmle in terms of (ref). That is, the worst-case MSE performance gap of close-npmle relative to the oracle contracts at the best possible rate, meaning that close-npmle mimics the oracle as well as possible.
We prove (ref) imposing the location-scale model (ref). This is an optimistic assessment of the performance of close-npmle. While (ref) nests precision independence, it may still be misspecified. This subsection explores the worst-case behavior of close-npmle without (ref).
We do so by considering an idealized version of close-npmle. So long as $\theta_i \mid \sigma_i$ has two moments, $\eta_0(\cdot) = (m_0(\cdot), s_0(\cdot))$ are well-defined as conditional moments. We will assume that $m_0, s_0$ are known. Without (ref), $G_0$ is ill-defined, but we assume that we obtain some pseudo-true value $G_0^*$ that has zero mean and unit variance. Thus, for estimating $\tau_i = \frac{\theta_i - m_0(\sigma_i)}{s_0(\sigma_i)}$, whose distribution is $\tau_i \mid \sigma_i \sim G_i$, this idealized procedure uses some misspecified prior $G_0^* \neq G_i$, where $G_0^*$ agrees with $G_i$ in the first two moments. The worst-case performance of the procedure that uses $G_0^*$ depends on how far posterior means under $G_0^*$ differs from posterior means under $G_i$.
\Copy{closegaussconstant}{We show in (ref) that this difference is bounded uniformly for all $G_0^*$ satisfying an additional tail assumption. This result implies that the maximum risk of this procedure is at most a constant multiple of the minimax risk; here, the minimaxity is defined with respect to a game between an analyst and an adversary, where the analyst knows $m_0, s_0$ and hopes to estimate $\theta_{1:n}$, and the adversary chooses the shape of the distribution $\tau_i \mid \sigma_i$. In this game, the oracle version of close-gauss (ref) is a minimax procedure ((ref)). }
Specifically, let $\mathcal P(m_0, s_0)$ denote the set of distributions of $\theta_{1:n} \mid \sigma_{1:n}$ where $\E[\theta_i \mid \sigma_i] = m_0(\sigma_i)$ and $\var(\theta_i \mid \sigma_i) = s_0^2(\sigma_i)$. Let \[\mathcal G_0(\lambda, \epsilon) \equiv \br{G_0^*: \E_{G_0^*}[\tau]=0, \var_{G_0^*} (\tau) = 1, G_0^* (-z) \vee (1-G_0^*(z)) \le \lambda z^{-2-\epsilon} \text{ for all $z>0$}} \]be the set of mean-zero, variance-one distributions satisfying an additional tail condition indexed by $\lambda > 0, \epsilon > 0$.\footnote{By Markov's inequality, this condition is satisfied if $G_0^*$ has its $ (2+\epsilon)$\th moment bounded by $\lambda$. A previous version of this paper stated (ref) without this additional tail condition, regrettably due to a technical error that is corrected in this version. See (ref).}
(ref) shows that the worst-case behavior of an idealized version of close-npmle comes within a factor of the minimax risk. Thus, close-npmle is not arbitrarily unreasonable, even under misspecification. We caution that (ref) is a fairly weak guarantee, in that the decision rule that simply outputs the prior conditional mean ($\delta_i = m_0 (\sigma_i)$) also satisfies it. Nevertheless, even so, (ref) {does not} hold for an idealized version of independent-gauss.\footnote{\Copy {indepgaussconstant} {That is, it does not hold for the implementation of independent-gauss that plugs in known unconditional moments $m_0 = \frac{1}{n} \sum_{i=1}^n m_0 (\sigma_i)$ and $s_0^2 = \frac{1}{n} \sum_{i=1}^n (m_0 (\sigma_i) - m_0)^2 + s_0^2 (\sigma_i)$. To wit, take $s_0 (\sigma_i) \approx 0$. Then, the minimax risk as a function of $(s_0 (\cdot), m_0(\cdot))$ is approximately zero, but $m_0(\cdot)$ can be chosen such that the risk of \textsc{independent-gauss} is bounded away from zero. See (ref) for a formal statement.} }
So far, our regret guarantees are only about estimation in MSE ((ref)). We now turn to two decision problems that involve ranking or selection and show similar guarantees for close-npmle in terms of regret for these decision problems. These decision problems are likely more economically relevant for, e.g., replacing low value-added teachers, recommending high-mobility tracts, or treatment choice gilraine2020new,bergman2019creating,manski2004statistical,stoye2009minimax,kitagawa2018should,athey2021policy.
The oracle Bayes decision rules $\bm{\delta}^\star$ in (ref) depend solely on the vector of oracle Bayes posterior means $\theta_ {1:n}^*$. Therefore, for these problems, the natural empirical Bayes decision rules simply replace oracle Bayes posterior means ($\theta_i^*$) with empirical Bayes ones ($\hat\theta_i$). It stands to reason that as $\hat\theta_i$ is close to $\theta_i^*$ in squared error, even when $\hat\theta_i$ implies the wrong selection decision, this decision is not too costly for the empirical Bayesian. We formalize this intuition in the following theorem, showing that if $\hat\theta_i$ are close to $\theta_i^*$ in MSE, then decisions plugging in $\hat\theta_i$ are also close to their oracle counterparts in terms of Bayes risk.
To specialize, let $\mathrm{UMRegret}_n$ denote $\mathrm{BayesRegret}_n$ for the loss function in (ref) and let $\mathrm{TopRegret}_{n}^{(m)}$ denote $\mathrm{BayesRegret}_n$ for (ref).
(ref) shows a sense in which (ref) are easier than (ref): The regret of the latter dominates those of the former. As a result, if we use close-npmle under (ref), our convergence rates from (ref) also upper bound regret rates for these two decision problems. In particular, for $m/n \to c \in (0,1)$, both regret rates (ref) and (ref) are of the form $n^{-p/(2p+1)} (\log n)^{C} = o(1)$ under (ref). Thus, the performance of the empirical Bayes decision rule approximates that of the oracle at least as fast as $O(n^{-p/(2p+1)})$, up to log factors.
We close this section with a procedure that provides unbiased estimates of the loss of arbitrary decision rules for (ref). Practitioners can use this procedure to evaluate the gain of close-npmle relative to other alternatives---we do so extensively in (ref). The validity of this validation depends only on the Gaussianity (ref)---without assuming $(\theta_i, \sigma_i)$ are random nor assuming the location-scale model (ref).
For some $\omega > 0$ and an independent Gaussian noise $W_i \sim \Norm(0,1)$, consider adding to $Y_i$ and subtracting from $Y_i$ some scaled version of $W_i$: \[ Y_{i}^{(1)} = Y_i + \sqrt{\omega} \sigma_i W_i \quad Y_{i}^{(2)} = Y_i - \frac{1}{ \sqrt{\omega}} \sigma_i W_i. \] oliveira2021unbiased call $(Y_{i}^{(1)}, Y_{i}^{(2)})$ the coupled bootstrap draws. Observe that the two draws are conditionally independent under (ref): \[ \colvecb{2}{Y_{i}^{(1)}}{Y_{i}^{(2)}} \mid \theta_i, \sigma_i^2 \sim \Norm\pr{ \colvecb{2}{\theta_{i}}{\theta_i},
}. \addtocounter{equation}{1}\tag{\theequation} \] The conditional independence allows us to use $Y_{i}^{(2)}$ as an out-of-sample validation for decision rules computed based on $Y_i^{(1)}$. We denote their variances by $\sigma_{i, (1)}^2$ and $\sigma_{i, (2)}^2$.
The coupled bootstrap can be thought of as approximating sample-splitting the micro-data without needing access. We could imagine splitting the micro-data into training and testing sets, and think of $Y_i^{(1)}$ as training-set estimates and $Y_{i}^{ (2)}$ as testing-set estimates. We might compute decisions based on $Y_i^ {(1)}$ and evaluate them honestly with fresh data $Y_i^{(2)}$. The coupled bootstrap precisely emulates this sample-splitting procedure.\footnote{To see this, suppose $Y_{i} = \frac{1}{n_i} \sum_{j=1}^ {n_i} Y_ {ij}$ is a sample mean of i.i.d. micro-data $Y_{ij}: j = 1,\ldots, n_i$. Suppose we split $Y_{ij}$ into two sets, with proportions $\frac{1} {\omega + 1}$ and $\frac{\omega}{\omega + 1}$, respectively. Let $Y_i^{(1)}$ and $Y_i^{(2)}$ be the sample means on each respective set. Then the central limit theorem motivates that, approximately, (ref) holds for $Y_i^{ (1)}$ and $Y_i^{(2)}$. For instance, coupled bootstrap with a value of $\omega = 1/9$ is statistically equivalent to splitting the micro-data with a 90-10 train-test split.}
The following proposition formalizes how to use coupled bootstrap to provide unbiased estimators for the loss of a generic decision rule.\footnote{oliveira2021unbiased state the unbiased estimation result for the mean-squared error estimation problem. They connect the coupled bootstrap estimator to Stein's unbiased risk estimate. Our calculation for other loss functions extends their unbiased estimation result. (ref) can also be easily generalized to other loss functions that admit unbiased estimators voinov2012unbiased.}
(ref) allows for an out-of-sample evaluation of decision rules, as well as uncertainty quantification around the estimate of loss, solely imposing the Gaussian model. This is a useful property in practice for comparing different empirical Bayes methods, especially if one is worried about the misspecification of (ref) or if one is unwilling to evaluate risk integrating over random $\theta_i$.
How does close-npmle perform in the field? We now consider two empirical exercises related to chetty2018opportunity and bergman2019creating. Using Census micro-data, chetty2018opportunity estimate a suite of tract-level children's outcomes in adulthood and publish an “Opportunity Atlas” of the estimates and the corresponding standard errors.\footnote{ Like prior work that uses this data andrews2021inference, we do not have access to the variance-covariance matrix of these estimates. Correlations across estimates are due to small proportion of movers between tracts and are anticipated to be small. } Taking these estimates, bergman2019creating conducted a program called {Creating Moves to Opportunity}. bergman2019creating provided assistance to treated low-income individuals to move to Census tracts with estimated posterior means in the top third. We view bergman2019creating's objectives as top-{ $m$} selection, for $m$ equal to one third of the number of tracts in Seattle and King County, WA.
The Opportunity Atlas published by chetty2018opportunity also includes tract-level covariates, a complication that we have so far abstracted away from. In the ensuing empirical exercises, following bergman2019creating, the estimates are residualized against the covariates as a preprocessing step fay1979estimates.\footnote{Alternatively, (ref) shows that flexibly modeling $\E[\theta_i \mid \sigma_i, X_i] = m_0(\sigma_i, X_i)$ and $\var(\theta_i \mid \sigma_i, X_i) = s_0^2(\sigma_i, X_i)$, as in (ref), induces substantial additional benefits, relative to simply projecting out the covariates linearly. Here, including $\sigma_i$ in the modeling remains important---modeling $\theta_i \mid X_i$ flexibly does not fully capture these benefits.} We now let $\tilde Y_i$ denote the raw Opportunity Atlas estimates for a pre-residualized parameter $\vartheta_i$ and let $(Y_i, \theta_i)$ be their residualized counterparts against a vector of tract-level covariates $X_i$, with regression coefficient $\beta$.\footnote{Precisely speaking, let $X_i$ be a vector of tract-level covariates. Let $(\tilde Y_i, \sigma_i)$ be the raw Opportunity Atlas estimates of a parameter $\vartheta_i$. Let $\beta$ be some vector of coefficients, typically estimated by weighted least-squares of $Y_i$ on $X_i$. Let $Y_i = \tilde Y_i - X_i'\beta$ and $\theta_i = \vartheta_i - X_i'\beta$ be the residuals. Since $\beta$ is precisely estimated, we ignore its estimation noise. Then, the residualized objects $(Y_i, \theta_i)$ obey the Gaussian sequence model $Y_i \mid \theta_i, \sigma_i \sim \Norm(\theta_i, \sigma_i^2).$ } We can apply the empirical Bayes procedures in this paper to $ (Y_i, \sigma_i^2)$ and obtain an estimated posterior for $\theta_i$. This estimated posterior for the residualized parameter $\theta_i$ then implies an estimated posterior for the original parameter $\vartheta_i = \theta_i + X_i' \beta$, by adding back the fitted values $X_i'\beta$. When there are no covariates, $\vartheta_i = \theta_i$ and $Y_i = \tilde Y_i$.
The covariates we use are included in the publicly available data from chetty2018opportunity and cross-referenced with their labels in (ref). They include: poverty rate in 2010, share of Black individuals in 2010, mean household income in 2000, log wage growth for high school graduates, fraction with college or post-graduate degrees in 2010, mean parent family income rank, mean parent family income rank for Black individuals, number of all and Black children under 18 with parents whose household income is below median in 2000 (in both levels and logs).
We consider 15 measures of economic mobility $\vartheta_i$. Each $\vartheta_i$ is the population mean of some outcome for individuals of some demographic subgroup growing up in tract $i$, whose parents are at the 25\th income percentile.\footnote{\Copy{fnalpha}{Since all measures of economic mobility have bounded support, as either percentile ranks or percentage rates, (ref) is automatically satisfied for $\theta_i$ with $\alpha = 2$, at least when there are no covariates.}} We will consider three types of outcomes:
for the following five demographic subgroups: all individuals (pooled), white individuals, white men, Black individuals, and Black men. Under these shorthands, the outcome in (ref) is top-{ 20} probability (Black), while bergman2019creating consider mean rank pooled.
The remainder of this section compares several methods on two exercises. In the first exercise, a calibrated simulation, we compare MSE performance of various methods to that of the oracle posterior. The second exercise is an empirical application to a scale-up of the exercise in bergman2019creating. It uses the coupled bootstrap ((ref)) to evaluate whether close-npmle selects more economically mobile tracts than alternatives.
\Copy{simdgp}{ We draw from a data-generating process estimated from the data. This data-generating process does not impose the location-scale assumption. On the data $(Y_i, \sigma_i)$, we estimate $\hat m(\cdot), \hat s^2(\cdot)$ via local linear regression. We then transform to obtain $\hat Z_i = \frac{Y_i - \hat m(\sigma_i)}{\hat s(\sigma_i)}$ and $\hat\nu_i = \frac{\sigma_i}{\hat s(\sigma_i)}$. We partition $\sigma_i$ into vingtiles. For the data $(\hat Z_i, \hat\nu_i)$ whose $\sigma_i$ falls in a given vingtile $v \in \br{1,2,3,4,5}$, we estimate a vingtile-specific $\hat G_{n,v}$ via npmle. We then normalize this estimated npmle to have mean zero and variance one, by affinely transforming the estimated distribution. Finally, to generate synthetic data, for a $\sigma_i$ corresponding to the $v(\sigma_i)$\th vingtile, we draw $\tau_i^* \mid \sigma_i \sim \hat G_{n,v(\sigma_i)}^{\text{normalized}}$, and set $\theta_i^* = \tau_i^*\hat s (\sigma_i) + \hat m(\sigma_i)$, $Y_i^* \mid \theta_i^*, \sigma_i \sim \Norm(\theta^*_i, \sigma_i^2)$ and $\tilde Y_i^* = Y_i^* + X_i'\beta$. Additional details for the sampling process and simulation setup are documented in (ref). }
On the simulated data, we then implement various empirical Bayes strategies. We consider the feasible procedures: naive, independent-gauss, independent-npmle, close-gauss (parametric), close-\textsc{gauss}, and \textsc{close}-\textsc{npmle}, as well as the infeasible \textsc{oracle}. Here,
(ref) plots the results from this calibrated simulation, focusing on MSE performance. For each method and each target variable, we display a relative measure of MSE gain. For each method, we calculate its MSE gain over naive, normalized by the MSE gain of oracle over naive. If we think of the oracle--naive difference as the total size of the “statistical pie,” then (ref) shows how much of this pie each method captures.
The first five columns show the relative mean-squared error performance {without} residualizing against covariates, applying empirical Bayes methods directly on $ (\tilde Y_i, \sigma_i)$. We see that methods which assume precision independence perform worse than methods based on close.\footnote{It may be surprising that independent-gauss can perform worse than naive even on MSE, since Gaussian empirical Bayes can be thought of as optimizing among a class of linear shrinkage estimators that include naive. We note that, as in bergman2019creating, when we estimate the prior mean and prior variance, we weight the data with precision weights proportional to $1/\sigma_i^2$. When the independence between $\theta$ and $\sigma$ holds, these precision weights typically improve efficiency. However, the weighting does mean that the resulting posterior means are no longer optimal, even asymptotically, among the class of linear shrinkage rules under precision dependence. To take an extreme example, if a particular observation has $\sigma_i \approx 0$, then that observation is highly influential for the prior mean estimate. If $\E[\theta_i \mid \sigma_i]$ is very different for that observation than the other observations, then the estimated prior mean is a bad target for shrinkage. } Across the 15 variables, the median proportion of possible gains captured by independent-gauss is only 31%. This value is 50% for \textsc{independent-npmle}, and 86% for \textsc{close}-\textsc{npmle}. Among the first five columns, \textsc{close}-\textsc{npmle} uniformly dominates all three other methods. This indicates that the standard error $\sigma_i$ is highly predictive of $\theta_i$, and using that information can be very helpful in the absence of additional covariates.
The next five columns show performance when the methods do have access to covariate information. For mean rank, after covariate residualization, the dependence between $\theta_i$ and $\sigma_i$ does not appear to substantially affect shrinkage decisions. independent-npmle and close-methods perform similarly, capturing almost all of the available gains. For the other two outcome variables, top-{ 20} probability and incarceration, the dependence between $\theta_i $ and $ \sigma_i$ is stronger, and close-based methods display substantial improvements over methods that assume precision independence. Among \textsc{close}-methods, those that are more flexible appear to reap a small benefit, though simple parametric models for $(m_0, s_0, G_0)$ remain competitive and significantly improve upon methods that assume precision independence. The most flexible method, \textsc{close}-\textsc{npmle}, achieves near-oracle performance across the different definitions of $\theta_i$ and again uniformly dominates all other feasible methods.\footnote{(ref) contains an alternative data-generating process in which the $\theta_i \mid \sigma_i$ distribution is Weibull, which has thicker tails and higher skewness. Under such a scenario, \textsc{npmle}-based methods more substantially outperform methods assuming Gaussian priors.}
Our second empirical exercise uses the coupled bootstrap described in (ref) for the policy problem in bergman2019creating. Viewing the policy problem in bergman2019creating as top-{ $m$} selection, can close-npmle make better selections?
Specifically, we imagine scaling up bergman2019creating's exercise and perform empirical Bayes procedures for all Census tracts in the largest 20 Commuting Zones (CZs). We then select the top third of tracts within each CZ, according to empirical Bayesian posterior means for $\vartheta_i$. Additionally, to faithfully mimic bergman2019creating, here we perform all empirical Bayes procedures within CZ. Throughout, we choose $\omega$ to emulate a 90-10 train-test split on the micro-data. See (ref) for details on the policy exercise setup.
(ref) shows the estimated performance of various methods. According to these estimates, close-npmle generally improves over independent-gauss.\footnote{close-npmle is worse by an estimated $0.006$ percentile ranks for mean rank \textsc{pooled} and worse by $0.04$ percentile ranks for \textsc{mean rank} for white men. In either case, the estimated disimprovement is small.} Strikingly, \textsc{independent-gauss} with covariates underperforms \textsc{naive} for four of the 15 variables, and \textsc{independent-gauss} without covariates underperforms for nearly all variables.
For the mean rank variables, using close-npmle generates substantial gains for mobility measures for Black individuals (0.63 percentile ranks for Black men and 0.43 percentile ranks for Black individuals). To put these gains in dollar terms, at the income level for experiment participants in bergman2019creating, an incremental percentile rank amounts to about \$1,000 per annum. Thus, the estimated gain in terms of mean income rank is roughly \$400--600. For the other two outcomes, top-{ 20} probability and incarceration, the gains are even more sizable. These gains are as high as 2--3 percentage points on average. Among close-methods, we again find that \textsc{close}-\textsc{npmle} generally performs the best, though by small margins.\footnote{Interestingly, the best performing method for \textsc{mean rank} (\textsc{pooled}) and \textsc{mean rank} (white men) is \textsc{close}-\textsc{gauss} (parametric), and the best performing method for \textsc{mean rank} (Black) and \textsc{mean rank} (Black men) is \textsc{close}-\textsc{npmle}, but without residualizing against covariates.} \Copy{empirical}{While \textsc{close}-\textsc{npmle} is a simple default that works uniformly well, in this case, simple parametric models that allow for dependence also appear competitive.}
We can think of the performance gap between independent-gauss and naive as the value of basic empirical Bayes. If practitioners find using the standard empirical Bayes method a worthwhile investment over screening on the raw estimates directly, perhaps they reveal that the value of basic empirical Bayes is economically significant. Across the 15 measures, the improvement of close-npmle over independent-gauss is on median 260% of the value of basic empirical Bayes, where the median is attained by \textsc{mean rank} for Black individuals. Thus, the additional gain of \textsc{close}-\textsc{npmle} over \textsc{independent-gauss} is substantial compared to the value of basic empirical Bayes. If the latter is economically significant, then it is similarly worthwhile to use \textsc{close}-\textsc{npmle} instead.
This paper studies empirical Bayes methods in the heteroskedastic Gaussian location model. We argue that precision independence---the assumption that the precision of estimates does not predict the true parameter---is often empirically rejected. Empirical Bayes methods that rely on precision independence can generate worse posterior mean estimates. Screening decisions based on these estimates can suffer as a result. They may even be worse than the selection decisions made with the unshrunk estimates directly.
Instead of treating $\theta_i$ as independent from $\sigma_i$, we model its conditional distribution as a location-scale family in $\sigma$-dependent location and scale parameters. This assumption leads naturally to a family of empirical Bayes strategies that we call close. The close-framework naturally subsumes and generalizes several existing proposals for accommodating precision dependence. We prove that close-npmle attains minimax-optimal rates in Bayes regret, extending previous theoretical results. That is, it approximates infeasible oracle Bayes posterior means as competently as statistically possible. Additionally, we show that an idealized version of close-npmle is robust, with finite worst-case Bayes risk. Finally, we further connect our main theoretical results to ranking-type decision problems in bergman2019creating.
Simulation and validation exercises demonstrate that close-npmle generates sizable gains relative to the standard parametric empirical Bayes shrinkage method. Across calibrated simulations, close-npmle attains close-to-oracle mean-squared error performance. In a hypothetical, scaled-up version of bergman2019creating, across a wide range of economic mobility measures, close-npmle consistently selects more mobile tracts than does the standard empirical Bayes method. The gains in the average economic mobility among selected tracts, relative to the standard empirical Bayes procedure, are often comparable to---or even multiples of---the value of basic empirical Bayes.