Empirical Bayes Estimates for Large-Scale

Prediction Problems
Bradley E FRON

Classical prediction methods, such as Fisher’s linear discriminant function, were designed for small-scale problems in which the number
of predictors N is much smaller than the number of observations n. Modern scientific devices often reverse this situation. A microarray
analysis, for example, might include n = 100 subjects measured on N = 10,000 genes, each of which is a potential predictor. This article
proposes an empirical Bayes approach to large-scale prediction, where the optimum Bayes prediction rule is estimated employing the data
from all of the predictors. Microarray examples are used to illustrate the method. The results demonstrate a close connection with the
shrunken centroids algorithm of Tibshirani et al. (2002), a frequentist regularization approach to large-scale prediction, and also with false
discovery rate theory.
KEY WORDS: Correlated predictors; Effect size estimation; Empirical Bayes; Local false discovery rate; Microarray prediction;
Shrunken centroid.

1. INTRODUCTION with  and Fn−2 the standard normal and tn−2 cumulative dis-
tribution functions (cdf’s), so that under the classical null hy-
An important class of prediction problems begins with the
observation of n independent vectors, pothesis, zi has a standard normal distribution,

(xj , yj ), j = 1, 2, . . . , n. (1.1) H0 : zi ∼ N (0, 1). (1.4)

Here xj is a N-vector of predictors and yj is a real-valued re- Figure 1 shows a histogram of all 6,033 z-values. The theo-
sponse, taken to be dichotomous in most of what follows. For retical N (0, 1) null distribution fits the center of the histogram
example, xj might include age, height, weight, sex, and other reasonably well, which makes sense because presumably most
data for person j, whereas yj indicates whether or not that person of the N genes have nothing to do with prostate cancer. How-
later developed cancer. Given a newly observed N-vector X, we ever, the histogram’s heavy tails suggest some “nonnull” genes
want to predict its corresponding Y value. Our task is to use the that express themselves differently in sick and healthy subjects,
“training data” (1.1) to construct an effective prediction rule. and these are the genes that should be useful for prediction. Just
Classic prediction methods, such as Fisher’s linear discrim- how to fashion a prediction rule from them is the subject of this
inant function, were fashioned for problems where N is much article. (Note that the zi ’s need not necessarily be obtained from
smaller than n, that is, where the number of predictors is less t-tests. Each of the N z-value calculations might involve a sep-
than the number of training cases. Current high-throughput sci- arate linear regression model, incorporating covariates such as
entific technology tends to produce just the opposite situation, age and weight.)
with N  n; modern equipment may permit thousands of mea- Large-scale prediction problems suffer from a surfeit of pos-
surements on a single individual, but recruiting new subjects sible predictors—6,033 of them in this case, most of which are
remains as difficult as ever.
useless. Even the genuinely nonnull cases appear in exagger-
Microarrays offer the prototypical example. Here xj is a vec-
ated form. Selection bias, the fact that we can only identify
tor of genetic expression measurements on subject j, one for
interesting possible predictors at the extremes of the N cases,
each of N genes, where N is typically several thousand. In the
means that an observed value at say, zi = 4 probably corre-
prostate cancer data (Singh et al. 2002) that we use for moti-
vation, N = 6,033 genes are measured in each of n = 102 men, sponds to a true effect considerably nearer the null hypothesis.
n1 = 50 healthy controls, and n2 = 52 prostate cancer patients. In this article empirical Bayes methods are used both to se-
Given a new microarray measuring the same 6,033 genes, we lect useful predictors and to undo selection bias in the evalu-
want to predict whether or not that man has prostate cancer. ation of their predictive power. This approach was suggested
Let ti be the two-sample t-statistic comparing sick versus by the “shrunken centroids” method of Tibshirani et al. (2002),
healthy subjects for gene i, described in Section 2.
 A simple model is introduced in Section 2 that would lead
x̄i2 − x̄i1
ti = c 0 (c0 = n1 n2 /n), (1.2) to an optimum prediction rule if the parameter values were
σ̂i known. Section 3 discusses Bayes estimation of the optimum
where x̄i1 and x̄i2 are the mean expression levels on gene i for rule, using a model of Brown (1971) and Stein (1981) to aid
the healthy and sick subjects and σ̂i is the usual pooled estimate the calculations (and showing a connection with the theory of
of standard deviation. For easier discussion later, we transform local false discovery rates [FDRs]). An empirical Bayes algo-
the ti ’s to a normal scale, rithm for approximating the Bayes solution is developed in Sec-
zi = −1 (Fn−2 (ti )), (1.3) tion 4. Section 5 modifies the empirical Bayes algorithm to

Bradley Efron is Proffesor, Department of Statistics, Stanford University, © 2009 American Statistical Association
Stanford, CA 94305 (E-mail: brad@[Link]). This work was sup- Journal of the American Statistical Association
ported in part by National Institutes of Health grant 8R01 EB002784 and Na- September 2009, Vol. 104, No. 487, Theory and Methods
tional Science Foundation grant DMS0505673. DOI: 10.1198/jasa.2009.tm08523

1015
1016 Journal of the American Statistical Association, September 2009

Figure 1. The 6,033 z-values from the prostate cancer study (Singh et al. 2002). A standard N (0, 1) density fits the histogram center, whereas
the heavy tails indicate the presence of nonnull genes that may be useful for prediction.

allow for correlation among the predictors. A different prob- two classes; nonnull cases, particularly those with large values
lem is considered in Section 6: the estimation of effect sizes for of |δi |, are promising ingredients for effective prediction.
those cases found to be nonnull, where our empirical Bayes ap- Let
proach provides an alternative to the false coverage rate theory
of Benjamini and Yekutieli (2005). Most of the article concerns Wi ≡ (Xi − μi )/σi , i = 1, 2, . . . , N, (2.2)
dichotomous responses yj , but the results are extended to gen-
eral response variables (e.g., survival times) in Section 7. Sec- be the standardized versions of Xi in (2.1). The optimal predic-
tion 8 concludes with remarks that expand on some of the tech- tion rule is based on the weighted sum
nical points and ideas (identified as Remark A, B, etc. through-

N
out the text). S= δi Wi ∼ N (±δ2 /2c0 , δ2 ), (2.3)
A healthy literature on large-scale prediction has grown up i=1
around innovative computer-intensive techniques, such as sup-
 2
port vector machines, lasso and ridge regression regularization δ2 = N 1 δi , with “±” indicating the two classes as in (2.1).
methods, the singular value decomposition, and sparse data rep- We predict
resentation. Chapter 18 of Hastie, Tibshirani, and Friedman
(2008) provides a nice overview of these techniques. A main “healthy” if S < 0,
goal here, besides presenting some new methodology, is to trace (2.4)
the inferential connections between Bayesian theory, regular- “sick” if S > 0.
ization methods like shrunken centroids, FDRs, and large-scale
prediction. Prediction error rates of the first and second kinds, confusing
healthy with sick or vice versa, both equal
2. A SIMPLE MODEL
α = (−δ/2c0 ). (2.5)
Motivation for our empirical Bayes prediction rules comes
from a simple idealized probability model for a vector of pre- Effective prediction requires a large δ vector. In what follows,
dictors X = (X1 , X2 , . . . , XN ). We assume that the individual prediction error is called simply “α”; see Remark B in Sec-
predictors Xi are independently normal, with (location, scale) tion 8.
parameters (μi , σi ), and with possibly different expectations in Rule (2.3), (2.4) is Fisher’s linear discriminant function ap-
the two subject classes, plied to situation (2.1) (Hastie, Tibshirani, and Friedman 2008),
   assuming equal prior probabilities for the two classes. Remark
Xi − μi ind δi “−” healthy class
∼N ± ,1 , (2.1) B in Section 8 discusses the case of unequal probabilities. Sec-
σi 2c0 “+” sick class,
tion 5 considers a more realistic version of (2.1) that allows for
with c0 = (n1 n2 /n)1/2 as in (1.2). (Here the classes have been correlations among the predictors Xi .
labeled “healthy” and “sick” in deference to the prostate exam- In practice, we need to estimate the parameters
ple. Section 7 discusses nondichotomous response variables.)
Null cases have δi = 0, indicating no difference between the (μi , σi , δi ), i = 1, 2, . . . , N, (2.6)
Efron: Empirical Bayes Estimates for Large-Scale Prediction Problems 1017

 Table 1. Shrunken centroids prediction for the prostate data (2.10),

entering into S = δi Wi . This is where the training data
(2.11) using R program pamr, CRAN
x = (xij ), i = 1, 2, . . . , N; j = 1, 2, . . . , n, and
(2.7) Shrinkage # nonzero CV error
y = (y1 , y2 , . . . , yn ), value λ genes rate
with yj equal +1 or −1 depending on the dichotomous clas- 0.00 6,033 0.34
sification of subject j, comes in. Ebay, the algorithm used 0.54 3,763 0.33
for the numerical calculations here, uses standard estimates for 1.08 1,931 0.23
(μi , σi ): 1.62 866 0.12
2.16 377 0.09
 
x̄i1 + x̄i2 SSi1 + SSi2 1/2 2.70 172 0.10
μ̂i = , σ̂i = , (2.8) 3.24 80 0.16
2 n−2
3.78 35 0.30
where x̄i1 and SSi1 are the mean and within-group sum of 4.32 4 0.41
squares for gene i measurements in the healthy subjects, and 4.86 1 0.48
likewise x̄i2 and SSi2 for the sick subjects. 5.29 0 0.52
If σi were known, then NOTE: The shrinkage parameter λ = 2.16 yields the smallest cross-validated error esti-
x̄i2 − x̄i1 mate, α̂CV = 0.09. The prediction statistic Ŝλ involves 377 of the 6,033 genes.
zi = c0 ∼ N (δi , 1) (2.9)
σi
would provide an obvious estimate of δi , say δ̄i = zi . With σi because cross-validation is involved in the choice of “best” λ,
unknown, we convert the t-statistic ti to the normal scale as in the estimated rate 0.09 may be downwardly biased. It would
(1.2), (1.3). Remark F considers this transformation more care- take a second level of cross-validation to correct this bias.
fully, but for now we will ignore it and use the approximation A small simulation study was run with N = 1,000, n1 =
zi ∼ N (δi , 1) for our actual z-values (1.3). ind
n2 = 10, and all xij ∼ N (0, 1). In this case δi = 0 for every i
Selection bias makes the δ̄i values overinflated estimates of in (2.1), so α = 0.50 at (2.5); but the minimum cross-validated
the true δi ’s. Suppose that for the prostate data, we decided to error rates observed in 100 repetitions of this set-up had median
use the genes with the 51 largest values of |δ̄i | for prediction. 0.30 with standard deviation ±0.16.
The vector of 51 δ̄i ’s has δ̄ = 27.3, suggesting α = 0.003 in This is an extreme example. Usually the downward bias is
(2.5) [using c0 = (50 · 52/102)1/2 = 5.05]. The empirical Bayes less severe, particularly when good prediction is possible. Nev-
calculations of Section 3 show that a more realistic estimate for ertheless, in what follows we will try to avoid such biases by
the actual 51-vector’s length is 19.8, giving α = 0.025. using rules in which the cross-validation calculations are not
The shrunken centroids algorithm of Tibshirani et al. (2002)
involved in the choice of tuning parameters.
counteracts selection bias by shrinking the estimates δ̄i = zi to-
ward zero according to a soft thresholding rule, 3. BAYESIAN PREDICTION
δ̂i = sign(zi ) · (|zi | − λ)+ . (2.10) Suppose that we had a Bayesian prior distribution for the pa-
In words, each value δ̄i = zi is shrunk toward zero by the rameters in model (2.1) that allowed us to calculate posterior
amount λ, under the restriction that shrinking never goes past expectations for the δi ’s, say
zero. A range of possible shrinkage parameters λ is tried, and
for each one a prediction rule like (2.3) is formed, using δ̃i = E{δi |z}. (3.1)

Ŝλ = δ̂i Ŵi [Ŵi = (Xi − μ̂i )/σ̂i ] (2.11) Bayes estimates are immune to selection bias; even if zi were
selected because it was the largest of the N z-values (zi = 5.29
for prediction as in (2.4). Cross-validation is then used to esti- for gene i = 610 in the prostate data), δ̃i would still be the cor-
mate αλ , the true error rate. (This description takes some liber- rect Bayes estimate for δi . We could,
ties with the details of the shrunken centroids procedure.)  for example, use the 50
largest values of |δ̃i | to form S̃ = δ̃i Ŵi , as in (2.3) or (2.11),
Note that only cases with |zi | > λ enter into the prediction while maintaining at least some confidence in the error rate es-
statistic Ŝλ . This is a favorable property; prediction is easier to timate α̃ = (−δ̃/2c0 ). (See Senn 2008 and Dawid 1994 for
implement and understand when the number of predictors is discussions of the “paradox” of Bayesian immunity to selection
small. effects, including its dangers.)
Table 1 shows a shrunken centroids analysis for the prostate
Brown (1971) and Stein (1981) developed a Bayesian model
data, carried out using pamr, a CRAN algorithm in the R lan-
that is especially convenient for calculating δ̃i in (3.1). For any
guage. Cross-validation suggests λ = 2.16 as the best shrink-
(δ, z) pair, we assume that δ has prior density g(δ),
age parameter [so, e.g., zi = 4 yields δ̂i = 1.84 in (2.11)] with
estimated error rate α̂CV = 0.09. 377 of the 6,033 genes are δ ∼ g(·) and z|δ ∼ N (δ, 1), (3.2)
involved in Ŝλ . Unlike the theoretical result (2.5), adding too
many predictors eventually decreases prediction accuracy in ac- so that z has marginal density
tual practice.  ∞ √
Looking at Table 1, it seems we should use λ = 2.16 in our f (z) = ϕ(z − δ)g(δ) dδ ϕ(z) = e−z
2 /2
/ 2π . (3.3)
prediction rule. There is a subtle danger lurking here, however; −∞
1018 Journal of the American Statistical Association, September 2009

Theorem 1. Under model (3.2), the posterior density of δ Suppose that we add to the Brown–Stein model (3.2) the as-
given z is sumption that the prior distribution of δ has a discrete atom of
probability at δ = 0 (see Remark C, Section 8),
g(δ|z) = eδz−ψ(z) e−δ
2 /2
g(δ) ,
p0 = Prob{δ = 0}. (3.7)
with ψ(z) = log(f (z)/ϕ(z)). (3.4)
Then Bayes theorem yields
Proof. According to the Bayes theorem,
fdr(z) ≡ Prob{δ = 0|z} = p0 ϕ(z)/f (z), (3.8)
g(δ|z) = ϕ(z − δ)g(δ)/f (z), (3.5)
where fdr(z) is the “local false discovery rate” (Efron 2008).
which reduces immediately to (3.4). Comparing this with (3.4), (3.6) gives
(3.4) represents an exponential family with sufficient statis- Corollary 2. Under model (3.2), (3.7),
tic δ, natural (canonical) parameter z, and cumulant generating
d
function (cgf) ψ(z). Therefore, the conditional cumulants of δ E{δ|z} = − log(fdr(z)) and
given z can be obtained by differentiating ψ with respect to z. dz
(3.9)
d2
Corollary 1. var{δ|z} = − 2 log(fdr(z)).
dz
E{δ|z} = ψ  (z) and var{δ|z} = ψ  (z). (3.6)
[Brown and Stein used multivariate versions of (3.6), derived It seemingly makes sense that only genes with low false dis-
differently, in their exploration of high-dimensional estimation covery rates should be utilized in prediction rules. The corollary
theory.] shows that this is roughly true, but in a rather surprising man-
The advantage of Corollary 1 is that ψ(z), and the cumulants ner; large values of δ̃i = E{δi |zi } depend on the rate of change of
of δ given z, are obtained directly from the marginal density log(fdr(zi )), not on fdr(zi ) itself. Small values of fdr(zi ) usually
f (z) without requiring specific calculation of the prior g(δ), fi- correspond to large values of |δ̃i |, but this does not have to be
nessing the usual difficulties of deconvolution. the case. Usually log(fdr(zi )) is nearly constant around z = 0,
The algorithm Ebay described in Section 4 approximates where fdr(z) = 1. This forces both Ê and SD to be small, as
E{δ|z} and var{δ|z} by substituting a smoothed estimate ψ̂(z) seen in Figure 2 (see Remark I).
into (3.6). Figure 2 displays the Ebay output Ê{δ|z} for the
4. EMPIRICAL BAYES PREDICTION
prostate data and compares it with the shrunken centroids curve
(2.10) for λ = 2.16, the preferred choice in Table 1. Ê is better The Ebay algorithm that produced Figure 2 uses empirical
matched to the choice λ = 1.42 in (2.10), suggesting that less Bayes methods to construct effective prediction rules; that is,
shrinkage is better here. it uses z, the vector of all N z-values, to estimate the Bayes

Figure 2. The heavy curve is Ê{δ|z} for the prostate data, Ebay algorithm Section 4, compared with the best shrunken centroids curve (2.10),
λ = 2.16. Also shown is SD = var{δ|z}1/2 . At z = 4, Ê = 2.49, shrunken centroid = 1.84, SD = 0.98. Remark I explains the slight positive slope
of Ê{δ|z} for z in (−2, 1.5).
Efron: Empirical Bayes Estimates for Large-Scale Prediction Problems 1019

prediction rule (2.3), (2.4). A schematic description of Ebay’s There are, potentially, many reasons why the nominal error
operation goes as follows: rate 0.025 might be overoptimistic: (μ̂i , σ̂i ) in (2.8) does not
equal (μi , σi ), the Xi ’s are not normally distributed, the Xi ’s are
(1) A target error rate α0 is selected (default α0 = 0.025).
not independent (see Section 5), and the empirical Bayes esti-
(2) An estimate f̂ (z) for the marginal density f (z), (3.3), is mates δ̂i differ from the actual Bayes estimates (3.1). This last
obtained using Poisson regression on z; see Remark D. point can cause particular trouble at the extremes of the z scale,
(3) The estimated cumulative generating function ψ̂(z) = just where |δ̂(z)| is largest but there are the fewest zi ’s for es-
log(f̂ (z)/ϕ(z)), (3.4), is numerically differentiated to timating δ̂. Figure 3 concerns the following artificial situation,
give using notation similar to that for the prostate data and model
δ̂i = ψ̂  (zi ) = Ê{δi |zi }, (4.1) (2.1):

as in (3.6). • N = 5,000, n1 = n2 = 20,

(4) Letting δ̂ I be the vector of I largest δ̂i ’s (in absolute ind
• δi ∼ N (1.5, 1) for i = 1, 2, . . . , 250,
value), I is chosen to be the smallest integer such that (4.4)
the nominal error rate (−δ̂ I /2c0 ), (2.5), is less than • δi = 0 for i = 251, 252, . . . , 5,000,
α0 ; that is, I is the minimum choice yielding ind
• xij ∼ N (±δi /2c0 , 1) for all i and j, c0 = 202 /40.
−1
δ̂ I  ≥ 2c0  (1 − α0 ). (4.2)
This results in
(5) The empirical Bayes prediction rule is based on the sign ind
of zi ∼ N (δi , 1) (4.5)
  Xi − μ̂i  at (2.9), with δi ∼ N (1.5, 1) for the first 250 genes and 0 other-
Ŝ = δ̂i , (4.3) wise.
σ̂i
I Figure 3 compares Ê{δ|z} from the Ebay algorithm with the
(μ̂i , σ̂i ) as in (2.8). true curve E{δ|z}. The estimates are reasonably accurate up to
(6) Repeated 10-fold cross-validation is used to provide an z = 4 but degenerate beyond that. Remark E of Section 8 de-
unbiased estimate of the rule’s prediction error; see Re- rives a delta-method formula for the standard error of Ê{δ|z}
mark G. that predicts this behavior. Modifying (4.4), (4.5) so that the
zi ’s were correlated, with a root mean squared correlation co-
Table 2 shows a portion of Ebay’s output for the prostate efficient of 0.1, increased the variability of Ê{δ|z} by roughly
data. Its prediction rule utilizes genes with the 51 largest val- 50%.
ues of |δ̂i |, at which point (4.2) is first satisfied (compared with An option in Ebay allows for truncation of the δ̂ estimation
377 genes for the apparently best shrunken centroids rule in Ta- procedure at some number “ktrunc ” of observations in from the
ble 1). An unbiased error estimate, based on 20 randomized 10- extremes. With ktrunc = 5, for instance, δ̂i for the five largest zi
fold cross-validation runs, was 0.092, the same as the minimum values is set equal to max{δ̂i : i ≤ N − 5}, and similarly at the
error seen in Table 1; see Remark G. negative end of the z scale.
Figure 4 shows the actual misclassification error probabilities
Table 2. Ebay prediction rule for the prostate data α for 200 simulations from model (4.4), each time using the
Ebay prediction rule with nominal error rate α0 = 0.025. As
Step Index z-value δ̂ α̂ α̂cor the truncation parameter ktrunc increases from 0 to 15, the ac-
1 610 5.29 4.30 0.335 0.335 tual prediction errors α decrease toward the target value 0.025.
2 1,720 4.83 3.78 0.285 0.281 Table 3 displays the means and standard deviations for the data
3 3,64 −4.42 −3.70 0.250 0.250 given in Figure 4.
4 3,940 −4.33 −3.64 0.222 0.222 Truncation had a less dramatic effect on the prostate data;
5 4,546 −4.29 −3.58 0.199 0.215 for ktrunc = 0, 5, 10, 15, the cross-validated error estimates were
6 4,331 −4.14 −3.40 0.182 0.189 0.092, 0.085, 0.070, and 0.077. Lowering the target rate from
7 332 4.47 3.34 0.167 0.181 α0 = 0.025 to 0.01 gave corresponding error estimates of 0.070,
8 914 4.40 3.24 0.154 0.166 0.062, 0.061, 0.058. Correlation among the predictors is part of
9 1,068 4.25 3.06 0.144 0.148 the problem here; see Section 5.
10 4,088 −3.88 −3.05 0.135 0.149 Our original error estimate α̂ = 0.092 is “honest,” that is,
.. .. .. .. .. ..
. . . . . . nearly unbiased for the Ebay rule produced with (α0 , ktrunc ) =
45 4,154 −3.38 −2.26 0.029 0.050 (0.025, 0). So are the α̂ estimates for the other (α0 , ktrunc ) com-
46 2 3.57 2.25 0.028 0.050 binations. But choosing the combination with the smallest α̂
47 2,370 3.56 2.24 0.028 0.049 again raises the possibility of overoptimism, as discussed at the
48 3,282 3.56 2.23 0.027 0.048 end of Section 2.
49 3,505 −3.33 −2.18 0.026 0.046 More elaborate “honest” selection criteria, beyond the cur-
50 905 3.51 2.18 0.025 0.047 rent capabilities of Ebay, might involve minimizing a linear
51 4,040 −3.33 −2.17 0.025 0.048 combination of nominal error rate and number of predictors,
NOTE: The rule uses the genes with the 51 largest |δ̂i | values, α̂ = (−δ̂/2c0 ) = say
0.025. The cross-validation error rate is 0.092 ± 0.004. The column α̂cor is explained in
Section 5. (−δ̂ I /2c0 ) + C · I (4.6)
1020 Journal of the American Statistical Association, September 2009

Figure 3. True curve E{δ|z} (heavy), compared with Ê{δ|z} = ψ̂(z) , 50 simulations from model (4.4). Estimates Ê are reasonably accurate
for z < 4 but fall apart for larger z-values.

over all choices of I; accounting for correlation as in Sec- Some “snooping” into the cross-validation estimates seems
tion 5, adjusting for nonnormality; using theoretical or data- inevitable in applications. Nevertheless, I believe that keeping
based techniques to choose the truncation parameter, and so on. snooping to a minimum is good practice for honest prediction

Figure 4. Actual prediction errors α of Ebay rule with nominal α0 = 0.025; 200 simulations from model (4.4). As truncation parameter
increases from 0 (rightmost histogram) to 15 (leftmost), the actual errors decrease toward nominal α0 .
Efron: Empirical Bayes Estimates for Large-Scale Prediction Problems 1021

Table 3. Means and standard deviations for actual prediction errors α Table 4. Ebay output for the Michigan lung cancer study
in simulation experiment for Figure 4
Step Index z-value δ̂ α̂ α̂cor
ktrunc
1 3,144 4.62 3.683 0.3290 0.329
15 10 5 0 2 2,446 4.17 3.104 0.2813 0.307
Mean 0.032 0.038 0.051 0.066 3 4,873 4.17 3.104 0.2455 0.256
SD 0.019 0.018 0.022 0.025 4 1,234 3.90 2.686 0.2234 0.225
5 621 3.77 2.458 0.2072 0.213
6 676 3.70 2.323 0.1942 0.228
7 2,155 3.69 2.313 0.1824 0.230
assessment, and that empirical Bayes methods, perhaps further 8 3,103 3.60 2.140 0.1731 0.236
refined, can be sufficiently accurate to allow for a nearly honest 9 1,715 3.58 2.103 0.1647 0.240
practical methodology. 10 452 3.54 2.028 0.1574 0.243
.. .. .. .. .. ..
5. CORRELATION CORRECTIONS . . . . . .
193 3,055 2.47 0.499 0.0519 0.359
The assumption of casewise independence in model (2.1) is 194 1,655 −2.21 −0.497 0.0518 0.359
likely to be untrue (perhaps spectacularly so) in many appli- 195 2,455 2.47 0.496 0.0517 0.359
cations. Suppose that the vector W of standardized predictors 196 3,916 2.47 0.496 0.0516 0.359
Wi = (Xi − μi )/σi , (2.2), actually has covariance matrix . 197 4,764 2.47 0.495 0.0515 0.359
Then both error probabilities in (2.5) become 198 1,022 −2.20 −0.492 0.0514 0.359
 199 1,787 −2.19 −0.490 0.0513 0.360
0 = δ/2c0 , 200 901 −2.18 −0.486 0.0512 0.360
α = (− 0 · η), where (5.1)
η = (δ t δ/δ t δ)1/2 . NOTE: The correlation error estimates α̂cor are much more pessimistic, as confirmed by
cross-validation.
Here 0 is the independence value, while η is a correction fac-
tor, usually <1, that increases the error rate α; see Remark K in
Section 8. prediction rule was made before the cross-validation calcula-
If we can estimate , then we can estimate correction fac- tions.
tor η, Sample correlation matrixes tend toward overdispersion
when n is small compared with the number of variates. Ebay
ˆ δ̂)1/2 .
η̂ = (δ̂ t δ̂/δ̂ t  (5.2) includes an option for empirical Bayes shrinkage of the ele-
ˆ see Remark H.
ments of ;
According to (2.1), cov(W) =  has diagonal elements 1 in
both classes, so the off-diagonal elements ρii are correlations. 6. EFFECT SIZE ESTIMATION
Notice that these need to be estimated only for the I cases
Current developments in large-scale simultaneous inference
selected by the Ebay algorithm at (4.2), not for all N cases.
have focused on hypothesis testing, where the goal is to iden-
For the prostate data, we need to estimate a 51 × 51 correla-
tify a small number of nonnull cases among a large number
tion matrix , from the 51 × 102 data submatrix xI of the full
of potential candidates (see Dudoit, Shaffer, and Boldrick 2003
6,033 × 102 matrix x whose rows are indexed by the first col- for a nice review). Benjamini and Yekutieli (2005) addressed a
umn of Table 2. more ambitious goal: to assess the effect sizes for the nonnull
The last column of Table 2 in Section 4 shows α̂cor , obtained cases, that is, to estimate how far away they lie from the null
from (5.1), (5.2), with ˆ the usual sample correlation matrix.
hypothesis. The empirical Bayes theory of Section 4 provides
Correlation degrades the nominal error probability from 0.025 an alternative approach to effect size estimation.
to 0.048 (closer to the cross-validation estimate 0.092). Much We begin with assumptions (3.2), (3.7), that
of the degradation is due to three large correlations,
zi ∼ N (δi , 1), i = 1, 2, . . . , N, (6.1)
r34,19 = 0.97, r36,15 = 0.65, r42,28 = 0.92, (5.3)
and that proportion p0 of the effects δi equal 0,
the subscripts referring to the steps in Table 2. p0 = Prob{δi = 0}, (6.2)
Table 4 concerns a microarray study having more severe cor-
relation problems, the Michigan lung cancer study discussed these being the uninteresting null cases. The local false discov-
by Subramanian et al. (2005). The data include N = 5,217 ery rate fdr(z) = p0 ϕ(z)/f (z) (3.8) is the Bayes posterior prob-
genes, n = 86 subjects, n1 = 62 “good outcomes,” and n2 = 24 ability Pr{δi = 0|zi }. If fdr(zi ), an estimate of fdr(zi ), is suitably
“poor outcomes.” Here the Ebay algorithm stopped after 200 small, then case i can be reported as “probably nonnull,” and
steps, without α̂ reaching the target value of α0 = 0.025. The we would like to put some sort of confidence limits on the ef-
correlation-corrected errors α̂cor were much more pessimistic, fect size δi . The prior g(δ) in (3.2) is now of the mixed form
actually increasing after the first six steps, eventually to α̂cor = g(δ) = p0 I0 (δ) + (1 − p0 )g1 (δ), (6.3)
0.360. A cross-validation error rate of 0.37 confirms the pes-
simism. Restricting Ebay to use at most I = 10 predictors re- where I0 (δ) is a delta function at 0 and g1 (δ) indicates the den-
duced the cross-validated error rate to 0.29, as suggested by Ta- sity of the nonnull cases (see Remark C). Then the mixture den-
sity f (z) (3.3) becomes
ble 4 (an example of the kind of “snooping” disparaged at the
end of Section 4, unless the decision to use the I = 10 Ebay f (z) = p0 ϕ(z) + (1 − p0 )f1 (z),
1022 Journal of the American Statistical Association, September 2009

where Corollary 3. Under model (6.1), (6.2),

 ∞
f1 (z) = ϕ(z − δ)g1 (δ) dδ. (6.4) E1 {δ|z} = E{δ|z}/[1 − fdr(z)]
−∞
and
Theorem 2. Under model (6.1), (6.2), the posterior density  
1 fdr(z)
of effect size δ given z, and given that δ = 0, is var1 {δ|z} = var(δ|z) − 2
E{δ|z} .
1 − fdr(z) 1 − fdr(z)
g1 (δ|z) = eδz−ψ1 (z) e−δ
2 /2
g1 (δ) , (6.10)

where Note. Because δ = 0 with probability fdr(z), we have

 
1 − fdr(z) 1 − p0 E{δ j |z} = [1 − fdr(z)] · E1 {δ j |z}. (6.11)
ψ1 (z) = log . (6.5)
fdr(z) p0
Using (6.11) with j = 1 and 2 leads to a quick verification of
Proof. Bayes rule says that g1 (δ|z) = ϕ(z − δ)g1 (δ)/f1 (z), (6.10).
yielding
Our prediction algorithm in Section 4 requires only the esti-
−δ 2 /2 mation of E{δ|z}. Effect size estimation is more difficult, requir-
g1 (δ|z) = e δz−log{f1 (z)/ϕ(z)}
e g1 (δ) . (6.6)
ing var{δ|z} and fdr(z) as well. The plug-in estimate of var1 {δ|z}
An equivalent form of (3.8) is in (6.10) may be particularly unstable, in which case we can
conservatively replace it with var{δ|z}, as shown next.
1 − fdr(z) = Pr{δ = 0|z} = (1 − p0 )f1 (z)/f (z), (6.7) Rearranging (6.10) yields
from which we obtain, using (6.4), var1 1 − fdr(z)Q(z)
= ,
var 1 − fdr(z)
f1 (z) p0 1 − fdr(z)
= . (6.8) E2
ϕ(z) 1 − p0 fdr(z) where Q(z) = (6.12)
(1 − fdr(z)) · var
Combining (6.8) and (6.6) verifies Theorem 2.
(with var = var{δ|z}, etc.), so that
As in (3.6), the conditional moments of a nonnull δ (one for
which δ = 0) given z are obtained by differentiating ψ1 (z), var1 ≤ var if var ≤ E2 /(1 − fdr(z)). (6.13)

E1 {δ|z} = ψ1 (z) and var1 {δ|z} = ψ1 (z), (6.9) Because var is usually near 1, this last condition is satisfied
whenever δ̂ = E{δ|z} is sufficiently large to be interesting—in
where the subscript “1” indicates conditioning on δ = 0. Some the case of the prostate data, for z ≥ 2.
calculation gives E1 and var1 in terms of E{δ|z} and var{δ|z} Figure 5 demonstrates effect size estimation for the prostate
in (3.6), data. The Poisson generalized linear model (GLM) estimate of

Figure 5. Effect size estimation for the prostate data. The band is an approximate 68% interval for δ given z and given δ = 0, also showing
fdr(z), the estimated probability δ = 0 given z. At z = 4, fdr(z) = 0.048, with interval [1.58, 3.64] if δ is nonnull.
Efron: Empirical Bayes Estimates for Large-Scale Prediction Problems 1023

Figure 6. Approximate effect size limits (6.14) for 25 replications of simulation model (4.4). The heavy straight lines are actual 68% Bayes
posterior limits for nonnull cases.

f (z) described in Remark D provides estimates of fdr(z), Ê{δ|z}, tion extends the empirical Bayes prediction methodology to
and var{δ|z}, as in Figure 2. The curved band in Figure 5 follows general univariate responses.
Let Y be a univariate response of interest, for example, a sur-
Ê{δ|z}/[1 − fdr(z)] ± var{δ|z}1/2 (6.14)
vival time that we wish to predict from X = (X1 , X2 , . . . , XN )
as in Corollary 3, showing approximate 68% intervals for δ as in Section 2. For convenience, we assume that Y has been
given z and given δ = 0, made more conservative by replac- standardized to have mean 0 and variance 1, denoted by
ing var1 with var. At z = 4, for example, we estimate that ei-
Y ∼ (0, 1), (7.1)
ther δ = 0 with probability fdr(4) = 0.048 or, if δ = 0, it lies
in the interval [1.58, 3.64) with estimated posterior probability although this plays no role in the actual methodology.
exceeding 0.68. [Remember that δ, as defined in (2.1), is the We suppose that Y influences the standardized variable Wi =
number of standard deviations separating the two class means, (Xi − μi )/σi (2.1) through linear regression,
multiplied by c0 .]
Wi = βi Y + i , i = 1, 2, . . . , N, (7.2)
Benjamini and Yekutieli’s (2005) false coverage rate algo-
rithm provides conservative frequentist confidence bounds on var(i ) = 1, where the vector of errors  is uncorrelated with
the cases declared nonnull by an FDR testing procedure, as- Y. In the dichotomous situation of (2.1), Y = −1 or 1 and βi =
suming independence of the zi ’s. There is a heavy price to pay, δi /2c0 . Effective prediction of Y depends on discovering those
however; the bounds tend to be very wide. For z = 4 in the Xi ’s with large values of |βi |; see Remark J in Section 8.
prostate example, their 68% interval is [1.36, 6.64) (using their The joint distribution of Y and W has mean vector and co-
definition 1, with q = 0.32). Part of the problem, as discussed variance matrix
in section 7 of Efron (2008), is that the Bejamini–Yekutieli pro-      
Y O 1 βt
cedure does not split off an atom of probability at δ = 0, even ∼ , , (7.3)
W O β ββ t + 
though splitting seems natural in the hypothesis testing frame-
work of (3.7) or (6.3).  indicating the covariance matrix of . The best linear predic-
The approximate 68% nonnull limits (6.14) were calcu- tor of Y from W is
lated for 25 replications of simulation model (4.4). They ap- Y † = β t (ββ t + )−1 W
pear in Figure 6,√along with the true Bayesian posterior limits
(z + 1.5)/2 ± 1/ 2. Using var instead of var1 in (6.14) makes 1
= β t  −1 W, (7.4)
the intervals too wide, but their overall performance is accept- 1+ 2

able as rough estimates of effect size. where 2 is the squared Mahalanobis distance
7. OTHER RESPONSE VARIABLES 2
= β t  −1 β. (7.5)
The development so far has concerned dichotomous response If  is the identity, as assumed in (2.1), then Y† = constant ·
variables, healthy versus sick in the prostate example. This sec- β t W, similar to (2.3).
1024 Journal of the American Statistical Association, September 2009

Combining (7.4) with (7.2) produces a simple expression for 4. We continue increasing I until either corI reaches some
the conditional mean and variance of Y † given Y, target value or I reaches a preselected upper bound, and
 2 2  use
Y † |Y ∼ Y, , (7.6)  
1  δ̂i Xi − μ̂i
I
1+ 2 1+ 2
Ŷ =
†
, (7.14)
from which (7.1) gives 1 + ˆ 2I i=1 2c0 σ̂i

cor(Y, Y † ) = / 1+ 2. (7.7) from (7.4), to predict Y from X.

Effective prediction of Y from W requires a large value of = Steps 3 and 4 assume uncorrelated Xi ’s (i.e.,  the identity
(β t β)1/2 . [In the context of Section 2, where  = I and β = matrix), but correlation can be incorporated as in Section 5.
δ/2c0 , we have = δ/2c0 , so the error probability α equals These steps were carried out for an ongoing lung cancer mi-
(− ) at (2.5).] croarray study involving n = 100 patients each measured on
To bring empirical Bayes methods to bear on the estimation N = 16,000 genes. All patients received the same new drug.
of Y † , we need to estimate posterior expectations for the regres- The response variable “Y” was a categorical assessment of
sion coefficients βi from the training data (1.1): the N ×n matrix improvement, adjusted for two covariates, running from −2
x and the n-vector of responses y = (y1 , y2 , . . . , yn )t . Let xti in- (worst) to +2 (best).
dicate the ith row of x. Applying model (7.2) independently to Figure 7 shows Ê{δ|z}, calculated by steps 1 and 2. It seems
each column of x gives a linear model for the rows, clear that any power of the microarray expression measure-
ments to predict Y must come from those genes having zi less
xi = μi 1n + σi (βi y +  i ), (7.8) than −2. Table 5 shows this to be true. Predictive power is mod-
est here, with a theoretical correlation of only 0.48 after I = 50
where 1n is a vector of n 1’s and the components of  i =
(i1 , i2 , . . . , in )t are independent and identically distributed, steps, asymptoting to 0.57 at I = 16,000.
with mean 0 and variance 1. Ordinary least squares applied to 8. REMARKS
(7.8) provides familiar estimates of μi , σi , and βi . In the di-
chotomous setting of (2.1), μ̂i and σ̂i are as given in (2.8), while The following remarks expand on some of the questions and
2c0 β̂i equals δ̄i = zi in (2.9). technical points raised earlier.
If we assume that the errors i are normally distributed, then
Remark A (Centroids interpretation).
 Prediction rule (4.3),
the t-statistic “ti ” for testing βi = 0 in (7.8) has a noncentral
which depends on the sign of Ŝ = δ̂i Ŵi , Ŵi = (Xi − μ̂i )/σ̂i ,
t distribution with n − 2 degrees of freedom and noncentrality
can be stated in more conventional centroid terminology. Let-
parameter proportional to βi ,
ting
 

n
ti ∼ tn−2 (δi ) δi ≡ 2c0 βi with c0 =
2
(yi − ȳ) /4 . (7.9)
2 D1 = Ŵ + δ̂/2c0  and D2 = Ŵ − δ̂/2c0 , (8.1)
1
we predict “healthy” if D1 < D2 and “sick” if D2 < D1 ; so
In usual practice, (7.9) remains a reasonable approximation as δ̂/2c0 and −δ̂/2c0 are the standardized centroids. An alterna-
long as the ij distribution does not have heavy tails. With di- tive statement refers to the hyperplane L̂ passing through the
chotomous yi , c0 = (n1 n2 /n)1/2 as before. origin of N-space orthogonal to the line segment connecting
We can transform ti to a z-value via δ̂/2c0 with −δ̂/2c0 : we predict healthy or sick is predicted de-
zi = −1 (Fn−2 (ti )), (7.10) pending on which side of L̂ the point Ŵ falls.

with  and Fn−2 as the standard normal and central tn−2 cdf’s. Remark B (Unequal prior probabilities). Prediction rule (4.3)
If n is large, then (7.9) gives tacitly assumes that our dichotomous response variable has
equal prior probabilities on the two categories irrespective of
zi ∼
˙ N (δi , 1) (7.11) the observed frequencies n1 and n2 in the training set. Suppose
as in (2.9). Remark F improves on approximation (7.11), but we that the prior probabilities are actually π1 and π2 . Starting with
take it as given here. model (2.1), calculations involving Fisher’s linear discriminant
We can now proceed as in Section 4: function imply the following change from Remark A: The pre-
diction boundary L̂ is translated to intersect the orthogonal line
1. z = (z1 , z2 , . . . , zN )t provides f̂ (z), an estimate of the mar- segment at a directed distance,
ginal density of the z-values (Remark D) and ψ̂(z) =  
f̂ (z)/ϕ(z). c0 π1
log , (8.2)
2. We then calculate δ̂ π2

δ̂i = ψ̂  (zi ) = Ê{δi |zi } for i = 1, 2, . . . , N. (7.12) from the origin. [The definition of Ŵ is still (X − μ̂)/σ̂ , with
(μ̂i , σ̂i ) as given in (2.8).]
3. δ̂ I , the vector of I largest δ̂i ’s in absolute value, gives
Remark C [The prior density g(δ)]. In the Brown–Stein
ˆ I = δ̂ I /2c0 and corI = ˆ I / 1 + ˆ 2I . model (3.2), the prior density g(δ) can be extended to a gen-
(7.13) eral probability distribution G(δ) incorporating discrete atoms
Efron: Empirical Bayes Estimates for Large-Scale Prediction Problems 1025

Figure 7. Lung cancer microarray study, N = 16,000 genes, n = 100 patients, ordered categorical response variable Y. The heavy curve
represents Ê{δ|z}, (7.14). Dashes indicate those z-values exceeding 3 in absolute value.

of probability as in (3.7). Theorem 1’s statement is almost un- are K = 90 bins, each of width 0.1, ranging from −4.5 to 4.5.
changed: The counts

dG(δ|z) = eδz−ψ(z) e−δ

2 /2
dG(δ). (8.3) ck = #{zi in bin k}, k = 1, 2, . . . , K, (8.4)
are the heights of the histogram bars. Let b indicate the K-
The factor e−δ /2 guarantees that the exponential family has
2
vector of bin midpoints. Then the estimate f̂ = (f̂1 , f̂2 , . . . , f̂K )t
natural parameter space including all values of z, justifying of f (z) at the points in b is obtained by Poisson regression of
Corollary 2 for all z. The same considerations apply to Theo- the counts on a natural spline function of the midpoints,
rem 2 as well.
f̂ = glm(c ∼ ns(b, df), poisson)$fit (8.5)
Remark D [Estimating f (z)]. Ebay estimates f (z), the mix- in R notation; the default degrees of freedom (df) equals 7 in
ture density (3.3), by means of a Poisson GLM applied to Ebay; and f̂ is the discretized maximum likelihood estimate of
binned counts of the N z-values. In Figure 1, for example, there f (z) in the seven-parameter exponential family defined by the
natural spline basis.
Table 5. Predictive analysis of lung cancer data Estimate (8.5) is the same one employed by locfdr, the lo-
cal FDR algorithm described by Efron (2008). Applied to the
Step Index z-value δ̂ β̂ ˆI corI prostate data, locfdr estimated p̂0 = 0.93 for the proportion of
1 12,404 −4.27 −2.44 −0.20 0.04 0.04 null genes (3.7), assuming that f (z) is the correct null density.
2 6,342 −3.98 −2.10 −0.17 0.07 0.07 Remark E (Accuracy formula for Ê{δ|z}). A closed-form
3 2,516 −3.92 −2.02 −0.16 0.10 0.10
delta method expression for the variance of δ̂i = Ê{δi |zi } can be
4 488 −3.89 −1.99 −0.16 0.12 0.12
derived if we are willing to assume that the zi ’s are independent
5 8,471 −3.84 −1.93 −0.16 0.15 0.15
6 25 −3.84 −1.92 −0.16 0.17 0.17 of one another. Let M be the K × m structure matrix ns(b, df)
7 2,872 −3.82 −1.90 −0.15 0.20 0.19 in (8.5), K = 90 and m = 8, diag(c) the K × K diagonal matrix
8 300 −3.78 −1.85 −0.15 0.22 0.21 with diagonal entries the bin counts ck , and G = Mt diag(c)M.
9 545 −3.78 −1.85 −0.15 0.24 0.23 Section 5 of Efron (2007) uses the relationship
−3.78 −1.84 −0.15
10 12,448 0.26 0.26
dˆ = MG−1 Mt dc (8.6)
.. .. .. .. .. .. ..
. . . . . . . for the derivative matrix of the K-vector ˆ = log(f̂) with respect
45 10,905 −2.83 −0.60 −0.05 0.54 0.47
to a continuized version of c.
46 390 −2.83 −0.59 −0.05 0.54 0.47
47 1,498 −2.82 −0.59 −0.05 0.54 0.48
Let D be the (K − 2) × K matrix whose kth row is
48 10,317 −2.81 −0.57 −0.05 0.54 0.48 (0, 0, . . . , 0, −1, 0, 1, 0, 0, . . .)/d0 , (8.7)
49 7,894 −2.79 −0.56 −0.05 0.55 0.48
50 13,263 −2.79 −0.55 −0.04 0.55 0.48 with −1 in the kth place, so Dˆ = ˆ  , the numerical derivative
ˆ This gives
of .
NOTE: The right column shows corI , (7.13), for the lung cancer data; I = 1–50. Final
value, cor16,000 = 0.57. dˆ  = DMG−1 Mt dc. (8.8)
1026 Journal of the American Statistical Association, September 2009

Table 6. Delta method standard errors for δ̂(z) = Ê{δ|z}, Suppose that the Brown–Stein model (3.2) is modified to
formula (8.10), for the prostate data have z|δ ∼ N (δ, σ 2 ). Then it is easy to show that
z −4 −3 −2 −1 0 1 2 3 4 E{δ|z} = z + σ 2  (z) and
SD 0.41 0.12 0.09 0.06 0.04 0.05 0.09 0.10 0.33 (8.16)
var{δ|z} = σ 2 + σ 4  (z),
where (z) is the log of the marginal density f (z). The afore-
The Poisson estimate cov(c) = diag(c) for the covariance ma- mentioned empirical Bayes estimate ζ̂i is given by
trix of c then yields cov(ˆ  ) = DMG−1 Mt Dt . But because
ζ̂i = zi + σ̂i2 ˆ (zi ), (8.17)
d
ψ  (z) = log{f (z)/ϕ(z)} = z +  (z), (8.9)
dz ˆ = log(f̂ (z)) and σ̂ 2 = σ 2 (zi ), where the variance function
(z) i
σ 2 (·) in (8.14) is calculated numerically. None of this gives an-
we have δ̂(k) ≡ ψ̂  (z = bk ) = bk + ˆk in (4.1), implying that swers significantly different than those derived by using (4.1)
directly, but the transformation effect becomes more important
cov(δ̂) = cov(ˆ  ) = DMG−1 Mt Dt . (8.10)
when n is smaller.
Table 6 gives estimated standard errors for δ̂ [square roots of Remark G (Cross-validation procedure). Both Ebay and
the diagonal elements in (8.10)] calculated for the prostate data. the shrunken centroids procedure default to 10-fold cross-
As in Figure 3, we can see an explosive increase in variability validation replicated R times. Each replication randomly splits
as |z| increases to 4. the N cases into 10 folds, with correctly proportional numbers
of “healthy” and “sick” in each fold. As usual, the prediction
Remark F (Transforming t-values to z-values). The ith row rule is refit 10 times with the cases of each fold withheld from
of x comprises n independent observations the training set in turn, the cross-validated rate α̂CV being the
ind overall proportion of errors on the withheld cases averaged over
xij ∼ N (μi ± σi δi /2c0 , σi2 ) for j = 1, 2, . . . , n (8.11) all R replications. The R replications also provide a standard er-
ror for α̂CV .
in the notation of Section 2, with n1 “−” values and n2 “+” val-
It is useful to remember that α̂CV is not an estimate of er-
ues. The corresponding two-sample t-statistic ti follows a non-
ror for the specific prediction rule selected by Ebay or pamr,
central t distribution with n − 2 degrees of freedom and non-
unlike the actual prediction errors in Figure 4, which were com-
centrality parameter δi ,
puted from knowledge of the simulation structure (4.4). Rather,
x̄2i − x̄1i it is the expected error rate for rules selected according to the
ti = c 0 ∼ tn−2 (δi ). (8.12) same recipe, as emphasized by Efron (1983). In this sense it
σ̂i
differs from the ideal Bayesian estimate α̃ = (−δ̃/2c0 ) fol-
Earlier we treated ti as zi ∼ N (δi , 1), but Ebay actually uses
lowing (3.1) or its empirical Bayes version α̂ = (−δ̂/2c0 ),
transformations that improve the accuracy of Corollary 3. both of which apply directly to the prediction rule at hand.
Let
Remark H (Empirical Bayes estimation of ). The his-
zi = −1 (Fn−2 (ti )), (8.13) togram of off-diagonal elements rii of a sample correlation ma-
trix usually will be more dispersed than the corresponding his-
as in (7.10), so if δi = 0, then zi ∼ N (0, 1). If δi = 0, then zi is
togram of true correlations ρii , because sampling error adds a
still surprisingly close to normal,
component of variance to the rii values. Ebay includes an em-
zi ∼
˙ N (ζi , σ 2 (ζi )) ζi = −1 (Fn−2 (δi )) , (8.14) pirical Bayes shrinkage option to account for overdispersion in
the estimation of , (5.2).
with σ (ζi ) < 1. For example, with δi = 4 and n = 102, zi from Let
(8.13) has (mean, standard deviation, skewness, kurtosis) equal √  
n−4 1 + ρii
to (3.845, 0.931, −0.046, 0.010). A plot of (8.14) superim- νii = log and
posed on (8.12) barely differentiates the two curves. 2 1 − ρii
√   (8.18)
The computation of δ̂i , (4.1) in the Ebay algorithm is actu- n−4 1 + rii
ally carried out using (8.14): vii = log
2 1 − rii
• The vector t = (t1 , t2 , . . . , tN )t is converted component- denote Fisher’s transform of ρii and rii , where the usual con-
wise to z = (z1 , z2 , . . . , zN ), as in (8.13). stant n − 3 has been reduced to n − 4 because two separate
• An estimate f̂ (z) is constructed from z as in Remark D. means are subtracted for the healthy and sick subjects sepa-
• A modified version of Corollary 2, described below, pro- rately. A standard normal theory approximation (Johnson and
vides empirical Bayes estimates ζ̂i . Kotz 1970, chap. 32, section 4), says that
• Finally, transformation (8.14) is inverted to give
vii ∼
˙ N (νii , 1), (8.19)
−1
δ̂i = Fn−2 ((ζ̂i )), (8.15)
implying that the histogram of the vii values will have variance
after which Ebay proceeds as in steps 4–6 in Section 4. about one unit greater than that for the true νii ’s.
Efron: Empirical Bayes Estimates for Large-Scale Prediction Problems 1027

Suppose that the ensemble of true νii values has (mean, vari- to consider δ̂ fixed in (5.2), then δ̂ t  δ̂ is a linear function of
ance) say (M, A), and that vii ∼ (νii , 1) as in (8.19), so that the ’s elements, estimated almost unbiasedly by δ̂ t  ˆ δ̂. The esti-
vii ensemble ∼ ˙ (M, A + 1). Then mation of  would be more crucial if we were attempting to
√ √ implement the general linear discriminant function rather than
ν̃ii = M(1 − C) + Cvii [C = A/(A + 1)] (8.20)
the simplified version (2.3), (2.4).
is the linear function of vii having (mean, variance) ∼ ˙ (M, A). Remark I (Overdispersed z-values). The z-value histogram
Ebay first obtains robust estimates of M and A + 1 from the for the prostate data in Figure 1 is slightly wider than N (0, 1)
set of values {vii }, and then substitutes M̂ and Ĉ = Â/(Â + 1) near z = 0: a fit to the center of the histogram gave
into (8.20) to give estimates ν̃ii . To protect genuine outliers like z∼ ˙ N (0, 1.062 ) [using the locdfr algorithm (Efron 2008)]. This
those in (5.3), Efron and Morris’ (1972) limited translation rule discrepancy is reflected in Figure 2 by the slight upward slope
is enforced: ν̃ii is not allowed to shrink further than one unit of Ê{δ|z} for z between −2 and 1.5. Theorem 1 and Corollary 2
away from vii . Finally, ν̃ii gives ρ̃ii by inverting transformation in Section 3 depend on the assumption z ∼ N (δ, 1). If actually
(8.18). [˜ may no longer be a correlation matrix, but that is not
z ∼ N (δ, σ 2 ), with σ 2 > 1, then the formula for E{δ|z} must
required for use in (5.2).] be modified as in (8.17). We can compensate for overdispersion
A small simulation experiment was run, comparing  ˜ with
by using the values z̃i = zi /1.06 rather than zi in the Ebay al-
the usual (unshrunk) estimate . ˆ The experiment began with gorithm. Doing so flattens Ê{δ|z} to 0 between −2 and 1.5 in
model (4.4), modified to instill correlation among the 5,000 en- Figure 2 and shrinks it slightly toward 0 for larger |z|.
tries in any one column of X. The root mean square of true Figure 8 concerns a leukemia microarray study from Golub
pairwise correlations was set equal to 0.10, about triple that for et al. (1999) where overdispersion is more severe. Here N =
the prostate study and half that for the Michigan lung cancer 7,129 genes were measured on n = 72 subjects in two subtypes,
study of Table 4. Each of 200 replications yielded δ̂ I as in (4.2), n1 = 45 and n2 = 12. Two-sample t-tests gave z-values of zi as
the I × I sample correlation matrix , ˆ and its empirical Bayes in (1.2), (1.3). The histogram of zi ’s corresponding to Figure 1
counterpart .˜ has z ∼ ˙ N (0.9, 1.682 ) near its center.
The corresponding estimates (5.2), Now the curve Ê{δ|z} based on the standardized values z̃i =
(zi − 0.09)/1.68 is much less optimistic than that based on the
ˆ δ̂ I )1/2
η̂ = (δ̂ tI δ̂ I /δ̂ tI  and original zi ’s, especially taking into account the decreased size of
(8.21)
˜ δ̂ I )1/2
η̃ = (δ̂ tI δ̂ I /δ̂ tI  the z̃i ’s. Prediction appears to be extremely easy with the zi ’s;
many genes have |δ̂i | values, (4.1), exceeding 6. However, |δ̂i |
are compared with tops out below 4 for the z̃i ’s. Ebay required only I = 10 genes
to reach target error α0 = 0.01 using the zi ’s, (4.2), compared
ηtrue = (δ̂ tI δ̂ I /δ̂ tI  δ̂ I )1/2 (8.22)
with I = 34 for the z̃i ’s.
in Table 7; η̃ is seen to offer only minor improvement over η̂. Which prediction rule is better? The answer depends on the
Robust estimates of standard deviation for η̃ − ηtrue compared reason for the overdispersion of zi ’s. If in fact zi ∼ N (δi , 1)
with η̂ − ηtrue were a little more decisive: 0.074, compared with and the appearance of overdispersion results from most of the
0.085. Root mean squared errors for estimating all of the el- δi ’s lying far from 0, then the I = 10 rule should perform well.
ements of  strongly favored  ˜ over ,
ˆ rms = 6.30 versus But overdispersion may indicate ephemeral effects due to un-
r
ms = 9.83. observed covariates in an observational study, that will not help
The α̂cor values in Table 2 and Table 4 were based on , ˆ with future predictions, in which case the z̃i analysis is more
Ebay’s default option. Using  ˜ gave smaller estimates of the realistic.
correlation effect in both cases. The choice is not crucial here, Remark J [Model (7.1), (7.2)]. The predictor variable Wi ap-
since the current version of Ebay does not involve α̂cor in con- pears as the response in (7.2), which may seem less natural
structing the prediction rule, but both methods convey useful there than in (2.1). This allows us to express each row xi of
information on the effects of correlation among the predictors. the predictor matrix x as a separate linear regression in y, (7.8),
Regularized estimation of correlation matrixes is a major facilitating the empirical Bayes estimation of parameters βi in
subject in its own right (see Warton 2008), and other methods (7.9)–(7.14). Notice that the correlation structure (7.3) implied
might further improve on . ˆ However,  ˆ performs relatively by (7.1), (7.2) leads directly to (7.4), where now W assumes its
well in our context for two reasons: The dimension “I” of  proper role as a predictor vector.
tends to not be too large, and, more importantly, we need only
estimate the function η, (5.1), not all of . If we are willing Remark K (Naive Bayes prediction). Rule (2.3)–(2.4) is
“naive Bayes” when applied in the correlated framework of
Section 5; that is, it ignores the possible decrease in prediction
Table 7. Estimates η̂ and η̃, (8.21), compared with true correlation
error available from using the complete form of Fisher’s linear
correction factor ηtrue , (8.22); 200 replications of correlated
discriminant function. Such gains are likely to be more hypo-
model (4.4)
thetical than genuine. The theory and simulations in Bickel and
ηtrue η̂ η̃ r
ms r
ms Levina (2004) and Dudoit, Fridlyand, and Speed (2002) show
that our naive Bayes prediction rules outperform more sophis-
Mean 0.597 0.588 0.598 9.83 6.30
ticated predictors in large-scale situations.
SD 0.138 0.171 0.164 9.13 5.74
NOTE: Here “rms values” are the root mean square errors for estimating the elements
of . [Received September 2008. Revised December 2008.]
1028 Journal of the American Statistical Association, September 2009

Figure 8. The solid curve is Ê{δ|z} for leukemia data (Golub et al. 1999); the dashed curve is Ê{δ|z} based on standardized values
z̃i = (zi − 0.09)/1.68. The top row of dashes indicates the 40 most extreme zi values; the lower row, the 40 most extreme z̃i values.

REFERENCES Golub, T. R., Slonim, D. K., Tamayo, P., Huard, C., Gaasenbeek, M., Mesirov,
J. P., Coller, H., Loh, M. L., Downing, J. R., Caligiuri, M. A., Bloomfield,
Benjamini, Y., and Yekutieli, D. (2005), “False Discovery Rate-Adjusted Mul- C. D., and Lander, E. S. (1999), “Molecular Classification of Cancer: Class
tiple Confidence Intervals for Selected Parameters” (with discussion), Jour- Discovery and Class Prediction by Gene Expression Monitoring,” Science,
nal of the American Statistical Association, 100, 71–93. 286, 531–537.
Bickel, P. J., and Levina, E. (2004), “Some Theory of Fisher’s Linear Dis- Hastie, T., Tibshirani, R., and Friedman, J. (2008), The Elements of Statistical
criminant Function, ‘Naive Bayes,’ and Some Alternatives When There Are Learning (2nd ed.), New York: Springer-Verlag.
Many More Variables Than Observations,” Bernoulli, 10, 989–1010. Johnson, N. L., and Kotz, S. (1970), Distributions in Statistics. Continuous Uni-
Brown, L. D. (1971), “Admissible Estimators, Recurrent Diffusions, and Insol- variate Distributions, Vol. 2, Boston, MA: Houghton.
uble Boundary Value Problems,” The Annals of Mathematical Statistics, 42, Senn, S. (2008), “A Note Concerning a Selection “Paradox” of Dawid’s,”
855–903. The American Statistician, 62, 206–210.
Dawid, A. P. (1994), “Selection Paradoxes of Bayesian Inference,” in Multi- Singh, D., Febbo, P. G., Ross, K., Jackson, D. G., Manola, J., Ladd, C., Tamayo,
variate Analysis and Its Applications (Hong Kong, 1992), Hayward, CA: P., Renshaw, A. A., D’Amico, A. V., Richie, J. P., Lander, E. S., Loda, M.,
IMS, pp. 211–220. Kantoff, P. W., Golub, T. R., and Sellers, W. R. (2002), “Gene Expression
Dudoit, S., Fridlyand, J., and Speed, T. P. (2002), “Comparison of Discrimi- Correlates of Clinical Prostate Cancer Behavior,” Cancer Cell, 1, 203–209.
nation Methods for the Classification of Tumors Using Gene Expression Stein, C. M. (1981), “Estimation of the Mean of a Multivariate Normal Distri-
Data,” Journal of the American Statistical Association, 97, 77–87. bution,” The Annals of Statistics, 9, 1135–1151.
Dudoit, S., Shaffer, J. P., and Boldrick, J. C. (2003), “Multiple Hypothesis Test- Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L.,
ing in Microarray Experiments,” Statistical Science, 18, 71–103. Gillette, M. A., Paulovich, A., Pomeroy, S. L., Golub, T. R., Lander, E. S.,
Efron, B. (1983), “Estimating the Error Rate of a Prediction Rule: Improvement and Mesirov, J. P. (2005), “Gene Set Enrichment Analysis: A Knowledge-
on Cross-Validation,” Journal of the American Statistical Association, 78, Based Approach for Interpreting Genome-Wide Expression Profiles,” Pro-
316–331. ceedings of the National Academy of Sciences of the USA, 102, 15545–
(2007), “Size, Power and False Discovery Rates,” The Annals of Sta- 15550.
tistics, 35, 1351–1377. Tibshirani, R., Hastie, T., Narasimhan, B., and Chu, G. (2002), “Diagnosis of
(2008), “Microarrays, Empirical Bayes, and the Two-Groups Model” Multiple Cancer Types by Shrunken Centroids of Gene Expression,” Pro-
(with discussion), Statistical Science, 23, 1–47. ceedings of the National Academy of Sciences of the USA, 99, 6567–6572.
Efron, B., and Morris, C. (1972), “Limiting the Risk of Bayes and Empirical Warton, D. I. (2008), “Penalized Normal Likelihood and Ridge Regularization
Bayes Estimators. II. The Empirical Bayes Case,” Journal of the American of Correlation and Covariance Matrices,” Journal of the American Statisti-
Statistical Association, 67, 130–139. cal Association, 103, 340–349.