The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
105,787 characters
Empirical Bayes When Estimation Precision Predicts Parameters
{\onehalfspacing
\begin{abstract}
Gaussian empirical Bayes methods usually maintain a \emph{precision independence}
assumption: The unknown parameters of interest are independent from the known standard errors of the
estimates.
This assumption is often theoretically questionable and empirically rejected. This paper
proposes to model the conditional distribution of the parameter given the standard errors
as a flexibly parametrized location-scale family of distributions, leading to a family of
methods that we call \textsc{close}. The \textsc{close}{} framework unifies and generalizes several
proposals under precision dependence. We argue that the most flexible member of
the \textsc{close}{} family is a minimalist and computationally efficient default for accounting
for precision dependence. We analyze this method and show that it is competitive in terms
of the regret of subsequent decisions rules.
Empirically, using \textsc{close}{} leads to sizable gains for selecting high-mobility Census
tracts.
\vspace{1em}
\noindent \textsc{JEL codes.} C10, C11, C44
\noindent \textsc{Keywords.} Empirical Bayes, $g$-modeling, regret, heteroskedasticity,
nonparametric maximum likelihood, Opportunity Atlas, Creating Moves to Opportunity
\end{abstract}
\maketitle
}
\newpage
\onehalfspacing
\section{Introduction}
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 \citep
{angrist2017leveraging,mountjoy2021returns,chandra2016health,doyle2017evaluating,hull2018estimating,einav2022producing,abaluck2021mortality},
place-based effects \citep
{chyn2021neighborhoods,finkelstein2021place,chetty2018opportunity,chetty2018impacts,diamond2021standard,baum2019microgeography,aloni2023one},
discrimination \citep
{kline2022systemic,kline2023discrimination,rambachan2021identifying,egan2022harry,arnold2022measuring,montiel2021empirical},
meta-analysis \citep
{azevedo2020b,meager2022aggregating,andrews2019identification,elliott2022detecting,wernerfelt2022estimating,dellavigna2022rcts,abadie2023estimating},
and correlated random effects in panel data \citep
{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 \emph{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 \citep{bergman2019creating} uses empirical Bayes
methods to shrink raw economic mobility estimates $(Y_i, \sigma_i)$ of low-income
children, curated by \citet{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.
\citet{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 \emph
{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
\emph{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 \emph{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 \citet
{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
\eqref{eq:gaussian_heteroskedastic_location}.
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 \emph{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,
\emph{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. \Cref {sec:model} 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 \underline{c}onditional \underline{lo}cation-\underline{s}cale
\underline{e}mpirical Bayes methods---which we call \textsc{close}---by estimating the
hyperparameters $(m_0 (\sigma), s_0(\sigma),G_0)$. The \textsc{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 \textsc{close}{} framework unifies and generalizes several proposals in the literature \citep
[among others,][]
{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 \emph{nonparametric maximum likelihood} (\textsc{npmle}) to estimate $G_0$
\citep{kiefer1956consistency,jiang2009general,koenker2014convex}. We refer to this
variant as \textsc{close}-\textsc{npmle}. We view \textsc{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 \textsc{close}-\textsc{npmle}{} in
\cref{sec:regret}. Our main result (\cref{cor:maintext,thm:minimaxlower}) establishes
that, under the
\textsc{close}{} assumptions, \textsc{close}-\textsc{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 \emph{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 \citep
{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 \textsc{close}-\textsc{npmle}{} to the \textsc{close}{} assumption, we study a population version of \textsc{close}-\textsc{npmle}
{} under misspecification of the location-scale model. \Cref{thm:worstcaserisk} 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 \citet{bergman2019creating}. \Cref
{thm:mserelevance} shows that the Bayes regret in squared error dominates the Bayes
regret for these other decision problems. Coupled with
\cref{cor:maintext}, this implies that \textsc{close}-\textsc{npmle}{} has good performance for these
ranking-related problems as well.
To illustrate our method, \cref{sec:empirical} applies \textsc{close}{} to two empirical
exercises \citep{chetty2018opportunity,bergman2019creating}. The
first exercise is a simulation calibrated to the Opportunity Atlas, the
dataset published by \citet{chetty2018opportunity}. For all 15
measures of economic mobility that we consider,
\textsc{close}-\textsc{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 \citep{bergman2019creating}, using an out-of-sample
validation procedure based on the coupled bootstrap that we introduce
\citep{oliveira2021unbiased}. \citet {bergman2019creating} use empirical Bayes procedures
to select high-mobility Census tracts in Seattle. In an exercise that mimics theirs, we
find that \textsc{close}-\textsc{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
\textsc{close}-\textsc{npmle}{} over the standard method are on median 2.6 times the \emph{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 \textsc{close}{} is likewise meaningful.
\section{Model and proposed method}
\label{sec:model}
\subsection{Empirical Bayes assumptions} 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
\citep{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$ \eqref{eq:gaussian_heteroskedastic_location}. The Gaussian model \eqref
{eq:gaussian_heteroskedastic_location} 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 \citep[cf. Example 2 in][]
{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
\eqref{eq:gaussian_heteroskedastic_location} 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
(\cref{sec:empirical}). 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 \cref{rmk:independence}.
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
\eqref{eq:oracle_bayes}:\footnote{To
emphasize the distinction between the true
expectation with respect to the data-generating process
\eqref{eq:eb_sampling} 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,
\eqref{eq:empirical_bayes_rule} generates empirical Bayes posterior means $\mathbf{E}_{\hat P}
[\theta_i \mid Y_i, \sigma_i]$, often referred to as shrinkage estimates
\citep{james1992estimation,efron1973stein}.
To simplify the estimation of $P_0$, popular empirical
Bayes methods often assume \emph{precision independence}: $\theta_i \indep \sigma_i$, or,
equivalently, $G_{(1)} =
\cdots = G_{(n)}$ in \eqref{eq:eb_sampling} 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)$ \citep{morris1983parametric}. State-of-the-art empirical Bayes methods
relax the parametric assumptions on $G_{(0)}$ and estimate $G_{(0)}$ with
\emph{nonparametric maximum likelihood}, or \textsc{npmle}{}
\citep{jiang2020general,gilraine2020new,soloff2021multivariate}. Henceforth, we refer to
these methods as \textsc{independent-gauss}{} and \textsc{independent-npmle}{}, respectively. The
``\textsc{independent}'' here emphasizes precision independence.
\subsection{Precision independence and its violation}
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 \citep{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.
\begin{figure}[!htb]
\centering
\includegraphics[width=0.9\textwidth]{final_assets/example_raw.pdf}
\begin{proof}[Notes]
All tracts within the largest 20 Commuting Zones are shown. Due to the regression
specification in \citet{chetty2018opportunity}, point estimates of $\theta_i \in [0,1]$
do not always lie within $[0,1]$. The orange line plots nonparametric regression
estimates of the conditional mean $\E[Y \mid
\sigma] = \E[\theta \mid \sigma] \equiv m_0(\sigma)$, estimated via local linear
regression implemented by \citet{calonico2019nprobust}. The orange shading shows a
95\% uniform confidence band, constructed by the max-$t$ confidence set over 50
equally spaced evaluation points. See \cref{sec:nuisance_estimation} for details on
estimating conditional moments of $\theta_i$ given $\sigma_i$.
\end{proof}
\caption{Scatter plot of $Y_i$ against $\log_{10}(\sigma_i)$ in
\citet{chetty2018opportunity}}
\label{fig:kfr_top20_black_raw_data}
\end{figure}
As this economic intuition predicts, precision independence is readily rejected for this
measure of economic mobility. \Cref {fig:kfr_top20_black_raw_data} 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. \Cref
{fig:kfr_top20_black_raw_data} shows the contrary---tracts with more imprecisely
estimated $Y_i$ indeed tend to have higher $\theta_i$.
\begin{figure}[!htb]
\centering
\includegraphics[width=0.9\textwidth]
{final_assets/example_eb_posterior_means.pdf}
\begin{proof}[Notes] The top left panel shows posterior mean estimates with \textsc{independent-gauss}.
The top right panel shows the same with \textsc{independent-npmle}. The bottom panel displays
posterior mean estimates from our preferred procedure, \textsc{close}-\textsc{npmle}. In the top panels,
the estimates for the unconditional mean and variance of $\theta_i$ are weighted by
the precision $1/\sigma_i^2$, following \citet{bergman2019creating}.
\end{proof}
\caption{Posterior mean estimates under precision independence}
\label{fig:kfr_top20_black_pooled_raw_grand_mean_shrink}
\end{figure}
What happens if we apply empirical Bayes methods that assume precision independence here?
\Cref{fig:kfr_top20_black_pooled_raw_grand_mean_shrink} overlays empirical Bayes posterior
means on the scatterplot. In the top left panel,
\textsc{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 \textsc{independent-gauss}{} {systematically undershoot}
$\theta_i$ for tracts with imprecise estimates. Similarly, the top right panel of
\cref{fig:kfr_top20_black_pooled_raw_grand_mean_shrink} shows that \textsc{independent-npmle}{} suffers
from the same undershooting. In contrast, the bottom panel of
\cref{fig:kfr_top20_black_pooled_raw_grand_mean_shrink} previews our preferred procedure,
\textsc{close}-\textsc{npmle}, which shrinks towards the conditional mean $\E[\theta_i \mid
\sigma_i]$, thus avoiding the undershooting.
\Copy{msevsranking}{Nonetheless, posterior means from \textsc{independent-gauss}{} or \textsc{independent-npmle}{} may
still be better predictors, on average, for $\theta_i$ in mean-squared error than the
noisy $Y_i$ \citep{james1992estimation,efron1975data}. However, the undershooting for
large $\sigma_i$ is particularly problematic if one hopes instead to \emph{select}
high-mobility Census tracts based on these posterior means, as do \citet
{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$ \citep[see, also, ][]
{mehta2019measuring}.
\begin{figure}[htb]
\centering
\includegraphics[width=0.9\textwidth]{final_assets/example_shrink_ranking.pdf}
\begin{proof}[Notes] This plot shows a subregion of \cref
{fig:kfr_top20_black_raw_data} and highlights two Census tracts. The two tracts are
those with $\log_{10}(\sigma_i) < -1.1$ with the highest $Y_i$, for which the selection
decisions from \textsc{independent-gauss}{} and \textsc{close}-\textsc{npmle}{} disagree. Like
\citet{bergman2019creating}, the selection decisions aim to select 1/3 of Census tracts
so as to maximize the average $\theta_i$ selected, by screening for the top 1/3 of
empirical Bayes posterior means (formally, see \cref {ex:topm}).
\end{proof}
\caption{Ranking decisions for two Census tracts}
\label{fig:example_shrink_ranking}
\end{figure}
To see this, \cref{fig:example_shrink_ranking} zooms into a subregion of
\cref{fig:kfr_top20_black_raw_data} 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, \textsc{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$---\textsc{independent-gauss}{} recommends tract $A$ over
$B$ instead. In contrast, our preferred method (\textsc{close}-\textsc{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.
\subsection{Conditional location-scale modeling of precision dependence}
\label{sub:close}
We propose the following \emph{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
\begin{align*}
\theta_i = m_0(\sigma_i) + s_0(\sigma_i) \tau_i &\quad\text{ where }\quad \tau_i \mid \sigma_i \iid
G_0 \quad \eta_0(\cdot) \equiv (m_0(\cdot), s_0(\cdot)) .
\addtocounter{equation}{1}\tag{\theequation}
\label{eq:location_scale}
\end{align*}
\eqref{eq:location_scale} 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 \emph{location} of the distribution and the function $s_0(\cdot)$ controls
the \emph{scaling}. The underlying \emph{shape} of the distribution is governed by $\tau_i
\sim G_0$ and
is restricted by \eqref{eq:location_scale} 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 \eqref{eq:empirical_bayes_rule} 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
\eqref{eq:location_scale} 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}
\]
\eqref{eq:location_scale_tau_form} makes clear that, first, the transformed triplet $(Z_i,
\tau_i, \nu_i)$ obeys an analogue of the Gaussian model
\eqref{eq:gaussian_heteroskedastic_location}, where $Z_i$ is a noisy Gaussian
signal on $\tau_i$ with variance $\nu_i^2$. Second, precision independence holds in
\eqref{eq:location_scale_tau_form}, 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 \textsc{close}:
\begin{center}
\begin{minipage}{0.95\textwidth}
\begin{enumerate}[label=$\boxed{\textbf{\textsc{close}--\textsc{step} {\small\arabic*}}}$,wide]
\item \label{item:close1} Estimate $m_0(\sigma), s_0^2(\sigma)$ according to
\eqref{eq:m_s_def}.
\item \label{item:close2} With the estimates $\hat\eta = (\hat m, \hat s)$, transform the
data according to
\eqref{eq:transformed_data}. Apply empirical Bayes methods under precision independence
to
estimate $G_0$ with some $\hat G_n$ on the transformed data $
(\hat Z_i, \hat \nu_i)$.
\item \label{item:close3} Having estimated $(\hat\eta, \hat G_n)$ and hence having
obtained $\hat P$, we then form empirical Bayes decision rules following \eqref
{eq:empirical_bayes_rule}.
\end{enumerate}
\end{minipage}
\end{center}
This framework produces a family of empirical Bayes strategies, since
\cref{item:close1,item:close2} 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 \textsc{close}-\textsc{npmle}. It uses
nonparametric regression for \cref{item:close1} and \textsc{npmle}{} for
\cref{item:close2}. We recommend this method as a flexible default and primarily analyze
it in \cref{sec:regret}. We conclude this section with several
self-contained discussions on implementations of
these two steps, the rationale for \eqref{eq:location_scale} and
\textsc{close}-\textsc{npmle}, and other miscellaneous issues.
\subsection{Discussions}
\label{sub:model_discussions}
\subsubsection{Implementation}
\label{sub:implementation}
For \cref{item:close1}, one can exploit \eqref{eq:m_s_def} 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 \citep[see][ for suggestions of the latter]
{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 \cref{rmk:practical}), 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 \cref{item:close2}, 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 \textsc{npmle}{} to estimate $G_0$ \citep
{koenker2019comment}. Formally, the
\textsc{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 \eqref{eq:npmle}
\citep{koenker2014convex}.\footnote{\citet{koenker2017rebayes} provide an efficient
software implementation for \eqref{eq:npmle}, 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 (\cref{as:npmle}).
}
\Copy{closegauss}{On the other hand, a default \emph{parametric} model for $G_0$ is to
simply assume that $G_0 \sim \Norm(0,1)$, which we refer to as \textsc{close}-\textsc{gauss}. This
approach amounts to using
\textsc{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 \eqref{eq:gaussian_cond_b}
enjoys strong robustness properties\footnote{\Cref{thm:worstcaserisk} shows that oracle versions of \textsc{close}-\textsc{npmle}
{} satisfy analogous but weaker robustness properties when the location-scale model
fails. } even without the location-scale model
\eqref{eq:location_scale} and the assumption that $G_0 \sim \Norm(0,1)$.
First,
\eqref{eq:gaussian_cond_b} is the optimal linear-in-$Y$ decision rule for estimating
$\theta_i$ in
squared error
\citep{weinstein2018group}; second,
\eqref{eq:gaussian_cond_b} 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 \cref{lemma:optimal_bayes_linear,lemma:minimax_close_gauss} for formal statements,
respectively). This method performs almost as well as \textsc{close}-\textsc{npmle}{} in our empirical
exercises.
\subsubsection{The location-scale assumption and \textsc{close}-\textsc{npmle}}
\label{subsub:rationale}
\begin{table}[tb]
\caption{Various existing methods fit into the \textsc{close}{} framework}
\label{tab:closetable}
\centering
\scriptsize
\begin{tabularx}{\textwidth}{m{0.24\textwidth}XX}
\toprule
& Step 1 & Step 2 \\ \midrule
\citet{weinstein2018group} & Partition-based nonparametric estimator for $m_0, s_0^2$ &
$G_0 \sim \Norm
(0,1)$ \\
\citet{george2017mortality} & Parametric models for $m_0, s_0^2$ & $G_0 \sim \Norm(0,1)$
\\
\citet{chamberlain1984panel} & Parametric models for $m_0,
s_0^2$ & $G_0 \sim \Norm(0,1)$ \\
\citet{efron2016empirical} & Constant $m_0, s_0^2$ & $G_0$ nonparametric (log-spline sieves)
\\
\citet{kline2023discrimination} & $m_0 = c_1 s(\cdot)$, $s_0 = c_2 s(\cdot)$ for
parametric $s(\cdot)$ & $G_0$ nonparametric (log-spline sieves) \\
\textsc{independent-npmle} & Constant $m_0, s_0^2$ & $G_0$ nonparametric \\
\textsc{independent-gauss} & Constant $m_0, s_0^2$ & $G_0 \sim \Norm(0,1)$ \\
\citet{ignatiadis2019covariate} & Nonparametric $m_0$, constant $s_0^2$ & $G_0 \sim \Norm
(0,1)$ \\
\citet{jiang2010empirical} & Constant $m_0$, $s_0(\sigma) = \sigma$
(see \cref
{rmk:alts_to_close}) & $G_0$ nonparametric
\\\bottomrule
\end{tabularx}
\end{table}
We argue that the location-scale assumption provides a unifying framework for a number of
existing methods, and \textsc{close}-\textsc{npmle}{} is a natural generalization of these methods within
this framework. We also briefly speculate how to generalize beyond \textsc{close}-\textsc{npmle}.
Several existing methods can be thought of as implementations of
\textsc{close}{} by making different choices in \cref{item:close1} and \cref{item:close2}. \Cref
{tab:closetable} summarizes how these methods fit into the \textsc{close}{} framework. Among these
methods, some choose nonparametric models for \cref{item:close1} and some choose
nonparametric models for \cref{item:close2}. For instance, \citet{weinstein2018group}
propose
\textsc{close}-\textsc{gauss}, with a partition-based nonparametric estimator for $m_0, s_0^2$.
\citet{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 \citep{efron2016empirical}.
\citet{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,
\textsc{close}-\textsc{npmle}{} may be particularly attractive due to its minimalism: The \textsc{npmle}{} is free of
tuning parameters \citep{koenker2019comment}, and tuning parameter choices for
nonparametric regression are relatively well-understood
\citep{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 \textsc{close}-\textsc{npmle}{} naturally generalizes the existing methods in
\cref{tab:closetable}, one might consider methods that do not impose
\eqref{eq:location_scale} 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 \citep
{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
\eqref{eq:location_scale} holds. See \cref{lemma:transform} for a precise statement.}
In this sense, these methods must depart substantially from those that impose precision
independence.}
A natural approach is to estimate \textsc{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
\textsc{independent-npmle}{} within each bin.\footnote{Our Monte Carlo exercise in \cref{sec:empirical}
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 \textsc{close}-\textsc{npmle}{} performs well relative to the oracle and thus
to this procedure (\cref{fig:mse_table}).} 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
\textsc{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 \textsc{close}-\textsc{npmle}.
\subsubsection{Additional remarks}
\begin{rmksq}[Negative $\hat s^2$ estimates]
\label{rmk:practical}
Analogue estimators for $s_0^2(\sigma_i) =
\var(Y_i \mid \sigma_i) - \sigma_i^2$ may take negative values.\footnote{The negative
estimated variance phenomenon is in part caused by estimation noise in $\var(Y_i \mid
\sigma_i)$. However, in our empirical application, there is some evidence that
observations with large estimated $\sigma_i$'s are underdispersed for the measures of
economic mobility in the Opportunity Atlas (see \cref{asub:variance_right_tail}).
\citet{armstrong2022robust} propose a Bayesian estimator for the conditional variance. }
In our experience, truncating $\hat s$ at zero does not seem to cause bad performance when
computing posterior means. Nevertheless, in \cref{sec:nuisance_estimation} and the
software implementation, we propose a heuristic but data-driven truncation rule that
produces strictly positive $\hat s$, borrowing from a statistics literature on estimating
non-centrality parameters for non-central $\chi^2$ distributions
\citep{kubokawa1993estimation}.
\end{rmksq}
\begin{rmksq}[Other transformations]
\label{rmk:alts_to_close} We summarize and compare \textsc{close}{} to two methodological
alternatives, deferring a detailed discussion on these and on several others to
\cref{sub:alt_methods}. First, \citet{jiang2010empirical} propose applying \textsc{npmle}{} on
the $t$-ratio $Z_i = Y_i / \sigma_i \sim \Norm(\theta_i/\sigma_i, 1)$; similar approaches
are used in \citet{efron2016empirical,kline2022systemic}. For estimating $\theta_i$, one
then uses $\hat\theta_i =
\sigma_i \cdot \mathbf{E}_{\hat G_n}[\theta_i/\sigma_i \mid Z_i]$. Interpreting $\hat\theta_i$
as an estimated posterior mean $\E_{P_0}[\theta_i \mid Y_i, \sigma_i]$ requires that
$\theta_i / \sigma_i \indep \sigma_i$---meaning that \eqref{eq:location_scale} holds
with $s_0(\sigma_i) = \sigma_i$ and constant $m_0(\cdot)$. Thus this $t$-ratio approach
can be viewed as a particular instance of \textsc{close}, if we wish to imbue it with an
empirical Bayesian interpretation.
Second, when $Y_i$ and $\theta_i$ are sample and population means of binary outcomes, the
estimated variance of $Y_i$ is mechanically correlated with $\theta_i$: $\sigma_i^2 =
\frac{Y_i (1-Y_i)}
{n_i}.$ A variance-stabilizing transform, e.g. $Z_i = \arcsin\sqrt{Y_i}$
\citep{10.1214/07-AOAS138}, results in \emph{approximately} Gaussian $Z_i \sim \Norm(\arcsin{
\sqrt{\theta_i}}, \frac{1}{4n_i})$ without the mechanical dependence. However, it is still
possible that $n_i$ predicts $\theta_i$, and when that happens, proper modeling of
$\theta_i \mid n_i$---e.g., via an analogue of \eqref{eq:location_scale}---can continue
to improve performance.
\end{rmksq}
\section{Theoretical results}
\label{sec:regret}
As a review, we observe $(Y_i, \sigma_i)_{i=1}^n$, where $(\theta_i,\sigma_i)$ satisfies
\eqref{eq:location_scale} and $(Y_i, \theta_i, \sigma_i)$ obeys
\eqref{eq:gaussian_heteroskedastic_location}. The procedure
\textsc{close}-\textsc{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
\cref{item:close1}. It then estimates $G_0$ via \textsc{npmle}{}
\eqref{eq:npmle} on $(\hat Z_i, \hat \nu_i)_{i=1}^n$. This section introduces a few
statistical guarantees on the performance of \textsc{close}-\textsc{npmle}{} in terms of \emph{regret}. To
unify presentation, we first review decision theory primitives and introduce regret.
Let $\bm{\delta}(Y_{1:n}, \sigma_{1:n})$ be
a \emph{decision rule} mapping the data $(Y_{1:n},
\sigma_{1:n})$ to \emph{actions}. Recall that $L(\bm{\delta},
\theta_{1:n})$ denotes a \emph{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 \emph{Bayes risk} of $\bm{\delta}$ under $P_0$. The oracle Bayes
decision rule $\bm{\delta}^\star$ \eqref{eq:oracle_bayes} is optimal in the sense that it
minimizes $R_{\mathrm{B}}$. Thus, a natural performance measure for the empirical Bayesian
\eqref{eq:empirical_bayes_rule} is the gap between the Bayes risks of $\bm{\delta}_ {\mathrm{EB}}$ and
$\bm{\delta}^\star$. We refer to this quantity as \emph{Bayes regret}:
\begin{align*}
\mathrm{BayesRegret}_n(\bm{\delta}_{\mathrm{EB}}) &= \E_{P_0}
[L
(\bm{\delta}_{\mathrm{EB}}, \theta_{1:n}) - L(\bm{\delta}^\star, \theta_{1:n}) \mid \sigma_{1:n}],
\addtocounter{equation}{1}\tag{\theequation} \label{eq:regret_def}
\end{align*}
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 \textsc{close}-\textsc{npmle}{} vanishes quickly as a
function of $n$.
\begin{rmksq}[Fixed vs. random $\theta$]
\label{rmk:james-stein}
Our results consider asymptotic optimality, in terms of \eqref{eq:regret_def}, of the empirical Bayes
decision rule when $\theta_i \mid \sigma_i$ is randomly sampled from $P_0$, following a
recent literature on nonparametric empirical Bayes
\citep{jiang2020general,soloff2021multivariate}. A separate literature considers instead
the frequentist risk $R_{\mathrm{F}}(\theta_{1:n}; \sigma_{1:n}) \equiv
\E\bk{
L(\bm{\delta}, \theta_{1:n}) \mid \theta_{1:n}, \sigma_{1:n} }$ under fixed $(\theta_{1:n},
\sigma_{1:n})$ \citep{robbins1956}. For instance,
\citet
{james1992estimation,bock1975minimax,10.1214/07-AOAS138,weinstein2018group} consider
shrinkage estimators that dominate $\delta_i = Y_i$ uniformly for all configurations of
$\theta_{1:n}$. \citet{xie2012sure,kwon2021optimal} consider choosing decision rules
within a restricted class that minimize an unbiased estimate of $R_{\mathrm{F}}$. In
particular, \citet{xie2012sure} can be thought of as implementing \textsc{independent-gauss}{} with
different ways of estimating the hyperparameters in
$\theta_i \mid\sigma_i \iid \Norm\pr{m_0, s_0^2}$, and \citet{weinstein2018group} can be
thought of as implementing \textsc{close}-\textsc{gauss}.
While these guarantees for $R_{\mathrm{F}}$ are preserved even if we further average the frequentist
risk over $\theta_{1:n} \mid
\sigma_{1:n} \sim P_0$, they are distinct from upper bounding
\eqref{eq:regret_def}.\footnote{For instance, the oracle Bayes rule for mean-squared error may not
dominate $\delta_i = Y_i$ in $R_{\mathrm{F}}$ uniformly for all $\theta_{1:n}$. Conversely,
decisions that merely dominate $\delta_i = Y_i$ may still be quite far from the oracle
Bayes rule.} In particular, they may leave much on the table if $R_{\mathrm{B}}$ is
targeted. Moreover, these guarantees in $R_{\mathrm{F}}$ are typically restricted to MSE. Our
example in
\cref{fig:example_shrink_ranking} shows that reasonable decisions for MSE may not be
reasonable for subsequent selection decisions. As a simple example, \citet
{bock1975minimax} considers spherical shrinkage rules of the form $\delta_{i} = c\pr
{\sum_j Y_j^2} Y_i$ for some function $c(\cdot)$. However, despite dominating
no-shrinkage in MSE, $\delta_i$ does not change the ranking of different units, and
hence does not improve on ranks over $Y_i$.
\end{rmksq}
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}$.
\subsection{Regret rate in squared error}
Our main result concerns the canonical statistical problem of
estimating the parameters $\theta_{1:n}$ under MSE.
\begin{probsq}[Squared-error estimation of $\theta_{1:n}$]
\label{ex:mse}
The action $\bm{\delta} =(\delta_1,\ldots,
\delta_n)$ collects estimates $\delta_i$ for $\theta_i$, evaluated with MSE:
$
L(\bm{\delta}, \theta_{1:n}) =
\frac{1}{n}
\sum_{i=1}^n (\delta_i - \theta_i)^2. $ The oracle Bayes decision rule $\bm{\delta}^\star
= (\theta_1^*,\ldots, \theta_n^*)$ here is the posterior mean under $P_0$, where
$\theta_i^* \equiv
\E_{P_0}\bk{\theta_i \mid Y_i, \sigma_i} $. The empirical Bayesian counterpart is $
\hat\theta_{i, \hat P} = \mathbf{E}_{\hat P}[\theta_i \mid Y_i, \sigma_i].$
\end{probsq}
\Copy{mseregretdef}{
For \cref{ex:mse}, define $\mathrm{MSERegret}_n$ as the excess loss of the empirical Bayes posterior means
relative to that of the oracle Bayes posterior means:
\begin{align*}
\mathrm{MSERegret}_n(G, \eta) &\equiv \frac{1}{n} \sum_{i=1}^n (\hat\theta_{i, G, \eta} - \theta_i)^2 -
\frac{1}{n} \sum_{i=1}^n \pr{\theta_i^* - \theta_i}^2,
\end{align*}
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 \eqref{eq:regret_def} for \textsc{close}-\textsc{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 \eqref{eq:mse_regret_and_mse} 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 (\cref{cor:maintext,thm:minimaxlower}) 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$.
\subsubsection{Assumptions for regret upper bound}
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, \citet{soloff2021multivariate}). This assumption also
accommodates for the fact that the \textsc{npmle}{} is approximated by a discrete distribution on
a grid.
\begin{restatable}{as}{asnpmle}
\label{as:npmle}
Let $\psi_i(Z_i,
\hat\eta, G) \equiv \log\pr{\int_{-\infty}^\infty
\varphi
\pr{\frac{\hat Z_i -\tau}{\hat \nu_i}} G(d\tau)}$ be the objective function in
\eqref{eq:npmle}, ignoring the factor $1/\hat\nu_i$ that does not involve $G$. We
assume that $\hat G_n$ satisfies \[
\frac{1}{n} \sum_{i=1}^n \psi_i(Z_i, \hat\eta, \hat G_n) \ge \sup_{H \in \mathcal P(\R)}
\frac{1}{n} \sum_{i=1}^n \psi_i(Z_i, \hat\eta, H) - \kappa_n
\addtocounter{equation}{1}\tag{\theequation} \label{eq:approx_mle}
\]
for tolerance $\kappa_n = \frac{2}{n} \log ({\frac{n}{
\sqrt{2\pi} e}})$. Moreover, we
require that $\hat G_n$ has support points within $[\min_i\hat Z_i,
\max_i \hat
Z_i]$. To ensure that $\kappa_n$ is positive, we assume that $n \ge 7 = \lceil \sqrt{2\pi}
e \rceil$.\footnote{The constants $\kappa_n \rateeq \frac{1}{n}\log(n)$ also feature in
\citet{jiang2020general} to ensure that the fitted likelihood is bounded away from zero.
The particular constants in $\kappa_n$ simplify expressions and are not
material to the result.}
\end{restatable}
We now state further assumptions on $\mathcal P_0$ beyond
\eqref{eq:location_scale}. First, we assume that $G_0$ is sufficiently thin-tailed such
that its moments grow slowly.\footnote{An equivalent statement to \cref
{as:moments} 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
\citep[see the definitions in][]{vershynin2018high}. \Cref{as:moments} 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 \citep{soloff2021multivariate,jiang2009general,jiang2020general}. } \Copy{alpha}
{The thickness
of its tail is parametrized by $\alpha \in (0,2]$, which subsequently affects the
log factors in \cref{cor:maintext}.}
\begin{restatable}{as}{moments}
\label{as:moments}
The distribution $G_0$ has zero mean, unit variance, and admits simultaneous moment
control: For some $\alpha \in (0,2]$ and $A_0 > 0$ such
that for all $p > 0$, $\pr{\E_{\tau \sim G_0}[|\tau|^p]}^{1/p} \le A_0 p^
{1/\alpha}.
$
\end{restatable}
Next, \cref{as:variance_bounds} 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 \citet{jiang2020general} and
\citet{soloff2021multivariate}.
\begin{restatable}{as}{variancebounds}
\label{as:variance_bounds}
The variances $(\sigma_{1:n}, s_0)$ admit lower and upper
bounds: There are positive reals $\sigma_\ell, \sigma_u, s_{0\ell}, s_{0u} >0$ such that, for all $i$ and all
$\sigma \in (\sigma_\ell, \sigma_u)$, $\sigma_\ell < \sigma_i < \sigma_u$ and $s_{0\ell} < s_0(\sigma) < s_
{0u}$.
\end{restatable}
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]$ \citep[Section 2.7.1,][]{vaart1996weak}.
\begin{restatable}{as}{holder}
\label{as:holder} Assume that
\begin{enumerate}
\item The true conditional moments are H\"older-smooth: $m_0, s_0 \in C_{A_1}^p([\sigma_\ell,\sigma_u])$.
\end{enumerate}
Additionally, let $\beta_0 > 0$ be a constant. Assume
that the estimators for $m_{0}$ and $s_0$, $\hat\eta = (\hat m, \hat s)$, satisfy:
\begin{enumerate}[resume]
\item For all
sufficiently large $C_{1,\H} > 0$ and all $n$, \[
\P\pr{\norm{
\hat\eta - \eta_0
}_\infty > C_{1,\H} n^{-
\frac{p}
{2p+1}} (\log n)^{\beta_0}} < \frac{1}{n^2}
\] where $\norm{\eta}_\infty \equiv \max(\norm{m}_\infty, \norm{s}_\infty)$ for $\eta
=(m,s)$.
\item $\hat\eta$ takes values in $\mathcal V$ almost surely: $\P
\pr{\hat m\in
\mathcal V, \hat s \in \mathcal V} = 1$, where $\mathcal{V}$ is a set of
functions supported on $[\sigma_\ell, \sigma_u]$ that (i) is uniformly bounded $\sup_{f \in
\mathcal V}
\norm{f}_\infty \le C_{A_1}$ and (ii) admits the metric entropy bound $\log N (\epsilon,
\mathcal V, \norm{\cdot}_\infty) \le C_{A_1,p, \sigma_\ell,\sigma_u} (1/\epsilon)^ {1/p}$.
\item The conditional variance estimator respects the conditional variance bounds in
\cref{as:variance_bounds}: $\P\pr{\frac{s_{0\ell}}{2} < \hat s < 2s_{0u}} = 1$.
\end{enumerate}
\end{restatable}
\Cref{as:holder} 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{
\label{fn:a4}\Cref{as:holder}(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$
\citep{tsybakov2008introduction,stone1980optimal}. Since the data is assumed to be thin-tailed in
\cref{as:moments}, such estimators also attain the stronger requirement in
\cref{as:holder}(2).
For \cref{as:holder}(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,
\cref{as:holder}(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 \cref{sec:nuisance_estimation} that a local linear regression estimator (with
$\hat s$ suitably truncated) satisfies weaker conditions than \cref{as:holder}(2)--(4)
that are nonetheless sufficient for the conclusion of \cref{cor:maintext}.
}
\Cref{as:moments,as:holder,as:variance_bounds} 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$.
\subsubsection{MSE regret results}
Our main result is a non-asymptotic upper bound for \eqref{eq:mse_regret_and_mse}: The MSE
regret of
\textsc{close}-\textsc{npmle}{}
converges to zero no slower than
$n^{-\frac{2p}{2p+1}}(\log n)^{C}$.
\begin{restatable}{theorem}{cormaintext}
\label{cor:maintext}
Under \cref{as:holder,as:moments,as:variance_bounds,as:npmle}, there exists a
constant $C_{0, \mathcal{H}} > 0$ such that the following upper bound holds:\[
\mathrm{BayesRegret}_n = \E\bk{
\mathrm{MSERegret}_n(\hat G_n, \hat\eta) } \le C_{0, \mathcal{H}} n^{-\frac{2p}{2p+1}} (\log n)^{\frac{2+\alpha}{\alpha} + 3 + 2\beta_0}.
\addtocounter{equation}{1}\tag{\theequation} \label{eq:regret_rate_final}
\]
\end{restatable}
Second, we give a corresponding minimax lower bound on the regret, which shows that
\cref{cor:maintext} cannot be improved by more than logarithmic factors.
\begin{restatable}{theorem}{thmminimaxlower}
\label{thm:minimaxlower}
Fix a set of valid hyperparameters $\mathcal{H}$. Let $\mathcal P (\mathcal{H}, \sigma_{1:n})$ be the set of
distributions $P_0$ on support points $\sigma_{1:n}$ which satisfy
\eqref{eq:location_scale}
and \cref{as:moments,as:holder,as:variance_bounds} corresponding to
$\mathcal{H}$.\footnote{This result additionally takes the supremum over the
support points $\sigma_{1:n}$. This is because the nonparametric regression problem would
be ``too easy'' for certain configurations of $\sigma_{1:n}$. For instance, when
$\sigma_{1:n}$ only takes $m \ll n$ unique values, nonparametric regression is possible at
rate $\sqrt{m/n}$. For the proof, it suffices to consider $\sigma_{1:n}$ being equally
spaced in $
[\sigma_\ell, \sigma_u]$. } For a given $P_0$, let
$\theta_i^* = \E_{P_0}[\theta_i \mid Y_i,
\sigma_i]$ denote the oracle posterior means. Then there exists a constant $c_
{\mathcal{H}} > 0$ such that
\[
\inf_{\hat\theta_{1:n}} \sup_{\substack{\sigma_{1:n} \in (\sigma_\ell, \sigma_u)\\ P_0 \in \mathcal
P(\mathcal
H, \sigma_{1:n})}} \E_{P_0} \bk{
\frac{1}{n} \sum_{i=1}^n (\hat\theta_{i} - \theta_i)^2 - \frac{1}{n} \sum_{i=1}^n
(\theta_{i}^* - \theta_i)^2
} \ge c_\mathcal{H} n^{-\frac{2p}{2p+1}},
\]
where the infimum is taken over all (possibly randomized) estimators of $\theta_{1:n}$.
\end{restatable}
\Cref{cor:maintext} continues a recent statistics literature on empirical Bayes methods
via \textsc{npmle}, by characterizing the effect of an estimated first-step parameter $\hat\eta$.
Our theory hews closely to---and extends---the results in \citet{jiang2020general} and
\citet{soloff2021multivariate}, which themselves extend earlier results in the
homoskedastic setting \citep{jiang2009general,saha2020nonparametric}. In particular,
\citet{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
\cref{cor:maintext} is deferred to the Online Appendix, but its main ideas are outlined in
\cref{asec:proof_main}.
\Cref{thm:minimaxlower} shows that the rate \eqref{eq:regret_rate_final} is optimal up to
logarithmic factors. These logarithmic factors partly reflect inefficiencies in the proof
of \cref{cor:maintext}, but in any case the gap is not large. We prove
\cref{thm:minimaxlower} 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$ \citep
{tsybakov2008introduction} then imply lower bounds for estimation of the oracle
posterior means $\theta_i^*$
\citep[see][ for a similar argument in a related setting] {ignatiadis2019covariate}.
We additionally note that these regret upper bounds readily extend to the case where
covariates are present and the location-scale assumption \eqref{eq:location_scale} 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
\cref{as:holder}. 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, \cref{cor:maintext,thm:minimaxlower} are statistical optimality guarantees
for \textsc{close}-\textsc{npmle}{} in terms of \cref{ex:mse}. That is, the worst-case MSE performance gap
of
\textsc{close}-\textsc{npmle}{} relative to the oracle contracts at the best possible rate, meaning that
\textsc{close}-\textsc{npmle}{} mimics the oracle as well
as possible.
\subsection{Robustness to the location-scale assumption \eqref{eq:location_scale}}
We prove \cref{cor:maintext,thm:minimaxlower} imposing the location-scale model
\eqref{eq:location_scale}. This is an optimistic assessment of the performance of
\textsc{close}-\textsc{npmle}. While \eqref{eq:location_scale} nests precision independence, it may still be
misspecified. This subsection explores the worst-case behavior of \textsc{close}-\textsc{npmle}{} without
\eqref{eq:location_scale}.
We do so by considering an idealized version of \textsc{close}-\textsc{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
\eqref{eq:location_scale}, $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 \cref{asec:max_gauss} 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 \textsc{close}-\textsc{gauss}{} \eqref
{eq:gaussian_cond_b} is a minimax procedure (\cref{lemma:minimax_close_gauss}). }
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
\cref{thm:worstcaserisk} without this additional tail condition, regrettably due to a
technical error that is corrected in this version. See \cref{asec:max_gauss}.}
\begin{restatable}{theorem}{worstcaserisk}
\label{thm:worstcaserisk}
Under the preceding setup and \eqref{eq:eb_sampling}, but not \eqref{eq:location_scale},
let $\hat \theta_ {i, G_0^*, \eta_0}$ denote the posterior mean for $\theta_i$ under a
prior $G_0^*$
for $\tau$. Let $\bar\rho =
\max_i s_0^2 (\sigma_i) / \sigma_i^2 < \infty$ be the maximal conditional
signal-to-noise ratio. Then, for some $0 < C_{\bar\rho, \lambda,
\epsilon}
<
\infty$ that solely depends on $\bar\rho, \lambda, \epsilon$, \[
\frac{\sup_{G_0^* \in \mathcal G_0(\lambda, \epsilon)}
\sup_{P_0 \in \mathcal P(m_0, s_0)} \E_{P_0}\bk{\frac{1}{n} \sum_{i=1}^n (\hat\theta_{i,
G_0^*, \eta_0} - \theta_i)^2}}{\inf_{\hat\theta_{1:n}} \sup_{P_0 \in \mathcal P(m_0, s_0)} \E_{P_0} \bk{\frac{1}{n}
\sum_{i=1}^n (\hat\theta_i - \theta_i)^2}} \le C_
{\bar\rho, \lambda, \epsilon}, \addtocounter{equation}{1}\tag{\theequation}
\label{eq:multiple_of}
\]
where the infimum in the denominator is over all (possibly randomized)
estimators of $\theta_i$ given $(Y_i, \sigma_i)_{i=1}^n$ and $\eta_0(\cdot)$.
\end{restatable}
\Cref{thm:worstcaserisk} shows that the worst-case behavior of an idealized version of
\textsc{close}-\textsc{npmle}{} comes within a factor of the minimax risk. Thus, \textsc{close}-\textsc{npmle}{} is not
arbitrarily unreasonable, even under misspecification. We caution that
\eqref{eq:multiple_of} 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,
\eqref{eq:multiple_of}
{does not} hold for an idealized version of \textsc{independent-gauss}.\footnote{\Copy
{indepgaussconstant} {That is, it does not hold for the implementation of \textsc{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 \cref{lemma:indepgauss_not_bounded} for a formal statement.}
}
\subsection{Other decision objectives and relation to squared-error loss}
\label{sub:other_decisions}
So far, our regret guarantees are only about estimation in MSE (\cref{ex:mse}). We now
turn to two decision problems that involve ranking or selection and show similar
guarantees for \textsc{close}-\textsc{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
\citep{gilraine2020new,bergman2019creating,manski2004statistical,stoye2009minimax,kitagawa2018should,athey2021policy}.
\begin{probsq}[\textsc{utility maximization by selection}]
\label{ex:utilmax}
Suppose $\bm{\delta} = (\delta_1,\ldots,\delta_n)$ consists of binary selection decisions
$\delta_i \in \br{0,1}$. For each population, selecting that population has net benefit
$\theta_i$. The decision maker wishes to maximize utility (i.e., negative loss): $
-L(\bm{\delta}, \theta_{1:n}) =
\frac{1}{n} \sum_{i=1}^n \delta_i \theta_i. $ The oracle Bayes rule selects all whose
posterior mean net benefit $\theta_i$ is nonnegative: $
\delta_i^\star = \one\pr{\theta_{i, P_0}^* \ge 0}.
$
One natural empirical Bayes decision rule replaces $\theta_{i, P_0}^*$ with $\theta_{i,
\hat P}^*$, following \eqref{eq:empirical_bayes_rule}.
\end{probsq}
\begin{probsq}[\textsc{top-{\small $m$} selection}]
\label{ex:topm}
Similar to
\textsc{utility maximization by selection}, suppose $\bm{\delta}$ consists of binary selection decisions, with the additional
constraint that exactly $m$ populations are chosen: $\sum_i \delta_i = m$. The decision
maker's utility is the average $\theta_i$ of the
selected set: \[ -L(\bm{\delta}, \theta_{1:n}) = \frac{1}{m} \sum_{i=1}^n \delta_i \theta_i.
\addtocounter{equation}{1}\tag{\theequation}
\label{eq:topm}
\]
The oracle Bayesian selects the populations corresponding to the $m$ largest posterior
means
$\theta_{i, P_0}^*$: $
\delta_i^\star = \one\pr{
\theta_{i, P_0}^* \text{ is among the top-$m$ of $\theta_{1:n, P_0}^*$}
}.
$
Again, the empirical Bayes recipe \eqref{eq:empirical_bayes_rule} replaces $P_0$ with the
estimate $\hat P$.
\end{probsq}
\begin{rmksq}
\label{rmk:mover}
The utility function \eqref{eq:topm} rationalizes the widespread practice of screening
based on empirical Bayes posterior means
\citep{gilraine2020new,chetty2014measuring,kane2008estimating,hanushek2011economic,bergman2019creating}.
In \citet{bergman2019creating}, for instance, where housing voucher holders are
incentivized to move to Census tracts selected according to economic mobility,
\eqref{eq:topm} represents the expected economic mobility of a mover were they to move
randomly to one of the selected tracts. Our theoretical results can accommodate slightly
less restrictive mover behavior (\cref{rmk:nonuniform}).
\end{rmksq}
The oracle Bayes decision rules $\bm{\delta}^\star$ in \cref {ex:utilmax,ex:topm}
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
\cref{ex:utilmax} and let $\mathrm{TopRegret}_{n}^{(m)}$
denote $\mathrm{BayesRegret}_n$ for \cref{ex:topm}.
\begin{restatable}{theorem}{mserelevance}
\label{thm:mserelevance} Suppose \eqref{eq:eb_sampling} holds but \eqref
{eq:location_scale} does not necessarily hold.
Let $\hat\delta_i$ be the plug-in decisions with any vector of estimates $\hat \theta_i$.
Then,
\begin{enumerate}
\item For \textsc{utility maximization by selection},
\[
\E[\mathrm{UMRegret}_n(\hat\bm{\delta})] \le\pr { \E\bk{\frac{1}{n} \sum_{i=1}^n (\hat\theta_i - \theta_i^*)^2}}^
{1/2} . \addtocounter{equation}{1}\tag{\theequation} \label{eq:utilmax_bound}
\]
\item For \textsc{top-{\small $m$} selection}, \[
\E[\mathrm{TopRegret}_{n}^{(m)}(\hat\bm{\delta})] \le 2\sqrt{\frac{n}{m}} \pr{\E\bk{\frac{1}{n} \sum_{i=1}^n (\hat\theta_i -
\theta_i^*)^2}}^{1/2}. \addtocounter{equation}{1}\tag{\theequation} \label{eq:topm_regret_bound}
\]
\end{enumerate}
\end{restatable}
\Cref{thm:mserelevance} shows a sense in which \cref{ex:utilmax,ex:topm} are easier than
\cref{ex:mse}: The regret of the latter dominates those of the former. As a result, if we
use \textsc{close}-\textsc{npmle}{} under \eqref{eq:location_scale}, our convergence rates from
\cref{cor:maintext} also
upper bound regret rates for these two decision problems. In particular, for $m/n \to c
\in (0,1)$, both regret rates
\eqref{eq:utilmax_bound} and \eqref{eq:topm_regret_bound} are of the form $n^{-p/(2p+1)}
(\log n)^{C} = o(1)$ under \cref{cor:maintext}. 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.
\begin{rmksq}[Tightness of \cref{thm:mserelevance}]
We suspect that the actual performance of
\textsc{close}-\textsc{npmle}{} for \cref{ex:utilmax,ex:topm} may be better than predicted by
\cref{thm:mserelevance}. The proof of \cref{thm:mserelevance} exploits the fact that when
the empirical Bayesian makes a selection mistake, the size of the mistake is not large if
the square-error regret is low. It does not exploit the fact that if squared error regret
is low, then the empirical Bayesian may be unlikely to make mistakes in the first place.
\footnote{Upper and lower bounds are derived in related but distinct settings by
\citet{audibert2007fast,bonvini2023minimax}; some upper bounds, under possibly stronger
assumptions, appear better than implied by \cref{thm:mserelevance}.
We speculate that the bound for $\textsc{utility maximization by selection}$ can be tightened by verifying a margin
condition, using Proposition 2 in \citet{bonvini2023minimax}.
Relatedly,
\citet{liang2000empirical} shows upper and lower bounds for \cref{ex:utilmax} of the form
$O ((\log n)^{1.5} /n)$ in a homoskedastic setting, assuming the oracle posterior means
fall on both sides of zero.
} Nevertheless, \cref{thm:mserelevance} is
competitive with recent results. For instance, in nonparametric settings, the rate in
\cref{thm:mserelevance} is more favorable than the upper bound derived in
\citet{coeyhung}, who also study \cref{ex:topm}.
\end{rmksq}
\subsection{Validating performance by coupled bootstrap}
\label{sub:coupled_bootstrap}
We close this section with a procedure that provides unbiased estimates of the loss of
\emph{arbitrary} decision rules for \cref{ex:mse,ex:topm,ex:utilmax}. Practitioners can
use this procedure to evaluate the gain of \textsc{close}-\textsc{npmle}{} relative to other
alternatives---we do so extensively in \cref{sec:empirical}. The validity of this
validation depends only on the Gaussianity
\eqref{eq:gaussian_heteroskedastic_location}---without assuming $(\theta_i,
\sigma_i)$ are random nor assuming the location-scale model \eqref{eq:location_scale}.
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.
\]
\citet{oliveira2021unbiased} call $(Y_{i}^{(1)}, Y_{i}^{(2)})$ the \emph{coupled
bootstrap} draws. Observe that the two draws are conditionally independent under
\eqref{eq:gaussian_heteroskedastic_location}: \[
\colvecb{2}{Y_{i}^{(1)}}{Y_{i}^{(2)}} \mid \theta_i, \sigma_i^2 \sim \Norm\pr{
\colvecb{2}{\theta_{i}}{\theta_i}, \begin{bmatrix}
(1 + \omega) \sigma_i^2 & 0\\ 0 & (1+\omega^{-1}) \sigma_i^2
\end{bmatrix}
}. \addtocounter{equation}{1}\tag{\theequation} \label{eq:coupled_bootstrap}
\]
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, \eqref{eq:coupled_bootstrap} 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{\citet
{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. \Cref{prop:unbiased} can also be easily generalized to other loss functions that
admit unbiased estimators \citep[Effectively, the loss is a function of a Gaussian
location $\theta_i$. For unbiased estimation of functions of Gaussian parameters, see
Table A1 in][]{voinov2012unbiased}.}
\begin{table}[htb]
\caption{Unbiased estimators for loss of decision rules and associated conditional
variance expressions (\cref{prop:unbiased})}
\label{tab:unbiased_estimation}
\centering
\scriptsize
\begin{tabularx}{\textwidth}{m{0.2\textwidth}YY}
\toprule
Problem & Unbiased estimator of loss, $T\pr{Y_{1:n}^{(2)}, \bm{\delta}}$ & $\var\pr{T\pr{Y_
{1:n}^{(2)}, \bm{\delta}} \mid \mathcal F}$
\\ \midrule
\Cref{ex:mse} & $\frac{1}{n} \sum_{i=1}^n \pr{Y^{(2)}_i - \delta_i(Y_{1:n}^{(1)})}^2 -
\sigma_{i, (2)}^2$ & $\frac{1}{n^2}\sum_{i=1}^n \var\pr{(Y_i^{(2)} - \delta_i(Y_{1:n}^{
(1)}))^2
\mid
\mathcal F}$ \\
\Cref{ex:utilmax} & $-\frac{1}{n}\sum_{i=1}^n \delta_i(Y_{1:n}^{(1)}) Y^{(2)}_i $ & $
\frac{1}{n^2}
\sum_{i=1}^n \delta_i(Y_{1:n}^{(1)}) \sigma_{i, (2)}^2$ \\
\Cref{ex:topm} & $-\frac{1}{m}\sum_{i=1}^n \delta_i(Y_{1:n}^{(1)}) Y^{(2)}_i $ & $
\frac{1}{m^2}
\sum_{i=1}^n \delta_i(Y_{1:n}^{(1)}) \sigma_{i, (2)}^2$ \\ \bottomrule
\end{tabularx}
\end{table}
\begin{restatable}{prop}{unbiased}
\label{prop:unbiased}
Suppose $(Y_i,\sigma_i)$ obey \eqref{eq:gaussian_heteroskedastic_location}. Fix some
$\omega > 0$ and let $Y_{1:n}^{(1)}, Y_{1:n}^{(2)}$ be the coupled bootstrap draws. For
some decision problem, let $\bm{\delta}(Y_{1:n}^{(1)})$ be some decision rule using only data
$\pr{Y_ {i}^{(1)},
\sigma_{i, (1)}^2}_{i=1}^n$. Let $\mathcal F = \pr{
\theta_{1:n}, Y_{1:n}^{(1)}, \sigma_{1:n, (1)}, \sigma_{1:n, (2)}}$, for
\cref{ex:mse,ex:utilmax,ex:topm}, the estimators $T(Y_{1:n}^{(2)}, \bm{\delta})$ displayed
in \cref{tab:unbiased_estimation} are unbiased for the corresponding loss: \[
\E \bk{T(Y_{1:n}^{(2)}, \bm{\delta}(Y_{1:n}^{(1)})) \mid \mathcal F } = L\pr{\bm{\delta}(Y_{1:n}^{
(1)}), \theta_{1:n}}.
\]
Moreover, their conditional variances are equal to those displayed in
\cref{tab:unbiased_estimation}.
\end{restatable}
\Cref{prop:unbiased} 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 \eqref{eq:location_scale} or
if one is unwilling to evaluate risk integrating over random $\theta_i$.
\section{Empirical illustration}
\label{sec:empirical}
How does \textsc{close}-\textsc{npmle}{} perform in the field? We now consider two empirical exercises
related to \citet{chetty2018opportunity} and \citet{bergman2019creating}. Using
Census micro-data, \citet {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{
\label{fn:correlation}Like prior work that uses this data
\citep [see, e.g., footnote 28 in] [] {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, \citet{bergman2019creating} conducted a program
called {Creating Moves to Opportunity}.
\citet{bergman2019creating} provided assistance to treated low-income individuals to move
to Census tracts with estimated posterior means in the top third. We view \citet
{bergman2019creating}'s objectives as \textsc{top-{\small $m$} selection}, for $m$ equal to one third of the number of
tracts in Seattle and King County, WA.
The Opportunity Atlas published by \citet{chetty2018opportunity} also includes
tract-level covariates, a complication that we have so far abstracted away from. In the
ensuing empirical exercises, following \citet{bergman2019creating}, the estimates are
residualized against the covariates as a preprocessing
step \citep{fay1979estimates}.\footnote{\label{fn:covariate_additive}Alternatively, \cref{asub:covariate_additive}
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
\eqref{eq:location-scale-covariates}, 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
\citet{chetty2018opportunity} and cross-referenced with their labels in
\cref{fn:covariates}. 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 \emph{some} outcome for individuals of \emph{some} demographic
subgroup growing up in tract $i$, whose parents are at the 25\th{} income
percentile.\footnote{\label{fn:alpha}\Copy{fnalpha}{Since all measures of economic
mobility have bounded support, as
either percentile ranks or percentage rates, \cref{as:moments} is automatically satisfied
for $\theta_i$ with $\alpha = 2$, at least when there are no covariates.}} We will
consider three types of outcomes:
\begin{enumerate*}[label=(\roman*)]
\item percentile rank of adult income (\textsc{mean rank}),
\item an indicator for whether the individual has incomes in the top 20
percentiles (\textsc{top-{\small 20} probability}), and
\item an indicator for whether the individual is incarcerated (\textsc{incarceration})
\end{enumerate*}
for the following five demographic subgroups: all individuals (\textsc{pooled}), white
individuals, white men, Black individuals, and Black men. Under these shorthands, the
outcome in \cref{sec:model} is
\textsc{top-{\small 20} probability}{} (Black), while \citet{bergman2019creating} consider \textsc{mean rank}{}
\textsc{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 \citet{bergman2019creating}. It uses the coupled bootstrap
(\cref{sub:coupled_bootstrap}) to evaluate whether
\textsc{close}-\textsc{npmle}{} selects more economically mobile tracts than alternatives.
\subsection{Calibrated simulation}
\label{sub:calibrated}
\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 \textsc{npmle}. We then normalize this
estimated \textsc{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 \cref{asub:sim_setup}.
}
On the simulated data, we then implement various empirical Bayes strategies. We consider
the feasible procedures: \textsc{naive},
\textsc{independent-gauss}, \textsc{independent-npmle}, \textsc{close}-\textsc{gauss}{} (parametric),
\textsc{close}-\textsc{gauss}, and
\textsc{close}-\textsc{npmle}, as well as the infeasible \textsc{oracle}.
Here,
\begin{itemize}
\item \textsc{naive}{} sets $\hat\theta_i = Y_i$.
\item \textsc{independent-gauss}{} weighs the estimation of the hyperparameters $(m_0, s_0)$ with
$1/\sigma_i^2$, following \citet{bergman2019creating}.
\item \textsc{close}-\textsc{gauss}{} (parametric) implements
\textsc{close}-\textsc{gauss}, where \cref{item:close1} models the conditional moments parametrically as
$m_0 (\sigma_i; a) = a_1 + a_2 \log \sigma_i$
and $s_0^2(\sigma_i; b) = \exp(b_1 + b_2 \log \sigma_i)$, and estimates $m_0, s_0$ via
least-squares.\footnote{That is, we fit $a_1, a_2$ via minimizing $\sum_i (Y_i - a_1 - a_2
\log \sigma_i)^2$. We then fit $b_1, b_2$ via minimizing $
\sum_{i} \br{(Y_i - \hat m(\sigma_i))^2 - \sigma_i^2 - \exp(b_1 + b_2 \log (\sigma_i))}^2
$.
We thank an anonymous referee for this suggestion.
}
\item The conditional moments $\eta_0 = (m_0(\cdot), s_0(\cdot))$ in \textsc{close}-\textsc{gauss}{} and
\textsc{close}-\textsc{npmle}{} are estimated via local linear
regression, where bandwidth is selected via plug-in IMSE-optimal bandwidth, as
implemented in
\citet{calonico2019nprobust}.\footnote{\label{fn:implementation}Specifically, $\hat m = \hat \E
[Y_i \mid \log\sigma_i]$ and $\hat s^2(\sigma_i) = \max(
\hat \E[(Y_i - \hat m(\sigma_i))^2 \mid
\log
\sigma_i] - \sigma_i^2, \tilde{s}^2(\sigma_i)),$ where $\hat \E[\cdot \mid \log
\sigma_i]$ implements local linear regression and $\tilde{s}(\sigma_i)$ implements a
data-driven truncation of $\hat s^2$, detailed in \cref
{sec:nuisance_estimation}. Replacing the truncation point $\tilde s(\sigma_i)$ with zero
(that is, we exclude the observations with $\hat s(\sigma_i) = 0$ from estimating $\hat
G_n$, and treat these observations as having empirical Bayes posterior degenerate at
$\hat m(\sigma_i)$) does not appear to qualitatively affect our results. }
\item Since we know
the ground truth data-generating process, we can also compute the
\textsc{oracle}{} procedure that uses posterior means under the true $P_0$.
\item None of the feasible procedures have access to $\beta$, which they must estimate
in the same way using weighted least squares with weight $1/\sigma_i^2$, following
\citet{bergman2019creating}.
\end{itemize}
\begin{figure}[htb]
\centering
\includegraphics[width=\textwidth]{final_assets/mse_table_calibrated.pdf}
\begin{proof}[Notes]
Each column is an empirical Bayes strategy that we consider, and each row is a different
definition of $\vartheta_i$. The table shows relative performance, defined as the
squared error improvement over \textsc{naive}, normalized as a percentage of the improvement of
\textsc{oracle}{} over \textsc{naive}. The last row shows the column median. Results are averaged over
1,000 Monte Carlo draws.
\end{proof}
\caption{Relative squared error Bayes risk for various empirical Bayes posterior means}
\label{fig:mse_table}
\end{figure}
\Cref{fig:mse_table} 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 \textsc{naive}{}, normalized by the
MSE gain of \textsc{oracle}{} over \textsc{naive}. If we think of the
\textsc{oracle}--\textsc{naive}{} difference as the total size of the ``statistical pie,'' then
\cref{fig:mse_table} 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
\textsc{close}.\footnote{\label
{foot:weighted_means}It may be surprising that \textsc{independent-gauss}{} can perform worse than
\textsc{naive}{} even on MSE, since Gaussian empirical Bayes can be thought of as optimizing
among a class of linear shrinkage estimators that include
\textsc{naive}. We note that, as in
\citet{bergman2019creating}, when we estimate the prior mean and prior variance, we
\emph{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
\textsc{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 \textsc{mean rank}, after covariate residualization, the dependence between
$\theta_i$ and $\sigma_i$ does not appear to substantially affect shrinkage decisions.
\textsc{independent-npmle}{} and
\textsc{close}-methods perform similarly, capturing almost all of the available gains. For the other two outcome variables, \textsc{top-{\small 20} probability} {} and
\textsc{incarceration}, the dependence between $\theta_i $ and $
\sigma_i$ is stronger, and \textsc{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{\Cref
{asub:weibull} 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.}
\subsection{Validation exercise via coupled bootstrap}
\begin{figure}[!htb]
\centering
\includegraphics[width=\textwidth]{final_assets/rank_table_additional.pdf}
\begin{proof}[Notes]
Each column is an empirical Bayes strategy that we consider, and each row is a different
definition of $\vartheta_i$. The table shows coupled-bootstrap estimates of average
$\vartheta_i$ among the Census tracts selected by each method---in terms of either
percentage points or percentile ranks---over 1,000 draws of coupled bootstrap. All
decision rules are estimated separately within CZs and select the top third of Census
tracts within each CZ. The color scheme within each row treats the performance of
\textsc{naive}{} as zero (grey) and \textsc{close}-\textsc{npmle}{} as one (dark green), and is hence not
comparable across rows. On this scale, a method that overperforms \textsc{naive}{} is colored
green; otherwise it is colored magenta. The best performer for each row is additionally
marked with an orange star.
\end{proof}
\caption{Performance of decision rules in top-$m$ selection exercise}
\label{tab:selection_validation_within_cz}
\end{figure}
Our second empirical exercise uses the coupled bootstrap described in
\cref{sub:coupled_bootstrap} for the policy problem in \citet {bergman2019creating}.
Viewing the policy problem in \citet{bergman2019creating} as \textsc{top-{\small $m$} selection}, can \textsc{close}-\textsc{npmle}{} make
better selections?
Specifically, we imagine scaling up \citet{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 \emph{within} each CZ, according to
empirical Bayesian posterior means for $\vartheta_i$. Additionally, to faithfully mimic
\citet{bergman2019creating}, here we perform all empirical Bayes procedures \emph
{within CZ}.
Throughout, we choose $\omega$ to emulate a 90-10 train-test split on the micro-data. See
\cref{asub:sim_setup} for details on the policy exercise setup.
\Cref{tab:selection_validation_within_cz} shows the estimated performance of various
methods.
According to these estimates, \textsc{close}-\textsc{npmle} {} generally improves over
\textsc{independent-gauss}.\footnote{\textsc{close}-\textsc{npmle}{} is worse by an estimated
$0.006$ percentile ranks for \textsc{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 \textsc{mean rank}{} variables, using \textsc{close}-\textsc{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 \citet{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, \textsc{top-{\small 20} probability}{} and \textsc{incarceration}, the gains are even more
sizable.
These gains are as high as 2--3 percentage points on average. Among \textsc{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 \textsc{independent-gauss}{} and \textsc{naive} {} as the
\emph{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 \textsc{close}-\textsc{npmle}{} over
\textsc{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.
\section{Conclusion}
\label{sec:conclusion}
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 \textsc{close}. The \textsc{close}-framework naturally subsumes and generalizes several
existing proposals for accommodating precision dependence. We prove that
\textsc{close}-\textsc{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
\textsc{close}-\textsc{npmle}{} is robust, with finite worst-case Bayes risk. Finally, we further connect
our main theoretical
results to ranking-type decision problems in
\citet{bergman2019creating}.
Simulation and validation exercises demonstrate that \textsc{close}-\textsc{npmle}{} generates sizable gains
relative to the standard parametric empirical Bayes shrinkage method. Across calibrated
simulations, \textsc{close}-\textsc{npmle}{} attains close-to-oracle mean-squared error performance. In a
hypothetical, scaled-up version of \citet{bergman2019creating}, across a wide range of
economic mobility measures, \textsc{close}-\textsc{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.