Rank conditional coverage and confidence intervals in high dimensional problems

TL;DR

This paper introduces rank conditional coverage (RCC), a new criterion for confidence intervals in high-dimensional settings, and proposes bootstrap-based methods to achieve accurate coverage conditioned on parameter ranking.

Contribution

It proposes the RCC criterion for confidence intervals in high-dimensional problems and develops bootstrap methods to ensure coverage aligns with the nominal level based on parameter rank.

Findings

Bootstrap intervals achieve coverage close to the nominal level

RCC provides a more accurate assessment of interval performance in high dimensions

The R package rcc implements the proposed methods

Abstract

Confidence interval procedures used in low dimensional settings are often inappropriate for high dimensional applications. When a large number of parameters are estimated, marginal confidence intervals associated with the most significant estimates have very low coverage rates: They are too small and centered at biased estimates. The problem of forming confidence intervals in high dimensional settings has previously been studied through the lens of selection adjustment. In this framework, the goal is to control the proportion of non-covering intervals formed for selected parameters. In this paper we approach the problem by considering the relationship between rank and coverage probability. Marginal confidence intervals have very low coverage rates for significant parameters and high rates for parameters with more boring estimates. Many selection adjusted intervals display the same pattern. This connection motivates us to propose a new coverage criterion for confidence intervals in multiple testing/covering problems --- the rank conditional coverage (RCC). This is the expected coverage rate of an interval given the significance ranking for the associated estimator. We propose interval construction via bootstrapping which produces small intervals and have a rank conditional coverage close to the nominal level. These methods are implemented in the R package rcc.

Figures

Equations

$$\text{mbox}ConditionalCoverageProbability=\frac{1}{p}\sum _{j=1}^{p}P[{\theta }_j\in C{I}_j|j\in S]$$

$$\text{mbox}FalseCoverageStatementRate=$$

$$Q=$$

$${\hat{\theta }}_{s(1)}\ge \cdots \ge {\hat{\theta }}_{s(p)}.$$

$${\text{Rank Conditional Coverage}}_i=$$

$$=$$

$$p({\theta }_i|{\hat{\theta }}_i)\propto p({\hat{\theta }}_i|{\theta }_i)p({\theta }_i).$$

$$\displaystyle\Big{(}\hat{\theta}-H^{-1}(1-\alpha/2),\ \hat{\theta}-H^{-1}(\alpha/2)\Big{)}$$

$${\delta }_{[i]}={\hat{\theta }}_{s(i)}-{\theta }_{s(i)}$$

$$\displaystyle CI_{s(i)}^{\text{exact}}=\Big{(}\hat{\theta}_{(i)}-H_{[i]}^{-1}(1-\alpha/2),\ \hat{\theta}_{(i)}-H^{-1}_{[i]}(\alpha/2)\Big{)}$$

$$P[{\theta }_{s(i)}\in C{I}_{s(i)}^{\text{exact}}]=1-\alpha$$

$$P[a\le {\theta }_{s(i)}\le b]=$$

$$=$$

$${\tilde{\delta }}_{k,[i]}={\vartheta }_{k,s_k(i)}-{\theta }_{s_k(i)}$$

$$CI_{s(i)}^{\text{oracle}}=\Big{(}\hat{\theta}_{s(i)}-\tilde{H}^{-1}_{[i]}(1-\alpha/2),\ \hat{\theta}_{s(i)}-\tilde{H}^{-1}_{[i]}(\alpha/2)\Big{)}$$

$$\displaystyle CI_{s(i)}^{\text{boot}}=\Big{(}\hat{\theta}_{s(i)}-\hat{H}^{-1}_{[i]}(1-\alpha/2),\ \hat{\theta}_{s(i)}-\hat{H}^{-1}_{[i]}(\alpha/2)\Big{)}$$

$${\hat{\delta }}_{k,[i]}={\vartheta }_{k,s_k(i)}-{\hat{\theta }}_{s_k(i)}$$

$$\displaystyle CI_{s(i)}^{\text{np-boot}}=\Big{(}\hat{\theta}_{s(i)}-\hat{H}^{-1}_{emp[i]}(1-\alpha/2),\ \hat{\theta}_{s(i)}-\hat{H}^{-1}_{emp[i]}(\alpha/2)\Big{)}$$

$${\hat{\delta }}_{k,[i]}={\vartheta }_{k,s_k(i)}-{\hat{\theta }}_{s_k(i)}$$

$${\mathbf{x}}_i\sim N_{1000}(0,\Sigma )$$

$$y_i={\mathbf{x}}_i\mathbf{\beta }+{ϵ}_i{ϵ}_i\sim N(0,1)$$

$${\beta }_j=$$

$$\left(E[y_i|tr{t}_i=1,w_i\le c_j]-E[y_i|tr{t}_i=0,w_i\le c_j]\right).$$

$$y_i={\beta }_0+{\beta }_{1}trt+{\beta }_{2,j}{1}_{w_i>c_j}+{\beta }_j{1}_{w_i>c_j}+{ϵ}_i,$$

$$E[y_i|tr{t}_i,w_i]=\{\begin{array}{ll}0& w_i<0.5\\ (w_i-\frac{1}{2})\cdot tr{t}_i& w_i\ge 0.5\end{array}.$$

$$FCR=$$

$$=$$

$${\hat{\delta }}_{k,[i]}={\vartheta }_{k,(i)}-{\tilde{\theta }}_{s_k(i)}$$

$$|{\vartheta }_{(1)}|\ge \cdots \ge |{\vartheta }_{(p)}|$$

$${\hat{\delta }}_{k,[i]}=\text{sign}({\vartheta }_{k,(i)})({\vartheta }_{k,(i)}-{\theta }_{s_k(i)})$$

$$\displaystyle\begin{cases}\Big{(}\hat{\theta}_{(i)}-\hat{H}^{-1}_{[i]}(1-\alpha/2),\ \hat{\theta}_{(i)}-\hat{H}^{-1}_{[i]}(\alpha/2)\Big{)}\qquad&\hat{\theta}_{s(i)}>0\\ \Big{(}\hat{\theta}_{(i)}+\hat{H}^{-1}_{[i]}(\alpha/2),\ \hat{\theta}_{(i)}+\hat{H}^{-1}_{[i]}(1-\alpha/2)\Big{)}\qquad&\hat{\theta}_{s(i)}\leq 0\end{cases}$$

Full text

In this paper we approach the problem by considering the relationship between rank and coverage probability. Marginal confidence intervals have very low coverage rates for significant parameters and high rates for parameters with more boring estimates. Many selection adjusted intervals display the same pattern. This connection motivates us to propose a new coverage criterion for confidence intervals in multiple testing/covering problems — the rank conditional coverage (RCC). This is the expected coverage rate of an interval given the significance ranking for the associated estimator. We propose interval construction via bootstrapping which produces small intervals and have a rank conditional coverage close to the nominal level. These methods are implemented in the R package rcc.

Keywords: bootstrap/resampling; Multiple comparisons; Selective inference; Winner’s curse;

1 Introduction

In many fields including genomics, proteomics, biomedical science, and neurology it is now common to conduct “high dimensional” studies in which thousands or millions of parameters are estimated. Often, one of the main goals of these studies is to select a small subset of features for description and future investigation. Estimates of selected parameters are often reported unadjusted and, when confidence intervals are not entirely omitted, either marginal or Bonferroni-corrected intervals are given. Many previous authors have demonstrated the undesirable features of these practices (Efron (2011), Sun and Bull (2005), and Simon and Simon (2013) among others). In particular, features selected for reporting are usually chosen because they have the largest or most significant estimates; unfortunately these most extreme estimates are also highly biased. This is known informally as the “winner’s curse”: Large statistics tend to come from large parameters, but they also tend to be large by chance.

There is a related phenomenon for confidence intervals: Marginal confidence intervals almost always fail to cover parameters associated with the most significant estimates because they fail to account for the bias of these estimates. This results in overly short intervals that are too far from zero. Often this problem is recognized by investigators, but the most commonly used alternative, the Bonferroni correction, yields enormous, uninformative intervals. Problems with the Bonferroni correction have been described by Benjamini and Yekutieli (2005), Zhong and Prentice (2008), Weinstein et al. (2013) and several others. This method does not recenter the intervals and inflates their size symmetrically in order to control the family-wise error rate. Intuitively, we know that the largest estimates are more often too large than too small so most of the upper extension provided by the Bonferroni confidence intervals is unnecessary. Furthermore, the family-wise error rate is a much more conservative criterion than is typically desired.

We will discuss several alternative coverage criteria for high dimensional settings and introduce rank conditional coverage (RCC). Intuitively, RCC is the expected coverage probability of an interval given the ranking of its corresponding estimate. The RCC captures the idea that, when many parameters are estimated, the rank of an estimate provides information about its bias and the coverage probability of the associated confidence interval. This criterion applies to all high dimensional studies regardless of whether or not selection is performed but is particularly relevant when rank is used as a selection criterion. One advantage to obtaining confidence intervals that control the RCC, rather than selection adjusted confidence intervals, is that these intervals are not dependent on a particular selection procedure — i.e. the interval calculated for a particular estimate will not change if the selection threshold is moved. We discuss several approaches for obtaining useful confidence intervals — these include a parametric and non-parametric bootstrap approach — and explore their behavior and RCC in a few common settings.

1.1 Coverage criteria after selection

Consider an analysis in which we would like to obtain estimates and confidence intervals for many parameters $\theta_{1},\dots,\theta_{p}$. Suppose we have point estimates for each parameter $\hat{\theta}_{1},\dots,\hat{\theta}_{p}$ and a marginally-valid procedure for constructing confidence intervals. By marginally-valid we mean that for each parameter $\theta_{j}$, the coverage probability of the $\alpha$ level confidence interval $CI_{j}$ satisfies $P[\theta_{j}\in CI_{j}(\alpha)]=1-\alpha$.

In the high dimensional setting, the marginal confidence interval will control the average coverage — i.e. if we construct 90% marginal confidence intervals, we can expect them to cover 90% of the parameters. Typically, however, the entire set of parameter estimates is not of interest and only the most significant estimates are reported. Benjamini and Yekutieli (2005) demonstrate that, for common selection rules, the expected rate of coverage within a selected subset of parameters will be much lower than the desired level. The marginal confidence interval achieves its average coverage by under-covering parameters associated with the most extreme or significant estimates and over-covering more boring parameters.

One might instead want to control the conditional coverage probability: the average coverage rate of selected parameters. More formally, Let $S$ be a set of indices for selected parameters (note that $S$ will be stochastic). The conditional coverage probability is defined as

[TABLE]

Unfortunately, it is not possible to construct confidence intervals guaranteed to control this criterion under many common selection procedures. For example, if selection is based on excluding 0 from the confidence interval but $\theta_{1}=\dots=\theta_{p}=0$, the conditional coverage probability will be 0 for any set of confidence intervals.

As an alternative, Benjamini and Yekutieli (2005) propose the concept of the false coverage statement rate (FCR). The FCR measures the probability of making a false coverage statement. This concept is very similar to conditional coverage probability. There are two differences: 1) it averages coverage over the selected set, and 2) if nothing is selected, it is counted as “no false statements” being made. Benjamini and Yekutieli (2005) also restrict application of FCR to intervals which are constructed after selection.

[TABLE]

To control FCR, Benjamini and Yekutieli (2005) give a procedure which, like the Bonferroni procedure, symmetrically inflates the size of marginally-valid intervals for all parameters. When selection is based on parameter estimates exceeding a threshold, these are equivalent to coverage $(1-\frac{|S|\alpha}{p})$ marginal intervals. This uniform inflation can be excessive in some cases. For example, if one parameter is very large, it will nearly always be selected. Thus, there is no need to inflate that interval at all (a $1-\alpha$ marginal interval will still control FCR). Additionally, the fact that highest ranked estimates are more often too large than too small suggests that confidence intervals should have longer tails extending towards the bulk of the estimates than extending away from the bulk.

These issues are partially accounted for in more recent literature: Zhong and Prentice (2008), Weinstein et al. (2013), and Reid et al. (2014) all propose FCR controlling intervals which return to the marginal interval for very large parameter estimates. Zhong and Prentice (2008) and Weinstein et al. (2013) both condition on selecting all estimates larger than a (possibly data-dependent) cutoff. Zhong and Prentice (2008) use a likelihood-based approach to obtain asymptotically correct FCR while Weinstein et al. (2013) calculate exact intervals under the assumption that parameter estimates are independent with a known symmetric unimodal distribution. Reid et al. (2014) condition on the identity of the selected set and construct exact intervals for finite sample sizes assuming that parameter estimates are drawn from independent Gaussian distributions. All three of these intervals are asymmetric about the original point estimate.

Despite these advances, controlling FCR remains an unsatisfying solution to confidence interval construction for high dimensional problems. FCR controlling methods achieve the correct coverage rate within a subset in the same way that the marginal intervals achieve the correct rate in the larger set — with under-coverage of the (more interesting) highest ranked parameters and over-coverage of (less interesting) more moderately ranked parameters. We illustrate this pattern through a simple example in Section 1.5.

1.2 Rank conditional coverage

In Section 1.1 we observed that, for unadjusted confidence intervals as well as the selection adjusted alternatives, the rank of an estimate is informative about the probability that the associated confidence interval covers its target (see also Section 1.5). We find this phenomenon undesirable since it means that parameters associated with top ranked estimates are covered at a much lower rate than parameters associated with less significant estimates. Additionally, this observation indicates that there is an opportunity to use more information and construct better intervals.

We first introduce the concept of rank conditional coverage (RCC) as a way to quantify the relationship between rank and coverage probability. In the majority of cases, the most interesting ranking of parameters is based on either the size of an associated test statistic or a $p$-value. In general, we assume that we have some ranking function $s$ where $s(i)$ gives the index of the $i$th ranked estimate. For example, if we are ranking simply based on the size of estimates, then

[TABLE]

In this paper, we will use the convention that a smaller rank indicates that an estimate is more significant, so the most significant estimate will have rank 1. In our examples, we focus on simple common rankings but the RCC could be defined for any scheme, and in fact the ranking scheme need not give a rank to every estimate. For example, if the parameter estimates can be grouped into highly correlated subsets such as LD blocks, we might choose the most significant estimate in each block and rank only this selected set. This type of ranking scheme is discussed at greater length along with simulation results in Section 5.3 of the Appendix.

We define the RCC at rank $i$ of a set of confidence intervals $CI_{1}\dots CI_{p}$ as

[TABLE]

This quantifies how often the interval formed around the $k$-th ranked estimate contains its target parameter. This is an appealing criterion, since we have a strong interest in ensuring that intervals around our most promising candidate features contain their targets. Something to note here is that we are not conditioning on which specific features achieve a given rank. Rather, we are averaging over all features (weighted by their probability of achieving that rank). While FCR summarizes the average coverage of a confidence interval procedure applied to a set of selected parameters in a single number, RCC gives a separate estimate of coverage probability for each rank and is not directly related to a selection procedure.

1.3 Implications of Controlling RCC

Intervals that control RCC do not provide guarantees for particular parameters. For example, suppose $\theta_{1}$ is of special interest. If we use an RCC controlling method with $\alpha=0.1$, we cannot say that, if the experiment were repeated many times, $\theta_{1}$ would be contained in $CI_{1}$ in 90% of experiments. Thus, if there is particular prior interest on one or a few parameters, RCC is not the correct criterion to control.

A correct statement that could be made about intervals controlling the RCC is that, if the experiment were repeated many times, we expect the parameter corresponding to the top ranked estimate to be contained in its interval in 90% of experiments. A similar statement could be made for any rank. While this property may seem less intuitive on its surface, it has important implications when parameters are selected based on rank or significance.

For example, suppose that a researcher publishes the results of many genome-wide association studies, each time reporting the most significant effect size estimates. If these estimates were paired with confidence intervals controlling the RCC at 90%, the researcher could expect that 90% of the published intervals for 1st ranked estimates (or 2nd etc.) averaging over studies contain their parameters. Most followup studies are conducted specifically for the most promising parameters, so this is precisely the type of guarantee needed to ensure these followup studies are worthwhile.

This guarantee is stronger than that made by the FCR. In fact, it is straightforward to see that confidence intervals which control RCC for every rank also control FCR for selection rules that choose the most significant parameters based on the same ranking used to define RCC (see Theorem 5.1 of the Appendix).

One major strength of the RCC over the FCR is that RCC controlling confidence intervals can be divorced from the selection procedure. For example, if FCR controlling intervals are published for the top 10 parameter estimates but we are only able to follow up on the top 5 then we will need to recompute new, wider intervals in order to guarantee coverage within the smaller set. Using RCC controlling intervals, the same intervals remain valid regardless of how many parameters are selected.

1.4 Relationship of RCC to Empirical Bayes Approaches

The observation that motivates the RCC is that, in a study estimating many parameters, the full set of estimates can provide information about the true underlying parameter values. This is the same idea that motivates empirical Bayes (EB) approaches to simultaneous inference problems. In a Bayesian paradigm, we are interested in estimating the posterior distributions of $\theta_{1},\dots,\theta_{p}$ which, assuming conditional independence of the estimates and using Bayes rule, we can express as

[TABLE]

The idea of EB approaches such as those of Efron (2008) and Stephens (2016) is to assume a theoretical distribution for $\hat{\theta}_{i}|\theta_{i}$ and use the large number of parameter estimates to estimate the prior $p(\theta_{i})$. For example, in the ashr method proposed by Stephens (2016), $p(\theta_{i})$ is assumed to be unimodal and centered at zero and, in one of several proposed variations, is estimated as a mixture of normal distributions. Provided the EB modeling assumptions hold, we can expect that, averaging over many realizations, the $1-\alpha$ EB credible intervals contain the true parameter $1-\alpha\%$ of the time and are immune to selection bias. That is, in a Bayesian system where a single realization of an experiment also includes resampling the parameter values, EB credible intervals should control the RCC.

The bootstrapping approaches we describe in Section 2 differ from EB methods in that they require fewer modeling assumptions and are derived from a frequentist perspective. We tend to view the parameters $\theta_{i}$ as fixed and use the large number of estimates to learn about the distribution of the bias $\hat{\theta}_{s(i)}-\theta_{s(i)}$. This method requires no assumptions about the form of $p(\theta_{i})$. In the non-parametric version, we are also able to avoid assumptions about the form of $p(\hat{\theta}_{i}|\theta_{i})$, provided we have access to the individual level data used to produce the original estimates. The flexibility of the non-parametric method does come at the expense of increased computational effort. The parameteric bootstrap can be fairly efficient and in the example in Section 1.5 is 8 times faster than ashr. The non-parametric bootstrap can be quite costly since we must be willing to repeat the entire analysis hundreds of times. However, in cases in which the parameter estimates are not independent or the theoretical distributions of test statistics are poor approximations, the non-parametric bootstrap is the most appropriate choice.

1.5 Example

Consider 1000 independent estimates $Z_{i}$ is drawn from a $N(\theta_{i},1)$ distribution and ranked according to their absolute value. Figure 1 shows the average RCC over 100 simulations for the top 20% of statistics for several different configurations of parameter values:

All parameters equal to zero. 2. 2.

All parameters small and non-zero: $\theta_{i}$ generated from a $N(0,1)$ distribution but fixed for all simulations. 3. 3.

A few large non-zero parameters: $\theta_{i}=3$ for $i=1\dots 100$ and $\theta_{i}=0$ for $i>100$. 4. 4.

A few small non-zero parameters: $\theta_{i}$ drawn from a $N(0,1)$ distribution but fixed over all simulations for $i=1\dots 100$ and $\theta_{i}=0$ for $i>100$.

The standard marginal confidence intervals are $CI_{i}=Z_{i}\pm\Phi^{-1}(1-\alpha/2)$. In configuration 1, this interval has an RCC of $\sim 0$% for the 65 most extreme observations but an RCC of $\sim 100$% for statistics closer to the median giving an overall average of 90% coverage.

We see a similar pattern in the intervals constructed by Weinstein et al. (2013) and Reid et al. (2014). Both provide intervals only for a selected subset of parameters (we selected the most extreme 10% of observations). Both methods control FCR but do so by under-covering parameters associated with the most significant statistics and over-covering parameters with more moderate statistics. In settings 2 and 3, the intervals of Reid et al. (2014) have poor RCC for both the most and least extreme parameters considered. The credible intervals generated by the ashr method of Stephens (2016) do control the RCC in this setting and can sometimes be very small. ashr is able to achieve a very small average interval width in settings 1, 3, and 4 because it attempts to shrink parameter estimates to zero. If the posterior probability that the parameter is equal to zero is larger than the desired level, the resulting credible interval will simply be $[0,0]$.

Figure 1 also shows the results of the parametric bootstrapping method described in Section 2. This method provides an RCC close to the nominal level for all ranks and, in most cases, a shorter interval. We see under-coverage for a small set of parameters in scenarios 3 and 4 in which most parameters are zero and a handful take on larger values. In all four scenarios, the bootstrapped confidence intervals are also shorter than either of the selection adjusted methods and sometimes shorter than the ashr credible intervals.

The deviations from the nominal level of RCC are a result of using estimated mean values to generate bootstrap samples. Hypothetically, if these values were known, we could produce an “oracle” estimate which, with enough Monte-Carlo samples, would achieve exactly the desired RCC for all order statistics. The oracle is shown in Figure 1 and provides the motivation for the methods described in Section 2 where it is discussed in more detail. An in depth walk-through of these results as well as code for replicating Figure 1 is available in https://github.com/jean997/rccSims/walkthroughs/compare_cis.pdf.

We explore the performance of these methods in two simulation studies designed to mimic common high dimensional analyses in Section 3. In the first, we simulate an outcome and many features and estimate marginal effect sizes via single-variable linear regression. This is a common approach to genome-wide association studies. In the second, we estimate the difference in average treatment effect for subsets of individuals defined the level of a biomarker. These estimates are highly correlated making the non-parametric bootstrap the most appropriate method.

2 Parametric and Non-Parametric Bootstrapping to Build Confidence Intervals

2.1 Rank conditional confidence intervals

First consider estimating a single parameter of a single distribution $\theta=T(F)$. Let $\delta=\hat{\theta}-\theta$ and $H(x)=P[\delta\leq x]$ be the cdf of $\delta$. When $H$ is known, a pivotal exact $1-\alpha$ confidence interval can be constructed as

[TABLE]

In the high dimensional setting, we are attempting to estimate $p$ parameters $\theta_{i}=T_{i}(F)$. We can construct a rank conditional analog of the classical pivotal interval in (5). Define $\hat{\theta}_{s(i)}$ as in section 1.2 where $s(i)$ gives the index of the $i$th ranked parameter estimate. We define the bias of the estimates at each rank

[TABLE]

where the subscript $[i]$ indicates ranked based indexing. Let $H_{[i]}=P[\delta_{[i]}\leq x]$ be the cdf of $\delta_{[i]}$. Were $H_{[i]}$ known, an exact $1-\alpha$ confidence interval for $\theta_{s(i)}$ could be constructed as

[TABLE]

We note that the rank conditional intervals are not pivotal because the distribution of $\delta_{[i]}$ depends on $\theta_{1},\dots,\theta_{p}$. This makes them more difficult to obtain when $H_{[i]}$ are unknown but doesn’t impact the coverage probability of (7).

Lemma 2.1**.**

The intervals in (7) have exact $1-\alpha$ coverage:

[TABLE]

Proof.

This proof is identical to the proof for the classical interval in (5) given by Wasserman (2005) among others. Let $a=\hat{\theta}_{s(i)}-H_{[i]}^{-1}(1-\alpha/2)$ and $b=\hat{\theta}_{s(i)}-H_{[i]}^{-1}(\alpha/2)$.

[TABLE]

∎

2.2 Generating oracle intervals with Monte Carlo sampling

Construction of the intervals in (7) requires knowledge of the quantiles of the cdfs $H_{[i]}$. Working directly with $H_{[i]}$ may be difficult. If, instead, we can easily sample from the joint distribution $G$ of $\hat{\boldsymbol{\theta}}=(\hat{\theta}_{1},\ldots,\hat{\theta}_{p})^{\top}$, then the quantiles of $H_{[i]}$ can be computed via Monte-Carlo. We now describe this oracle Monte-Carlo procedure which is detailed in algorithm 1 and illustrated in Figure 2.

First, we draw $K$ independent $p$-vectors $\boldsymbol{\vartheta}_{1}\dots\boldsymbol{\vartheta}_{K}$ from $G$. Let $\vartheta_{k,s_{k}(i)}$ be the $i$th rankedelement of $\boldsymbol{\vartheta}_{k}$ and $s_{k}$ be the ranking permutation function for $\boldsymbol{\vartheta}_{k}$. Define the observed bias in sample $k$ at rank $i$ as

[TABLE]

In the left panels of Figures 2(a) and 2(b), the true distribution $G$ is shown in the background and the wider distributions of three sets of Monte-Carlo replicates are overlayed. In Figure 2(a), solid lines mark $\vartheta_{k,(1)}$ for $k=1,2,3$ and dashed lines mark the corresponding true parameter value $\theta_{s_{k}(1)}$. The difference between these lines is the bias $\tilde{\delta}^{k,[1]}$. In Figure 2(b), the same lines are shown for the 100th rank.

Next, for each rank, we use the sample quantiles of $\left\{\tilde{\delta}_{1,[i]}\dots\tilde{\delta}_{K,[i]}\right\}$ to estimate the quantiles of $H_{[i]}$. We denote these sample quantiles as $\tilde{H}^{-1}_{[i]}(\cdot)$. In the middle panels of Figures 2(a) and 2(b), we show the distributions of $\left\{\tilde{\delta}_{1,[1]}\dots\tilde{\delta}_{K,[1]}\right\}$ and $\left\{\tilde{\delta}_{1,[100]}\dots\tilde{\delta}_{K,[100]}\right\}$ and the 0.05 and 0.95 quantiles are marked by vertical dashed lines. Note that the distribution of bias for the largest statistic is right shifted compared to the distribution of bias for the 100th statistic.

Substituting these estimates into the interval in (7) gives the oracle confidence interval

[TABLE]

This pivot is illustrated in the right panels of Figures 2(a) and 2(b). For both of the ranks shown, the oracle Monte-Carlo intervals are shorter than the marginal interval and contain the true parameter value.

These intervals are called oracle confidence intervals because they use knowledge of $G$ and $\theta_{1},\dots,\theta_{p}$. We know that $\hat{H}^{-1}_{[i]}(x)\rightarrow H^{-1}_{[i]}(x)$ as the number of Monte-Carlo samples increases, so using (9) we can achieve the correct $1-\alpha$ confidence level (within any $\epsilon$ tolerance). This can be seen in Figure 1 where the oracle intervals have very close to the target 90% rank conditional coverage at all ranks and are shorter than other methods. The following sections describe bootstrapping methods for estimating $H_{[i]}^{-1}$ when the distribution of $\hat{\boldsymbol{\theta}}$ is unknown.

The intervals, $CI_{s(i)}^{\text{oracle}}$, in (9) do not necessarily contain $\hat{\theta}_{s(i)}$. This is particularly true if the point-estimates, $\hat{\theta}_{s(i)}$, have not been adjusted for multiplicity/selection bias — for example, if each $\hat{\theta}_{s(i)}$ is a maximum likelihood estimate. In Figure 2(a), both the $1-\alpha/2$ and $\alpha/2$ quantiles of the observed-bias are positive so the confidence interval lies completely below $\hat{\theta}_{s(1)}$. A more natural point estimate to pair with this confidence interval is the de-biased estimate proposed by Simon and Simon (2013), $\hat{\theta}_{s(i)}-\frac{1}{K}\sum_{k=1}^{K}\tilde{\delta}_{k,[i]}$, which will generally lie within the confidence interval.

The following bootstrap methods for confidence interval construction are extensions of point estimation methods proposed by Simon and Simon (2013) and Tan et al. (2014). Those former are based on estimating the mean of $\delta_{[i]}$ while we estimate quantiles. We consider two bootstrapping strategies — a parametric bootstrap, useful when the distribution and covariance of the parameter estimates are known or can be approximated well, and a non-parametric bootstrap method which can be applied to any set of estimates based on iid samples, but is more computationally costly. Both of these are general strategies, where the specifics of the algorithm may vary depending on the specific application.

2.3 Parametric Bootstrap

The parametric bootstrap parallels the Monte Carlo algorithm in Algorithm 1 replacing $G$ and $\theta_{1},\dots,\theta_{p}$ with estimates based on the data.

We assume that $G$ is a member of a parametric family of distributions and estimate its parameters. Though, in principal we could use any family of distributions, in most cases we will wish to assume that $\hat{\theta}_{i}\sim N(\theta_{i},\sigma_{i}^{2})$ where $\sigma_{i}^{2}$ is either known or can be estimated and $\hat{\theta}_{i}$ are independent. This type of parametric bootstrap is best suited for scenarios in which the estimator has an asymptotically normal distribution with known variance such as linear regression. $\hat{G}$ can be chosen by replacing $\theta_{i}$ with an estimate such as $\hat{\theta}_{i}$ itself or de-biased estimates of Simon and Simon (2013) or Tan et al. (2014). The latter choices will involve two stages of bootstrapping — the first to generate a de-biased mean estimate and the second to generate confidence intervals.

The quantiles estimated through Monte Carlo simulation in (7) are replaced by bootstrapped quantiles $\hat{H}^{-1}_{[i]}(x)$ obtained by sampling $p$-vectors from $\hat{G}$ rather than from $G$. This gives the bootstrap intervals

[TABLE]

This procedure is described in algorithm 2 for the case when $G=N(\boldsymbol{\theta},\mathbf{I})$ where $\boldsymbol{\theta}=(\theta_{1},\dots,\theta_{p})$ and $\mathbf{I}$ is the $p\times p$ identity matrix. Many variations on this procedure are possible. For example, it may be easier to specify a distribution for a transformation of $\theta_{i}$.When the family assumptions on $G$ are correct or asymptotically correct, a better RCC is obtained when the estimates used to calculate the bias are closer to the true parameter values. Algorithm 4 of the Appendix shows the additional steps necessary when using a debiased mean estimate. If the ranking scheme is based on the absolute value of the parameter estimates (for example, using the magnitude of a $t$-statistic) it is necessary to reflect the the pivot across zero for negative parameter estimates. Algorithm 5 of the Appendix gives the parametric bootstrap procedure for absolute value based rankings. Algorithm 2 and the two variations described in the Apppendix are implemented in the par_bs_ci function of the R package rcc.

2.4 Non-Parametric Bootstrap

The parametric bootstrap can be applied when $\hat{\boldsymbol{\theta}}$ has a distribution that is well approximated by a member of a parametric family. It is particularly convenient for statistics which are asymptotically normal and either independent or with a covariance that can be estimated well. Many high dimensional problems possess complex dependence structures which are not easy to estimate. Furthermore, sometimes our estimators do not have a known distribution. In these cases, the parametric bootstrap, the selection adjusted FCR controlling methods, and EB methods that assume conditional independence between estimates are unsuitable. In general, it is not possible to estimate a general $G$ without making any structural assumptions about the true parameters. However, it is possible to generate bootstrap samples non-parametrically if individual data are available.

The non-parametric bootstrap is based on sampling from the data used to compute $\hat{\boldsymbol{\theta}}$ and computing new estimates using the re-sampled data. This is implicitly sampling from a distribution $\hat{G}$ without requiring an analytical form.We assume the data consist of $n$ independent data vectors $\mathbf{y}_{1},\ldots,\mathbf{y}_{n}$. These may be vectors of genotypes, biometric, or image data for $n$ individuals. They may be a mix of data types and include covariates. We assume only that there is a procedure which takes $\mathbf{y}_{1},\ldots,\mathbf{y}_{n}$ as inputs and generates estimates $\hat{\boldsymbol{\theta}}$ and statistics indicating the significance of each estimate.

A bootstrap $p$-vector can be generated by sampling $n$ data vectors from $\mathbf{y}_{1},\ldots,\mathbf{y}_{n}$ without replacement and applying the original estimation procedure. From this point confidence intervals may be constructed identically to the parametric case.

More formally, if $\mathbf{y}_{1},\ldots,\mathbf{y}_{n}$ are iid draws from $\Pi$, and $\hat{\boldsymbol{\theta}}\equiv\hat{\boldsymbol{\theta}}\left(\mathbf{y}_{1},\ldots,\mathbf{y}_{n}\right)$ is a function of those observations, then $G\equiv G\left(\Pi\right)$ is directly a function of $\Pi$. To estimate $G$, we can use the estimate induced by the empirical distribution of the $\mathbf{y}_{i}$: $\hat{G}_{emp}\equiv G\left(\Pi_{n}\right)$. From here we can estimate the quantiles in (7) by $\hat{H}^{-1}_{emp[i]}(x)$ obtained by sampling repeatedly from $\hat{G}_{emp}$. This leads us to the non-parametric bootstrap intervals:

[TABLE]

The specifics of this procedure are shown in Algorithm 3 and implemented in the nonpar_bs_ci function of the R package rcc.Non-parametric bootstrapping can potentially be very time consuming. If the original analysis was computationally expensive it may be infeasible to repeat it many times to obtain confidence intervals.

3 Simulations

3.1 Linear Regression with Correlated Features

In this set of simulations, we explore how correlation among parameter estimates effects the rank conditional coverage rates of different methods of confidence interval construction. Code replicating these results as well can be found at https://github.com/jean997/rccSims/walkthroughs/linreg_sims.pdf.

We consider a common analysis procedure used in genetic and genomic studies. In these studies, researchers measure far more features (such as gene expression levels) than there are samples. They therefore focus on estimating the marginal association between each feature and an outcome. We consider a setting wherein the features occur in correlated blocks leading to correlated parameter estimates.

In each simulation, we simulate 1000 normally distributed features features for 100 samples. Let $x_{i,j}$ denote the value of the $j$th feature for the $i$th individual and $\mathbf{x}_{i}=(x_{i,1},\dots,x_{i,1000})^{T}$. The features are simulated as

[TABLE]

where the covariance matrix, $\Sigma$, is block diagonal with 100 $10\times 10$ blocks. The diagonal elements of each block are equal to 1 and the off diagonal elements are equal to $\rho$. The outcome for individual $i$ is simulated as

[TABLE]

where elements of $\boldsymbol{\beta}$ (the vector conditional effect sizes) are equal 0 at all but 100 elements. In each block the effect size for the 5th feature is drawn from a $N(0,1)$ distribution while the effects for the other features are 0. These effects are fixed over all simulations.

In this analysis we estimate the marginal rather than conditional effect sizes, $\beta^{(marg)}~{}=~{}\Sigma\boldsymbol{\beta}$. We estimate $\beta_{j}^{(marg)}$ through a univariate linear regression of $\mathbf{x}_{j}=(x_{1,j},\dots,x_{n,j})^{T}$ on $\mathbf{y}$. This is a standard analysis strategy for many genomic studies such as genome wide association studies and gene expression studies.

We consider four levels of correlation between the features by setting $\rho$ equal to $0,0.3,0.8$ and $-0.1$. Rank conditional coverage and interval widths averaged over 400 simulations for each scenario are shown in Figure 3. In these results, parameter ranking is based on the absolute value of the $t$-statistic $\hat{\beta}_{j}^{(marg)}/\hat{se}(\hat{\beta}_{j})$. In Section 5.3 of the Appendix we consider ranking the parameters by first selecting the parameter with the most significant estimate in each block and then ranking only these 100 selected parameters based on the absolute value of the $t$-statistic.

We find that both the parametric and non-parametric confidence intervals perform well in all four settings and are quite similar, even though the parametric bootstrap assumes independence between the estimates. None of the other methods provides an RCC close to the nominal level except for ashr in the highest correlation scenario. The ashr method does poorly in the other scenarios because the marginal effects are not sparse. Since ashr attempts to shrink parameters to zero, the credible intervals will often be to close to zero and too small in settings when the true parameter values are not sparse.

3.2 Treatment Effects in Nested Subgroups

In Section 3.1, we found that the parametric bootstrap performed well even when the assumption of independence between estimates was violated. Here we provide an example of how the parametric bootstrap can fail when estimates are very highly correlated. Code replicating these results can be found at https://github.com/jean997/rccSims/walkthroughs/biomarker_sims.pdf.

This example is motivated by the use of biomarkers in clinical trials. Suppose we have conducted a clinical trial in which participants are randomized into two groups. For participant $i$, we record the treatment group, $trt_{i}$, an outcome $y_{i}$ and the value of a biomarker, $w_{i}$. We expect that the treatment will have a greater effect in individuals with higher values of the biomarker but don’t know the exact relationship between the biomarker and the treatment effect. In an exploratory analysis, we define a series of cut-points $c_{1},\dots,c_{p}$. For each cut-point we estimate the difference in treatment effects for participants with biomarker measurements above and below the cut-point:

[TABLE]

We estimate $\beta_{j}$ as the OLS estimate fitting the regression

[TABLE]

where $1_{w_{i}>c_{j}}$ is an indicator that $w_{j}>c_{j}$. We then rank these estimates by the absolute value of their $t$-statistics in order to select a cut-point that gives the most significant difference in treatment effect between groups. This cut-point might be used to design future clinical trials.

In each simulation, we generate data for 200 study participants, 100 randomized to the treatment arm and 100 randomized to the control arm. We simulate the value of the biomarker as uniformly distributed between 0 and 1. The true relationship between the biomarker, the treatment, and the outcome given by

[TABLE]

The observed outcome for individual $i$ ($i\in 1,\dots,200$) is $y_{i}=E[y_{i}|w_{i},trt_{i}]+\epsilon_{i}$ where $\epsilon_{i}\sim N(0,0.25)$.

We chose 100 cut-points evenly spaced between 0.1 and 0.9. Rank conditional coverage and interval width averaged over 400 simulations are shown in Figure 4. In this scenario, parameter estimates are very highly correlated. This results in very poor performance for the parametric bootstrap which assumes independence between estimates. Interestingly, the standard marginal intervals do well despite making the same assumption. The non-parametric bootstrap also controls the RCC though it has slight under-coverage for the least significant parameters. Unlike the marginal intervals, the non-parametric bootstrap controls the RCC by modeling the correlation structure between parameter estimates and also performs well in the simulations in Section 3.1 making it a more reliable choice.

4 Discussion

Interval estimation when the number of parameters is large is a challenging problem often ignored in large scale studies. Out of caution, these studies are often limited to hypothesis testing but this limitation is unnecessary in many cases. We have shown that the full set of parameter estimates contains information and can be used to correct bias and generate useful confidence intervals. We have also introduced a more granular, informative concept of coverage which can be applied to confidence intervals constructed for numerous parameters.

Rank conditional coverage is an important criterion to consider in evaluating confidence intervals for large parameter sets. As a finer grained criterion, it reveals problems that are ignored by the FCR. In many cases, using an FCR controlling procedure (even after selecting top parameter-estimates) results in very low coverage probabilities for the very largest parameters.

In our simulations we found that rank conditional coverage is a more difficult criterion to control than the false coverage statement rate of Benjamini and Yekutieli (2005). The two proposed bootstrapping methods almost always outperformed other methods and produced smaller intervals than all methods except the ashr method of Stephens (2016).

SUPPLEMENTARY MATERIAL

Appendix:

Contains a proof that controlling the RCC guarantees control of the FCR, two variations of Algorithm 2, and additional simulation results referenced in the text.

R-package rcc:

R-package implementing Algorithm 2, Appendix Algorithms 4, and 5, and Algorithm 3. (GNU zipped tar file)

R-package rcc-sims:

R-package replicating the simulations shown in Section 1.5, and 3. (GNU zipped tar file)

5 Appendix

5.1 Proof that controlling RCC guarantees FCR control

Theorem 5.1**.**

If a set of confidence intervals $CI_{1}\dots CI_{p}$ controls RCC for a particular ranking scheme at level $\alpha$, then these intervals control FCR for any selection procedure which selects the top $r$ using the same ranking scheme.

Proof.

Let $S(\mathbf{t})$ be a selection procedure which chooses the $r$ most significant parameters based on the vector of statistics $\mathbf{t}=(t_{1}\dots t_{p})$. The number of selected parameters may be data dependent. Using the identity from Benjamini and Yekutieli (2005)

[TABLE]

∎

5.2 Variations of the Bootstrapping Procedure in Algorithm 2 of the Main Text

Here we present two variations of the Parametric bootstrap presented in Algorithm 2 of the main text. In Algorithm 4, we show the additional steps necessary for bootstrapping from de-biased parameter estimates. In Algorithm 4 we show a variation of the algorithm for absolute value based ranking schemes. We define an absolute value based ranking scheme as any ranking scheme that is invariant to changes in the sign of any or all parameter estimates. For example, ranking estimates based on the absolute value of the $t$-statistic $\hat{\theta}_{i}/\hat{se}(\hat{\theta})_{i}$ is an absolute value based ranking scheme.

5.3 Simulations Results Using a Block-Based Ranking Scheme

Here we consider the Simulations described in Section 3.1 of the main text, but instead of simply ranking by test statistic, we use a ranking scheme that incorporates the block correlation structure of the features. In Section 3.1 of the main text, we describe how 1000 features are simulated in 100 blocks each containing 10 correlated features. Consider a ranking scheme in which we first select the most significant feature in each block and then rank these selected features.

This is similar to the way that genome-wide association study results are often presented. In these studies, we often discover that many variants in a small genomic region are associated with the trait. Usually this arises because there are many variants that are correlated with a single causal variant so researches will typically describe only the top variant in the region.

Rank conditional coverage and interval widths averaged over 400 simulations for each scenario are shown in Figure 5. The methods of Reid et al. (2014) and Weinstein et al. (2013) are not shown because these are based on selection rules based on the test statistic alone. These results show a similar pattern as those presented in Section 3.1 of the main text — the parametric and non-parametric boostrap confidence intervals have close to the nominal RCC and are shorter than the naive intervals. Both ash and the naive intervals have much lower RCC than the nominal level. In the scenario with the largest correlation between features ($\rho=0.8$), we find undercoverage of the least significant parameters by the parametric bootstrap. This is caused by correlation between the parameter estimates which are assumed to be independent in this implementation of the parametric bootstrap.