Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
122,536 characters · 25 sections · 0 citation commands
Invidious Comparisons: Ranking and Selection as Compound Decisions
\pagestyle{myheadings} \markboth{\sc Invidious Comparisons}{\sc Gu and Koenker}
In the wake of Wald's seminal monograph on statistical decision theory there was a growing awareness that the Neyman-Pearson testing apparatus was inadequate for many important statistical tasks. Ranking and selection problems featured prominently in this perception. Motivated by a suggestion of Harold Hotelling, \citeasnoun{Bahadur50} studied selection of the best of several Gaussian populations. Assuming that sample means were observed for each of $K$ populations with means, $\theta_k$ and common variance, the problem of selecting the best population, $\theta^* = \max_i \{\theta_1, \ldots , \theta_K \}$, was formulated as choosing weights $z_1, \cdots , z_K$ to minimize, \[ L(\theta, z) = \theta^* - \sum_{k=1}^K z_k \theta_k/\sum_{k=1}^K z_k. \] Bahadur showed that among “impartial decision rules,” i.e. permutation equivariant rules, it was uniformly optimal to select only the population with the largest sample mean, that is to choose $z_i^* = 1$ if $\bar X_i = \max \{ \bar X_1, \cdots , \bar X_K \} $ and $z_i^* = 0$ otherwise, thereby clearly demonstrating that procedures that did preliminary tests of equality of means and then chose $z_i > 0$ for several or even all of the populations when tests failed to reject were inadmissible. This finding was reinforced in \citeasnoun{BR50} who focused on the two-sample setting but relaxed the common variance assumption. In related work, \citeasnoun{Bechhofer} and \citeasnoun{Gupta56} sought to optimize the number of selected populations as well as their identities, see \citeasnoun{GuptaPan} and \citeasnoun{BKS} for extensive reviews of subsequent developments.
\citeasnoun{GoelRubin} pioneered the hierarchical Bayesian approach to selection that has been adopted by numerous authors in the ensuing decades, early on by \citeasnoun{BergerDeeley} and \citeasnoun{LairdLouis}. \citeasnoun{Portnoy} showed that rankings based on best linear predictors were optimal in Gaussian multivariate variance components models, but cautioned that departures from normality could easily disrupt this optimality. A notable feature of the hierarchical model paradigm is the recognition that sample observations may exhibit heterogeneous precision; this is typically accounted for by assuming known variances for observed sample means. As ranking and selection methods became increasingly relevant in genomic applications there has been renewed interest in loss functions and linkages to the burgeoning literature on multiple testing. Our perspective is informed by recent developments in the nonparametric estimation of mixture models and its relevance for a variety of compound decision problems. This approach seeks to reduce the reliance on Gaussian distributional assumptions that pervades the earlier literature. As we have argued elsewhere, \citeasnoun{ProbRob}, and \citeasnoun{Minimalist} nonparametric empirical Bayes methods offer powerful complementary methods to more conventional parametric hierarchical Bayes for multiple testing and compound decision problems. Our primary objective in this paper is to elaborate this assertion for ranking and selection applications. Throughout we try to draw parallels and contrasts with the literature on multiple testing. We will restrict our attention to settings where we observe a scalar estimate of an unobserved latent quality measure accompanied by some measure of its precision, thereby evading more complex multivariate settings, as in \citeasnoun{accuracy} who employ quantile regression methods.
An important motivation for revived interest in ranking and selection problems in econometrics has been the influential work of \nocite{Chettya,Chettyb,ChettyMob} Chetty and his collaborators on teacher evaluation and geographic mobility in the U.S. This has stimulated the important recent work of \citeasnoun{MRSW} proposing new resampling methods for constructing confidence sets for ranking and selection for a finite population. \citeasnoun{AKPM} propose an innovative approach to the construction of confidence intervals for classical, linear shrinkage, empirical Bayes estimators of the type used by Chetty. Recent work by \citeasnoun{AKM} and \citeasnoun{GuoHe} propose new confidence interval constructions for highly ranked individuals or treatments influenced by recent contributions to the “inference after model selection” literature. In contrast to these inferential approaches we focus instead on the complementary perspective of compound decision making, constructing decision rules for selecting the best, or worst, populations subject to control of the expected number of elements selected and among those selected, the expected proportion of false discoveries. Rather than treating each selection decision in isolation, the compound decision framework tries to exploit their common structure to produce improved collective performance. Our approach is thus more closely aligned to that of \citeasnoun{KlineWalters} who study decision rules for assessing employer discrimination from experiments involving fictitious job applications using closely related GMM methods for binomial mixture models.
\citeasnoun{GGM} study teacher value added estimation employing nonparametric maximum likelihood methods for estimating Gaussian mixture models as we advocate below. Their analysis of data from both North Caroline and Los Angeles illustrates the advantages of more flexible mixture models for latent value added. In contrast to the present work, they focus on Bayes rules for posterior means that are often used to study teachers' influence on students' future outcomes. These more flexible nonparametric empirical Bayes methods improve upon traditional linear shrinkage rules especially in the tails of the distribution where policy attention is usually focused. This is a valuable, complementary perspective to the ranking and selection objectives of the present work.
Before proceeding it is important to acknowledge that despite its universal appeal and application there is something inherently futile about many ranking and selection problems as intimated by our title. If the latent measure of true quality is Gaussian, as assumed in virtually all of the econometric applications of the selection problem, and we wish to select the top ten percent of individuals given that their true quality is contaminated by Gaussian noise, accurate selection can be very challenging when the signal to noise ratio is low. We will see that conventional linear shrinkage as embodied in the classical James-Stein formula can improve performance considerably over naive maximum likelihood (fixed effects) procedures, and some further improvement is possible by carefully tailoring the decision rules for tail probability loss. However, we find that even oracle decision rules that incorporate complete knowledge of the precise distributional features of the problem may not be able achieve better than about even odds that selected individuals have latent ability above the selection thresholds when measurement error is comparable in magnitude to Gaussian variability in latent ability. When the latent distribution of ability is heavier tailed then selection becomes somewhat easier, and more refined selection rules are more advantageous, but as we will show the selection problem still remains quite challenging.
Thus, a secondary objective of the paper is to add another cautionary voice to those who have already questioned the reliability of existing ranking and selection methods. A critical overview of the role of ranking and selection in public policy applications is provided by \citeasnoun{LeagueTables}. It is widely acknowledged that league tables as currently employed can be a pernicious influence on policy, a viewpoint underscored in \citeasnoun{GelmanPrice}. While much of this criticism can be attributed to inadequate data collection and inherently low signal to noise ratios, we believe that there is also room for methodological improvements.
Section (ref) provides a brief overview of compound decision theory and describes nonparametric methods for estimation of Gaussian mixture models. Section (ref) introduces a basic framework for our approach to ranking and selection in a setting with homogeneous precision of the observed measurements. In Section (ref) we introduce heterogeneous precision of known form, and Section (ref) considers settings in which the joint distribution of the observed measurements and their precision determines the form of the ranking and selection rules. Optimal ranking and selection rules are derived in each of these sections under the assumption that the form of the mixing distribution of the unobserved, latent quality of the observations is known. Section (ref) introduces feasible ranking and selection rules and conditions under which they attain the same asymptotic performance as the optimal rules. Section (ref) then compares several {\it{feasible}} ranking and selection methods, some that ignore the compound decision structure of the problem, some that employ parametric empirical Bayes methods and some that rely on nonparametric empirical Bayes methods. Finally, Section (ref) describes an empirical application on evaluating the performance of medical dialysis centers in the United States. Proofs of all formal results are collected in Appendix A.
\citeasnoun{Robbins.51} posed a challenge to the nascent minimax decision theory of \citeasnoun{Wald50}: Suppose we observe independent Gaussian realizations, $Y_i \sim \mathcal N(\theta_i , 1), \; i = 1, \cdots, n$ with means $\theta_i$ taking either the value $+1$ or $-1$. We are asked to estimate the $n$-vector $\theta = (\theta_1, \cdots , \theta_n)$ subject to mean absolute error loss, \[ L(\hat \theta , \theta) = n^{-1} \sum_{i=1}^n | \hat \theta_i - \theta_i |. \] When $n = 1$ Robbins shows that the minimax decision rule is $\delta (y) = \mbox{sgn} (y)$; in the least favorable variant of the problem malevolent nature chooses $\pm 1$ with equal probability, and the optimal response is to estimate $\theta_i = +1$ when $Y_i$ is positive, and $\theta_i = -1$ otherwise. Robbins goes on to show that when $n > 1$, this rule remains minimax, each coordinate is treated independently as if viewed in complete isolation. This is also the maximum likelihood estimator, and may be viewed in econometrics terms as a classical fixed-effects estimator. But is it at all reasonable?
Doesn't our sample convey information about the relative frequency of $\pm 1$ that might potentially contradict the pessimistic presumption of the minimax rule? If we happened to know the unconditional probability, $p = \mathbb{P} (\theta_i = 1)$, then the conditional probability that $\theta = 1$ given $Y_i = y$, is given by, \[ \mathbb{P}(\theta = 1 | y) = \frac{p \varphi(y-1)}{p \varphi(y-1) + (1-p) \varphi(y + 1)}, \] where $\varphi$ denotes the standard Gaussian density. We should guess $\hat \theta_i = 1$ if this probability exceeds 1/2, giving us the revised decision rule, \[ \delta_p (y) = \mbox{sgn} ( y - \tfrac{1}{2} \log ( (1 - p)/p)). \] Each observed, $y_i$, is modified by a simple logistic perturbation before computing the sign. Our observed random sample, $\bm{y} = (y_1, \cdots , y_n)$, is informative about $p$. We have the log likelihood, \[ \ell_n (p | y) = \sum_{i=1}^n \log ( p \varphi(y_i - 1) + (1-p) \varphi(y_i + 1)), \] which could be augmented by a prior of some form, if desired, to obtain a posterior mean for $p$ and a plug-in Bayes rule for estimating each of the $\theta_i$'s. The Bayes risk of this procedure is substantially less than the minimax risk when $p \neq 1/2$ and is asymptotically equivalent to the minimax risk when $p = 1/2$. This is the first principle of compound decision theory: borrowing strength across an entire ensemble of related decision problems yields improved collective performance.
What happens when we relax the restriction on the support of the $\theta$'s and allow support on the whole real line? We now have a general Gaussian mixture setting where the observed $Y_i$'s have marginal density given by the convolution, $f = \varphi * G$, that is, \[ f(y) = \int \varphi(y - \theta) dG(\theta), \] and instead of merely needing to estimate one probability we need an estimate of an entire distribution function, $G$. \citeasnoun{KW}, anticipated by an abstract of \citeasnoun{R50}, established that the nonparametric maximum likelihood estimator (NPMLE), \[ \hat G = \mbox{argmin}_{G \in \mathcal G} \{ - \sum_{i=1}^n \log f(y_i) \; | \; f(y_i) = \int \varphi(y_i - \theta)dG(\theta) \} \] where $\mathcal G$ is the space of probability measures on $\mathbb{R}$, is a consistent estimator of $G$. This is an infinite dimensional convex optimization problem with a strictly convex objective subject to linear constraints. See \citeasnoun{Lindsay.95} and \citeasnoun{KM} for further details on the geometry and computational aspects of the NPMLE problem. \citeasnoun{HS} pioneered this approach in econometrics to argue that more flexible models of heterogeneity were needed to get reliable estimates of duration dependence in survival models.
A powerful consequence of the seemingly innocuous condition that $G$ must be non-decreasing is that $\hat G$ must be atomic, a discrete distribution with fewer than $n$ atoms. A secondary consequence is that the NPMLE is “self-regularizing,” that is, the number, locations and mass of the atoms are all determined jointly by the optimization without any recourse to auxiliary tuning parameters. This is all a consequence of the classical Carath\'eodory theorem, but until quite recently little was known about the precise growth rate of the number of atoms characterizing the solutions, although empirical experience suggested it was quite slow. \citeasnoun{PW20} have recently established that for $G$ with sub-Gaussian tails the cardinality of its support, i.e. the number of atoms, of $\hat G$, grows like $\mathcal O(\log n)$. Thus, without any further penalization, maximum likelihood automatically selects a highly parsimonious $\hat G$. This is in sharp contrast to the notorious difficulties with maximum likelihood for finite dimensional mixture models, or with Gaussian deconvolution employing Fourier methods.
Having seen that the upper bound on the complexity of the NPMLE $\hat G$ was only $\mathcal O (\log n)$, one might wonder whether $\mathcal O (\log n)$ mixtures are “complex enough” to adequately represent the process that generated our observed data. \citeasnoun{PW20} also address this concern: they note that for any sub-Gaussian $G$, there exists a discrete distribution, $G_k$, with $k = \mathcal O ( \log n)$ atoms, such that for $f_k = \varphi * G_k$, the total variation distance, $TV(f, f_k) = o (1/n)$, and consequently there is no statistical justification for considering estimators of $G$ whose complexity grows more rapidly than $\mathcal O (\log n)$. This observation is related to recent literature on generative adversarial networks, e.g. \citeasnoun{AIMM20}, that target models and estimators that, when simulated, successfully mimic observed data.
Other nonparametric maximum likelihood estimators for $G$ are potentially also of interest. \citeasnoun{E16} has proposed an elegant log-spline sieve approach that yields smooth estimates of $G$; this has advantages especially from an inferential perspective, at the cost of reintroducing the task of selecting tuning parameters. An early proposal of \citeasnoun{LL91} merged parametric empirical Bayes estimation of $G$ with an EM step that pulled the parametric estimate back toward the NPMLE.
Given an estimate, $\hat G$, it is straightforward to compute posterior distributions for each sample observation, or for that matter, for out-of-sample observations. In effect, we have estimated the prior, as in \citeasnoun{Robbins.51} binary means problem, but we have ignored the variability of $\hat G$ when we adopt plug-in procedures that use it. This may account for the improved performance of smoothed estimates of $G$ in certain inferential problems, as conjectured in \citeasnoun{K20} and studied in more detail in \citeasnoun{JZ21}. In the sequel we will compare ranking and selection procedures based on various functionals of these posterior distributions. A leading example is the posterior mean, but ranking and selection problems suggest other functionals of potential interest.
If we are asked to estimate the $\theta_i$'s subject to quadratic loss, and assuming standard Gaussian noise, the Bayes rule is given by the posterior mean,
\citeasnoun{efron.11} refers to this as Tweedie's formula, it appears in \citeasnoun{Robbins.56} credited to M.C.K. Tweedie. Appendix A of \citeasnoun{Haydn} provides an elementary derivation. The nonlinear shrinkage term takes a particularly simple affine form when $G$ happens to be Gaussian, since in this case $f$ is itself also Gaussian and the formula reduces to well-known linear shrinkage variants of classical Stein rules.
We have focused in this brief overview on compound decision problems for Gaussian location mixtures and posterior means, however the NPMLE is adaptable to a wide variety of other mixture problems and other loss functions that imply other posterior functionals as we will see in the next section. \citeasnoun{E19} and the discussion thereof offers a broader perspective on related methods. Implementation of several NPMLE options are described in \citeasnoun{REBayesVig} and are available in the R package REBayes of \citeasnoun{REBayes}.
Suppose that you are given real-valued measurements, $y_i : i = 1, 2, \cdots , n$ of some attribute like test score performance for students or their teachers, survival rates for hospital surgical procedures, etc., and are told that the measurements are exchangeable and approximately Gaussian with unknown means $\theta_i$ and known variances $\sigma_i^2$ assumed provisionally to take the same value $\sigma^2$. Your task, should you decide to accept it, is to choose a group of size not to exceed $\alpha n$ of the elements with the largest $\theta_i$'s. One's first inclination might be to view each $y_i$ as the maximum likelihood estimate for the corresponding $\theta_i$, and select the $\alpha n$ largest observed values, but the compound decision framework suggests that it would be better to treat the problems as an ensemble. A second natural inclination might be to compute posterior means of the $\theta$'s with some linear or nonlinear shrinkage rule, rank them and select the $\alpha$ best, but we will see that this too may be questionable.
A natural alternative to ranking by the posterior means is to rank by posterior tail probabilities. Let $\theta_\alpha = G^{-1} (1 - \alpha)$, and define, $v_{\alpha}(y):= \mathbb{P}(\theta \geq \theta_\alpha | Y = y)$, then ranking by posterior tail probability gives the decision rule, \[ \delta(y) = \11\{v_\alpha(y) \geq \lambda_\alpha\} \] where the threshold $\lambda_\alpha$ is chosen so that $\mathbb{P}(v_\alpha(Y) \geq \lambda_\alpha) = \alpha$. This ranking criterion has been proposed by \citeasnoun{HN} motivated as a ranking device for a fixed quantile level $\alpha$. It can be interpreted in multiple testing terms: $1-v_{\alpha}(y)$ is the local false discovery rate of \citeasnoun{ETST} and \citeasnoun{storey02}, for testing the hypothesis $H_0: \theta < \theta_\alpha$ vs. $H_A: \theta \geq \theta_\alpha$. To see this, let $h_i$ be a binary random variable $h_i= \11\{\theta_i \geq \theta_\alpha\}$, the loss function for observation $i$ is \[ L(\delta_i, \theta_i) = \lambda \11\{h_i = 0, \delta_i = 1) + \11\{h_i = 1, \delta_i = 0\}, \] for a generic Lagrange multiplier, $\lambda$. The compound Bayes risk is, \[ \mathbb{E}[\sum_{i=1}^n L(\delta_i, \theta_i)] = n [ \alpha + \int \delta(y) [(1-\alpha) \lambda f_0(y) - \alpha f_1(y)] dy] \] where $f_0(y) = (1-\alpha)^{-1} \int_{-\infty}^{\theta_\alpha} \varphi(y|\theta,\sigma^2) dG(\theta)$ and $f_1(y) = \alpha^{-1} \int _{\theta_\alpha}^{+\infty} \varphi(y|\theta, \sigma^2) dG(\theta)$, $\varphi(y| \theta, \sigma^2) = \varphi((y - \theta)/\sigma)/\sigma$. The Bayes rule for a fixed $\lambda$ is \[ \delta(y_i) = \11\Big \{ v_\alpha(y_i) \geq \frac{\lambda}{1+\lambda}\Big \} \] where $v_\alpha(y) =\alpha f_1(y)/f(y) = \mathbb{P}(\theta \geq \theta_\alpha | Y = y)$, and $f(y) = (1 - \alpha) f_0(y) + \alpha f_1(y).$ Provided that $v_\alpha(y)$ is monotone in $y$ a unique $\lambda_\alpha$ can be found such that $\mathbb{P}(\delta(Y) = 1) = \mathbb{P}(v_\alpha(Y)\geq \lambda_\alpha/(1+\lambda_\alpha)) = \alpha$.
Any implementation of such a Bayes rule requires an estimate of the mixing distribution, $G$, or something essentially equivalent that would enable us to compute the local false discovery rates $v_\alpha (y)$ and the cut-off $\theta_\alpha$. The NPMLE, or perhaps a smoothed version of it, will provide a natural $\hat G$ for this task.
Rather than assessing loss by simply counting misclassifications we might consider weighting such misclassifications by the magnitude of $\theta$, for example, \[ L(\delta_i , \theta_i) = \sum_{i=1}^n (1 - \delta_i ) \11 (\theta_i \geq \theta_\alpha) \theta_i. \] This presumes, of course, that we have centered the distribution $G$ in some reasonable way, perhaps by forcing the mean or median to be zero. Minimizing with respect to $\delta$ subject to the constraint that $\mathbb{P}(\delta(Y) = 1) = \alpha$ leads to the Lagrangian, \[ \underset{\delta}{\min} \int \int (1-\delta(y)) \11\{\theta \geq \theta_\alpha\} \theta \varphi(y|\theta, \sigma^2) dG(\theta) dy + \lambda \Big[\int \int \delta(y) \varphi(y|\theta, \sigma^2) dG(\theta)dy - \alpha\Big] \] which is equivalent to
Ignoring the first term since it doesn't depend upon $\delta$, the Oracle Bayes rule, becomes, choose $\delta (y) = 1$ if, \[ \frac{\int \11\{\theta \geq \theta_\alpha \} \theta \varphi(y|\theta,\sigma^2) dG(\theta)}{\int \varphi(y|\theta,\sigma^2) dG(\theta)} \geq \lambda, \] with $\lambda$ chosen so that $\mathbb{P} (\delta(Y) = 1) = \alpha$. Such criteria are closely related to expected shortfall criteria appearing in the literature on risk assessment. Again, the NPMLE can be employed to construct feasible posterior ranking criteria.
Several other loss functions are considered by \citeasnoun{LLPR06} including some based on global alignment of the ranks. While intuitively appealing, such loss functions are considerably less tractable than those we consider in the remainder of the paper.
Although our loss functions yield distinct criteria for ranking, their decision rules lead to the same selections when the precision of the measurements is homogeneous. When variances are homogeneous there is a global cut-off, $\eta_\alpha$, and a decision rule, $\delta_\alpha(Y) = 1(Y \geq \eta_\alpha)$, determining a common selection for all decision rules.
The marginal false discovery rate for selection in our Gaussian mixture setting is, \[ mFDR = \mathbb{P} (\theta < \theta_\alpha | \delta_\alpha (Y) = 1) = \alpha^{-1} \int_{- \infty}^{\theta_\alpha} \Phi ((\theta - \eta_\alpha)/\sigma) dG(\theta), \] The marginal false non-discovery rate is, \[ mFNR = \mathbb{P} (\theta \geq \theta_\alpha | \delta_\alpha (Y) = 0) = (1-\alpha)^{-1} \int_{\theta_\alpha}^\infty \Phi ((\eta_\alpha - \theta)/\sigma) dG(\theta), \]
Figure (ref) shows the false discovery rate, and false non-discovery rate for a range of capacity constraints, $\alpha$, when the mixing distribution, $G$, is standard Gaussian and $\sigma^2 = 1$. In this low signal to noise ratio case, the cut-off value $\eta_\alpha$ is the $(1-\alpha)$ quantile of $\mathcal{N}(0, 2)$, and it is very difficult to distinguish the meritorious from the merely lucky. For selecting individuals at the top $\alpha$ quantile, the false discovery rate is alarmingly high especially for smaller $\alpha$, implying that the selected set may consist of a very high proportion of false discoveries. When $\alpha = 0.10$ the proportion of selected observations with $\theta$ below the threshold $\theta_\alpha$ is slightly greater than 50 percent.
{\tiny{
}}
When the variance of the $\theta$'s, the signal-to-noise ratio, increases from 1 to 5, the selection problem becomes somewhat easier. This is reflected not only by the false discovery rate decreasing from above one half to about a third, but is also reflected in results in Table (ref) that show that the wrongly selected individuals have their true values of $\theta$ clustered closer to the thresholding value $\theta_\alpha$ measured in terms of standard deviation.
When $\sigma^2 = 1$, we have about 50% falsely selected, and the $\theta$'s among the correctly-selected and the wrongly selected individuals are roughly symmetrically distributed around the thresholding value. In such a case, even oracle decision rules that incorporate complete knowledge of the precise distributional features of the problem may not be able to achieve better than about even odds that selected individuals have latent ability above the selection thresholds when measurement error is comparable in magnitude to Gaussian variability in latent ability. As variances of alpha increases to 5 then selection becomes somewhat easier, only 1/3 among the selected are falsely selected, and most (80%) of these $\theta$ values are within 0.5 standard deviations away from the selection threshold.
It is perhaps worth stressing that at the margin, near the decision boundary, it will always be difficult to distinguish true from false discoveries, but FDR measures the proportion of all selections that are incorrect, not just those near the threshold. Other loss functions that penalize in a more continuous way may be considered to reflect information in Table (ref). For example, losses that weight the classification error by the magnitude of the discrepancy between the latent effect and the threshold could be considered. Such losses, however, make it more difficult to incorporate conventional forms of error control.
Thus far we have implicitly assumed that the size of the selected set is predetermined by the parameter $\alpha$. Having established a ranking based on a particular loss function, we simply select a subset of size $\lceil {\alpha n} \rceil$ consisting of the highest ranked observations. In the next subsection we begin to consider modifying this strategy by constraining the probability of false discoveries. This will allow the size of the selected set to adapt to the difficulty of the selection task.
Recognizing the risk of false “discoveries” among those selected, we will consider an expanded loss function,
where $h_i = \11\{\theta_i \geq \theta_\alpha\}$. If we set $\tau_1$ to zero, then minimizing the expected loss leads to the Bayes rule discussed in Section (ref). On the other hand, if we set $\tau_2$ to zero, then minimizing expected loss leads to a decision rule that is equivalent to a multiple testing problem with null hypothesis $H_{0i}: \theta_i \leq \theta_\alpha$; the goal is to minimize the expected number of over-looked discoveries subject to the constraint that the marginal FDR rate is controlled at level $\gamma$, that is, $\mathbb{E}[\sum_{i=1}^n (1-h_i) \delta_i]/ \mathbb{E}[\sum_{i=1}^n \delta_i] \leq \gamma.$ When $\tau_1 = 0$, the risk can be expressed as, \[ \mathbb{E}_{\bm{\theta}|Y}\Big [L(\bm{\delta}, \bm{\theta})\Big ] = \sum_{i=1}^n (1-\delta_i) v_\alpha(Y_i) + \tau_2 \Big( \sum_{i=1}^n \delta_i - \alpha n\Big) \] where $v_\alpha(y_i) = \mathbb{P}(\theta_i \geq \theta_\alpha|Y_i = y_i)$. Taking another expectation over $Y$, and minimizing over both $\bm{\delta}$ and $\tau_2$, leads to the decision rule, \[ \delta_i^* =
\] The Lagrange multiplier is chosen so that the constraint $\mathbb{P}(\delta_i = 1) \leq \alpha$ holds with equality: \[ \tau_2^* = \text{min}\{\tau_2: \mathbb{P}(v_\alpha(y_i) \geq \tau_2) \leq \alpha\} \] Each selection improves the objective function by $v_\alpha(y_i)$, but incurs a cost of $\tau_2$. Since all selections incur the same cost, we may rank according to $v_\alpha(y_i)$, selecting units until the capacity constraint $\alpha n$ is achieved. Selection of the last unit may need to be randomized to exactly satisfy the constraint, as we note below.
When $\tau_2 = 0$ the focus shifts to the marginal FDR, the ratio of the expected number of false discoveries, to the expected number of selections. This is slightly different from the original FDR as defined in \citeasnoun{BH}. However, when $n$ is large the two concepts are asymptotically equivalent as shown by \citeasnoun{GenoWass}. Our objective becomes, \[ \mathbb{E}_{\bm{\theta}|\bm{Y}}\Big[ L(\bm{\delta}, \bm{\theta})\Big] = \sum_{i=1}^n (1-\delta_i) v_\alpha(Y_i) + \tau_1 \Big(\sum_{i=1}^n \{ \delta_i (1- v_\alpha(Y_i)) - \gamma \delta_i \}\Big). \] Taking expectations again over $Y$ and minimizing over both $\bm{\delta}$ and $\tau_1$ yields, \[ \delta_i^* =
\] and the Lagrange multiplier takes a value $\tau_1^*$ to make the marginal FDR constraint hold with equality.
When both constraints are incorporated we must balance the power gain from more selections and the cost that occurs from both the capacity constraint and FDR control. The Bayes rule solves, \[ \min_{\bm{\delta}} \mathbb{E}\Big [\sum_{i=1}^{n} (1-\delta_i) v_{\alpha}(y_i) \Big] + \tau_1 \Big ( \mathbb{E}\Big[ \sum_{i=1}^n \Big \{(1- v_\alpha(y_i)) \delta_i - \gamma \delta_i \Big\}\Big] \Big) + \tau_2 \Big( \mathbb{E}\Big[ \sum_{i=1}^n \delta_i \Big] - \alpha n \Big). \] Given the discrete nature of the decision function, this problem appears to take the form of a classical knapsack problem, however following the approach of \citeasnoun{Basu} we will consider a relaxed version of the problem in which units are selected sequentially until one or the other constraint would be violated, with the final selection randomized to satisfy the constraint exactly.
We have seen that when variances are homogeneous, the optimal selection rule thresholds on $Y$, so it is clear then that any ranking that is based on a monotone transformation of $Y$ will lead to an equivalent selected set. We should also stress that we have focused on a null hypothesis that depends on $\alpha$, while the multiple testing literature, for example \citeasnoun{ETST}, \citeasnoun{SunCai} and \citeasnoun{Basu}, typically focuses on the null hypothesis of $H_{0i}: \theta_i = 0$. When variances are homogeneous, it doesn't matter whether we use an $\alpha$ dependent null or the conventional zero null, because the transformation based on the conventional null, $\mathbb{P}(\theta > 0|Y=y)$, is also a monotone function of $Y$, and therefore yields an equivalent decision rule. However, when variances are heterogeneous, this invariance no longer holds; different transformations of the pair $(y, \sigma)$ lead to distinct decision rules that lead to distinct performance, and using the conventional null hypothesis is no longer advisable for the ranking and selection problem as we will show in the next section.
The homogeneous variance assumption of the preceding section is unsustainable in most applications. Batting averages are accompanied by a number of “at bats” and mean test score performances are accompanied by student sample sizes. In this section we will consider the expanded model, \[ Y_i \sim \mathcal{N}(\theta_i, \sigma_i^2), \quad \text{and } \quad \theta_i \sim G, \quad \sigma_i \sim H, \quad \sigma_i \perp\!\!\!\perp \theta_i \] We will assume that we observe $\sigma_i$, an assumption that will be relaxed in the next section.
With the same alternative hypothesis as $H_{A}: \theta \geq \theta_\alpha$, it is natural to consider the posterior tail probability again, now as a function of the pair, $(y_i, \sigma_i)$, \[ v_\alpha(y_i,\sigma_i) = \mathbb{P}(\theta_i \geq \theta_\alpha | y_i, \sigma_i) = \frac{\int_{\theta_\alpha}^{+\infty} \varphi(y_i | \theta,\sigma_i^2) dG(\theta)}{\int_{-\infty}^{+\infty} \varphi(y_i | \theta,\sigma_i^2)dG(\theta)}. \]
Solving the same decision problem with the loss function specified in ((ref)), we have the conditional risk, \[ \mathbb{E}_{\bm{\theta}|\bm{Y}, \bm{\sigma}}\Big[ L(\bm{\delta}, \bm{\theta})\Big] = \sum_{i=1}^n (1- \delta_i) v_\alpha(Y_i, \sigma_i) + \tau_1 \Big( \sum_{i=1}^n \{\delta_i(1- v_\alpha(Y_i, \sigma_i)) - \gamma \delta_i \}\Big) + \tau_2 \Big( \sum_{i=1}^n \delta_i - \alpha n\Big). \] Taking another expectation with respect to the joint distribution of the $(Y_i, \sigma_i)$'s, the Bayes rule solves \[ \min_{\bm{\delta}} \mathbb{E}\Big [\sum_{i=1}^{n} (1-\delta_i) v_{\alpha}(y_i, \sigma_i ) \Big] + \tau_1 \Big ( \mathbb{E}\Big[ \sum_{i=1}^n \Big \{(1- v_\alpha(y_i, \sigma_i )) \delta_i - \gamma \delta_i \Big\}\Big] \Big) + \tau_2 \Big( \mathbb{E}\Big[ \sum_{i=1}^n \delta_i \Big] - \alpha n \Big) \] The optimal selection rule can again be characterized as a thresholding rule on $v_\alpha(y_i, \sigma_i)$.
The posterior tail probability criterion is motivated by viewing the ranking and selection problems as hypothesis testing while allowing the null hypothesis to be $\alpha$ dependent. The particular construction of the null hypothesis turns out to be critical for the ranking exercise. In this subsection we present a simple example to illustrate that tail probability based on the conventional null hypothesis of zero effect does not lead to a powerful ranking device. Consider data generated from a three component normal mixture model,
Instead of transforming the data by $v_\alpha$, we consider the transformation, \[ \texttt{T}(y_i,\sigma_i^2) = \mathbb{P}(\theta_i > 0| y_i ,\sigma_i^2) = \frac{\int_0^{+\infty} \varphi(y_i | \theta,\sigma_i^2) dG(\theta)}{\int_{-\infty}^{+\infty} \varphi(y_i | \theta,\sigma_i^2) dG(\theta)} \] and rank individuals accordingly. This transformation corresponds to the procedure proposed in \citeasnoun{SM}, and is motivated for multiple testing problems under the conventional null hypothesis $H_{0}: \theta \leq 0$. The decision rule $\delta_i^{\texttt{T}} = \11\{\texttt{T}(y_i,\sigma_i) \geq \lambda\}$ then chooses the cutoff value $\lambda$ that respects both the capacity constraint and the FDR control constraint for selecting the top $\alpha$ proportion.
Figure (ref) compares the selection region for the two ranking procedures with $\alpha = 5\%$ and marginal FDR control at level $10\%$. The solid black line corresponds to the selection boundary using ranking based on transformation $v_\alpha$ and the dashed red line corresponds to the selection boundary using ranking based on the transformation $\texttt{T}$. The black highlighted area below the black selection boundary corresponds to a region where the ranking method based on $\texttt{T}$ will select but the ranking method based on $v_\alpha$ does not. On the other hand, the blue highlighted area corresponds to a region selected by $v_\alpha$, but not for $\texttt{T}$. The transformation $\texttt{T}$ ranks those in the black region higher than those in the blue region because although they have a relatively smaller mean effect $y$, their associated variances are also smaller indicating stronger evidence that such individuals have a positive $\theta$ than those located in the blue area. However, our task is to find individuals with true effects, $\theta_i$, in the upper tail. For $\alpha = 5\%$, we aim to select all individuals with $\theta=5$, individuals in the black region present strong evidence that their true effect can not be too large because their observed effect $y$ is small and their associated variance is also small, while those in the blue region, although their observed mean effects are associated with larger variances offer reasonable evidence that their associated true effect $\theta$ may be large. This evidence is not apparent in transformation $\texttt{T}$, but is captured in the transformation $v_\alpha$.
Indeed, the average power of ranking based on the two different transformation $v_\alpha$ and $\texttt{T}$ differ significantly. Defining the power of the selection rule as $\beta(\delta) := \mathbb{P}(\theta_i \geq \theta_\alpha, \delta_i = 1)/\mathbb{P}(\theta_i \geq \theta_\alpha)$, the proportion of true top $\alpha$ cases selected based on decision rule $\bm{\delta}$, then $\beta(\bm{\delta^{\texttt{T}}}) = 39\%$ and $\beta(\bm{\delta^*}) = 69\%$. Thus, although much of the literature relies on ranking and selection rules based on some form of posterior means and conventional hypothesis testing apparatus we would caution that such methods can be quite misleading and inefficient.
If we were to relax the capacity constraint to allow a larger proportion, $\alpha_1 > \alpha_0$ be be selected, while maintaining our initial false discovery control, we would expect that members selected under the more stringent capacity constraint should remain selected under the relaxed constraint. We now discuss sufficient conditions under which we obtain this nestedness of the selection sets when using the posterior tail probability rule. This is a natural condition in applications like our analysis of ranking and selection of dialysis centers especially because we would like to assign “letter grades” to several several subgroups of the centers.
The optimal Bayes rule defines the selection set for each pair of $(\alpha, \gamma)$ as \[ \Omega_{\alpha, \gamma} := \{(y, \sigma): v_\alpha(y, \sigma) \geq \lambda^*(\alpha,\gamma)\} \] and when $\sigma$ is known, $v_\alpha(y,\sigma)$ is monotone in $y$ as shown in Lemma (ref) for each fixed $\sigma$, hence the selection set can also be represented as \[ \Omega_{\alpha, \gamma} =\{(y,\sigma): y \geq t_{\alpha}(\lambda^*(\alpha,\gamma),\sigma)\} \] It is also convenient for later discussion to define
which are respectively the selection sets when the false discovery rate constraint or the capacity constraint is binding. It is easy to see that $\Omega_{\alpha,\gamma} = \Omega_{\alpha,\gamma}^{FDR} \cap \Omega_{\alpha}^C$.
In this section we consider several examples beginning with the simplest classical case in which the $\theta_i$ constitute a random sample from the standard Gaussian distribution. This Gaussian assumption on the form of the mixing distribution $G$ underlies almost all of the empirical Bayes literature in applied economics; it is precisely what justifies the linear shrinkage rules that are typically employed.
Assuming that the $\sigma_i$'s are known, up to a common scale parameter, may be plausible in some applications such as baseball batting averages, but it is frequently more plausible to adopt the view that we are simply confronted with estimates of scale available perhaps from longitudinal data. In such cases we need to consider the pairs, $(y_i , S_i)$ as potentially jointly dependent random variables arising from the longitudinal model, \[ Y_{it} = \theta_i + \sigma_i \epsilon_{it}, \quad \epsilon_{it} \sim_{iid} \mathcal{N}(0,1), \quad (\theta_i, \sigma_i^2) \sim G, \] with sufficient statistics, $Y_i = T_i^{-1} \sum_{t=1}^{T_i} Y_{it}$ and $S_i = (T_i-1)^{-1} \sum_{t=1}^{T_i} (Y_{it} - Y_i)^2$, for $(\theta_i,\sigma_i^2)$. Conditional on $(\theta_i,\sigma_i^2)$, we have $Y_i |\theta_i,\sigma_i^2 \sim \mathcal{N}(\theta_i, \sigma_i^2/T_i)$ and $S_i |\sigma_i^2$ is distributed as Gamma with shape parameter $r_i = (T_i-1)/2$, scale parameter, $\sigma_i^2/r_i$, and density function denoted as $\Gamma(S_i | r_i,\sigma_i^2/r_i)$.
Given the loss function ((ref)) and defining $\theta_\alpha$ as $\alpha = \mathbb{P}(\theta_i \geq \theta_\alpha) = \int \int_{\theta_\alpha}^{+\infty} dG(\theta, \sigma^2)$, the conditional risk is,
with
Taking expectations with respect to $(Y, S)$, the Bayes rule solves, \[ \underset{\bm{\delta}}{\min} \mathbb{E} \Big [ \sum_{i=1}^n (1-\delta_i) v_\alpha(y_i,s_i)\Big] + \tau_1 \Big( \mathbb{E}\Big[ \sum_{i=1}^n \Big\{(1-v_\alpha(y_i,s_i))\delta_i - \gamma \delta_i\Big\}\Big]\Big) + \tau_2 \Big (\mathbb{E}\Big[ \sum_{i=1}^n \delta_i\Big] - \alpha n\Big). \]
Before characterizing the Bayes rule any further, we should observe that when variances $\sigma^2$ are not directly observed, the tail probability $v_\alpha(Y,S)$ may no longer have the monotonicity property we have described above.
Based on the Bayes rule, the selected set is defined as \[ \Omega_{\alpha, \gamma} = \{(Y,S): v_\alpha(Y,S) \geq \lambda^*(\alpha,\gamma)\}. \]
Suppose we have balanced panel data $y_{i1}, \dots, y_{iT} \sim \mathcal{N}(\theta,\sigma^2)$ with sample means $Y_i = \frac{1}{T}\sum_t y_{it}$ and sample variances $S_i = \frac{1}{T-1}\sum_t (y_{it} - Y_i)^2$. Further, suppose that $G(\theta,\sigma^2)$ takes the normal-inverse-chi-squared form, $\mbox{NIX} (\theta_0, \kappa_0, \nu_0, \sigma_0^2) = \mathcal{N}(\theta | \theta_0, \sigma^2/\kappa_0) \chi^{-2}(\sigma^2 | \nu_0, \sigma_0^2)$. Integrating out $\sigma^2$, the marginal distribution of $\theta$ becomes a Student $t$ distribution, \[ \frac{\theta - \theta_0}{\sigma_0/\sqrt{\kappa_0}} \sim t_{\nu_0} \] where $t_{\nu_0}$ is the $t$-distribution with degree of freedom $\nu_0$. Therefore, the $1-\alpha$ quantile of $\theta$, denoted $\theta_\alpha$ is simply, \[ \theta_\alpha =\theta_0 + \frac{\sigma_0}{\sqrt{\kappa_0}} F^{-1}_{t_{\nu_0}}(1-\alpha) \] where $F^{-1}_{t_{\nu_0}}$ denotes the quantile function of $t_{\nu_0}$.
Conjugacy of the distribution $G$ implies that the posterior distribution of $(\theta, \sigma^2|Y,S)$ follows $NIX(\theta_T, \kappa_T, \nu_T, \sigma_T^2) = \mathcal{N}(\theta | \theta_T, \sigma_T^2/\kappa_T) \chi^{-2}(\sigma^2| \nu_T, \sigma_T^2)$ with
Integrating out $\sigma^2$, the marginal posterior of $\theta$ again follows a $t$-distribution, \[ \frac{\theta - \theta_T}{\sigma_T/\sqrt{\kappa_T}} \sim t_{\nu_T}. \] It is thus clear that the posterior mean of $\theta$, is simply a linear function of $Y$ and independent of $S$, \[ \mathbb{E}[\theta | Y,S] = \theta_T = \frac{\kappa_0 \theta_0 + T Y}{\kappa_T}, \] and the posterior tail probability is given by, \[ v_{\alpha}(Y,S) = \mathbb{P}(\theta \geq \theta_\alpha|Y,S) = \mathbb{P}\Big(\frac{\theta - \theta_T}{\sigma_T/\sqrt{\kappa_T}} \geq \frac{\theta_\alpha - \theta_T}{\sigma_T/\sqrt{\kappa_T}} | Y,S\Big) = 1- F_{t_{\nu_T}}\Big( \frac{\theta_\alpha - \theta_T}{\sigma_T/\sqrt{\kappa_T}}\Big). \]
To illustrate this case, suppose $\theta_0=0$, $\kappa_0 = 1$ , $\sigma_0^2 = 1$ and $\nu_0 = 6$ and $T=9$, it can be verified that $v_\alpha(Y,S)$ is in fact a monotone function of $Y$ for each fixed $S$ and any $\alpha > 0$, hence in this example we can invert the function $v_\alpha(y,s)$ to obtain the level curves. The left panel of Figure (ref) shows the level curves for $v_\alpha(Y,S)$ and $\mathbb{E}(\theta|Y,S)$ for $\alpha = 5\%$. It is clear that the posterior mean is a constant function of $S$, while the posterior tail probability exhibits more exotic behaviour with respect to $S$, especially for more extreme values of $Y$. If we fix $S = s_0$, then $v_\alpha(Y,s_0)$ is an increasing function of $Y$. On the other hand, fixing $Y = y_0$, for small $y_0$ implies that $v_\alpha(y_0,S)$ is a increasing function of $S$, while for $y_0$ large, $v_\alpha(y_0,S)$ becomes a decreasing function of $S$.
A capacity constraint of size $\alpha$ implies the thresholding rule, \[ \mathbb{P}(v_\alpha(Y,S)\geq \lambda_2^*) = \alpha, \] while FDR control at level $\gamma$ leads to a cutoff value $\lambda_1^*$ defined as \[ \gamma = \mathbb{E}[(1-v_\alpha(Y,S) )\11\{v_\alpha(Y,S)\geq \lambda_1^*\}]/\mathbb{P}(v_\alpha(Y,S) \geq \lambda_1^*). \] The larger of the two thresholds, denoted $\lambda^* = \max \{\lambda_1^*, \lambda_2^*\}$ defines the selection region based on posterior tail probability ranking $\Omega_{\alpha,\gamma}= \{(Y,S): v_\alpha(Y,S) \geq \lambda^*\}$. For $\alpha = 5\%$ and $\gamma = 10\%$, the selection region based on the tail probability rule is $ \{(Y,S): v_\alpha(Y,S)\geq 0.72\}$. The posterior mean ranking is defined as $\{(Y,S): \mathbb{E}[\theta|Y,S]\geq 2.2\}$. These selection boundaries are depicted as the red dashed line and black solid line respectively in the right panel of Figure (ref). In this case, the FDR constraint binds. If only the capacity constraint were in place, we would have a cutoff for tail probability at $0.40$ and the cutoff for the posterior mean at $1.84$. Figure (ref) further shows the comparison of the selected set based on a sample realization from the model.
In Appendix (ref) we consider a more complex bivariate discrete example that illustrates somewhat more exotic behavior of the decision boundaries and compares performance of several different ranking and selection rules.
We have assumed that the only scale heterogeneity is driven by $\sigma_i$ in the above model, but often there may be more heteroskedasticity that should be allowed in $\epsilon_{it}$. Here we consider a variant where \[ Y_{it} = \theta_i + \sigma_i \epsilon_{it}, \quad \epsilon_{it} \sim \mathcal{N}(0, 1/w_{it}), \quad (\theta_i, \sigma_i^2) \sim G. \] We will assume that $w_{it} \sim H$ are known quantities and are independent from $(\theta_i, \sigma_i^2)$. Denoting $w_i = \sum_{t=1}^{T_i} w_{it}$, the sufficient statistics now take the form $Y_i = \sum_{t=1}^{T_i} w_{it} Y_{it}/w_i$ and $S_i = (T_i-1)^{-1} \sum_{t=1}^{T_i} (Y_{it} - Y_i)^2$. In \citeasnoun{Bayesball} we have illustrated this formulation for predicting baseball batting averages; in that setting “at bats” for player $i$ in year $t$ are given by the $w_{it}$, but there is still some player specific heterogeneity in the $\sigma_i$'s representing consistency of batting performance. It is easy to show that $Y_i |\theta_i, \sigma_i^2 \sim \mathcal{N}(\theta_i, \sigma_i^2/w_i)$ and $S_i|\sigma_i^2$ follows Gamma distribution with shape parameter $r_i = (T_i-1)/2$ and scale parameter $\sigma_i^2/r_i$. The decision rules now become a function of the tuple $(Y_i, S_i, w_i)$, for instance the tail probability can be specified as \[ v_\alpha(y_i, s_i, w_i) = \mathbb{P}(\theta_i \geq \theta_\alpha| y_i, s_i, w_i) = \frac{\int_{\theta_\alpha}^{+\infty} f(y_i | \theta, \sigma^2/w_i) \Gamma( s_i | r_i, \sigma^2/r_i) dG(\theta, \sigma^2) } {\int_{-\infty}^{+\infty} f(y_i | \theta, \sigma^2/w_i) \Gamma( s_i | r_i, \sigma^2/r_i) dG(\theta, \sigma^2) }. \] and the posterior mean takes the form \[ \mathbb{E}[\theta_i | y_i,s_i,w_i] = \int \theta f(y_i, s_i|\theta,\sigma^2, w_i) dG(\theta,\sigma^2). \] The threshold values under either the capacity or the FDR constraint can be worked out in a similar fashion. For any ranking statistics $\delta(Y_i,S_i,w_i)$ together with a decision rule $\11\{\delta(Y_i,S_i,w_i) \geq \lambda\}$, the capacity constraint requires choosing a thresholding value $\lambda_2^*(\alpha)$ such that \[ \alpha = \int \int \11\{\delta(y, s, w) \geq \lambda_2^*(\alpha)\} f(y, s| \theta, \sigma^2, w) dG(\theta, \sigma^2) dH(w), \] while the thresholding value to control in addition the FDR rate under size $\gamma$ requires solving for $\lambda_1^*(\alpha,\gamma)$ such that \[ \gamma = \frac{\mathbb{P}(\delta(y,s,w) \geq \lambda_1^*(\alpha, \gamma); \theta < \theta_\alpha)}{\mathbb{P}(\delta(y,s,w) \geq \lambda_1^*(\alpha, \gamma))} \] which can be further represented as \[ \gamma = \frac{\int \int \11\{\delta(y, s, w) \geq \lambda_1^*(\alpha, \gamma)\} (1-\alpha) f_0(y, s|\theta, \sigma^2, w) dG(\theta, \sigma^2) dH(w)} {\int \int \11\{\delta(t, s, w) \geq \lambda_1^*(\alpha, \gamma)\} f(y, s|\theta, \sigma^2, w) dG(\theta, \sigma^2) dH(w)} \] where $f_0(y, s|\theta, \sigma^2, w)$ is the density of $(y, s)$ under the null hypothesis $\theta < \theta_\alpha$.
We can again consider selection region as those plotted in Figure (ref) and Figure (ref) to appreciate how different decision criteria determines the selection. As soon as the ranking statistics depends on $w$, the selection region of the thresholding rule $\11\{\delta(y,s,w) \geq \lambda^*\}$ will also depend on the magnitude of $w$.
The previous sections propose Bayes rules for minimizing the expected number of missed discoveries subject to both capacity and FDR constraints under several modeling environments. In each of these environments, the Bayes rule takes the form $\delta^* = 1\{v_\alpha \geq \lambda^*\}$ where $v_\alpha$ is defined as the posterior probability of $\theta \geq \theta_\alpha$ conditional on the data. The thresholding value $\lambda^*$ is defined to satisfy both the capacity and FDR constraints. The Bayes rule involves several unknown quantities, in particular the $v_\alpha$'s and the optimal thresholding value, $\lambda^*$, that require knowledge on the distribution of $\theta_i$ or the joint distribution of $(\theta_i, \sigma_i^2)$ when the variances are latent variables. For estimating this distribution of the latent variables, we propose a plug-in procedure that is very much in the spirit of empirical Bayes methods pioneered by \citeasnoun{Robbins.56}. In this Section we also establish that the resulting feasible rules achieve asymptotic validity and asymptotically attain the same performance as the infeasible Bayes rule.
We begin by discussing properties of the Oracle procedure assuming that $v_\alpha$ is known and we only need to estimate the optimal thresholding value. We establish asymptotic validity of this Oracle procedure and then propose a plug-in method for both $v_\alpha$ and the thresholding value thereby establishing the asymptotic validity of the empirical rule. Before presenting the formal results, we introduce regularity conditions that will be required. We distinguish two cases depending on whether the $\sigma_i^2$'s are observed.
Whether $\sigma_i^2$ is observed or estimated, the optimal thresholding value can be defined in a unified manner by $\lambda^* = \max \{\lambda_1^*, \lambda_2^*\}$ with
where $H_v$ denotes the cumulative distribution of either $v_\alpha(y_i, \sigma_i)$ or $v_\alpha(y_i, s_i)$, induced by the marginal distribution of the data, either as the pair $\{y_i, \sigma_i\}$ when variances are observed or the pair $\{y_i, s_i\}$ otherwise. Hence $\lambda_1^*$ is the $1-\alpha$ quantile of $H_v$.
The function $Q(t)$ is defined as $Q(t) = \mathbb{E}[(1-v_\alpha) 1\{v_\alpha \geq t\}] / \mathbb{E}[1\{v_{\alpha} \geq t\}]$. Its formulation recalls Proposition (ref) and the existence of $\lambda_2^*$ is guaranteed as long as $\alpha < 1- \gamma$. The thresholding value is also equivalent to those defined in Proposition (ref) and Proposition (ref). In particular, the thresholding values $t_1^*$ and $t_2^*$ in Proposition (ref) are cast in terms of $Y$ directly and it is easy to see $\lambda_j^* = v_\alpha(t_j^*)$ for $j = 1,2$ when variances are homogeneous. In a similar spirit, the explicit formulae for $\lambda_1^*$ and $\lambda_2^*$ in Proposition (ref) are a result of invoking the monotonicity of $v_\alpha(y, \sigma)$ with respect to $y$ for each fixed value of $\sigma$. The function $Q(t)$ is the mFDR of the procedure $\delta= 1\{v_\alpha \geq t\}$ for any $\alpha\in (0,1)$, and is monotonically decreasing in $t$. Monotonicity of $Q(t)$ is crucial to justify this thresholding procedure insuring that either the capacity constraint or the mFDR constraint must be binding. \citeasnoun{cao2013optimal} have observed that a sufficient condition for monotonicity for a broad class of multiple testing procedures is that the ratio of densities under the null and alternative of the test statistics employed for ranking be monotone and they discuss the consequences of the violation of this condition. For the posterior tail probability criterion this monotone likelihood ratio condition, as we will see, can be verified directly.
Recall that mFDR is defined as $\sum_{i=1}^n \mathbb{P}[\delta_i = 1, \theta_i < \theta_\alpha] / \sum_{i=1}^n \mathbb{P}(\delta_i = 1)$. It suffices to show that $\mathbb{P}[\delta_i = 1, \theta_i < \theta_\alpha] = \mathbb{E}[(1-v_{\alpha,i}) \delta_i]$. Since $v_{\alpha,i} = P[\theta_i \geq \theta_\alpha |D_i] = \alpha f_1(D_i) / f(D_i)$ where $D_i$ is the individual data being either $\{y_i, \sigma_i\}$ or $\{y_i, s_i\}$ depending on the model and $f_1$ is the marginal density of the data when $\theta_i \geq \theta_\alpha$ and $f$ is the marginal density of $D_i$. Then it is easy to see that $\mathbb{P}[\delta_i = 1, \theta_i < \theta_\alpha] = (1-\alpha) \int 1\{v_{\alpha,i} \geq t \} f_0(D_i) dD_i = \int 1\{v_{\alpha, i} \geq t\} (1-v_{\alpha,i}) f(D_i) dD_i = \mathbb{E}[(1-v_{\alpha,i}) 1\{v_{\alpha,i} \geq t\}]$. Then $Q(t) = \int_{t}^1 (1- v) h_v dv/\int_{t}^1 h_v dv$ where $h_v$ is the density function of $v_\alpha$. Monotonicity of $Q(t)$ can be easily verified by showing that the derivative with respect to $t$ of the right hand side quantity is nonpositive.
The only unknown quantity in the Oracle procedure is the thresholding value and we now discuss how to estimate it to achieve asymptotic validity. $H_v$ and $Q$ can be estimated by the following quantities:
and the associated thresholding values are then defined as $\lambda_n = \max\{ \lambda_{1n}, \lambda_{2n}\}$, with
In practise the posterior tail probability also involves the unknown quantity $\theta_\alpha = G^{-1} (1 - \alpha)$ that needs to be estimated. We propose a plug-in estimator in the spirit of the empirical Bayes method: estimating $G$ by its nonparametric maximum likelihood estimator $\hat G_n$ and estimating $\theta_\alpha$ as its $1-\alpha$ quantile.
Consistency of the nonparametric maximum likelihood estimator, $\hat G_n$, was first proven by \citeasnoun{KW} using Wald type arguments. A Hellinger risk bound for the associated marginal density estimate and adaptivity of $\hat G_n$ and a self-regularization property have been recently established in \citeasnoun{saha2020nonparametric} and \citeasnoun{PW20}. In particular, the following established result, stated here as an assumption, is crucial for establishing the asymptotic validity of the adaptive procedure.
When variances are homogeneous or when variances are unknown but we have longitudinal data so that we have a mixture model for the pair $\{Y_i, S_i\}$, the Hellinger convergence is established in \citeasnoun{van1993hellinger}. When variances are heterogeneous but known, the Hellinger bound for marginal density has been established recently in \citeasnoun{jiang2020general}.
The plug-in estimators for the posterior tail probability, $v_\alpha(y_i, \sigma_i)$ when variances are known or $v_\alpha(y_i, s_i)$ when variances are unknown is then defined respectively as
where $f$ is the density function for $(y_i, s_i)$ which is a product of Gaussian and gamma densities. Abbreviating the estimated posterior tail probability by $\hat v_{\alpha,i}$, we mimic the Oracle procedure and estimate the thresholding value by $\hat \lambda_n = \max \{\hat \lambda_{1n}, \hat \lambda_{2n}\}$, where,
It is clear that given the Lagrangian formulation of the compound decision problem, it can be viewed equivalently as a constrained optimization problem. See also the discussion in Remark (ref). We seek to maximize power defined as $\beta(t) := \mathbb{P}(\theta_i \geq \theta_\alpha, \delta_i = 1)/\alpha $ subject to two constraints: the first is the marginal FDR rate and the other is the selected proportion. For each fixed pair of $\{\alpha, \gamma\}$, the Bayes rule achieves the best power among all thresholding procedures that respect the two constraints. The next theorem establishes that our feasible, adaptive procedure achieves the same power as the oracle rule asymptotically. It is supported by the simulation evidence presented in the next section. In practice we suggest convolution smoothing of the discrete $\hat G$ as in \citeasnoun{JZ21} with a bandwidth slowly tending to zero with $n$. As they show, the smoothed mixing distribution is also consistent, hence fulfilling Assumption (ref) and therefore all our adaptivity results.
In this Section we describe two small simulation exercises designed to illustrate performance of several competing methods for ranking and selection. As a benchmark for evaluating performance we consider several Oracle methods that presume knowledge of the true distribution, $G$, generating the $\theta$'s as well as several feasible methods that rely on estimation of $G$. These are contrasted with more traditional methods that are based on linear shrinkage rules of the Stein type. The linear shrinkage rule is the posterior mean of $\theta$ under the assumption that $G$ follows a Gaussian distribution with unknown mean and variance parameters. This is the commonly used estimator for ranking and selection in applied work, notably Chetty, Friedman and Rockoff (2014a, 2014b) for teacher evaluation and \citeasnoun{chetty2018impacts} for studying intergenerational mobility.
Typically the linear shrinkage estimator is used in the context of heterogeneous known variances, this will be the model we focus on in our simulation experiments. The linear shrinkage formula defined in (ref) easily adapts to the heterogeneous variances case and leads to the Jame-Stein shrinkage rule with heterogeneous known variances. \citeasnoun{EfronMorris} introduced some further modifications. As we have already demonstrated, when variances are heterogeneous, the linear shrinkage estimator provides a different ranking than the posterior tail probability rules. Further complications arise when we seek procedures that also control false discovery. To estimate the false discovery rate for different thresholding values we requires knowledge of $G$. If the Gaussian assumption on $G$ underlying the linear shrinkage rules is misplaced, it may lead to an inaccurate estimates of FDR, and consequently to procedures that fail to control for false discovery.
Performance will be evaluated primarily on the basis of power, which we define as the proportion of individuals whose true $\theta_i$ exceeds the cutoff $\theta_\alpha = G^{-1} (1 - \alpha)$, who are actually selected. This is the sample counterpart of $\mathbb{P}(\delta_i = 1, \theta_i \geq \theta_\alpha)/\mathbb{P}(\theta_i \geq \theta_\alpha)$. FDR is calculated as the sample counterpart of $\mathbb{P}(\delta_i = 1, \theta_i < \theta_\alpha)/\mathbb{P}(\delta_i = 1)$, that is the proportion of selected individuals whose true $\theta_i$ falls below the threshold. While our selection rules are {\it intended} to constrain FDR below the $\gamma$ threshold, as in other testing problems they are not always successful in this objective in finite samples so empirical power comparisons must be interpreted cautiously in view of this. Nonetheless, asymptotic validity is assured by the results in Section (ref). We compare performance for three distinct $\alpha$ levels, $\{ 0.05, 0.10, 0.15 \}$ and three $\gamma$ levels $\{ 0.05, 0.10, 0.15 \}$.
Our first simulation setting focuses on the effect of tail behavior of the distribution on performance of competing rules. For these simulations we take $G$ to be a discrete approximation to Student $t$ distributions with degrees of freedom in the set $\{ 1,2,3,5,10\}$, and supported on the interval $[-20,20]$. The scale parameters of the Gaussian noise contribution are independent and uniformly distributed on the interval $[0.5,1.5]$. We report power performance for several alternative ranking and selection rules:
The KWs rule uses $\tilde G = \hat G * K_h$, with biweight kernel $K$ and bandwidth $h$ equal to half the mean absolute deviation from the median of $\hat G$. The Efron rule uses his suggested default of a natural spline basis with five degrees of freedom and penalty parameter 0.1.
We illustrate the results in Figure (ref), where we plot empirical power against degrees of freedom of the $t$ distribution for a selected set of values for the capacity constraint, $\alpha \in \{ 0.05, 0.10, 0.15 \}$ and FDR constraint, $\gamma \in \{ 0.05, 0.10, 0.15 \}$ as indicated at the top of each panel of the figure. The most striking conclusion from this exercise is the dramatic decrease in power as we move toward the Gaussian distribution. At the Cauchy, $t_1$, power is quite respectable for all choices of $\alpha$ and $\gamma$, but power declines rapidly as the degrees of freedom increases, reenforcing our earlier conclusion that the Gaussian case is extremely difficult. We would stress, in view of this finding, that classical linear shrinkage procedures designed for the Gaussian setting are poorly adapted to heavy tailed settings in which the reliability of selection procedures is potentially greatest.
Careful examination of this figure also reveals that there is a slight advantage to the posterior tail probability rules over the posterior mean procedures, both for the Oracle rules and for our feasible procedures. There is surprisingly little sacrifice in power in moving from the Oracle methods to the Efron or Kiefer-Wolfowitz rules. The Efron and Morris selection rule is very competitive in the almost Gaussian, $t_{10}$ setting but sacrifices considerable power in the lower degrees of freedom settings due to the misspecification of the distribution $G$ and consequent inaccurate estimation of the false discovery rate.
Our second simulation setting is based on a discrete approximation of the data structure employed in \citeasnoun{GGM} to study teacher value-addded methods. Several longitudinal waves of student test scores from the Los Angeles Unified School District were combined in this study. Here we abstract from many features of the full longitudinal structure of this data, and focus instead on comparing performance of several selection methods. We maintain our standard known variance model in which we observe $Y_i \sim \mathcal N ( \theta_i, \sigma_i^2)$ with $\theta_i$'s drawn iidly from a distribution $\tilde G$ estimated by \citeasnoun{GGM}. This distribution was estimated from the full longitudinal LA sample using the nonparametric maximum likelihood estimator of Kiefer and Wolfowitz and then smoothed slightly by convolution with a biweight kernel and illustrated in the left panel of Figure (ref). Variances, in keeping with our hypothesis in Section (ref), are drawn from a distribution with density illustrated in the right panel of Figure (ref). We focus on selection from the left tail of the resulting distribution since it is those teachers whose jobs are endangered by recent policy recommendations in the literature. (see for instance \citeasnoun{hanushek2011economic}).
We draw samples of size 10,000 from the foregoing distribution and compute performance measures based on 100 replications. The fitted densities for this simulation exercise are based on a sample of roughly 11,000 teachers, so the simulation sample size is chosen to be commensurate with this. In Table (ref) we report power, FDR and the proportion selected by ten selection rules. The Oracle rules, OTP and OPM, based ranking by the tail probability and posterior mean criteria can be considered benchmarks for the remaining feasible procedures. Only the Oracle procedures can be considered reliable from the perspective of adhering to the capacity and FDR constraints. Consequently, some caution is required in the interpretation of the power comparisons since feasible procedures can exhibit good power at the expense of violating these constraints. This is analogous to the common difficulty in interpreting power in testing problems when different procedures have differing size. When FDR is constrained to 5%, even the Oracle is only able to select about half of the deserving individuals; OTP is consistently preferable to OPM as expected and power performance improves somewhat as the capacity constraint is relaxed. Among the feasible $G$-modeling selection procedures, the Efron rules have good power performance, but fail to meet the FDR constraints. We conjecture that somewhat less aggressive smoothing than the default, $df = 5, c_0 = 0.1$ might help to rectify this. In contrast the smoothed Kiefer-Wolfowitz rules are somewhat overly conservative in meeting the FDR constraints and might benefit from somewhat more aggressive smoothing.
\nocite{Chettya,Chettyb} Among the other procedures the linear posterior mean rule, LPM, as employed by Chetty, Friedman and Rockoff (2014a, 2014b), and the linear posterior mean rule, EM, of \citeasnoun{EfronMorris} behave identically and exhibit somewhat erratic FDR control due to the misspecified Gaussian assumption on $G$; this leads to weaker power performance. As a further comparison, when the linear shrinkage rules are implemented without any FDR constraint, denoted LPM* and EM* in the table, as they typically would be used in practice, the false discovery proportion is considerably higher than the targeted $\gamma$. We also report the performance of MLE and P-value rules, implemented without FDR control; again both yield a higher FDR rate, making it difficult to evaluate their power performance.
{\tiny{
}}
Motivated by important prior work on ranking and selection by Lin, Louis, Paddock and Ridgeway (2006, 2009) \nocite{LLPR06,LLPR09} illustrated by applications to ranking U.S. dialysis centers, we have chosen to maintain this focus to illustrate our own approach. Kidney disease is a growing medical problem in the U.S and considerable effort has been devoted to data collection and evaluation of the relative performance of the more than 6000 dialysis centers serving the afflicted population. Centers are evaluated on multiple criteria, but the primary focus of center ranking is their standardized mortality rate, or SMR, the ratio of observed deaths to expected deaths for center patients. Allocating patients to centers is itself a complex task since patients may move from one center to another in the course of a year. Centers also vary considerably in the mix of patients they serve. Predictions from an estimated Cox proportional hazard model that attempts to account for this heterogeneity are employed to estimate expected deaths for each center.
Our analysis focuses exclusively on the SMR evaluation of centers using longitudinal data from 2004-18 as reported in \citeasnoun{DataSource}. We restrict attention to 3230 centers that have consistently reported SMR data over this sample period. Observed deaths, denoted $y_{it}$ for center $i$ in year $t$ are conventionally modeled as Poisson, \[ y_{it} \sim \mbox{Pois} (\rho_i \mu_{it}) \] where $\mu_{it}$ is center $i$'s expected deaths as predicted by the Cox model in year $t$ and $\rho_i$ is the center's unobserved mortality rate. We view $\mu_{it}$ as the effective sample size for the center, after adjustment for patient characteristics of the center. Center characteristics are expliciitly excluded from the Cox model. The classical variance stabilizing transformation for the Poisson brings us back to the Gaussian model, \[ z_{it} = \sqrt{y_{it}/\mu_{it}} \sim \mathcal N (\theta_i, 1/w_{it}), \] where $\theta_i = \sqrt{\rho_i}$ and $w_{it} = 4 \mu_{it}$. Exchangeability of the centers yields a mixture model in which the parameter $\theta_i$, is effectively assumed to be drawn iidly from a distribution, $G$. The predictions of expected mortality, $\mu_{it}$, are assumed to be sufficiently accurate that we treat $w_{it}$ as known, and independent of $\theta_i \sim G$.
Over short time horizons like 3 years we assume that centers have a fixed draw of $\theta_i$ from $G$, and thus we have sufficient statistics for $\theta_i$, as, \[ T_i = \sum_{t \in \mathcal{T}} w_{it} z_{it}/w_i \sim \mathcal N (\theta_i, 1 /w_{i}), \] where the set $\mathcal{T}$ is the corresponding three year window and $w_i = \sum_t w_{it}$. Given these ingredients it is straightforward to construct a likelihood for the mixing distribution, $G$, and proceed with estimation of it.
Our objective is then to select centers based on the posterior distributions of their $\theta_i$'s. For example, the posterior tail probability of center $i$ is given by, \[ v_\alpha(t_i, w_i) = \mathbb{P}(\theta_i \geq \theta_\alpha| t_i, w_i) = \frac{\int_{\theta_\alpha}^{+\infty} f(t_i | \theta, w_i) dG(\theta) } {\int_{-\infty}^{+\infty} f (t_i | \theta, w_i) dG(\theta) }, \] where $f$ is the density function of $T_i$ conditional on $\theta_i$ and $w_i$. The capacity constraint requires choosing a thresholding value $\lambda_2^*(\alpha)$ such that \[ \alpha = \int \int \11\{v_\alpha(t, w) \geq \lambda_2^*(\alpha)\} \varphi (t| \theta, w) dG(\theta) dH(w), \] which can be approximated by $\frac{1}{n}\sum_i \11\{v_\alpha(t_i, w_i)\geq \lambda_2^*(\alpha) \}$, and inverted to obtain the threshold. Based on the discussion in Section (ref), for the FDR constraint we chose a thresholding value $\lambda_1^*(\alpha, \gamma)$ such that
where $H$ is the marginal distribution of the observed portion of the variance effect. The numerator can be approximated by $\frac{1}{n}\sum_i (1- v_\alpha(t_i, w_i)) \11\{v_\alpha(t_i,w_i) \geq \lambda_1^*(\alpha,\gamma)\}$ and the denominator can be approximated by, $\frac{1}{n} \sum_i \11\{v_\alpha(t_i,w_i) \geq \lambda_1^*(\alpha,\gamma)\}$.
The posterior mean ranking, in contrast, is based on, \[ \delta(t_i,w_i) = \mathbb{E}[\theta_i | t_i,w_i] = \int \theta f(t_i|\theta, w_i) dG(\theta). \] For the capacity constraint we choose a thresholding value $C_2^*(\alpha)$ such that \[ \alpha = \int \int \11\{\delta(t,w) \geq C_2^*(\alpha)\} f(t| \theta, w) dG(\theta ) dH(w). \] For FDR constraint, we pick a thresholding value $C_1^*(\alpha, \gamma)$ such that \[ \gamma = \frac{\mathbb{P}(\delta(t,w) \geq C_1^*(\alpha, \gamma); \theta < \theta_\alpha)}{\mathbb{P}(\delta(t, w) \geq C_1^*(\alpha, \gamma))}. \] The right hand side of the FDR constraint can be approximated by \[\frac{1}{n} \sum_i \11\{\delta(t_i,s_i,w_i)\geq C_1^*(\alpha,\gamma)\} (1-v_\alpha(t_i,w_i))/ \frac{1}{n} \sum_i \11\{\delta(t_i,w_i)\geq C_1^*(\alpha,\gamma)\}, \] while the right hand side of the capacity constraint can be approximated by \[ \frac{1}{n} \sum_i \11\{\delta(t_i,w_i)\geq C_2^*(\alpha)\}, \] so $C_2^*(\alpha)$ is simply the empirical quantile of the $\delta(t_i, w_i)$.
We will compare the foregoing ranking and selection rules with more naive rules based upon the Poisson and Gaussian MLEs, $\sum_{t \in \mathcal{T}} y_{it} / \sum_{t\in \mathcal{T}} \mu_{it}$, and $T_i$, respectively, a variant of the much maligned P-value, as well as a linear shrinkage procedure. For these rules we do not attempt to control for FDR since this is how they are typically implemented in practice.
To help appreciate the difficulty of the selection task, Table (ref) reports estimated FDR rates for several selection rules under a range of capacity constraints $\alpha$ for both right and left tail selection based on the data from 2004 - 2006. Right tail selection corresponds to identifying centers whose mortality rate is higher than expected; left tail selection to centers with mortality lower than expected. To estimate FDR we require an estimate of the distribution of distribution, $G$. For this purpose we use the smoothed version of the Kiefer-Wolfowitz NPMLE introduced in Section (ref). The biweight bandwidth for the smoothing was chosen as the mean absolute deviation from the median of the discrete NPMLE, $\hat G$. The assessment of FDR reported in Table (ref) reflects the considerable uncertainty associated with the selected set of centers deemed by the capacity constraint to be in the upper (or lower) $\alpha$ quantile based upon our estimate of the distribution, $G$, of unobserved quality.
The MLE rule ranks centers based on their Gaussian MLE, $T_i$, while the Poisson-MLE rule ranks on $\sum_t y_{it} /\sum_t \mu_{it}$, which is the MLE of $\rho_i$ from the Poisson model. P-value ranks centers based on the variance stabilizing transformation from the Poisson model under the null hypothesis $\rho_i = 1$ and $\rho_i > 1$ as the alternative hypothesis for right selection and $\rho_i < 1$ for the left selection. All these rules ignore the compound decision perspective of the problem entirely.
Among the compound decision rules, we consider the linear (James-Stein) shrinkage rule, $\hat \mu_\theta + (T_i - \hat \mu_\theta) \hat \sigma_\theta^2 / (\hat \sigma_\theta^2 + 1/w_i)$ which is the posterior mean of $\theta_i$ based on the model $T_i \sim \mathcal{N}(\theta_i, 1/w_i)$ assuming that the latent variable $\theta_i$ follows a Gaussian distribution with mean $\mu_\theta$ and variance $\sigma_\theta^2$. We also consider the \citeasnoun{EfronMorris} estimator which is a slight modification of the James-Stein estimator.
Finally, PM and TP are the posterior mean of $\theta$ and posterior tail probability of $\theta \geq \theta_\alpha$, for right selection, and $\theta \leq \theta_\alpha$ for left selection based on our estimated $\hat G$. For both left and right tail selection, as $\alpha$ increases, the FDR rate decreases, indicating the selection task becomes easier. All rules that account for the compound decision perspective of the problem have slightly lower FDRs than those that consider each center individually.
{\tiny{
}}
The \citeasnoun{Stars} assigns ratings of five stars down to one star to centers in the proportions $\{0.22, 0.30, 0.35, 0.09, 0.04\}$ respectively. We will abbreviate these ratings to the conventional academic scale of A-F. To illustrate the conflict between the selection criteria we plot in Figure (ref) the centers selected for the grade A (Five star, which consists 22% of the centers that suppose to have their true mortality rate being the lowest) category with and without FDR control. Centers are characterized by pairs, $(T_i, w_i)$, consisting of their weighted mean standarized mortality, $T_i$, and their estimate of the precision, $w_i$, of these mortality estimates. In each plot the solid curves represent the decision boundaries of the selection rule under comparison. Centers with low mortality and relatively high precision appear toward the northwest in each figure.
Panel (a) of the figure compares the posterior tail probability selection with the MLE, or fixed effect, selection. The selection boundary for the MLE is the (red) vertical line, since the MLE ignores the precision of the estimates entirely. The selection boundary for the tail probability rule is indicated by the (blue) curve. A few centers with high precision excluded by the MLE rule are selected by the TP rule, and on the contrary a few centers with low precision are selected by the MLE rule, but excluded by the TP criterion. Panel (b) imposes FDR control with $\gamma = 0.20$ on the TP selection with an estimated thresholding value implied by the FDR constraint using the smoothed NPMLE. The MLE selection is the same as in Panel (a) without the FDR control. We see that under TP rule with FDR control, the number of selected centers is reduced considerably. Instead of selecting 711 centers allowed by the capacity constraint, it selects only 230 centers. In comparison, the MLE rule under capacity constraint has an estimated FDR rate at 0.431. Panel (c) compares selected centers by the TP rule with those selected by a James-Stein linear shrinkage rule. Now the TP rule tolerates a few more low precision centers, while it is the James-Stein rule that demands higher precision to be selected. Finally, in Panel (d) we again subject the TP rule to FDR control of 20 percent, while the James-Stein rule continues to adhere only to the capacity constraint. The TP boundary scales back substantially, suggesting that a large proportion of the extra selections made by James-Stein linear shrinkage rules are likely to be false discoveries. In fact, the estimated FDR rate of the James-Stein rule under just capacity constraint is also 0.431, the same as that of the MLE rule.
Given the longtitudinal structure of the Dialysis data, it would be possible to consider the models in Section (ref) that allow for unobserved variance heterogeneity. We refrain from doing so partly due to space considerations and because we are reluctant to assume stationarity of random effects over longer time horizons.
Given the longitudinal nature of the data, it is natural to ask, “How stable are rankings over time, and isn't there some temporal dependence in the observed data that should be accounted for?” Perhaps surprisingly, the year-to-year dependence in the observed mortality is quite weak. In Figure (ref) we plot a histogram of estimated AR(1) coefficients for the 3230 centers; it is roughly centered at zero and slightly skewed to the left. We do not draw the conclusion from this that there is no temporal dependence in the observed $y_{it}$, but only that there is considerable heterogeneity in the nature of this dependence with roughly as many centers exhibiting negative serial dependence as those with positive dependence. Our approach of considering brief, 3-5 year, windows of presumed stability in center performance is consistent with the procedures of the official ranking agency. In each of these windows we can compute a ranking according to one of the criteria introduced above, and it is of interest to see how much stability there is in these rankings.
To address this question we consider rankings based on the posterior tail probability criterion for three year windows. In each of the 5 3-year windows we assign centers letter grades, A-F, with proportions $\{0.22, 0.30, 0.35, 0.09, 0.04 \}$ respectively. Table (ref) reports the estimated transition matrix between these categories, so entry $i,j$ in the matrix represents the estimated probability of a center in state $i$ moving to state $j$ in the next period.
It is obviously difficult to maintain an “A” rating for more than a couple of periods, but centers with poor performance are also likely to move into the middle of the rankings. Although, as we have seen, there is no guarantee that the posterior tail probability criterion yields a nested ranking, nestedness does hold in this particular application. Posterior mean ranking yields similar transition behavior. The high degree of mobility between rating categories reenforces our conclusion that ranking and selection into rating categories is subject to considerable uncertainty.
Robbins's compound decision framework is well suited to ranking and selection problems, and nonparametric maximum likelihood estimation of mixture models offers a powerful tool for implementing empirical Bayes rules for such problems. Posterior tail probability selection rules perform better than posterior mean rules when precision is heterogeneous. Ranking and selection is especially difficult in Gaussian settings where classical linear shrinkage methods are most appropriate. Nonparametric empirical Bayes methods can substantially improve upon selection methods based on linear shrinkage and traditional $p$-values when the latent mixing distribution is not Gaussian both in terms of power and false discovery rate.