Robust Empirical Bayes Confidence Intervals

Robust Empirical Bayes Confidence Intervals

Introduction

Empirical researchers in economics are often interested in estimating effects for many individuals or units, such as estimating teacher quality for teachers in a given geographic area. In such problems, it is common to shrink unbiased but noisy preliminary estimates of these effects toward baseline values, say the average effect for teachers with the same experience. In addition to estimating teacher quality KaSt08,JaLe08,cfr14i, shrinkage techniques are used in a wide range of applications including estimating school quality ahpw17, hospital quality hull20, the effects of neighborhoods on intergenerational mobility chetty_impacts_2018, and patient risk scores across regional health care markets fghw17.

The shrinkage estimators used in these applications can be motivated by an \ac{EB} approach. One imposes a working assumption that the individual effects are drawn from a normal distribution (or, more generally, a known family of distributions). The \ac{MSE} optimal point estimator then has the form of a Bayesian posterior mean, treating this distribution as a prior distribution. Rather than specifying the unknown parameters in the prior distribution ex ante, the \ac{EB} estimator replaces them with consistent estimates, just as in random effects models. This approach is attractive because one does not need to assume that the effects are in fact normally distributed, or even take a “Bayesian” or “random effects” view: the \ac{EB} estimators have lower \ac{MSE} (averaged across units) than the unshrunk unbiased estimators, even when the individual effects are treated as nonrandom JaSt61.

In spite of the popularity of \ac{EB} methods, it is currently not known how to provide uncertainty assessments to accompany the point estimates without imposing strong parametric assumptions on the effect distribution. Indeed, Hansen2016 describes inference in shrinkage settings as an open problem in econometrics. The natural \ac{EB} version of a \ac{CI} takes the form of a Bayesian credible interval, again using the postulated effect distribution as a prior morris83. If the distribution is correctly specified, this parametric \acf{EBCI} will cover 95%, say, of the true effect parameters, under repeated sampling of the observed data and of the effect parameters. We refer to this notion of coverage as “\ac{EB} coverage”, following the terminology in morris83. Unfortunately, we show that, in the context of a normal means model, the parametric \ac{EBCI} with nominal level 95% can have actual \ac{EB} coverage as low as 74% for certain non-normal effect distributions. The potential undercoverage is increasing in the degree of shrinkage, and we derive a simple “rule of thumb” for gauging the potential coverage distortion.

To allow easy uncertainty assessment in \ac{EB} applications that is reliable irrespective of the degree of shrinkage, we construct novel robust \acp{EBCI} that take a simple form and control \ac{EB} coverage regardless of the true effect distribution. Our baseline model is an (approximate) normal means problem $Y_i \sim N(\theta_{i}, \sigma_i^2)$, $i=1, \dotsc, n$. In applications, $Y_i$ represents a preliminary estimate of the effect $\theta_i$ for unit $i$. Like the parametric \ac{EBCI} that assumes a normal distribution for $\theta_i$, the robust \ac{EBCI} we propose is centered at the normality-based \ac{EB} point estimate $\hat{\theta}_i$ that shrinks $Y_{i}$ toward some baseline value, but it uses a larger critical value to account for bias due to shrinkage.\footnote{ Our methods are implemented in the Stata package ebreg, R package ebci, and Matlab package ebci_matlab, which are available at SSC, CRAN, and GitHub, respectively.} \ac{EB} coverage is controlled in the class of all distributions for $\theta_{i}$ that satisfy certain moment bounds, which we estimate consistently from the data (similarly to the parametric \ac{EBCI}, which uses the second moment). We show that the baseline implementation of our robust \ac{EBCI} is “adaptive”: its length is close to that of the parametric \ac{EBCI} when the $\theta_{i}$'s are in fact normally distributed. Thus, little efficiency is lost from using the robust \ac{EBCI} in place of the non-robust parametric one.\footnote{ If the $\theta_{i}$'s are not normally distributed, our robust \acp{EBCI} are valid but may leave room for greater efficiency improvement, as we discuss in (ref).}

In addition to controlling \ac{EB} coverage, the robust \acp{EBCI} with level $1-\alpha$ have a frequentist average coverage property: If the means $\theta_{1}, \dotsc, \theta_{n}$ are treated as fixed, the coverage probability, averaged across the $n$ parameters $\theta_i$, is at least $1-\alpha$. In fact, under mild conditions, at least a fraction $1-\alpha$ of the $n$ \acp{EBCI} will contain their respective parameters (with high probability as $n \to\infty$). This weakening of the usual requirement of coverage for each parameter $\theta_i$ allows our robust \ac{EBCI} to be shorter than the usual \ac{CI} centered at the unshrunk estimate $Y_{i}$, and often substantially so.\footnote{Relaxing the usual notion of coverage in some way is necessary to obtain intervals that reflect the efficiency improvement of the empirical Bayes approach. In particular, the results in pratt61 imply that for \acp{CI} with coverage 95%, one cannot achieve expected length improvements greater than 15% relative to the usual unshrunk \acp{CI}, even if one happens to optimize length for the true parameter vector $(\theta_{1}, \dotsc, \theta_{n})$. See, for example, Corollary 3.3 in armstrong_optimal_2018 and the discussion following it.} Intuitively, the average coverage criterion only requires us to guard against the average coverage distortion induced by the biases of the individual shrinkage estimators $\hat{\theta}_i$, and the data is quite informative about whether most of these biases are large, even though individual biases are difficult to estimate. To complement the frequentist properties, our \acp{EBCI} can be viewed as Bayesian credible sets that are robust to the prior on $\theta_i$, in terms of ex ante coverage.

The average coverage criterion has the same motivation as the usual frequentist justification of the \ac{EB} point estimator: the \ac{EB} point estimator achieves lower \ac{MSE} on average across units at the expense of potentially worse performance for some individual units Efron2010. Thus, researchers who use \ac{EB} estimators instead of the unshrunk $Y_{i}$'s prioritize favorable group performance over protecting individual performance. Our average coverage intervals make an analogous tradeoff: they guarantee coverage and achieve short length on average across units at the expense of giving up on a coverage guarantee for every individual unit. We examine this tradeoff in more detail in (ref).

We caution, however, that the average coverage criterion is typically inappropriate in applications where shrinkage point estimation is unattractive. This includes settings where one is interested in the magnitude or the identity of the largest $\theta_{i}$, or the true effect for the largest observed $Y_{i}$ (as in, for example, HuFi19, or Andrews2021).\footnote{As we show in (ref), our methods do extend to settings where we keep a subset of units $i$ that exceed a given cutoff. However, we do not allow this cutoff to diverge with the sample size, such as when one focuses on the unit $i$ with the single largest observed $Y_{i}$.} It also includes settings where a particular effect, say $\theta_{1}$, is of primary interest, or, more generally, settings where the effects are not exchangeable, and their ordering is relevant GrRi19. Our methods are also not applicable if one is interested in functionals of the random effects distribution (as in bonhomme2020, or ignatiadis2019), rather than in the effects themselves. Finally, the justification for our methods is asymptotic in the number of parameters $n$. In our Monte Carlo simulations, we find that our \acp{EBCI} have close to nominal coverage over a range of \acp{DGP} once $n$ is greater than $100$.

We illustrate our results by computing \acp{EBCI} for the causal effects of growing up in different U.S. neighborhoods (specifically commuting zones) on intergenerational mobility. We follow chetty_impacts_2018, who apply \ac{EB} shrinkage to initial fixed effects estimates. Depending on the specification, we find that the robust \acp{EBCI} are on average 12--25% as long as the unshrunk \acp{CI}.

Our underlying ideas extend to other linear and non-linear shrinkage settings with possibly non-Gaussian data. For example, our techniques allow for the construction of robust \acp{EBCI} that contain (nonlinear) soft thresholding estimators, as well as average coverage confidence bands for nonparametric regression functions.

The average coverage criterion was originally introduced in the literature on nonparametric regression (Wahba1983; Nychka1988; Wasserman2006). Cai2015 construct adaptive average coverage confidence bands. These procedures are challenging to implement in our \ac{EB} setting, and lack a clear finite-sample justification, unlike our procedure. Liu2019 construct forecast intervals in a dynamic panel data model that guarantee average coverage in a Bayesian sense (for a fixed prior). We discuss alternative approaches to inference in \ac{EB} settings in (ref).

The rest of this paper is organized as follows. (ref) illustrates our methods in the context of a simple homoskedastic Gaussian model. (ref) presents our recommended baseline procedure and discusses practical implementation issues. (ref) presents our main results on the coverage and efficiency of the robust \ac{EBCI}, and on the coverage distortions of the parametric \ac{EBCI}; we also verify the finite-sample coverage accuracy of the robust \ac{EBCI} through extensive simulations. (ref) compares our \ac{EBCI} with other inference approaches. (ref) discusses extensions of the basic framework. (ref) contains an empirical application to inference on neighborhood effects. (ref) give details on finite-sample corrections, computational details, and formal asymptotic coverage results. The Online Supplement contains proofs as well as further technical results. Applied readers are encouraged to focus on (ref).

Simple example

This section illustrates the construction of the robust \acp{EBCI} that we propose in a simplified setting with no covariates and with known, homoskedastic errors. (ref) relaxes these restrictions, and discusses other empirically relevant extensions of the basic framework, as well as implementation issues.

We observe $n$ estimates $Y_{i}$ of elements of the parameter vector $\theta=(\theta_{1}, \dotsc, \theta_{n})'$. Each estimate is normally distributed with common, known variance $\sigma^{2}$,

In many applications, the $Y_{i}$'s arise as preliminary least squares estimates of the parameters $\theta_{i}$. For instance, they may correspond to fixed effect estimates of teacher or school value added, neighborhood effects, or firm and worker effects. In such cases, $Y_{i}$ will only be approximately normal in large samples by the \ac{CLT}; we take this explicitly into account in the theory in (ref).

A popular approach to estimation that substantially improves upon the raw estimator $\hat{\theta}_{i}=Y_{i}$ under the compound \ac{MSE} $\sum_{i=1}^{n}E[(\hat{\theta}_{i}-\theta_{i})^{2}]$ is based on \acf{EB} shrinkage. In particular, suppose that the $\theta_{i}$'s are themselves normally distributed,

Our discussion below applies if (ref) is viewed as a subjective Bayesian prior distribution for a single parameter $\theta_i$, but for concreteness we will think of (ref) as a “random effects” sampling distribution for the $n$ mean parameters $\theta_1, \dotsc, \theta_n$. Under (ref), it is optimal to estimate $\theta_i$ using the posterior mean $\hat{\theta}_i=w_{EB}Y_i$, where $w_{EB}=1-\sigma^{2}/(\sigma^{2}+\mu_{2})$. To avoid having to specify the variance $\mu_{2}$, the \ac{EB} approach treats it as an unknown parameter, and replaces the marginal precision of $Y_{i}$, $1/(\sigma^{2}+\mu_{2})$, with a method of moments estimate $n/\sum_{i=1}^{n}Y_{i}^{2}$, or the degrees-of-freedom adjusted estimate $(n-2)/\sum_{i=1}^{n}Y_{i}^{2}$. The latter leads to the classic estimator of JaSt61, $\hat{w}_{EB}=1-\sigma^2(n-2)/\sum_{i=1}^{n}Y_{i}^{2}$.

One can also use (ref) to construct \acp{CI} for the $\theta_{i}$'s. In particular, since the marginal distribution of $w_{EB}Y_{i}-\theta_{i}$ is normal with mean zero and variance $(1-w_{EB})^{2}\mu_{2}+w_{EB}^{2}\sigma^{2}=w_{EB}\sigma^{2}$, this leads to the $1-\alpha$ \ac{CI}

where $z_{\alpha}$ is the $\alpha$ quantile of the standard normal distribution. Since the form of the interval is motivated by the parametric assumption (ref), we refer to it as a parametric \ac{EBCI}. With $\mu_{2}$ unknown, one can replace $w_{EB}$ by $\hat{w}_{EB}$.\footnote{Alternatively, to account for estimation error in $\hat{w}_{EB}$, morris83 suggests adjusting the variance estimate $\hat{w}_{EB}\sigma^{2}$ to $\hat{w}_{EB}\sigma^{2}+2Y_{i}^{2}(1-\hat{w}_{EB})^{2}/(n-2)$. The adjustment does not matter asymptotically.} This is asymptotically equivalent to (ref) as $n\to\infty$.

The coverage of the parametric \ac{EBCI} in (ref) is $1-\alpha$ under repeated sampling of $(Y_{i}, \theta_{i})$ according to (ref). To distinguish this notion of coverage from the case with fixed $\theta$, we refer to coverage under repeated sampling of $(Y_{i}, \theta_{i})$ as “empirical Bayes coverage”. This follows the definition of an \acf{EBCI} in morris83 and CaLo00. Unfortunately, this coverage property relies heavily on the parametric assumption (ref). We show in (ref) that the actual \ac{EB} coverage of the nominal $1-\alpha$ parametric \ac{EBCI} can be as low as $1-1/\max\{z_{1-\alpha/2},1\}$ for certain non-normal distributions of $\theta_{i}$ with variance $\mu_{2}$; for 95% \acp{EBCI}, this evaluates to 74%. This contrasts with existing results on estimation: although the \ac{EB} estimator is motivated by the parametric assumption (ref), it performs well even if this assumption is dropped, with low \ac{MSE} even if we treat $\theta$ as fixed.

This paper constructs an \ac{EBCI} with a similar robustness property: the interval will be close in length to the parametric \ac{EBCI} when (ref) holds, but its \ac{EB} coverage is at least $1-\alpha$ without any parametric assumptions on the distribution of $\theta_{i}$. To describe the construction, suppose that all that is known is that $\theta_{i}$ is sampled from a distribution with second moment given by $\mu_{2}$ (in practice, we can replace $\mu_2$ by the consistent estimate $n^{-1}\sum_{i=1}^{n}Y^{2}_{i}-\sigma^{2}$). Conditional on $\theta_{i}$, the estimator $w_{EB}Y_{i}$ has bias $(w_{EB}-1)\theta_{i}$ and variance $w_{EB}^{2}\sigma^{2}$, so that the $t$-statistic $(w_{EB}Y_{i}-\theta_{i})/w_{EB}\sigma$ is normally distributed with mean $b_{i}=(1-1/w_{EB})\theta_{i}/\sigma$ and variance $1$. Therefore, if we use a critical value $\chi$, the non-coverage of the \ac{CI} $w_{EB}Y_{i}\pm \chi w_{EB}\sigma$, conditional on $\theta_{i}$, will be given by the probability

where $Z$ denotes a standard normal random variable, and $\Phi$ denotes its cdf. Thus, by iterated expectations, under repeated sampling of $\theta_{i}$, the non-coverage is bounded by

where $E_{F}$ denotes expectation under $b\sim F$. Although this is an infinite-dimensional optimization problem over the space of distributions, it turns out that it admits a simple closed-form solution, which we give in (ref) in (ref). Moreover, because the optimization is a linear program, it can be solved even in the more general settings of applied relevance that we consider in (ref).

Set $\chi=\operatorname{cva}_{\alpha}(\sigma^{2}/\mu_{2})$, where $\operatorname{cva}_{\alpha}(t)=\rho^{-1}(t, \alpha)$, and the inverse is with respect to the second argument. Then the resulting interval

will maintain coverage $1-\alpha$ among all distributions of $\theta_{i}$ with $E[\theta_{i}^{2}]=\mu_{2}$ (recall that we estimate $\mu_{2}$ consistently from the data). For this reason, we refer to it as a robust \ac{EBCI}. (ref) in (ref) gives a plot of the critical values for $\alpha=0.05$. We show in (ref) below that by also imposing a constraint on the fourth moment of $\theta_{i}$, in addition to the second moment constraint, one can construct a robust \ac{EBCI} that “adapts” to the Gaussian case in the sense that its length will be close to that of the parametric \ac{EBCI} in (ref) if these moment constraints are compatible with a normal distribution.

Instead of considering \ac{EB} coverage, one may alternatively wish to assess uncertainty associated with the estimates $\hat{\theta}_i=w_{EB}Y_{i}$ when $\theta$ is treated as fixed. In this case, the \ac{EBCI} in (ref) has an average coverage guarantee that

provided that the moment constraint can be interpreted as a constraint on the empirical second moment on the $\theta_{i}$'s, $n^{-1}\sum_{i=1}^{n}\theta_{i}^{2}=\mu_{2}$. In other words, if we condition on $\theta$, then the coverage is at least $1-\alpha$ on average across the $n$ \acp{EBCI} for $\theta_{1}, \dotsc, \theta_{n}$. To see this, note that the average non-coverage of the intervals is bounded by (ref), except that the supremum is only taken over possible empirical distributions for $\theta_{1}, \dotsc, \theta_{n}$ satisfying the moment constraint. Since this supremum is necessarily smaller than $\rho(\sigma^{2}/\mu_{2}, \chi)$, it follows that the average coverage is at least $1-\alpha$. In fact, if the $Y_i$'s exhibit limited dependence across $i$, a stronger property holds: the probability that at least a fraction $1-\alpha$ of the $n$ \acp{EBCI} contain their respective true parameters converges to 1 as $n \to \infty$, cf.\ (ref) below.

The usual \acp{CI} $Y_{i}\pm z_{1-\alpha/2}\sigma$ also of course achieve average coverage $1-\alpha$. The robust \ac{EBCI} in (ref) will, however, be shorter, especially when $\mu_{2}$ is small relative to $\sigma^{2}$---see (ref) below. The reduction in length is achieved by weakening the requirement that each \ac{CI} covers its true parameter $1-\alpha$ percent of the time to the requirement that the coverage probability equal $1-\alpha$ on average across the \acp{CI}. It may seem surprising that we can construct a narrower \ac{CI} by centering it at the shrinkage estimates $w_{EB}Y_{i}$. The intuition for this is that the shrinkage reduces the variability of the estimates, at the expense of introducing bias in the estimates. The bias necessitates the use of a larger critical value $\operatorname{cva}_{\alpha}(\sigma^{2}/\mu_{2})$. Because under the average coverage criterion we only need to control the bias on average across $i$, rather than for each individual $\theta_{i}$, this increase in the critical value is smaller than the reduction in the standard error.

Practical implementation

We now describe how to compute a robust \ac{EBCI} that allows for heteroskedasticity, shrinks towards more general regression estimates rather than towards zero, and exploits higher moments of the bias to yield a narrower interval. In (ref), we describe the empirical Bayes model that motivates our baseline approach. (ref) describes the practical implementation of our baseline approach.

Motivating model and robust EBCI

In applied settings, the unshrunk estimates $Y_i$ will typically have heteroskedastic variances. Furthermore, rather than shrinking towards zero, it is common to shrink toward an estimate of $\theta_i$ based on some covariates $X_i$, such as a regression estimate $X_i'\hat{\delta}$. We now describe how to adapt the ideas in (ref) to such settings.

Consider a generalization of the model in (ref) that allows for heteroskedasticity and covariates,

The covariate vector $X_i$ may contain just the intercept, and it may also contain (functions of) $\sigma_{i}$. To construct an \ac{EB} estimator of $\theta_{i}$, consider the working assumption that the sampling distribution of the $\theta_i$'s is conditionally normal:

The hierarchical model (ref)--(ref) leads to the Bayes estimate $\hat{\theta}_{i}=\mu_{1, i}+w_{EB, i}(Y_i-\mu_{1, i})$, where $w_{EB, i}=\frac{\mu_{2}}{\mu_{2}+\sigma_{i}^{2}}$. This estimate shrinks the unrestricted estimate $Y_{i}$ of $\theta_{i}$ toward $\mu_{1,i}=X_{i}'\delta$. In contrast to (ref), the normality assumption (ref) typically cannot be justified simply by appealing to the \ac{CLT}; the linearity of the conditional mean $\mu_{1,i}=X_i'\delta$ may also be suspect. Our robust \ac{EBCI} will therefore be constructed so that it achieves valid \ac{EB} coverage even if assumption (ref) fails. To obtain a narrow robust \ac{EBCI}, we augment the second moment restriction used to compute the critical value in (ref) with restrictions on higher moments of the bias of $\hat{\theta}_{i}$. In our baseline specification, we add a restriction on the fourth moment.

In particular, we replace assumption (ref) with the much weaker requirement that the conditional second moment and kurtosis of $\varepsilon_i=\theta_{i}-X_{i}'\delta$ do not depend on $(X_{i}, \sigma_{i})$:

where $\delta$ is defined as the probability limit of the regression estimate $\hat{\delta}$.\footnote{Our framework can be modified to let $(X_i, \sigma_i)$ be fixed, in which case $\delta$ depends on $n$. See the discussion following (ref) below.} We discuss this requirement further in (ref) below, and we relax it in (ref) below.

We now apply analysis analogous to that in (ref). Let us suppose for simplicity that $\delta$, $\mu_{2}$, $\kappa$, and $\sigma_{i}$ are known; we relax this assumption in (ref) below, and in the theory in (ref). Denote the conditional bias of $\hat{\theta}_i$ normalized by the standard error by $b_{i} = (w_{EB, i}-1)\varepsilon_i/(w_{EB, i}\sigma_i) = -\sigma_{i}\varepsilon_{i}/\mu_{2}$. Under repeated sampling of $\theta_i$, the non-coverage of the \ac{CI} $\hat{\theta}_i \pm \chi w_{EB, i}\sigma$, conditional on $(X_{i}, \sigma_{i})$, depends on the distribution of the normalized bias $b_{i}$, as in (ref). Given the moments $\mu_{2}$ and $\kappa$, the maximal non-coverage is given by

where $b$ is distributed according to the distribution $F$. Here $m_{2,i} = E[b_{i}^{2} \mid X_{i}, \sigma_{i}] =\sigma^{2}_{i}/\mu_{2}$. Observe that the kurtosis of $b_{i}$ matches that of $\varepsilon_{i}$. (ref) shows that the infinite-dimensional linear program (ref) can be reduced to two nested univariate optimizations. We also show that the least favorable distribution---the distribution $F$ maximizing (ref)---is a discrete distribution with up to 4 support points (see (ref)).

Define the critical value $\operatorname{cva}_{\alpha}(m_{2,i}, \kappa)=\rho^{-1}(m_{2,i}, \kappa, \alpha)$, where the inverse is in the last argument. (ref) plots this function for $\alpha=0.05$ and selected values of $\kappa$. This leads to the robust \ac{EBCI}

which, by construction, has coverage at least $1-\alpha$ under repeated sampling of $(Y_i, \theta_i)$, conditional on $(X_{i}, \sigma_{i})$, so long as (ref) holds; it is not required that (ref) holds. Note that both the critical value and the \ac{CI} length are increasing in $\sigma_{i}$.

Baseline implementation

Our baseline implementation of the robust \ac{EBCI} plugs in consistent estimates of the unknown quantities in (ref), based on the data $\{Y_{i}, X_{i}, \hat{\sigma}_{i}\}_{i=1}^{n}$, where $\hat{\sigma}_{i}$ is a consistent estimate of $\sigma_{i}$ (such as the standard error of the preliminary estimate $Y_{i}$), and $X_{i}$ is a vector of covariates that are thought to help predict $\theta_{i}$.

We provide fast and stable software packages that automate these steps (see (ref)). We now discuss the assumptions needed for validity of the robust \ac{EBCI}.

Main results

This section provides formal statements of the coverage properties of the \acp{CI} presented in (ref). Furthermore, we show that the \acp{CI} presented in (ref) are highly efficient when the mean parameters are in fact normally distributed. Next, we calculate the maximal coverage distortion of the parametric \ac{EBCI}, and derive a rule of thumb for gauging the potential coverage distortion. Finally, we present a comprehensive simulation study of the finite-sample performance of the robust \ac{EBCI}. Applied readers interested primarily in implementation issues may skip ahead to the empirical application in (ref).

Coverage under baseline implementation

In order to state the formal result, let us first carefully define the notions of coverage that we consider. Consider intervals $CI_{1}, \dotsc, CI_{n}$ for elements of the parameter vector $\theta=(\theta_1, \dotsc, \theta_n)'$. The probability measure $P$ denotes the joint distribution of $\theta$ and $CI_1, \dotsc, CI_n$. Following morris83 and CaLo00, we say that the interval $CI_i$ is an (asymptotic) $1-\alpha$ \acf{EBCI} if

We say that the intervals $CI_i$ are (asymptotic) $1-\alpha$ \acp{ACI} under the parameter sequence $\theta_1, \dotsc, \theta_n$ if

The average coverage property (ref) is a property of the distribution of the data conditional on $\theta$ and therefore does not require that we view the $\theta_i$'s as random (as in a Bayesian or “random effects” analysis). To maintain consistent notation, we nonetheless use the conditional notation $P(\cdot\mid \theta)$ when considering average coverage. See (ref) for a formulation with $\theta$ treated as nonrandom.

Observe that under the exchangeability condition that $P(\theta_i\in CI_i)=P(\theta_j\in CI_j)$ for all $i, j$, if the \ac{ACI} property (ref) holds almost surely, then the \ac{EBCI} property (ref) holds, since then

We now provide coverage results for the baseline implementation described in (ref). To keep the statements in the main text as simple as possible, we (i) maintain the assumption that the unshrunk estimates $Y_i$ follow an exact normal distribution conditional on the parameter $\theta_i$, (ii) state the results only for the homoskedastic case where the variance $\sigma_i$ of the unshrunk estimate $Y_i$ does not vary across $i$, and (iii) consider only unconditional coverage statements of the form (ref) and (ref). In (ref), we allow the estimates $Y_i$ to be only approximately normally distributed and allow $\sigma_i$ to vary, and we verify that our assumptions hold in a linear fixed effects panel data model. We also formalize the statements about conditional coverage made in (ref).

(ref) follows immediately from (ref) in (ref). In order to cover both the \ac{EB} coverage condition (ref) and the average coverage condition (ref), (ref) considers a random sequence of parameters $\theta_1, \dotsc, \theta_n$, and shows average coverage conditional on these parameters. See (ref) for a formulation with $\theta$ treated as nonrandom.

The condition on the moments $\mu_2$ and $\kappa$ avoids degenerate cases such as when $\mu_{2}=0$, in which case the \ac{EB} point estimator $\hat{\theta}_i$ shrinks each preliminary estimate $Y_i$ all the way to $X_i'\hat{\delta}$. Note also that the theorem does not require that $\hat{\delta}$ be the \ac{OLS} estimate in a regression of $Y_{i}$ onto $X_{i}$, and that $\delta$ be the population analog; one can define $\delta$ in other ways, the theorem only requires that $\hat{\delta}$ be a consistent estimate of it. The definition of $\delta$ does, however, affect the plausibility of the moment independence assumption in (ref) needed for conditional coverage results stated in (ref).\footnote{The specification of $\mu_{1i}=X_{i}'\delta$ also affects the \ac{EBCI} width through its effect on $\mu_{2}$ and $\kappa$.}

Relative efficiency

The robust \ac{EBCI} in (ref), unlike the parametric \ac{EBCI} $\hat{\theta}_{i}\pm z_{1-\alpha/2}\sigma_{i}\sqrt{w_{EB, i}}$, does not rely on the normality assumption in (ref) for its validity. We now show that this robustness does not come at a high cost in terms of efficiency: if the normality assumption (ref) in fact holds, the efficiency loss is limited unless the signal-to-noise ratio $\mu_{2}/\sigma^{2}_{i}$ is very small.

There are two reasons for the inefficiency of the robust \ac{EBCI}. First, the robust \ac{EBCI} only makes use of the second and fourth moment of the conditional distribution of $\theta_{i}-X_{i}'\delta$, rather than its full distribution. Second, if we only have knowledge of these two moments, it is no longer optimal to center the \ac{EBCI} at the estimator $\hat{\theta}_{i}$: one may need to consider other, perhaps non-linear, shrinkage estimators, as we do below in (ref).

We decompose the sources of inefficiency by studying the relative length of the robust \ac{EBCI} relative to the \ac{EBCI} that picks the amount of shrinkage optimally. For the latter, we maintain assumption (ref), and consider a more general class of estimators $\tilde{\theta}(w_{i})=\mu_{1, i}+w_{i}(Y_{i}-\mu_{1, i})$. For tractability, we focus on fixed-length \acp{CI} based on linear shrinkage estimators, but allow the amount of shrinkage $w_{i}$ to be optimally determined. The normalized bias of $\tilde{\theta}(w_{i})$ is given by $b_{i}=(1/w_{i}-1)\varepsilon_{i}/\sigma_{i}$, which leads to the \ac{EBCI}

The half-length of this \ac{EBCI}, $\operatorname{cva}_{\alpha}((1-1/w_{i})^{2}\mu_{2}/\sigma^{2}_{i}, \kappa)w_{i}\sigma_{i}$, can be numerically minimized as a function of $w_{i}$ to find the \ac{EBCI} length-optimal shrinkage. Denote the minimizer by $w_{opt}(\mu_{2}/\sigma^{2}_{i}, \kappa, \alpha)$. Like $w_{EB, i}$, the optimal shrinkage depends on $\mu_{2}$ and $\sigma_{i}^{2}$ only through the signal-to-noise ratio $\mu_{2}/\sigma^{2}_{i}$. Numerically evaluating the minimizer shows that $w_{opt}(\cdot, \kappa, \alpha)\geq w_{EB, i}$ for $\kappa\geq 3$ and $\alpha\in\{0.05,0.1\}$. The resulting \ac{EBCI} is optimal among all fixed-length \acp{EBCI} centered at linear estimators under (ref), and we call it the optimal robust \ac{EBCI}.\footnote{Since the optimal robust \ac{EBCI} is always shorter than the robust \ac{EBCI} in (ref), the former is preferable on efficiency grounds. It may not contain the \ac{MSE}-optimal point estimator $\hat{\theta}_i$, however.}

(ref) plots the ratio of lengths of the optimal robust \ac{EBCI} and robust \ac{EBCI} relative to the parametric \ac{EBCI}, for $\alpha=0.05$. The figure shows that to maintain efficiency relative to the normal benchmark, it is important to impose the fourth moment constraint. If this constraint is imposed, the efficiency loss of the robust \ac{EBCI} is modest unless the signal-to-noise ratio is very small: if $w_{EB, i}\geq 0.1$ (which is equivalent to $\mu_{2}/\sigma^{2}_{i}\geq 1/9$), the efficiency loss is at most $11.4\%$ for $\alpha=0.05$; up to half of the efficiency loss is due to not using the optimal shrinkage. For $\alpha=0.1$ (not plotted), the results are very similar; in particular, if $w_{EB, i}\geq 0.1$, the efficiency loss is at most $12.9\%$.

When the signal-to-noise ratio is very small, so that $w_{EB, i}<0.1$, the efficiency loss of the robust \ac{EBCI} is higher (up to 39% for $\alpha=0.05$ or $0.1$). Using the optimal robust \ac{EBCI} ensures that the efficiency loss is below 20%, irrespective of the signal-to-noise ratio. On the other hand, when the signal-to-noise ratio is small, any of these \acp{CI} will be significantly tighter than the unshrunk \ac{CI} $Y_{i}\pm z_{1-\alpha/2}\sigma_{i}$. To illustrate this point, (ref) plots the efficiency of the robust \ac{EBCI} that imposes the second moment constraint only, relative to this unshrunk \ac{CI}. It can be seen from the figure that shrinkage methods allow us to tighten the \ac{CI} by 44% or more when $\mu_{2}/\sigma^{2}_{i}\leq 0.1$.

Undercoverage of parametric EBCI

The parametric \ac{EBCI} $\hat\theta_i \pm z_{1-\alpha/2}w_{EB, i}^{1/2}\sigma_{i}$ is an \ac{EB} version of a Bayesian credible interval that treats (ref) as a prior. We now assess its potential undercoverage when (ref) is violated.

Given knowledge of only the second moment $\mu_{2}$ of $\varepsilon_i=Y_i-X_i'\delta$, the maximal undercoverage of this interval is given by

since $w_{EB, i}= \mu_{2}/(\mu_{2}+\sigma_i^2)$. Here $\rho$ is the non-coverage function defined in (ref). (ref) plots the maximal non-coverage probability as a function of $w_{EB, i}$, for significance levels $\alpha=0.05$ and $\alpha=0.10$. The figure suggests a simple “rule of thumb”: if $w_{EB, i} \geq 0.3$, the maximal coverage distortion is less than 5 percentage points for these values of $\alpha$.

The following lemma confirms that the maximal non-coverage is decreasing in $w_{EB, i}$, as suggested by the figure. It also gives an expression for the maximal non-coverage across all values of $w_{EB, i}$ (which is achieved in the limit $w_{EB, i} \to 0$).

The maximal non-coverage probability $1/\max\{z_{1-\alpha/2}^2,1\}$ equals $0.260$ for $\alpha=0.05$ and $0.370$ for $\alpha=0.10$. For $\alpha>2\Phi(-1)\approx 0.317$, the maximal non-coverage probability is 1.

If we additionally impose knowledge of the kurtosis of $\varepsilon_i$, the maximal non-coverage of the parametric \ac{EBCI} can be similarly computed using (ref), as illustrated in the application in (ref).

Monte Carlo simulations

Here we show through simulations that the robust \ac{EBCI} achieves accurate average coverage in finite samples.

Design

The \ac{DGP} is a simple linear fixed effects panel data model. We first draw $\theta_i$, $i=1,\dotsc, n$, i.i.d.\ from a random effects distribution specified below. Then we simulate panel data from the model

where the errors $U_{it}$ are mean zero and i.i.d.\ across $(i, t)$ and independent of the $\theta_i$'s. The unshrunk estimator of $\theta_i$ is the sample average of $W_{it}$ for unit $i$, with standard error obtained from the usual unbiased variance estimator:

We draw $U_{it}$ from one of two distributions: (1) a normal distribution and (2) a (shifted) chi-squared distribution with 3 degrees of freedom. In case (1), $Y_i$ is exactly normal conditional on $\theta_i$, but $\hat{\sigma}_i^2$ does not exactly equal $\operatorname{var}(Y_i \mid \theta_i)$ for finite $T$. In case (2), $Y_i$ is non-normal and positively skewed (conditional on $\theta_i$) for finite $T$.

We consider six random effects distributions for $\theta_i$ (see \Cref*{sec:sim-details} for detailed definitions):

Given $T$, we scale the $\theta_i$ distribution to match one of four signal-to-noise ratios $\mu_2/\operatorname{var}(Y_i \mid \theta_i) \in \lbrace 0.1,0.5,1,2 \rbrace$, for a total of $6 \times 4 = 24$ \acp{DGP} for each distribution of $U_{it}$. We shrink towards the grand mean ($X_i=1$ for all $i$). We construct the robust \acp{EBCI} following the baseline implementation in (ref) (with $\omega_{i}=1/n$), as well as a version that does not impose constraints on the kurtosis.

As $T \to \infty$, we recover the idealized setting in (ref), with $(Y_{i}-\theta_{i})/\sqrt{\operatorname{var}(Y_i \mid \theta_i)}$ converging in distribution to a standard normal (conditional on $\theta_i$), and $\hat{\sigma}_i^2/\operatorname{var}(Y_i \mid \theta_i)$ converging in probability to 1, for each $i$.

Results

(ref) shows that the 95% robust \acp{EBCI} achieve good average coverage when the panel errors $U_{it}$ are normally distributed. This is true for all \acp{DGP}, panel dimensions $n$ and $T$, and whether we exploit one or both of the (estimated) moments $\mu_2$ and $\kappa$. When the time dimension $T$ equals 10, the maximal coverage distortion across all \acp{DGP} and all cross-sectional dimensions $n \in \lbrace 100,200,500\rbrace$ is 3.2 percentage points. For $T \geq 20$, the coverage distortion of the robust \acp{EBCI} is always below 2.1 percentage points.

(ref) shows that coverage distortions are somewhat larger when the panel errors $U_{it}$ are chi-squared distributed and $T$ is small. The robust \acp{EBCI} undercover by up to 7.2 percentage points when $T=10$ due to the pronounced non-normality of $Y_i$ given $\theta_i$. However, the distortion is at most 4.3 percentage points when $T=20$, and at most 2.4 percentage points when $T \geq 50$. The coverage distortion due to non-normality when $T$ is small is similar to the coverage distortion of the usual unshrunk \ac{CI} (not reported).

Importantly, in all cases considered in (ref), the worst-case coverage distortion of the parametric \ac{EBCI} substantially exceeds that of the corresponding robust \acp{EBCI}, sometimes by more than 10 percentage points. Nevertheless, the cost of robustness in terms of extra \ac{CI} length is modest and consistent with the theoretical results in (ref).

Both the estimation of the standard errors $\sigma_i$ and the estimation of the moments $\mu_2$ and $\kappa$ contribute to the finite-sample coverage distortions. The “ora” columns in (ref) exploit the oracle (true) values of $\mu_2$, $\kappa$, and $\sigma_i=\sqrt{\operatorname{var}(Y_i \mid \theta_i)}$, while the $T=\infty$ columns use oracle standard errors but not oracle moments. By comparing these columns, we see that estimation of $\mu_2$ and $\kappa$ is responsible for modest coverage distortions when $n=100$ or $200$. However, estimation of the standard errors $\sigma_i$ also contributes to the distortions, as can be seen by comparing the $T=10$ and $T=\infty$ columns.

In \Cref*{sec:sim-heterosk} we show that the robust \ac{EBCI} also has good coverage in a heteroskedastic design calibrated to the empirical application in (ref) below.

Comparison with other approaches

Here we compare our \ac{EBCI} procedure with other approaches to confidence interval construction in the normal means model. We also discuss other related inference problems.

Average coverage vs.\ alternative coverage concepts

The average coverage requirement in (ref) is less stringent than the usual (pointwise) notion of frequentist coverage that $P(\theta_i\in CI_i\mid \theta) \geq 1-\alpha$ for all $i$. An even stronger coverage requirement is that of simultaneous coverage: $P(\forall i\colon \theta_i\in CI_{i}\mid \theta) \geq 1-\alpha$. As outlined in (ref), under the pointwise coverage criterion, one cannot achieve substantial reductions in length relative to the unshrunk \ac{CI}. Under the simultaneous coverage criterion, it is likewise impossible to substantially improve upon the usual sup-$t$ confidence band based on the unshrunk estimates Cai2015. Thus, undercoverage for some $\theta_{i}$'s must be tolerated if one wants to use shrinkage to improve \ac{CI} length.

The fact that our \acp{EBCI} achieve improvements in average length at the expense of undercovering for certain units $i$ is analogous to well-known properties of \ac{EB} point estimators. We now show that the units $i$ for which our \ac{EBCI} undercovers are quantitatively similar to the units for which the shrinkage estimator $\hat{\theta}_{i}$ has higher \ac{MSE} than the unshrunk estimator $Y_i$. Let $\varepsilon_{i}=\theta_{i}-X_{i}'\delta$ be the “shrinkage error” defined in (ref). The pointwise coverage of our \ac{EBCI} is decreasing in the normalized shrinkage error $\abs{\varepsilon_{i}}/\sqrt{\mu_2}$, for a fixed signal-to-noise ratio $\mu_2/\sigma_i^2$.\footnote{The pointwise coverage (conditional on $X_{i}$) equals $1-r(\sqrt{1/w_{EB, i}-1}\cdot \abs{\varepsilon_{i}}/\sqrt{\mu_{2}}, \operatorname{cva}_{\alpha}(1/w_{EB, i}-1,\kappa))$, with $r$ defined in (ref) and $w_{EB, i}=\mu_{2}/(\mu_{2}+\sigma_{i}^{2})$.} Hence, the units $i$ for which our \ac{EBCI} undercovers are those whose covariate-predicted value $X_i'\delta$ fails to approximate their true effect $\theta_i$ well. The \ac{MSE} of the shrinkage estimator (for an individual unit $i$), normalized by the \ac{MSE} of the unshrunk estimator, is similarly increasing in $\abs{\varepsilon_{i}}/\sqrt{\mu_2}$.\footnote{The ratio of \acp{MSE} equals $E[(\hat{\theta}_i-\theta_i)^2 \mid \theta_i, X_i]/\sigma_{i}^{2}=w_{EB, i}^{2}+(1-w_{EB, i}) w_{EB, i}\cdot \abs{\varepsilon_{i}}/\sqrt{\mu_2}$.}

(ref) shows that the knife-edge value of $\abs{\varepsilon_{i}}/\sqrt{\mu_2}$ for which the pointwise coverage of our \ac{EBCI} equals $1-\alpha$ is quantitatively close to the value of $\abs{\varepsilon_i}/\sqrt{\mu_2}$ for which the MSE of the shrinkage estimator equals that of the unshrunk estimator. In other words, to the extent that one worries about undercoverage for certain types of $\theta_i$ values, one should simultaneously worry about the relative performance of the shrinkage point estimator for those same values.

We stress that the pointwise coverage depends on the unobservable shrinkage error $\varepsilon_{i}$, which cannot be gauged directly from the observables $(Y_{i}, X_{i})$. If one wishes to avoid systematic differences in coverage across units $i$ with different genders, say (i.e.,\ one is worried that $\varepsilon_{i}$ correlates with gender) one can simply add gender to the set of covariates $X_{i}$: the baseline procedure in (ref) ensures control of average coverage conditional on the covariates $X_i$. In (ref), we show how to adapt our \acp{EBCI} to settings where one focuses the analysis on a subset of units $i$ based on the values of their unshrunk estimates $Y_i$ (e.g., keeping only the estimates that exceed a given threshold).

From a Bayesian point of view, our robust \ac{EBCI} can be viewed as an uncertainty interval that is robust to the choice of prior distribution in the unconditional gamma-minimax sense: the coverage probability of this \ac{CI} is at least $1-\alpha$ when averaged over the distribution of the data and over the prior distribution for $\theta_{i}$, for any prior distribution that satisfies the moment bounds. This follows directly from the derivations in (ref), reinterpreting the random effects distribution for $\theta_i$ as a prior distribution. In contrast, conditional gamma-minimax credible intervals, discussed recently by Kitagawa2019, are too stringent in our setting. This notion requires that the posterior credibility of the interval be at least $1-\alpha$ regardless of the choice of prior, in any data sample, which would require reporting the entire parameter space (up to the moment bounds).

Finite-sample vs.\ asymptotic coverage

Our procedures are asymptotically valid as $n\to\infty$, as proved in (ref). These asymptotics do not capture the impact of estimation error in the “hyper-parameters” $\hat{\sigma}_i$, $\hat{\delta}$, $\hat{\mu}_2$, and $\hat{\kappa}$, or the impact of lack of exact normality of the $Y_{i}$'s, on the finite-sample performance of the \acp{EBCI}. As detailed in (ref) and (ref), we do apply a finite-sample adjustment to the moments $\hat{\mu}_2$ and $\hat{\kappa}$, which is motivated by the same heuristic arguments that morris83pebci,morris83 uses to motivate finite-sample adjustments to the parametric \ac{EBCI}.\footnote{An alternative approach would be to adapt the bootstrap adjustment proposed by CaLo00 in the context of parametric \ac{EBCI} construction efron2019. As with the morris83pebci,morris83 adjustment, we are not aware of a formal result justifying it.} The promising simulation results in (ref) notwithstanding, these adjustments do not ensure exact average coverage control in finite samples.\footnote{One could account for hyperparameter uncertainty by computing the critical value $\sup_{\tilde{\sigma}_i, \tilde{\mu}_2,\tilde{\kappa} \in \hat{\mathcal{C}}_i} \operatorname{cva}_{\alpha}(\tilde{\sigma}_i^{2}/\tilde{\mu}_{2}, \tilde{\kappa})$ over an initial confidence set $\hat{\mathcal{C}}_i$ for the hyper-parameters, coupled with a Bonferroni adjustment of the confidence level $1-\alpha$. This approach appears to be highly conservative in practice.}

Our results are thus analogous to standard results on coverage of Eicker-Huber-White \acp{CI} in cross-sectional \ac{OLS}: asymptotic validity follows by consistency of the \ac{OLS} variance estimate and asymptotic normality of the outcomes, while adjustments to account for finite-sample issues (such as the HC2 or HC3 variance estimators studied in MaWh85) are justified heuristically. Deriving \acp{EBCI} with finite-sample coverage guarantees is an interesting problem that we leave for future research; the problem appears to be challenging even in the context of constructing parametric \acp{EBCI}.

Local vs.\ global optimality

Our \acp{EBCI} are designed to provide uncertainty assessments to accompany linear shrinkage estimates that, as the Introduction argues, have been popular in applied work. Our procedure's global validity, as well as local near-optimality when the $\theta_{i}$'s are normal (cf. (ref)), is analogous to Eicker-Huber-White \acp{CI} for \ac{OLS} estimators: these \acp{CI} are optimal under normal homoskedastic regression errors, but remain valid when this assumption is dropped.

Similar to the Eicker-Huber-White \acp{CI}, our \acp{EBCI} are not globally efficient: when the $\theta_{i}$'s are not Gaussian, it is generally inefficient to restrict attention to \acp{CI} that are centered at a linear point estimator and have fixed width. While we expect our \acp{EBCI} to remain near-efficient under mild departures from normality, substantial efficiency gains may be possible if the effect size distribution is, for example, heavy-tailed or bimodal.\footnote{Indeed, if the true effect distribution puts mass $1/2$ on $\theta_{i}=K$ and $\theta_{i}=-K$, then, as $K$ gets large, our \acp{EBCI} become arbitrarily conservative relative to an oracle that reports the highest posterior density set under this prior.} (ref) shows how our method can be adapted to construct \acp{EBCI} that are locally near-optimal under non-normal baseline priors using non-linear shrinkage, such as soft thresholding. Since the distribution of $\theta_{i}$ is nonparametrically identified under the normal model (ref), it is in principle possible to construct \acp{EBCI} that are globally efficient using nonparametric methods. In the context of the homoskedastic model with no covariates in (ref), various approaches to nonparametric point estimation of the $\theta_{i}$'s have been proposed, including kernels BrGr09, splines efron2019, or nonparametric maximum likelihood KiWo56,JiZh09,KoMi14. An interesting problem for future research is to adapt these methods to \ac{EBCI} construction, while ensuring asymptotic validity, good finite-sample performance, and allowing for covariates, heteroskedasticity, and possible dependence across $i$.

Other inference problems

A number of alternative inference procedures have been proposed in the context of the normal means model. Efron2015 develops a formula for the frequentist standard error of \ac{EB} estimators, but this cannot be used to construct \acp{CI} without a corresponding estimate of the bias. There is a substantial literature on shrinkage confidence balls, i.e., confidence sets of the form $\{\theta\colon\sum_{i=1}^{n} (\theta_i-\hat\theta_i)^2\le \hat c\}$ (see casella_shrinkage_2012, for a review). While theoretically interesting, these sets can be difficult to visualize and report in practice.\footnote{Confidence balls can be translated into average coverage intervals using Chebyshev's inequality Wasserman2006. However, such intervals are very conservative compared to the ones we construct.}

Finally, while we focus on \ac{CI} length in our relative efficiency comparisons, our approach can be fruitfully applied when the goal of \ac{CI} construction is to discern non-null effects, rather than to construct short \acp{CI}. In particular, suppose one forms a test of the null hypothesis $H_{0,i}:\theta_i=\theta_0$ for some null value $\theta_0$ by rejecting when $\theta_0\notin CI_i$, where $CI_i$ is our robust \ac{EBCI} given in (ref). In \Cref*{sec:power-details}, we show that the test based on our \ac{EBCI} has higher average power than the usual $z$-test based on the unshrunk estimate when $X_i'\delta$ (the regression line towards which we shrink) is far enough from the null value $\theta_0$, and that these power gains can be substantial. Furthermore, such tests can be combined with corrections from the multiple testing literature to form procedures that asymptotically control the \ac{FDR}, a commonly used criterion for multiple testing.\footnote{In particular, storey_direct_2002 shows that the benjamini_controlling_1995 procedure asymptotically controls the \ac{FDR} so long as the $p$-values do not exhibit too much statistical dependence and the proportion of rejected null hypotheses does not converge too quickly to zero. While storey_direct_2002 assumes that the uncorrected tests control size in the classical sense, the argument goes through essentially unchanged so long as the tests invert \acp{CI} that satisfy ((ref)), which holds so long as the \acp{CI} do not exhibit too much statistical dependence, as discussed in (ref). We note, however, that this does not hold for modifications of the benjamini_controlling_1995 procedure that use initial estimates of the proportion of true null hypotheses.}

Extensions

We now discuss two extensions of our method: adapting our intervals to general, possibly non-linear shrinkage, and constructing intervals that achieve coverage conditional on $Y_{i}$ falling into a pre-specified interval.

General shrinkage

Our method can be generalized to cover general, possibly non-linear shrinkage based on possibly non-Gaussian data. Let $\mathcal{S}(y; \chi, \tilde{X}_{i})\subseteq \mathbb{R}$ be a family of candidate confidence sets for a parameter $\theta_{i}$, which depends on the data $Y_{i}=y$, a tuning parameter $\chi\in\mathbb{R}$ to be selected below, and covariates $\tilde{X}_{i}$ (that include any known nuisance parameters) that we treat as fixed. We assume that $\mathcal{S}$ is increasing in $\chi$, in the sense of set containment, and that the non-coverage probability conditional on $\theta$ satisfies

where $a_{i}$ is some function of $\theta_i$, $\tilde{X}^{(n)}=(\tilde{X}_{1}, \dotsc, \tilde{X}_{n})$, and $\tilde r$ is a known function (perhaps computed numerically or through simulation). Similarly to linear shrinkage in the normal means model, (ref) may only hold approximately if the set $\mathcal{S}$ depends on estimated parameters (such as standard error estimates or tuning parameters), or if we use a large-sample approximation to the distribution of $Y_{i}$. We assume that $a_{i}$ satisfies the moment constraints $E_{F}[g(a_{i})\mid \tilde{X}^{(n)}]=m$, where $g$ is a $p$-vector of moment functions, and the expectation is over the conditional distribution $F$ of $a_{i}$ conditional on $\tilde{X}^{(n)}$.\footnote{The moment functions $g$ need not be simple moments, and could incorporate constraints used for selection of hyper-parameters, such as constraints on the marginal data distribution or, if an unbiased risk criterion is used, the constraint that the derivative of the risk equals zero at the selected prior hyper-parameters.} To guarantee \ac{EB} coverage, we compute the maximal non-coverage

analogously to (ref). This is a linear program, which can be computed numerically to a high degree of precision even with several constraints; see (ref) for details. Given an estimate $\hat m$ of the moment vector $m$, we form a robust \ac{EBCI} as

Coverage after selection

In some applications, researchers may be primarily interested in parameters corresponding to those units $i$ whose initial estimates $Y_{i}$ fall in a given interval $[\iota_1,\iota_2]$, where $-\infty \leq \iota_1 < \iota_2 \leq \infty$. For example, in a teacher value added application, we may only be interested in the ability $\theta_i$ of those teachers $i$ whose fixed effect estimates $Y_i$ are positive, corresponding to setting $\iota_1=0$ and $\iota_2=\infty$. Because of the selection on outcomes, na\"{i}vely applying our baseline \ac{EBCI} procedure to the selected sample $\lbrace i\colon Y_i \in [\iota_1,\iota_2]\rbrace$ does not yield the desired average coverage across the selected units $i$. We now show how to correct for the selection bias in the simple homoskedastic model $Y_i \mid \theta_i \sim N(\theta_i, \sigma^2)$ without covariates from (ref) (reintroducing the extra model features in (ref) only complicates notation).

We seek a critical value $\chi$ such that the average coverage of the \ac{CI} $[\hat{\theta}_i \pm \chi w_{EB} \sigma]$ is at least $1-\alpha$ conditional on the sample selection, i.e.,

under repeated sampling of $(Y_i, \theta_i)$, regardless of the distribution for $\theta_i$ (we maintain focus on linear shrinkage for simplicity, but our approach extends to nonlinear shrinkage using the ideas in (ref)). Straightforward calculations show that the non-coverage, conditional on $\theta_{i}$ and on selection, equals

where $b_{i}=(1-1/w_{EB})\theta_{i}/\sigma$ as in (ref). Among all distributions for $\theta_i$ consistent with the conditional moment $\tilde{\mu}_{2,\iota_1,\iota_2}=E\left[\theta_i^2 \mid Y_i \in [\iota_1,\iota_2]\right]$, the worst-case non-coverage probability, conditional on selection, is given by

where $E_{F}$ denotes expectation under $\theta_{i}\sim F$. This is an infinite-dimensional linear program that can be solved numerically to a high degree of accuracy, cf.\ (ref). To achieve robust conditional coverage, we solve numerically for the $\chi$ such that $\tilde{\rho}_{\iota_1,\iota_2}(\tilde{\mu}_{2,\iota_1,\iota_2}, \chi) = \alpha$.

We can estimate the conditional second moment $\tilde{\mu}_{2,\iota_1,\iota_2}$ as follows. Denote the log marginal density of $Y_i$ by $\ell(y) \equiv \log \int \phi(y-\theta)\, d\Gamma_0(\theta)$, where $\Gamma_0$ is the true distribution of $\theta_i$. Tweedie's formulas efron2019 imply

Let $\hat{\ell}(y)$ be a kernel estimate of the log marginal density function of the data $Y_1, \dotsc, Y_n$. Then the estimate

will be consistent as $n\to\infty$ for $\tilde{\mu}_{2,\iota_1,\iota_2}$ in (ref) under mild regularity conditions.

The criterion (ref) can be viewed as the \ac{EB} analogue of the criterion $P(\theta_i \in CI_i \mid Y_i \in [\iota_1,\iota_2], \theta) \geq 1-\alpha$, which requires frequentist coverage conditional on the event $\{Y_i\in [\iota_1,\iota_2]\}$. The latter criterion has been considered in the recent “selective inference” literature benjamini_false_2005,lee_exact_2013,HuFi19,Andrews2021. In contrast to this literature, we cannot allow $\iota_1$ to be given by the maximum of the initial estimates Andrews2021, as we require $\iota_{1}$ and $\iota_{2}$ to converge in probability to distinct nonrandom limits. On the other hand, weakening the notion of frequentist coverage to \ac{EB} (or average) coverage allows for improvements in the length of the intervals, similar to the analysis in (ref) in the absence of selection.

Empirical application

We illustrate our methods using the data and model in chetty_impacts_2018, who are interested in the effect of neighborhoods on intergenerational mobility.

Framework

We follow chetty_impacts_2018 in using two definitions of a “neighborhood effect” $\theta_{i}$. The first focuses on effects for children growing up in low-income families, and defines $\theta_{i}$ as the effect of spending an additional year of childhood in \ac{CZ} $i$ on children's rank in the income distribution at age 26, for children with parents at the 25th percentile of the national income distribution. The second definition is analogous, except it focuses on children growing up in high-income families, and consequently conditions on children with parents at the 75th percentile. chetty_impacts_2018 argue that these definitions approximately capture the mean rank effects for children in below-median and above-median income families. Using de-identified tax returns for all children born between 1980 and 1986 who move across \acp{CZ} exactly once as children, chetty_impacts_2018 exploit variation in the age at which children move between \acp{CZ} to obtain preliminary fixed effect estimates $Y_{i}$ of $\theta_{i}$.

Since these preliminary estimates are measured with noise, to predict $\theta_{i}$, chetty_impacts_2018 shrink $Y_{i}$ towards average outcomes of permanent residents of \ac{CZ} $i$ (children with parents at the same percentile of the income distribution who spent all of their childhood in the \ac{CZ}). To give a sense of the accuracy of these forecasts, chetty_impacts_2018 report estimates of their unconditional \ac{MSE} (i.e.,\ treating $\theta_{i}$ as random), under the implicit assumption that the moment independence assumption in (ref) holds. Here we complement their analysis by constructing robust \acp{EBCI} associated with these forecasts.

Our sample consists of 595 U.S. \acp{CZ}, with population over 25,000 in the 2000 census: this is the sample for which chetty_impacts_2018 report baseline estimates $Y_{i}$ of the effects $\theta_{i}$. These baseline estimates are normalized so that their population-weighted mean is zero. We may therefore interpret $\theta_{i}$ as the effect relative to an “average” \ac{CZ}. We follow the baseline implementation from (ref) with standard errors $\hat{\sigma}_{i}$ reported by chetty_impacts_2018, and covariates $X_{i}$ corresponding to a constant and the average outcomes for permanent residents. In line with the original analysis, we use precision weights $\omega_i=1/\hat{\sigma}^{2}_{i}$ when constructing the estimates $\hat{\delta}$, $\hat{\mu}_{2}$ and $\hat{\kappa}$.

Results

Columns (1) and (2) in (ref) summarize the main estimation and efficiency results. The shrinkage magnitude and relative efficiency results are similar for children with parents at the 25th and 75th percentiles of the income distribution. In both columns, the estimate of the kurtosis $\kappa$ is large enough so that it does not affect the critical values or the form of the optimal shrinkage: specifications that only impose constraints on the second moment yield identical results.\footnote{The truncation in the $\hat{\kappa}$ formula in our baseline algorithm in (ref) binds in columns (1) and (2), although the non-truncated estimates 345.3 and 5024.9 are similarly large; using these non-truncated estimates yields identical results.} In line with this finding, a density plot of the $t$-statistics (reported as Figure S2 in akp20v2) exhibits a fat lower tail. As a robustness check, columns (3) and (4) show that results based on nonparametric moment estimates (see (ref)) are very similar to our baseline specification. Indeed, the $R^{2}$ gain in predicting $\hat{\varepsilon}_{i}^{2}-\hat{\sigma}^{2}_{i}$ using $\hat{\mu}_{2i}$ is less than $0.001$ in both specifications, indicating that there is little evidence in the data against the moment independence assumption.

The baseline robust 90% \acp{EBCI} are 75.2--87.7% shorter than the usual unshrunk \acp{CI} $Y_{i}\pm z_{1-\alpha/2}\hat{\sigma}_{i}$. To interpret these gains in dollar terms, for children with parents at the 25th percentile of the income distribution, a percentile gain corresponds to an annual income gain of \$818 chetty_impacts_2018. Thus, the average half-length of the baseline robust \acp{EBCI} in column (1) implies \acp{CI} of the form $\pm \$160$ on average, while the unshrunk \acp{CI} are of the form $\pm \$643$ on average. These large gains are a consequence of a low signal-to-noise ratio $\mu_{2}/\sigma_{i}^{2}$ in this application. Because the shrinkage magnitude is so large on average, the tail behavior of the bias matters, and since the kurtosis estimates suggests these tails are fat, it is important to use the robust critical value: the parametric \ac{EBCI} exhibits average potential size distortions of 12.7--17.8 percentage points. Indeed, for over 90% of the \acp{CI} in the specifications in columns (1) and (2), the shrinkage coefficient $w_{EB, i}$ falls below the “rule of thumb” threshold of 0.3 derived in (ref).

To visualize these results, (ref) plots the unshrunk 90% \acp{CI} based on the preliminary estimates, as well as robust \acp{EBCI} based on \ac{EB} estimates for cities in the state of New York for children with parents at the 25th percentile. While the \acp{EBCI} for large \acp{CZ} like New York City or Buffalo are similar to the unshrunk \acp{CI}, they are much tighter for smaller \acp{CZ} like Plattsburgh or Watertown, with point estimates that shrink the preliminary estimates $Y_{i}$ most of the way toward the regression line $X_{i}'\hat{\delta}$.

In summary, shrinkage allows us to considerably tighten the \acp{CI} based on preliminary estimates. This is true even though the \acp{CI} effectively only use second moment constraints---imposing kurtosis constraints does not affect the critical values in this application.