On Estimating Many Means, Selection Bias, and the Bootstrap

Introduction

Often, in modern applications, researchers are interested in testing and estimating effect sizes for many different features at once. In the simplest cases one is interested in estimating population means from a sample (often with the most extreme means as the most interesting).

In his revolutionary paper (Stein 1956), Charles Stein showed that, for estimating the means of three or more Gaussian random variables, one can do better (in terms of MSE) than simply using the sample means — that cleverly shrinking the effect sizes will strictly dominate the obvious estimate. Efron and Morris 1973 illustrate this on baseball batting averages and give some insight into the phenomenon. If the population means are very close together, then the sample means will be too spread out — the largest sample mean likely came from a large population mean, but also got lucky and won the “largest statistic” competition, so it is biased high (similarly true for the smallest sample mean, though biased low).

This regression-to-the-mean type idea is concerning for estimating effect sizes in high-throughput experiments. Generally scientists look at the most extreme effects in the sample and report the unadjusted sample estimates of these effect-sizes. As such, on attempts to confirm these effects, effect sizes (even of very statistically significant effects) are often much less extreme than originally reported.

This selection bias has been explored using ideas in empirical Bayes and compound decision theory (Robbins 1956, Robbins 1985 among others), though much of this focus was on asymptotically sub-minimax estimators, rather than selection bias in particular, and before the explosion of high dimensional data and modern computing, this was more a theoretical than practical pursuit.

With new important applications and cheap computing, these problems have gained popularity. Recently Efron 2011 gave an elegant approach to correct for this selection bias by applying empirical Bayes via Tweedie’s formula. Unlike the estimator of James and Stein 1961, which shrinks everything toward the overall mean (ignoring all information beyond the square sum of the statistics), Efron’s formulation gives locally adaptive shrinkage (shrinking based on the local shape of the histogram of statistics). It is particularly appealing as it requires no parametric assumptions on the prior. These ideas have been extended by (Wager 2013, Jiang and Zhang 2009, Brown and Greenshtein 2009, among others), with very efficient estimates of the marginal likelihood.

In this paper, we give a frequentist formulation of the selection bias problem. We show how this bias affects mean square error, and give intuition for classical frequentist shrinkage ideas. We also discuss connections between this “frequentist selection bias” and the more common Bayesian shrinkage in the literature.

Motivated by our definition of frequentist selection bias, we give a simple procedure based on a parametric bootstrap to estimate the frequentist bias. We compare this bootstrap approach to several flavors of empirical Bayes. In situations where empirical Bayes is applicable, the two perform comparably (though empirical Bayes is slightly stronger), however we also detail many situations where empirical Bayes solutions are intractable while our bootstrap shrinkage is simple and effective.

Selection Bias

Suppose we have a large number ($p$) of features, and for each feature ($i$) we have a sample estimate ($z_{i}$) of its mean ($\mu_{i}$) which is normally distributed with variance $1$

This is approximately the scenario we get from many t-tests (as in Tusher et al. 2001 and others). Now, clearly for a fixed $i$, $z_{i}$ is an unbiased estimate of $\mu_{i}$. Generally however we select the largest $z_{i}$ (or $|z_{i}|$) and would like to estimate its corresponding mean. Because we have selected an extreme statistic, if we use the unadjusted statistic as an estimate of the mean, we incur a selection bias.

To explore this bias let us first introduce some notation. Let $z_{[k]}$ denote the $k$-th order statistic, and $i(k)$ denote the index of the $k$-th order statistic (i.e., $z_{i(k)}=z_{(k)}$). Note that since the ordering of our statistics is stochastic, $i(k)$ is a stochastic index (the inverse rank of $z_{(k)}$). The bias we are interested in is

Note, both the order statistic $z_{}$ and the mean $\mu_{i(1)}$ are random variables (the mean is a random variable because the index $i(1)$ is stochastic).

For each rank, $k$, we can use the same idea and define our bias as

If we knew these biases (unrealistic in practice) then we could estimate the means of the extreme statistics by

Our risk decreases by the sum of the squared biases. For the remainder of the manuscript, we will refer to the estimates in Eq (3), with the true biases known, as our “oracle estimates.”

For ease of reading we will define the following notation. We will use $\boldsymbol{\mu}$ to refer to a vector of means, with $\mu_{i}$ to denote the $i$-th element of $\boldsymbol{\mu}$. Let $\beta\left(\boldsymbol{\mu}\right)$ denote the vector of biases for a given mean vector $\boldsymbol{\mu}$. More specifically

In practice we will never know the bias (2) and must estimate it. We propose a simple $2$ step method. We first estimate $\boldsymbol{\mu}$ by maximum likelihood, giving, in this case $\hat{\boldsymbol{\mu}}=z$. We then use the biases for this estimated model, as estimates of the bias for our original model

2 Calculating the Bias of the Estimated Model

Now, given known means $\mu_{i}$, $i=1,\ldots,p$ we need to calculate the bias. This is most tractable by monte-carlo. Though the monte-carlo is straightforward, we give it in full detail.

Simulate $z_{1}^{b},\ldots,z_{p}^{b}$ a $p$-vector of Gaussians with means $\mu_{1},\ldots,\mu_{p}$ and variance 1

Find the $k$-th order statistic $z_{(k)}^{b}$ and the index of its corresponding mean $i(k)^{b}$

As $B$ grows, this will give increasingly accurate calculations of the bias.

3 Second Order Bias

In many cases we can improve the estimate in (4). Heuristically the further spread out the means are, the smaller the bias is. Intuitively, our sample means are more spread out than the true means — thus our estimates of the bias, are themselves biased. Often, the bias estimate in the bootstrap sample will be smaller than the true bias.

If we knew this quantity, then we could update our estimate $\hat{\beta}\left(\boldsymbol{\mu}\right)$ to

Unfortunately, this quantity is never known in practice. To approximate it, we use the same trick as before

It is straightforward to calculate the estimate of second-order bias in (5) via monte-carlo, though, for the sake of brevity, we leave the details to the reader. For the remainder of the manuscript we will refer to

as the “first-order” and “second-order” bootstrap estimates of effect-size.

This second-order bias estimate $\widehat{\beta\beta}\left(\mu\right)$ will also be biased for the true second-order bias $\beta\beta\left(\mu\right)$. One might consider higher order de-biasing. There is a “bias-variance” tradeoff here, and in practice, while the second-order correction has been useful in some problems, we have not seen improvement past second-order corrections.

4 Simple Example

We will give a simple example illustrating the difference between the “local” shrinkage of our method and the “global” shrinkage of James-Stein estimation. In our example we simulate $1000$ features, $990$ of which have $\mu_{i}=0$; the remaining $10$ have $\mu_{i}=6$. We can see the results in Figure 1. James-Stein gives great shrinkage for the bulk of the features, however it completely shrinks away the interesting effects. In contrast, our resampling approach can take local behavior into account, and while it shrinks the bulk of the effect sizes to $0$, it correctly pushes the estimates for the interesting effect towards $6$.

5 Empirical Bayes

We will briefly review the empirical Bayes approach of Efron 2011. If we assume some prior density ($g$) on the means

and let $f(z)$ denote the marginal distribution of $z$

then the posterior expectation of $\mu$ given $z$ is

This result is known as Tweedie’s theorem (Efron 2011). We will term $f^{{}^{\prime}}(z)/f(z)$ as the Bayesian bias — it adjusts the naive estimate and accounts for selection bias (this will be made more clear in Section 2.6). This approach is elegant as (9) does not require $g$, or a direct estimate of $g$. Instead, one needs only estimate $f$: a smoothed histogram of the $z_{i}$ (and its derivative). This is tractable if the number of features is large — though estimating the derivative is still difficult, and the degree of smoothing can influence results.

Efron suggests using Lindsey’s method (Efron and Tibshirani 1996) for the density estimate. For Lindsey’s method, one bins the data, and uses poisson regression with a spline or polynomial basis and an offset for bin-size to estimate the density. More recently Wager 2013 and others have given nonparametric approaches which generally outperform Lindsey-based approaches. In our comparisons in Section 4, we include the estimate of Wager 2013, which we will refer to as nlpden.

6 Reconciling the differences

We have two approaches to selection bias (frequentist and Bayesian) that at first glance are very different however on deeper consideration they are actually quite similar. The Bayesian bias that we estimate in the Bayesian approach is

whereas in the frequentist approach the bias is

Though frequentist, (11) already has a Bayesian flavor as $i(k)$, the index of our mean, is stochastic. In the Bayesian framework (where $\mu$ has some prior distribution), we can show that (under some assumptions) these two biases are asymptotically the same.

Then, as $p\rightarrow\infty$, for any $t\in(0,1)$ we have

where $\lfloor\cdot\rfloor$ is the “floor” function and $F^{-1}\left(t\right)$ is the $t$-th quantile of the marginal distribution $dF(z)=\int\phi\left(z-\mu\right)dG(\mu)$.

The proof is given in the appendix. The assumption of bounded support for $G$ can easily be weakened, but it simplifies the proof.

This raises a similar question in the frequentist framework. Namely, if we denote the empirical distribution of the means by $G_{p}$:

and $G_{p}\rightarrow G$, for some well-behaved $G$, then can we treat $G$ like a prior distribution, and get the same asymptotic equivalence? We believe the answer is yes (and simulations support this). The mathematics for this problem becomes significantly more complex, and is beyond the scope of this manuscript. Questions of this flavor (trying to treat the empirical distribution of the means as a prior) have been explored in the compound decision theory literature (though none that we have seen use this formal definition of frequentist selection bias).

Extensions to General Problems

An important strength of the bootstrap approach is that it is not restricted to this Gaussian scenario — the bootstrap is flexible and can be applied to more complex problems intractable for empirical Bayes. Consider the more general scenario, where the $z$ vector has a joint distribution parametrized both by our parameters of interest $\boldsymbol{\mu}$, and some nuisance parameters $\boldsymbol{\Theta}$.

for some known functional form $F$. We can generalize our bias from before as

where $\hat{\boldsymbol{\mu}}$ and $\hat{\boldsymbol{\Theta}}$ are maximum likelihood estimates. We illustrate this on estimating $\rho^{2}$ values in regression with categorical variables in Section 4, but it can be further applied to a wide variety of problems involving non-gaussian distributed statistics, dependence, and more complicated model-based estimates.

Empirical Results

We give empirical results for the bootstrap method both for the simple gaussian scenario as well as the more complicated categorical variable regression problem. In the gaussian scenario, we compare the performance of the bootstrap to empirical Bayes methods, on real and simulated data. We see that both empirical Bayes and the bootstrap effectively reduce selection bias (in some cases quite drastically). The second order bootstrap outperforms the first order, and is comparable though somewhat outperformed by the best empirical Bayes estimates (nlpden). However, these empirical Bayes methods cannot handle the categorical variable regression problem, while the bootstrap approach still performs well there.

We know that the potential gain of these procedures is based on the bias of our rank estimates, which in turn is based on the spacing of the means. We consider $6$ simulated scenarios with varying mean spacings to explore the bias estimation of the procedures in different regimes.

All $\mu_{i}=0$ — hypothesis testing with a global null.

$500$ of the $\mu_{i}=0$, and $500$ of the $\mu_{i}=6$ — a simple mixture model.

$900$ of the $\mu_{i}=0$, and $100$ of the $\mu_{i}=6$ — hypothesis testing with a strong clustered set of alternatives

$900$ of the $\mu_{i}=0$ and $100$ of the $\mu_{i}\sim N(0,2)$ — hypothesis testing with a weak diffuse set of alternatives

$5$ clusters each with $200$ features. In the $j$-th cluster ($j=1,\ldots,5$) all $\mu_{i}=6*j$ — separated mixture model.

Performance can be see in Table 1. While competitive with empirical Bayes, the bootstrap is outperformed by the nonparametric estimates of nlpden. This table also illustrates the improvement from our second stage of bootstrap debiasing.

We also have plots showing the shrinkage from the bootstrap and empirical Bayes for scenarios $1$, $3$,$4$,$6$ (from left to right and top to bottom). We can see in Plot 2 that the bootstrap and empirical Bayes estimates are very similar. Though entirely unrealistic as an applied scenario, the results in scenario $6$ for both the bootstrap and nlpden are particularly neat.

2 Non-Gaussian Statistics

One of the major strengths of our procedure is the ability to deal with nongaussian statistics. We now give a simulated example; estimating $\rho^{2}$, the coefficient of determination, in the regression of a continuous variable on a $3$-level categorical variable. We simulated an $n\times p$ matrix $X$ of features, and an $n\times p$ matrix $Y$ of responses (each feature had its own, response). The entries in each column of $X$ were sampled with equal probability from the $3$ classes (though we ensured at least $2$ observations per class in each column). The $\rho^{2}$ values for each regression were selected based on various schemes, and each entry of $Y$ was independently simulated as

where $\epsilon_{ij}\sim N(0,1)$, and $\boldsymbol{\theta}_{j}=\{\theta_{j0},\theta_{j1},\theta_{j2}\}$ were chosen so that the regression had the prespecified $\rho^{2}$-value. We apply the general method (with nuisance parameters) detailed in Section 3. More specifically, we estimate our regression coefficients (and the variance of $\epsilon_{\cdot j}$), and calculate the bias from those estimated models.

For these examples we used $p=1000$ features with $n=50$ observations. The schemes used for choosing $\rho^{2}$ were:

$\rho^{2}\sim\operatorname{exponential}\left(10\right)$ (with values truncated at $0.99$)

$800$ of the $\rho^{2}\sim\operatorname{exponential}\left(20\right)$, and $200$ $\rho^{2}\sim N(0.55,1/20)$ (again with truncation at $0.99$)

Performance of the bootstrap methods on these examples can be seen in Table 2. While perhaps not as extreme as the improvement in our gaussian examples, the bootstrap shrinkage does still show substantial gain over the naive estimates. We illustrate the shrinkage plots for scenarios $2$ and $3$ in Figure 3. As we see from both Table 2 and Figure 3, our estimates are close to the oracle estimates.

3 Real Data

While simulated data can give insight, it can also be misleading as the efect-sizes for real data rarely follow our simplified simulation schemes. To compare the behavior of the bootstrap and empirical Bayes approaches on real data, we applied both methods to the prostate data of Singh et al. 2002. This data has $102$ samples, $52$ from prostate cancer tumors and $50$ from healthy tissue. For each sample, there are $6033$ gene expression measurements. For each of these $6033$ features we calculated a two sample $t$-statistic (which we transformed to be approximately normal under the null).

As shown in Figure 4 our bootstrap and empirical Bayes estimates are similar, though the empirical Bayes methods provide more shrinkage in the center. Unfortunately there is no gold standard for this data, so it is hard to compare the performance of the methods beyond this.

Discussion

Large scale estimation of effect-sizes has become of increased importance, in particular in biomedical and epidemiological studies with the development of whole genome assays like gene expression microarrays, genotyping arrays and next generation DNA sequencing. Typical studies attempt to identify univariate differences between disease cases and control. Methods for control of false discoveries are commonly used in such studies, but the effect sizes of the features selected as important are rarely adjusted for selection. This can result in a misleading assessment of the importance of the feature and an overestimation of its value for predictive purposes.

In light of this, shrinkage estimators are becoming more and more important. Classical global shrinkage (as in James-Stein type estimators) unfortunately tends to overshrink the estimated size of the most extreme (and therefore interesting) effects. Modern shrinkage estimators (such as our proposal and non-parametric empirical Bayes) leverage local information to give data-adaptive estimates, with much better performance for the effect-sizes of “interesting” effects.

In this paper we give a general formalism for frequentist selection bias. We show how this formalism connects to Bayesian ideas for removing selection bias. We propose a resampling estimator for estimating this frequentist selection bias which leads to accurate estimates of effect-size based on locally adaptive shrinkage. Unlike empirical Bayes methods which must be hand-crafted for each scenario, our approach is general and applicable in a wide array of problems. Our method is widely applicable, simple to apply, and performs well in practice.

Appendix

Here we will give a proof for Theorem 2.1.

To begin we note that sample quantiles converge to population quantiles, or

so, because $G$ has finite support, we have that: $\operatorname{E}\left[z_{[\lfloor pt\rfloor]}\right]\rightarrow q_{t}$.

Continuing to the second term, using the fact that $\mu_{i}$ are iid, we have

The last line follows through concentration of measure (combined with the finite support of $G$), because, as we have shown, $z_{[\lfloor pt\rfloor]}\rightarrow q_{t}$, thus completing the proof $\blacksquare$

References