Inference with approximate local false discovery rates

The two-group model introduced in Efron et al., (2001), is a popular mixture model-based procedure for large scale inference. Within this framework it is common to control false discovery by controlling the False Discovery Rate (FDR, Benjamini and Hochberg, 1995), or its closely related variant the marginal FDR (mFDR, Genovese and Wasserman, 2002).

Sun and Cai, (2007) developed an adaptive multiple testing rule for mFDR control in the two-group model and showed that such a procedure is optimal in the sense that it minimizes the false non-discovery rate while controlling the mFDR under independence of tests. This rule is based on thresholding the marginal local FDR (locFDR), which is the conditional probability of the hypothesis being null, given only the observed test statistic for that hypothesis. If we denote the true state for test $i$ by $h_{i}\in\{0,1\}$ and the observed test statistics by $Z_{i}$, then marginal $locFDR={\mathbb{P}}(h_{i}=0|Z_{i})$. In a series of papers (Efron,, 2007, 2010), Efron has investigated the effect of correlations on the outcome of multiple testing procedures, being concerned that approaches that ignore dependence may not correctly control the FDR, so it is necessary to accommodate these correlations in multiple testing methodology.

Putting power at the forefront, Xie et al., (2011) introduced the general two-group model with arbitrary dependence structure among z-scores. They showed that the optimal policy for mFDR control is thresholding the “oracle” locFDR statistics, $\mathbb{P}(h_{i}=0|Z_{1},\ldots,Z_{K})=\mathbb{P}(h_{i}=0|\boldsymbol{Z})$, with a fixed threshold. The oracle locFDR statistics explicitly use the correlation among the z-scores. Recently Heller and Rosset, (2021) showed that the oracle locFDRs are the optimal test statistics for controlling FDR as well, under general dependence. They showed that for some special dependence structures like block dependence or equicorrelated setting, it may be practical to calculate these oracle locFDRs for dependent data. However beyond these special cases there is no computationally feasible method for calculating the full locFDR under dependence. A standard approach is to approximate the oracle locFDR statistics by sub optimal marginal locFDR statistics which are much easier to compute. By ignoring the correlation among tests statistics, this approach leads to substantial loss of power, as demonstrated empirically in our analyses below.

As a real life example, in GWAS we are interested in finding significant SNPs associated with a phenotype among millions of SNPs in the whole genome. It is well understood that nearby SNPs on the genome are highly correlated due to linkage disequilibrium (LD), while SNPs which are far apart are very weakly correlated. Marginal analysis is a standard procedure for finding associations in GWAS, but it does not identify the SNPs that are associated with the phenotype given all other SNPs. For a joint analysis, one approach is to approximate the oracle locFDRs using a complete Bayesian treatment, which relies on methods like Markov Chain Monte Carlo (MCMC) for inference. For example, Zhu and Stephens, (2017) developed a Bayesian approach for finding association using GWAS summary statistics which takes into account the correlations among the SNPs. Their approach relies on publicly available estimates of effect sizes for $K$ SNPs $\hat{\beta}_{j}$ for $1\leq j\leq K$, calculated from marginal regressions without taking correlations into account, and builds a complete model taking into account the correlations between these estimates, using a known correlation (LD) matrix. They then implement a full Bayesian approach using MCMC to ultimately infer the posterior null probabilities $\mathbb{P}(\beta_{j}=0|\hat{\boldsymbol{\beta}})$ where $\hat{\boldsymbol{\beta}}=(\hat{\beta_{1}},\ldots,\hat{\beta_{K}})$. In practice, their procedure is very time consuming, as each MCMC iteration takes substantial computation, and the entire procedure has to be repeated several times to confirm convergence of the MCMC solution.

In this paper we propose a novel approach for mFDR control under dependence in the two-group model, which can be viewed as a compromise between the simplistic approach of relying on marginal locFDRs and the computationally prohibitive approach of calculating the full oracle locFDRs and obtaining optimal power. We propose to define the N-neighborhood for each test, for some small $N\in\mathbb{N},$ and only consider multiple testing procedures that base the decision for test $i$ only on its $N$ neighborhood. We enote this class of procedures (or decision policies) by $C_{N}.$ For example, in our motivating application to GWAS, LD structure implies physically close SNPs are correlated, and the range of correlation is roughly similar throughout the genome, so the neighborhood for each SNP is defined by the group of $2\times N$ SNPs surrounding it.

The rest of the paper is organized as follows.

In § 2 we define all the necessary notations and explain our goal in this paper. In § 3 we prove in Theorem 3.1 that within the reduced class $C_{N}$, the optimal test that maximizes power while controlling mFDR in the two-group model, relies on thresholding the posterior probabilities $T_{i,N}=\mathbb{P}(h_{i}=0|Z_{i-N},\dots,Z_{i},\dots,Z_{i+N}),$ which we term $locFDR_{N}.$We further show in Corollary 3.1.1 that the power (expected number of true discoveries) is an increasing function of $N$ for the resulting rules, making explicit the trade off between power and computation that the neighborhood approximation provides. In § 4 we show that the computation of the $locFDR_{N}$ statistics is linear in the number of hypothesis. Furthermore we show how we can estimate the $locFDR_{N}$ statistics assuming a reasonable model for GWAS. Finally we summarize our complete data driven algorithm. In § 5 we perform extensive simulations with different covariance structures and we observe that the $locFDR_{N}$ statistics achieves more power compared to other practical methods. In § 6 we provide a complete analysis pipeline for GWAS using summary statistics. Required input data are marginal SNP effect sizes ($\hat{\beta}_{j}$), standard error of the effect (s.e.($\hat{\beta}_{j}$) $\approx\hat{\sigma}_{j}$) and estimated matrix $\hat{R}$ of LD (correlations) among Z-scores ($Z_{j}={\hat{\beta}_{j}}/{\hat{\sigma}_{j}}$) of SNPs.

We assume the data are generated from the generalized two group model. The hypothesis states vector, $\boldsymbol{h}=(h_{1},\dots,h_{K})$, has entries sampled independently from the $Ber(\pi)$ distribution. Given the hypothesis states $\boldsymbol{h}$, the sequence of test statistic for testing the $K$ null hypothesis, $h_{i}=0$, $i=1,\ldots,K$, follows $\boldsymbol{Z}=(Z_{1},\dots,Z_{K})\sim g(\boldsymbol{z}|\boldsymbol{h})$. Based on $\boldsymbol{Z}$, the goal is to construct a multiple testing procedure, defined by the decision rule $\boldsymbol{D}(\boldsymbol{Z})=(D_{1}(\boldsymbol{Z}),\dots,D_{K}(\boldsymbol{Z})):\mathbb{R}^{K}\to\{0,1\}^{K}$, which has good power while controlling a meaningful error.

Let $T_{i}(\boldsymbol{Z})$ denote the oracle locFDR, which is the statistic for optimal decisions (Xie et al.,, 2011; Heller and Rosset,, 2021):

The expected number of the true and false positives are simple expressions of the oracle locFDRs:

A relevant error rate is the marginal false discovery rate (mFDR),

where $\quad\mathbb{E}_{\boldsymbol{Z},\boldsymbol{h}}(R(\boldsymbol{D}))=\mathbb{E}_{\boldsymbol{Z},\boldsymbol{h}}(V(\boldsymbol{D}))+TP(\boldsymbol{D})$ is the expected number of rejections.

Note that some expectations above are over the joint distribution of $(\boldsymbol{Z},\boldsymbol{h})$ and some are over $\boldsymbol{Z}$. From now on we will not denote the random variable with respect to which we take the expectation since it is easily understandable from the context. It has been shown in Xie et al., (2011) that the optimal rule over the class of all decision functions $\boldsymbol{D}:\mathbb{R}^{K}\to\{0,1\}^{K}$ for the above problem is thresholding the locFDR statistic with a fixed threshold. So the decision rule is $D_{i}(\boldsymbol{z})=\mathbb{I}(T_{i}(\boldsymbol{z})\leq c)$, where $c$ is a fixed constant depending on $\alpha$. The complexity required for naively calculating all the locFDR statistics is $O(K2^{K})$, which is typically infeasible. Hence the optimal procedure is impractical, except for some very special types of dependencies across the $z$-scores, e.g., block dependence with small block sizes (Heller and Rosset,, 2021). Due to the difficulty in finding the globally optimal solution, we suggest considering a smaller class of decision functions which depend only on $N$ neighbors, for $N\in\{0,\ldots,K\}$:

In this definition we assume that the tests are ordered and each test’s ”neighbors” are defined by the ones next to it in the order. We also assume that the neighborhood size $N$ is fixed for all tests. Both of these assumptions can be relaxed, but we make them here for simplicity of notation and presentation of the resulting algorithms.

Given a choice of $N$ and the resulting class $C_{N},$ we show in § 3 that the optimal testing procedure, which is limited to rules in $C_{N}$, relies on the localFDR after marginalization so that the test statistic is within this class:

Henceforth, we refer to $T_{i,N}$ as the $locFDR_{N}$. For $N=0$, $C_{0}$ includes only marginal rules and hence $locFDR_{0}$ is the marginal locFDR (Efron et al.,, 2001), while for $N=K$, $locFDR_{K}$ (Sun and Cai,, 2009; Xie et al.,, 2011; Heller and Rosset,, 2021) is the true locFDR. We can view increasing the size of $N$ as offering a trade-off between computation and power: for small $N$, calculating $locFDR_{N}$ and the resulting optimal rules within the limited class $C_{N}$ is computationally efficient but can suffer loss of power, and as $N$ increases computations become exponentially more complex but power improves. As we demonstrate below, in GWAS-like situations with a natural neighborhood structure and strong dependence only within neighborhood, the power increase can be rapid when $N$ is small, and we can obtain excellent power-computation tradeoff (Figure 1(a)). For other types of correlation structures we also get a nice trade-off (Figure 1(b), 1(c)) . Here our aim is not to find the optimal solution (which is not implementable), but a relaxation of the optimal solution which improves power significantly compared to existing marginal test statistics. We show in § 5 that the power gain using the $locFDR_{N}$ statistic, compared to existing methods, can be substantial.

In the following proposition we state a simple but useful result.

A similar result was given in Heller and Rosset, (2021) for the marginal locFDR statistics $\{T_{1,0},\dots,T_{K,0}\}$, but we expect more power by using the procedure based on $T_{i,N}$, $N>0$.

The optimal multiple testing procedure with mFDR control is a solution to the problem:

Now we concentrate on the class $\mathcal{C}_{N}$ and find an optimal solution within this restricted class, so the optimization problem is

Assume we have a standardized phenotype $(\boldsymbol{y}^{n\times 1})$ and standardized Genotype data $(X^{n\times p})$ available for $n$ individuals and $p$ different SNPs. We want to find the SNPs that are associated with the phenotype, while controlling the mFDR at level 0.05, using the $locFDR_{N}$ statistics for $N>0$. In our GWAS example we have a sufficient number of individuals so $n>p$. In order to apply our method we need z-scores and correlation among z-scores as input data. For finding the z-scores, we use the linear model,

and we know the least squares estimate of $\boldsymbol{\beta}$ is $\hat{\boldsymbol{\beta}}=(X^{T}X)^{-1}X^{T}\boldsymbol{y}$ with $\hat{\boldsymbol{\beta}}\sim\mathcal{N}(\boldsymbol{\beta},\sigma^{2}(X^{T}X)^{-1})$. Let $S=diag((X^{T}X)^{-1})$, and let the estimated variance of $\epsilon_{i}$ be $\hat{\sigma^{2}}=\frac{1}{n-p}{(\boldsymbol{y}-X\hat{\boldsymbol{\beta}})^{T}(\boldsymbol{y}-X\hat{\boldsymbol{\beta}})}$. Our z-scores vector is $\boldsymbol{z}=\frac{1}{\hat{\sigma}}S^{-1/2}\hat{\boldsymbol{\beta}}$, and the correlation among the z-scores is $\Sigma=S^{-1/2}(X^{T}X)^{-1}S^{-1/2}$. So, our input data is $(\boldsymbol{z},\Sigma)$ and our working model is the two group model (7).

We apply our method to some selected SNPs of chromosome 20, as described next. We have the genotype data of around 500K individuals and approximately 20K SNPs for each of them from the UK biobank (Sudlow et al.,, 2015). The individuals with missing height (a small number) were removed from our analysis. In the raw data for each SNP we also have many missing individuals. Moreover, SNPs that are adjacent to each other tend to be highly correlated. Therefore, we choose about 3.5K SNPs from chromosome 20 such that between any two of them there are at least 5 SNPs and none of these 3.5K SNPs has more than 10% missing percentage. We impute the number of major allele count of each missing individuals for a particular SNP by the mean of all the available individuals’ major allele count of that SNP. Next, we compute the $z$ scores vector and $\Sigma$, and use them in the “Est-S&C” data driven procedure. The estimated mixture components are

In Table 3 we show simulation results with parameters estimated from the data $(\boldsymbol{y},X)$ (essentially a parametric bootstrap experiment) after repeating the process $B$ times and aggregating the results. First row of the table shows the result where we assume the parameters are unknown and we estimate the parameters in each run (we take $B=200$ here). Second row of the table shows the result where we assume the parameters are known and we take $B=500$ here. Our procedures control the mFDR at the nominal 0.05 level.

The number of rejected SNPs for the real data at level $\alpha=0.05$ by our method and the competitors with 3.5K SNPs are given in Table 4. Though “S$\&$C” and “$\boldsymbol{T}_{0}$-rule” use the same statistics for inference, due to different estimation methods of these statistics and different way of determining cutoffs, they reject different number of SNPs. The “$\boldsymbol{T}_{2}$-rule” rejects $30$ more SNPs than the “S$\&$C” method. Manhattan plots of all the methods are given in Figure 2. In all the subplots in Figure 2, we indicate by a red rectangle a specific region of the Genome where $T_{i,N}$ statistics with $N>0$ finds more associations than “BH”, “ABH” and “$\boldsymbol{T}_{0}$-rule”. Note that in this specific region there are two associations ($1292$-th and $1405$-th SNP) which are found by all the methods. However between these two SNPs, there are two other close SNPs ($1350$-th and $1363$-th SNP) which are only identified by $T_{i,N}$ statistics with $N>0$. Interestingly, we observe from the correlation matrix of the $z$-scores that these two newly identified SNPs are not highly correlated.

In this paper we have studied practical multiple testing procedures in the two-group model under general dependence. We provide motivation and need for finding an optimal solution in a restricted decision space $\mathcal{C}_{N}$ and indeed show empirically that we do not lose much by concentrating on the smaller class $\mathcal{C}_{N}$. We provide practical approaches to estimate the parameters of a two group model under short range dependence structures which encompasses typical real applications. Using this approach we give a practical application in GWAS by finding more associations than existing practical method using summary statistics. Our developed method works in GWAS without resorting to extensive modelling assumptions of complete bayesian framework as done in Zhu and Stephens, (2017).

Though concentration on the class $\mathcal{C}_{N}$ of decision functions makes sense, it is natural to ask if there is better way of defining such a class which yields more practical and powerful solution. Our solution does not use the full correlation matrix - it only uses the correlations within $N$ main diagonals of the correlation matrix. In other words, the solution ignores all the correlation outside the $N$-band. One possible approach can be to define $i$-th locFDR statistic by conditioning on top $N$ highest correlated $z$-scores and possibly we can choose this $N$ depending on $i$ determined by strength of correlation of the test statistic with the other test statistics. It will be interesting to explore how we can modify our method to accommodate the whole correlation matrix so we do not lose much “information”.

Proposition 3.1 enables us to see how we can devise a method for mFDR control using the statistics $T_{i,N}$. If the entire joint distribution of the $z$-score vector is known, we can always define a policy of the form $\boldsymbol{D}^{t}=(D_{1}^{t}(\boldsymbol{Z}_{1,N}),\ldots,D_{K}^{t}(\boldsymbol{Z}_{K,N}))$, $D_{i}^{t}(\boldsymbol{Z}_{i,N})=\mathbf{I}(T_{i,N}\leq t)$ and determine

with $Err\in\{FDR,pFDR,FDX\}$. This procedure $\boldsymbol{D}^{t_{\alpha,N}}$ will control the chosen error by construction. However we may not be able to prove similar optimality result like Theorem 3.1 for the policy $\boldsymbol{D}^{t_{\alpha,N}}$ with these error rates. To find optimal solutions, we need to solve the generic problem with $Err\in\{FDR,pFDR,FDX\}$,

The difficulty with finding such a solution is that other error measures such as FDX, FDR, pFDR do not decompose in the same way as mFDR as shown in equation (A). It is a challenge to find similar solution for these error measures.

It is an interesting direction of future work to develop $T_{i,N}$ based method for $p>>n$ problems with properly chosen $z-$scores which takes into account the joint effect of covariates to some extent at least.

Sun and Cai, (2009) developed a method assuming hypothesis states follow the hidden markov model (HMM) structure and the z-scores given the hypothesis states are independent. So their method does not use correlation among the $z$-scores. Interestingly, proposition 3.1 holds for any joint distribution of $(\boldsymbol{Z},\boldsymbol{h})$ on $\mathbb{R}^{K}\times\{0,1\}^{K}$. Although it is normally assumed that $h_{i}\overset{i.i.d.}{\sim}Ber(\pi)$, we can also assume that $\boldsymbol{h}$ follows a HMM as done in Sun and Cai, (2009) and we can draw inference based on the statistics $T_{i,N}$ in these situations also. An interesting direction to extend our approach is thus to also allow an HMM structure of hypothesis states in addition to the correlation among $z$-scores. We expect that it will be possible to derive a practical solution for this setting.

In particular, our estimation method is such that we can easily extend to this HMM based model as well, thereby producing a multiple testing procedure which we can directly apply to real data sets.

Our approach for parameter estimation in the two-group model relies on being able to extract a subset of test statistics which are roughly independent. We believe this is a realistic assumption in most real-life applications, in particular in GWAS. However it remains a challenge to design and test approaches that also apply in cases of long-range strong dependence. Note also that our estimation method do not take into account the known joint dependence among the statistics. Though this estimation approach is robust, taking joint dependence into account can lead to efficient estimate of the model parameters and in turn a boost in power. It may be interesting to explore how we can accommodate the correlations among the test statistics into estimation.

One major observation from our simulation results is that taking even part of the (short range) dependence into account improves power. In fact, we observe large power gains even when there is very limited and not very strong dependence. Hence our recommendation is to use $locFDR_{N}$ statistics based rules, with $N>0$, whenever there is known dependence, or the dependence within neighborhoods can be well approximated.

SUPPLEMENTARY MATERIAL

Github Link of all relevant codes:

https://github.com/rajeslab/locFDR_N

Height data set:

We use Genotype ($X^{n\times p}$) and Height ($\boldsymbol{y}^{n\times 1}$) data of around 500K individuals provided by uk biobank. Unfortunately we can not make the data available publicly. For reproducibility, readers can use their own data in the notebook available in the above github link.

Acknowledgments

This study was supported by Israeli Science Foundation grant 2180/20. UK Biobank research has been conducted using the UK Biobank Resource under Application Number 56885.