Selection-adjusted inference: an application to confidence intervals for cis-eQTL effect sizes
Snigdha Panigrahi Junjie Zhu Chiara Sabatti
stat.AP



The goal of eQTL studies is to identify the genetic variants that influence the expression levels of the genes in an organism. High throughput technology has made such studies possible: in a given tissue sample, it enables us to quantify the expression levels of approximately 20,000 genes and to record the  alleles present at millions of genetic polymorphisms. While obtaining this data is relatively cheap once a specimen is at hand, obtaining human tissue remains a costly endeavor:  eQTL studies continue to be based on relatively small sample sizes, with this limitation particularly  serious for tissues as brain, liver, etc.—often the organs of most immediate medical relevance. 

Given the high dimensional nature of these datasets and the large number of hypotheses tested, the scientific community has adopted early on multiplicity adjustment procedures. These testing procedures  primarily control the  false discoveries rate for the identification of genetic variants with influence on the expression levels. In contrast, a problem that has not received much attention to date is that of providing estimates of the effect sizes associated with these variants, in a way that accounts for the considerable amount of selection. Yet, given the difficulty of procuring additional samples, this challenge is of practical importance.

We illustrate in this work how the recently developed  conditional inference approach can be deployed to obtain confidence intervals for the eQTL effect sizes with reliable coverage. In addition to interval estimates, we also provide a point estimate that approximately counters the effect of selection bias to calibrate the strength of discovered associations. The procedure we propose is based on a randomized hierarchical strategy. Such a strategy has a two-fold contribution: one, it reflects the selection steps typically adopted in state of the art investigations and two, it introduces the use of randomness instead of data splitting to maximize the use of available data. Analysis of the GTEx Liver dataset (v6) suggests that naively obtained confidence intervals would likely not cover the true values of effect sizes and that the number of local genetic polymorphisms influencing the  expression level of genes might be underestimated.




PDNet: Semantic Segmentation integrated with a Primal-Dual Network for Document binarization
    [
    
============================================================================================






§ INTRODUCTION




The goal of an eQTL (expression quantitative trait loci) study is to identify  the genetic variants  that regulate the expression of genes in different biological contexts and quantify their effects. Using statistical terminology,  the outcome variables (typically on the order of 20,000) are molecular measurements of the gene expression and the predictors are genotypes  for single nucleotide polymorphisms (SNP), typically on the order of  1,000,000). The variants that are discovered to regulate gene expression are referred as eVariants and careful estimation of their effect sizes is often deferred to follow-up studies.

One commonly studied sub-type of eVariants are those referred to as  cis-eQTL: the DNA variants in the neighborhood of a gene that influence its expression directly.  The majority of eQTL investigations have focused on detecting these cis-variants owing to their relative simple biological interpretation as well as the fact that restricting attention to this subset of gene and variant pairs reduces the number of tested hypotheses and leads to improved power. Still, even when concentrating on cis regulation, eQTL studies face a formidable multiplicity problem and also, a subsequent winner's curse during effect size estimation of discovered associations, with approximately 20,000 genes and an average of 7500 variants in each cis region.  


Since the first studies <cit.>, the eQTL research community has recognized the False Discovery Rate (FDR) as a relevant global error rate and adopted corresponding controlling strategies. As the density of SNP genotyping increased over time, it became apparent that naive application of FDR controlling strategies <cit.> might lead to excessive false discoveries. Many correlated variants may all be correctly identified as associated with changes in gene expression, even when there is really only one causal effect. This abundance of  discoveries, all corresponding to one true signal, artificially increases the denominator of the false discovery proportion, leaving room for some additional false findings. To address this difficulty, more recent works 
<cit.> adopt a hierarchical strategy, which controls the FDR of eGene discoveries. Firstly, for each gene one tests the null hypothesis of no association with any local variants and the p-values corresponding to these tests are passed to an FDR controlling procedure.  Secondly, for those genes for which this null is rejected (eGenes), scientists proceed to try to identify which among the cis variants have an effect (eVariants). To this end both marginal testing <cit.> and multivariate regression models <cit.> have been used to report the cis variants with significant association with eGenes. 


Once eVariants have been detected, the logical next step is to attempt to estimate their effect sizes; and given the scarcity of biological samples, it is tempting to do so using the data at hand. However, since these discoveries have been selected out of a large number of possible associations, naive estimators based on the same data used for selection would result in inaccurate estimates. Indeed, this is a situation similar to that of GWAS, where this problem of “winner's curse” has been noted before <cit.>, with some proposed solutions aiming at reducing bias (see also <cit.>). Other approaches to this general challenge include an empirical Bayesian approach in <cit.>, simultaneous inference methods in <cit.>, and  a differential privacy take on data-adaptivity in <cit.>.
One clear way out, of course, is offered by the classical concept of data-splitting, see <cit.>.  However in settings where the sample size is already small, reserving a hold-out data set for inference  is beyond affordability. This is often the case in human eQTL studies: for tissues other then the easily accessible blood and skin, the  relations between 20,000 genes and million of SNPs is typically studied with  a number of specimens in the hundreds at best. The same difficulties that lead to small sample sizes, make it unrealistic to simply defer the task of estimating effects to a new dataset.

In this work, we explore the potential for this problem of recent developments in the statistical literature: specifically,  the notion of conditional inference after selection and the power of randomization strategies. We offer a 
pipeline for the identification of eVariants and estimation of their effect sizes that mimics the hierarchical analysis, representing the state of the art in eQTL studies. 
Comparing the results of our pipeline with alternative strategies in simulations and real data analysis helps us understand the severity of the challenges of inference after selection in the context of eQTL. Our contribution supports the investigators in their choice of optimal use of the limited samples available in a single study, balancing the need to efficiently discover relations with that of making inference on the effect sizes of the discoveries. 



 §.§ Approach: a randomized conditional perspective

Our methods build upon a conditional inference perspective that was introduced in <cit.>.  The central idea is that, to counter selection bias, inference on the parameters (effect sizes, in the case of eQTL analysis) should be based on an adjusted likelihood, obtained by conditioning the data generative model upon the selection event. Conditioning has the effect of discarding the information in the data  used in selecting the eVariants, so that effect sizes are estimated on “unused data.” Other contributions that employ this perspective include <cit.>.
In addition, we capitalize on the observation that it is possible to generalize data-splitting in a manner that allows one to make a more efficient use of the information in the sample by introducing some randomization at the selection stage <cit.>. 

We provide a pipeline to construct confidence intervals and point estimates for the discovered effect sizes. 
The 100·(1-α)%  selection-adjusted confidence intervals are such that the probability with which each of them does not cover its target population parameter is at  most α, conditional on selection. We will report eVariants when the selection-adjusted confidence interval for their effect size does not cover zero. This guarantee is first described in <cit.>, where it is called selective false-coverage rate control. 
As a point estimator, we employ the maximum likelihood estimator (MLE) calculated from the conditional law, introduced in <cit.>. We call this estimator the selection-adjusted MLE: this serves as a quantification of the strengths of the discovered associations.


We conclude this introduction with examples of the challenges presented by selection and of how the conditional inference approach addresses them: this gives us the opportunity to introduce terminology as well as to discuss the concepts of randomization, target of inference and statistical model in this context. The rest of the paper deals with the more complicated hierarchical setting of eQTL research, which requires novel methodological results and  is organized as follows.  Section 2 describes the randomized selection pipeline that we employ as the eQTL identification strategy; section 3 presents our proposal for selection adjusted inference; section 4 contains the results of a simulation study and the concluding section 5 presents the analysis of data from GTEx. 




 §.§ Motivating examples


While the substantial contribution of the present work is to design a selection and inference pipeline that adapts to the hierarchical strategy typically adopted in eQTL studies, we start by considering a couple of “cartoon” examples that 
 illustrate the effect of selection bias and the overall sprit of our strategy, bypassing the complications associated to eGene selection. We focus on one gene, and imagine that the entire goal of the study is to find which of its cis variants influence its expression and with what effect sizes.  We consider two selection strategies: the first (1) consists in choosing the DNA variant that is most strongly associated with the gene (see <cit.>); the second (2) uses the LASSO to identify a set of variants. In both cases, we are interested in inferring the effect sizes of the selected variants.



Introducing some notation, let y∈^n be the response and X∈^n× p be a matrix collecting the values of predictors. Assume, without any loss of generality, that  X is both centered and scaled to have columns of norm 1. Strategy (1) identifies 
 the variant with the greatest marginal t-statistic. Let this variant correspond to the  j_0-th column of X, denoted as X_j_0 which satisfies

    |X_j_0^T y /σ | ≥ |X_j^T y /σ|  for  j∈{1,2, ⋯, p}∖ j_0;

σ being the noise-variance in the outcome variable, here assumed known.
For strategy (2), we consider a LASSO selection with a small ridge penalty ϵ >0  (for numerical stability) given by

    minimize_β1/2y-Xβ_2^2 + λβ_1 + ϵ/2β_2^2.

We choose the tuning parameter as λ= 𝔼[X^T Ψ_∞],  Ψ∼𝒩(0,σ^2 I), a theoretical value advocated in <cit.> and known to recover the true support asymptotically. 

We will be introducing randomization schemes for both strategies, which perturb the selection with the addition of some gaussian noise ω∈^p∼𝒩(0, τ^2 I_p). 
 Specifically, the randomized version of strategy (1) leads to the identification of variant j_0  when it satisfies:

    |X_j_0^T y /σ + ω_j_0| ≥ |X_j^T y /σ + ω_j|  for  j∈{1,2, ⋯, p}∖ j_0,

and the randomized version of strategy (2) solves the following  modified optimization problem   
    minimize_β1/2y-Xβ_2^2 -ω^T β + λβ_1 + ϵ/2β_2^2.


If we now consider the problem of inference following these selections, we note that we need to (a) specify a model for the data with respect to which evaluate the properties of estimators and (b) we need to formally identify the target parameters.
Of course, (a) is challenging in a context of model selection, where by definition we do not know what is the “true” model. Nevertheless it makes sense to work with what we might term the “full model,” where the mean of response variable Y  is parametrized by a ^n vector μ, without specifying a relation with X and Y∼𝒩(μ, σ^2 I): this is also the choice made in <cit.>. 

While this full model allows us to talk precisely about the distribution of Y, the target of inference (b) in both these examples  is not μ, but depends on the outcome of selection. This adaptive target can be described as the projection of μ on the space spanned by the selected variables and can be interpreted as the best linear parameter which is identified through selection. For the marginal example, the adaptive target is given by X_j_0^Tμ.  Denoting X_E as the selected sub-matrix, this target for the LASSO is given by (X_E^T X_E)^-1 X_E^T μ, the partial regression coefficients that are obtained by fitting a linear model to the selected set of variables E.

 We are now ready to explore via simulations the performance of three different inferential methods. The first (i) is  “Vanilla” inference that ignores selection and relies on the estimates for the adaptive targets that one would use if these were specified in advance: in the first example, this is simply the largest t-statistic, and for the second example, this is the least squared estimator with the corresponding intervals centered around these point estimates. Secondly (ii) we consider the  adjusted inference described in <cit.> (“Lee et al.”) which corrects for selection bias based on a screening without randomization (note that  <cit.> does not provide a selection-adjusted point estimate, but gives a recipe for confidence intervals). The third approach (iii) is in the spirit of the methods we develop in the rest of the paper: we condition on the outcome of randomized selection and provide adjusted confidence intervals and point estimates (“Proposed”).

In our simulations,  X is fixed and a response vector y∈^n is generated in each round from Y ∼𝒩(0, I_n), independently of X. The perturbation ω∈^p is generated from Ω∼𝒩(0, τ^2 I_p);  τ^2 = 0.5, independently of Y. We simulate 100 such instances and compare the three inferential procedures in terms of coverage, length of the confidence intervals and risk of the point estimates. The risk metric we use is the averaged squared error deviation of the point estimates from the respective adaptive targets, described above for the two selection strategies. The results are summarized in Figure <ref>.

As expected, the unadjusted intervals based on “Vanilla" inference fall way short of coverage. The adjusted intervals, both in the form described in <cit.> and in the present work, based on randomization,  achieve the target coverage. These two adjusted methods, however, differ in terms of interval length: while the “Proposed” intervals are 1.5 times longer than the unadjusted intervals (a price we pay for selection), they are much shorter than the exact intervals post non-randomized screenings in <cit.>. This is the advantage of randomization at the selection stage: there is more “left-over” information available at the time of inference. The risks of the “Vanilla" point-estimates post a non-randomized screening is also seen to be higher than the adjusted MLE, considered as a point estimate for effect sizes in this paper.
 
In the interest of clarity, we note that while the comparisons in Figure  <ref> are meaningful on average, the selections for the non randomized and randomized versions of strategies (1) and (2) are not the same,  so that while “Vanilla" and “Lee” are working with the same targets of inference in each realization of the simulations, “Proposed” might have slightly different targets.

In the rest of the paper, we will develop a pipeline for the analysis of eQTL data that takes advantage of randomization in the selection step to preserve more information for the inferential stage.




§ RANDOMIZED HIERARCHICAL SCREENING TO IDENTIFY INTERESTING EQTL EFFECTS




We start by introducing the notation that we will be using throughout. 
Our data will include measurements on the expression levels of G genes and V variants in n subjects.  We denote the outcome variables using Y: Y^(g)∈^n  is the  vector of gene-expression levels for gene  g∈{1,2, ⋯ ,G}. Let  V_g be the number of variants measured within a defined cis-window around gene g: we indicate the matrix of the corresponding  potential explanatory variables with X^(g)∈^n × V_g. We  assume, without  loss of generality, that X^(g) is centered and the columns are scaled to have norm 1. Denote the variance in gene-expression response Y^(g) by σ^2. 
In this section we are going to outline a randomized selection procedure that mimics the hierarchical strategy often adopted in the literature <cit.> and specify the derived “target of inference.” 




 §.§ A hierarchical randomized selection



Following the practice in eQTL studies, we want to first identify a set G of genes that appear likely to be cis-regulated, and then, for each of the genes g∈ G identify a set of potential e-variants E^(g) belonging to the cis region and appearing to have an effect on the expression of  gene g.
In order to both use the entire sample n to guide our selection and, at the same time, reserve enough information for the following inference on the effect size, we explore the potential of a selection strategy that includes randomization in both these stages. The “cartoon” examples with which we concluded the previous section illustrate the advantages of randomization. For a more comprehensive discussion, please see <cit.> and  <cit.>, which place this into the context of data re-use.

Stage-I: Randomized selection of eGenes   To discover promising genes, we test the global nulls that the expression of gene g is not influenced by any of the cis-variants V^(g): G collects the set of genes for which we reject this null controlling FDR at level q. The p-values for these global null hypotheses are calculated with randomization. 

Step 1.  Compute for each gene g ∈{1,2, ⋯ ,G}, a univariate t-test statistic based on marginal correlation of local variant X^(g)_j, j∈{1,2,⋯, V_g} with gene expression Y^(g) and added Gaussian randomization ω_j^(g). The perturbed t-statistic (z-statistic for a known σ) is given by

    T^(g)_j = X^(g)_j^T y^(g)/ σ + ω_j^(g),

where ω_j^(g) is a realization of a  Gaussian random variable Ω_j^(g)∼𝒩(0, γ^2) with γ^2 controlling the amount of perturbation. Further, perturbations Ω_j^(g) are independent for j ∈{1,2,⋯, V_g} and also, independent across all genes g ∈{1, 2, ⋯, G}.
For settings where σ is unknown, we use a marginal estimate of σ; see Remark <ref>.

Step 2.  Compute a global p-value based on Bonferroni p̃^(g) = V_g  p^(g)_(1) 
where

    p^(g)_j = 2·(1-Φ(|T^(g)_j|/√(1 +γ^2)))  for  j∈{1,2,⋯, V_g}, g∈{1,2, ⋯ ,G}.

Step 3.  We apply a Benjamini-Hochberg (BH-q) procedure at level q to the global p-values 
{p̃^(1), p̃^(2), ⋯, p̃^(G)}. Based on the rejection set of BH, we report K_0 eGenes with the number of rejections calculated as

    K_0 = max_k∈{1,2,⋯, G}{p̃_(j)≤j/Gq  for at most  k  many p-values}.
 This set of identified genes is denoted as G.
The analogous non-randomized selection procedure is based on t-statistics in (<ref>) with no added perturbation. 


Stage-II: Randomized selection of potential eVariants   The second stage identifies promising variants  E^(g) for each of the eGenes g in 𝒢. 
We use a randomized version of penalized regression <cit.>, where the ℓ_1 penalty induces sparsity (selection) and a small  ℓ_2 penalty  is used to regularize the problem. The set E^(g)  is identified by solving: 
    mimimize_β1/2y^(g) -X̃^(g)β_2^2 -ζ^(g)^T β +λβ_1+ϵ/2β_2^2,

where X̃^(g)∈^n × p_g,  p_g≤ V_g indicate a pruned set of cis variants.
The value of ϵ is small, with the only goal of ensuring the existence of a well-defined solution, while λ  is set to be the theoretical value considered in Section <ref>. In the same spirit as the randomized screening of eGenes, ζ^(g) is based on a Gaussian variable Z ∼𝒩(0, τ^2 I).  


Note that E^(g) indicates the set of active variants for the problem (<ref>) for gene g: this is the result of our selection step. However, the final results of our analysis will consist of a set  𝒬^(g)⊂ E^(g), corresponding to the   subset of variants with significant selection-adjusted p-values at a chosen level of selective Type-I control α. 

We noted that problem (<ref>)
is defined on pruned set of variants. The variants in the cis region of a gene can be highly correlated in the sample.  This makes it hard for any multi-regression analysis like the LASSO to recover the set of variables truly associated
and hence, pruning becomes essential; prior works <cit.> have recognized these challenges. We give details of an unsupervised pruning of variants using a hierarchical clustering scheme, based on the empirical correlations between variants as a distance measure in Appendix <ref>. 





In real data analysis, of course, the variance of gene expressions is not known: we use a plug-in estimate. For the marginal screening step, we consider a marginal estimate of variance to compute a t-statistic for the association between variant j and outcome. For the randomized LASSO screening, we use an estimate of σ from a refitted least squared regression post the LASSO. 




 §.§ Model and adaptive target



To proceed with inference post a hierarchical screening, we assume a full Gaussian model for gene expression, as in the motivating examples: 
    Y^(g) = μ^(g) + ϵ^(g), ϵ^(g)∼𝒩(0, σ^2 I).

The generative law in this framework parametrizes the Gaussian mean as μ^(g)∈^n, g ∈{1,2,⋯, G}. In addition, we assume that errors ϵ^(g) are independent across genes. 

Gene expression measurements and genotype data often are affected by  confounding factors, including gender, demography, platform etc. For the purposes of this paper we assume  that the data has undergone preprocessing that eliminates the effects of hidden confounding factors and measured covariates.  This allows us to assume that the noise is approximately independent across genes. Details on how these confounding factors are regressed out from the real data we analyze is provided in  Appendix <ref>.


Having described a model, we define the adaptive target of inference as the projection of the model parameters onto the space spanned by E^(g) for each gene g:

    b_E^(g) = (X̃_E^(g)^T X̃_E^(g))^-1X̃_E^(g)^T μ∈^E^(g).

Denoting the j-th component of the adaptive partial regression coefficient b_E^(g) as b_j; E^(g) = e_j^T b_E^(g), we note that  unadjusted inference for b_j; E^(g) would be based on the j-th least squared estimator

    b̂_j; E^(g) = e_j^T (X̃_E^(g)^T X̃_E^(g))^-1X̃_E^(g)^T Y^(g).

Under an independent Gaussian noise model, b̂_j; E^(g) is the MLE for b_j; E^(g) and naive confidence intervals for the same target are centered around b̂_j; E^(g) with a length of 2· z_1-α/2· (X̃_E^(g)^T X̃_E^(g))^-1_j,j;  z_1-α/2 being the standard normal quantile. In the next section we will describe alternative point and interval estimates for these quantities that take into account the selection we have described. 




§ SELECTION-ADJUSTED INFERENCE FOR EVARIANTS EFFECTS



We outline a recipe to provide selection-adjusted inference for the adaptive partial regression coefficients in (<ref>) in a coordinate-wise manner. That is, we provide the steps to get a tractable selection adjusted law: ℒ(b̂_j; E^(g)| g ∈𝒢,  Ê^(g)= E^(g)), the conditional law of the j-th component of the least squared estimator in (<ref>) conditioned upon the event that gene g was screened marginally in Stage-I and E^(g) was screened by randomized LASSO in Stage-II.

We start again by listing the notation and terminology we will use. We then give an outline of the steps that allow us to provide a practical selection adjustment for inference, describing them in detail in following sections. 



  * Data and randomization   At the core of the adjustment for selection is the calculation of a conditional likelihood. Because of the randomization we introduced, the relevant likelihood includes both the original data and the outcome of the randomization. The observed data y^(g) and randomizations (ω^(g)∈^V_g, ζ^(g)∈^p_g) are realizations of the random variables (Y^(g), Ω^(g), Z^(g)) respectively; (Ω^(g), Z^(g)) represent the random variables for Gaussian randomizations used in both stages of selection. The entire data across all genes is denoted by (y, ω, ζ) where 

    y = {y^(g) : g∈ G},  ω ={ω^(g): g ∈ G},  ζ= {ζ^(g): g ∈ G}.
 

  * Parameters in model   In providing coordinate-wise inference, we note that the target parameter in our case is 

    b_j;E^(g)= c_j^T μ;  c_j = X̃_E^(g) (X̃_E^(g)^T X̃_E^(g))^-1 e_j.
 Based on (<ref>), there are nuisance parameters that are given by
n_j^(g)= (μ - b_j;E^(g)c_j/c_j_2). In a nutshell, n_j^(g) are parameters that we are not interested in inferring about while providing inference for b_j;E^(g). 

  * Target and null statistic   A target statistic is a statistic that can be used to make inference on the target parameter b_j;E^(g). Specifically, we will consider the j-th coordinate of the least squared estimator, denoted by b̂_j;E^(g) as the target statistic for target parameter- b_j;E^(g). Note that b̂_j;E^(g) also, serves as the target statistic for naive (unadjusted) inference for b_j;E^(g). Except now, the law of this target statistic is no longer a Gaussian distribution centered around b_j;E^(g): the goal of the remaining section is to provide a tractable selection adjusted law for b̂_j;E^(g).

The null statistic according to model (<ref>) is a data vector that is orthogonal to b_j;E^(g) based on decomposition
y^(g)= b̂_j; E^(g)c_j/c_j_2+ 𝒰_j^(g). In model (<ref>), the mean of the null statistic is the nuisance parameter, defined as n_j^(g). Conditioning out the null statistics in the selection-adjusted law eliminates the nuisance parameters in the model.

  * Variance of target and null statistic   Finally, denote as σ_j; E^(g)^2 = (X̃_E^(g)^T X̃_E^(g))^-1_j,j, the variance of b̂_j;E^(g) and the variance of 𝒰_j^(g) as Σ_U^(g) under (<ref>).




 §.§ An outline of conditional inference after randomized selection in eQTL

 In order to successfully apply conditional inference,   we need to be able to efficiently compute with the selection-adjusted law. To achieve this goal in the eQTL inference problem that we have described, we take the following steps.



I. A selection-adjusted law across the genome   The randomized selection 𝒢̂(y, ω) of eGenes and of promising variants for each eGene Ê^(g)(y^(g),ζ^(g)) define a selection event 

    {(y, ω, ζ): 𝒢̂(y, ω) = 𝒢,  Ê^(g)(y^(g),ζ^(g))= E^(g) for  g ∈𝒢}.

Conditioning on this  event modifies the  law of data and randomizations as follows:

    ∏_g∈𝒢ℒ(y^(g); μ)×ℒ(ω^(g); γ)×ℒ(ζ^(g); τ)1_{g ∈𝒢, Ê^(g)= E^(g)}

where ℒ(y^(g); μ)×ℒ(ω^(g); γ)×ℒ(ζ^(g); τ) is the unadjusted law of (Y^(g), Ω^(g), Z^(g)). 

The above conditional law results from the fact that we infer about target (<ref>) for a gene g only if it has been screened in the set of eGenes and the set E^(g) of variants are chosen by a randomized version of LASSO. We point out to our readers that the selection of gene g as an eGene depends not only on the data specific to gene g, but also on data for all other genes. This leads to a complicated selection event, which simplifies to a great extent when we condition additionally on some more information beyond knowing the set of eGenes: 𝒢 and the respective active variants E^(g), chosen in the next stage. This extra information includes the signs of the t-statistics for the selected genes, the signs of Lasso coefficients, the total number of eGenes K_0 etc.
 
II. A simplified adjusted law across eGenes  

We obtain a simplified selection event {ℋ̂(y^(g), ω^(g), ζ^(g)) = ℋ} by conditioning on some additional information together with the selection of eGene g and variants E^(g) selected by LASSO. Such an event is easier to handle in the sense that it is described only in terms of data specific to gene g and thereby, allows a decoupling of the joint law under (<ref>) into an adjusted law for each gene in G. A full description of the conditioning event, a superset of the event {g ∈𝒢, Ê^(g)= E^(g)} is available in Section <ref>. 

Under (<ref>) and Gaussian randomizations, the simplified selection-adjusted density for each eGene g in terms of target and null statistic (b̂_j;E^(g), 𝒰_j^(g)), randomizations (ω^(g), ζ^(g)) and an observed selection ℋ, is given by

    exp(-(b̂_j;E^(g) -b_j;E^(g))^22σ_j;E^(g)^2)·exp(-Σ_U^-1/2(𝒰_j^(g)-n_j^(g))_2^2)
       ×exp(-ω^(g)_2^2/2γ^2)×exp(-ζ^(g)_2^2/2τ^2)× 1_{ℋ̂(y^(g), ω^(g), ζ^(g)) = ℋ}.



III. Selection-adjusted interval and point estimates  Since, we are interested in the selection-adjusted law of the target statistic b̂_j;E^(g), we marginalize over the randomizations in the joint law in (<ref>) and condition on nuisance statistics 𝒰_j^(g) in order to eliminate nuisance parameters n^(g)_j. This finally, leads to a selection-adjusted density for the target statistic b̂_j;E^(g), proportional to

    exp(-(b̂_j;E^(g) -b_j;E^(g))^2/2σ_j;E^(g)^2)𝒫_ℋ(b̂_j;E^(g)).

In this law, 𝒫_ℋ(t) = ℙ((b̂_j; E^(g), Ω^(g), Z^(g)) ∈ℋ|b̂_j;E^(g)=t) is the conditional probability of selection given b̂_j;E^(g)=t.  This is the density we will use to carry out inference on the effect sizes of eVariant.

Using the adjusted law of the target statistic to define 
    T(b̂_j;E^(g); b_j;E^(g), σ_j;E^(g)) = ∫_b̂_j;E^(g)^∞exp(-(t -b_j;E^(g))^2/2σ_j;E^(g)^2)  𝒫_ℋ(t) dt∫_-∞^∞exp(-(t -b_j;E^(g))^2/2σ_j;E^(g)^2)  𝒫_ℋ(t)dt,

an adjusted (two-sided) p-value can be computed as:

    p(b̂_j;E^(g); b_j;E^(g), σ_j;E^(g)) = 2·min (T(b̂_j;E^(g); b_j;E^(g), σ_j;E^(g)), 1- T(b̂_j;E^(g); b_j;E^(g), σ_j;E^(g))).
 
Confidence intervals for b_j;E^(g) with target coverage 100(1-α)% are constructed as

    {b ∈ : p(b̂_j;E^(g); b, σ_j;E^(g)) ≤α}.

And, a maximum likelihood estimator is obtained  solving

    [0.95]b_j;E^(g)minimize(b̂_j;E^(g) -b_j;E^(g))^2/2σ_j;E^(g)^2  + log∫exp(-(t -b_j;E^(g))_2^2/2σ_j;E^(g)^2)𝒫_ℋ(t)dt.



IV. Approximate and tractable inference   Even though, we derived a selection adjusted law in (<ref>), the term 𝒫_ℋ(b̂_j;E^(g)), the conditional selection probability in (<ref>), lacks a tractable closed form expression.
This contributes to the computational bottleneck in constructing intervals or solving a MLE problem based on the adjusted law of b̂_j;E^(g). To offer tractable inference based on (<ref>), we provide an approximation for 𝒫_ℋ(·).
 
 Details on the approximation are given in Section <ref>. Approximate inference in the form of interval and point estimates is based on plugging 𝒫̂_ℋ(·) in (<ref>) and (<ref>).


We now make precise the steps involved in obtaining the simplified conditional law (II) and the approximation of the marginal adjusted law (IV).



 §.§ A selection-adjusted law for b̂_j; E^(g)


With the unadjusted law of data vector and randomizations 

    {(b̂_j; E^(g), 𝒰_j^(g), Ω^(g), Z^(g));  g∈𝒢}
 based on a decomposition of the response vector y^(g) into target statistic, b̂_j; E^(g) and null statistics, 𝒰_j^(g) (see in list of terminologies)
as our starting point, we realize that the selection event in (<ref>) modifies this (unadjusted) density by truncating it as 

    ∏_g ∈𝒢exp(-(b̂_j;E^(g) -b_j;E^(g))^22σ_j;E^(g)^2-Σ_U^-1/2(𝒰_j^(g)-n_j)_2^2-ω^(g)_2^22γ^2-ζ^(g)_2^22τ^2)×
               × 1_{g∈𝒢̂(y, ω)=𝒢,  Ê(y^(g),ζ^(g))= E^(g)}.

By conditioning on some more additional information, we obtain s a simplified selection-adjusted law for b_j; E^(g) in each eGene g. At the core of this simplification is a decoupling of the selection event (with extra information) into an event exclusive to eGene g and an event that is dependent on data from all genes that are not screened in Stage-I. This is formally described in Theorem <ref>. The event with extra information that we condition on is given by 

    {g ∈𝒢̂,   K = K_0,    j_(1)^(g)=  j_0,  T_j_(2)^(g)= T_0, sign(T^(g)_j_0) = s_j_0,
            Ê(y^(g),ζ^(g))= E^(g), sign(β̂_E^(g)) = s_E^(g)}
 
where


  * K_0 is the number of eGenes or equivalently, the number of rejections determined by BH as

    K_0 = max_k∈{1,2,⋯, G}{p̃_(j)≤j/Gq  for at most  k  many p-values},

based on ordered (Bonferroni-adjusted) p-values
{p̃_(1), p̃_(2), ⋯, p̃_(G)}. 



  * j_(1)^(g)= j_0= min_j p^(g)_j for eGene g is the index in the set {1,2,⋯, V_g } associated with the largest t-statistic in (<ref>) or correspondingly, the smallest marginal p-value. 

  * T_j_(2)^(g)= T_0= max_j≠ j_(1)^(g) |T^(g)_j|
 is the second largest t-statistic in magnitude for gene g.

  * sign(T^(g)_j_0)= s_j_0 denotes the sign of the largest t-statistics corresponding to the index j_0 for each eGene. 

  * s_E^(g)= sign(β̂_E^(g)) denotes the signs of the active coefficients from the LASSO.




Separability result
  The selection event in (<ref>) 
decouples as 

    {ℋ̂(y^(g), ω^(g), ζ^(g))= ℋ}∩{𝒥̂(y^(g'), ω^(g'), g'∉𝒢) =𝒥} where

 {ℋ̂(y^(g), ω^(g), ζ^(g))= ℋ} is an event based on data and randomization associated exclusively with eGene g and 
{𝒥̂(y^(g'), ω^(g'), g'∉𝒢) =𝒥} is a function of data and randomization for genes that are not rejected by the BH procedure.

A proof for Theorem <ref> is included in Appendix <ref>. A consequence of the resulting separability is that the contribution from {𝒥̂(y^(g'), ω^(g'), g'∉𝒢) =𝒥} in the selection-adjusted law for b̂_j; E^(g) is only that of a constant, due to the assumed independent noise structure across genes. This results in a more simplified conditional law for (b̂_j; E^(g), 𝒰_j^(g), Ω^(g), Z^(g)) in (<ref>).

Finally, we obtain (<ref>) from the joint law in (<ref>) by marginalizing over the randomizations and conditioning on nuisance statistics 𝒰_j^(g) in order to eliminate nuisance parameters n^(g), as provided by the following Lemma.



The marginal selection-adjusted law of b̂_j; E^(g), conditioned on nuisance statistics 𝒰^(g)_j= u^(g)_j is an exponential family with a natural parameter b_j; E^(g) that is given by

    exp(-(b̂_j;E^(g) -b_j;E^(g))^2/2σ_j;E^(g)^2) ·𝒫_ℋ(𝒰^(g)_j)(b̂_j;E^(g));  with


    𝒫_ℋ(u^(g)_j)(t) = ℙ[(b̂_j; E^(g),Ω^(g), Z^(g)) ∈ℋ(u^(g)_j)  |b̂_j; E^(g)=t,  𝒰^(g)_j = u^(g)_j]

where conditional on 𝒰_j^(g), the selection-induced region in terms of the target statistic and randomizations is denoted by
{(b̂_j; E^(g), ω^(g), ζ^(g))∈ℋ(𝒰_j^(g))}.





 §.§ Approximate and tractable selection-adjusted law



A technical hurdle in providing inference about b_j;E^(g) is posed by the intractability of function 𝒫_ℋ(·) in the selection-modified exponential family of distributions. We make use of two key steps in obtaining an approximate selection-adjusted distribution based on a surrogate 𝒫̂_ℋ(·), as outlined below.



  * K.K.T. maps for hierarchical selection. The selection map associated with the identification of eGene g is given by

    ω^(g)_j_0 =  P^(g)_j  b̂_j; E^(g)  + η^(g) + q^(g)_j,

where j_0 corresponds to the index of the largest t-statistic for eGene g and P_j^(g), q_j^(g) are fixed coefficients of the linear map. The screening of eGene g imposes a simple thresholding constraint on η^(g)∈ of the form 
η^(g)>C where C is a threshold determined by the information we conditioned upon in Section <ref> and induces a selection region 𝒞_1 where

    [0.95]𝒞_1 = {η: η≥max( √(1+γ^2)·Φ^-1(1- K_0/2V_g Gq), |T_0^(g)|) }.




Similarly, the selection of set of active variants E^(g) with signs s_E^(g) by randomized LASSO can be written as

    ζ^(g) =  A^(g)_jb̂_j; E^(g)  + B^(g)_j o^(g) + c^(g)_j

where A^(g)_j, B^(g)_j , c^(g)_j are fixed matrices/ vector. Selective constraints on o^(g)∈^p_g that are variables associated with solving the LASSO objective take the form of sign and cube restrictions, which we call 𝒞_2 where

    𝒞_2 ={o: sign(o_E^(g)) = s_E^(g),  o_-E^(g)_∞≤λ}.

We provide comprehensive derivations of the maps in (<ref>) and (<ref>) associated with the hierarchical randomization strategy  and the corresponding selection regions induced by these mining strategies in (<ref>) and in (<ref>) in Appendix <ref>.


  * Change of variables formulae. We use the K.K.T. maps in (<ref>) and (<ref>) as change of variable formulae from (b̂_j; E^(g), ω^(g) _j_0, ζ^(g)) to (b̂_j; E^(g),η^(g), o^(g)). The motivation behind this is that the hierarchical selection induces complicated polyhedral geometry on (b̂_j; E^(g), ω^(g) , ζ^(g)), but the constraints on (η^(g), o^(g)) are simple threshold or sign restrictions, thereby, simplifying the computation of 𝒫_ℋ(·). 
 The selection-adjusted law of b̂_j; E^(g) induced by the change of variables formulae is provided in Lemma <ref> in Appendix <ref> as

    exp(-(b̂_j;E^(g) -b_j;E^(g))^22σ_j;E^(g)^2) ·𝒫_𝒞_1; 𝒞_2(b̂_j;E^(g)) with


    𝒫_𝒞_1; 𝒞_2(t)=ℙ[η^(g)∈𝒞_1, o^(g)∈𝒞_2  |b̂_j; E^(g)=t].

Recall that 𝒞_1 and 𝒞_2 are selection-induced regions described in terms of optimization variables η^(g) and o^(g) respectively.


  * A Chernoff-based approximation for 𝒫_𝒞_1; 𝒞_2(.). For selection regions 𝒞_1 and 𝒞_2, Theorem <ref> in Appendix <ref> obtains an expression for 𝒫_𝒞_1; 𝒞_2(·) as a function of the statistic of interest b̂_j;E^(g). 

Finally, an upper bound in Theorem <ref> gives an approximation for the logarithm of 𝒫_𝒞_1; 𝒞_2(t) and is computed as

    log C_1(t)-inf_sign(o_E) = s_E^(g){ A_E^(g), j  t + B_E^(g), j o_E + c_E^(g), j_2^2/2τ^2- log C_2(o_E, t)}

where 

    C_1(t)  = Φ̅({L+ s_j_0P^(g)_j  t + s_j_0q^(g)_j}/γ)  and


    [0.9]C_2(o_E, t ) = ∏_k ∉ E^(g){Φ({A_-E^(g),j^kt + B_-E^(g),j^ko_E+ λ}/τ) - Φ({A_-E^(g),j^k t + B_-E^(g),j^k o_E -  λ} /τ)}.

In the above, L is the lower threshold on η^(g); A_-E^(g),j^k represents the k-th row of A_-E^(g),j and B_-E^(g),j^k denotes the k-th row of B_-E^(g),j. 

In implementations of this scheme, an approximation for the logarithm of 𝒫_𝒞_1; 𝒞_2(·) is finally computed as 

    log C_1(t) -inf_o_E{ A_E^(g), j  t + B_E^(g), j o_E + c_E^(g), j_2^2/2τ^2- log C_2(o_E, t) + ℬ(o_E)}

where we use a barrier penalty function ℬ(o_E) to reflect the sign constraints on o_E in the Chernoff bound, (<ref>). Specifically, the barrier penalty that we use is ℬ(o) = ∑_k=1^E b(o_k) with

    b(o) = log(1+ 1/s o)     if sign(o) =s 
    ∞   otherwise.






































§ PERFORMANCE IN A CIS-EQTL  STUDY



In this section and the next, we examine the properties of the proposed pipeline 
by testing its performance on data collected in one eQTL study <cit.> and comparing its results with those of two other “benchmark” procedures.
We first rely on simulations of the outcome variables according to a known model based on genotype data in order to evaluate a number of performance metrics and then turn to the analysis of the real data.

The dataset we analyze is the collection of liver samples in V6p of the GTEx study <cit.>. It comprises  97 individuals, with genotypes for 7,207,740 variants (these variants are obtained as the output of the default imputation pipeline implemented in <cit.>) 
 and expression quantification for 21,819 genes. Details on data acquisition and pre-processing are provided in the Appendix Section <ref>.

In studying the performance of the proposed pipeline, we compare it with (1) a slightly modified version of the analysis strategy in the original GTEx paper, followed by a naive construction of the confidence interval for the effect sizes of the identified eVariant (GTEx+vanilla); (2) a strategy that uses “out-of-the-box” selective inference tools and employs no randomization during the selection steps (Bonf+Lasso+Lee et al.). Specifically,


  * GTEx+vanilla (G.V.)  The selection of eGenes is done controlling FDR at level 0.1 with the  BH procedure applied to  p-values for the global null obtained via the  Simes' combination rule <cit.>. We note that this differs from the original GTEx paper in two minor ways: the p-values for the global null for eGene discoveries are obtained via Simes rather than permutations in <cit.>, and the adopted FDR controlling procedure at this stage is  BH rather than Storey's procedure referred in the GTEx paper.

In the second stage, the eVariants are selected by an adaptive forward-backward selection:  a variable is added to the regression only if it p-value is lower than the largest p-value for global nulls in the set of discovered eGenes  (with details given in <cit.>). 

To estimate effect size, we fit the least squared estimator on the selected model for the point estimate and use the normal quantiles to construct intervals, thereby “naively” ignoring the selection occurred in the two screening stages.


  * Bonf+Lasso+Lee et al. (B.L.L.)  The eGene screening is conducted in a two-stage procedure similar to our proposed pipeline, with the exception that there is no randomization in the t-statistics or in the second stage selection of eVariants. The confidence intervals for effect sizes are constructed using the adjustments for the LASSO as detailed in <cit.>. We note that this approach (based on exact evaluation of the selection conditional probability) only accounts for the eVariant selection the second stage, but not for eGene selection in the first stage. Further, <cit.> does not give a selection-adjusted point estimate, and so we report the unadjusted least squared estimator on the selected variants as a point-estimate. 


Note that the proposed pipeline differs from G.V. in both stages of screening: the selection of eGenes and the screening of potential eVariants. These changes (swapping the construction of Simes global p-values with Bonferroni p-values in eGene screening and swapping the adaptive forward-backward selection with a LASSO screening of variants associated with eGenes) allow an easier explanation of the mining scheme and provides a more natural way to introduce randomization. Comparing our proposed procedure with these two alternatives, therefore, enables us to study how its performance relates to that of (1) a state-of-the-art  method for the identification of eVariants that does not take into account at all the effect of selection at the inferential stage, and to that of (2) an approach to correct for selection that relies on out-of-the box tools, without fully capturing the hierarchical identification of eVariants and not capitalizing on the possible power increases due to randomization.



 §.§ Simulation study


We start with a simulation study to explore the performance of the different procedures:
using the genotypes from the GTEx liver sample, we generate artificial gene expression values and investigate how the three approaches reconstruct the effect sizes of eVariants.



  §.§.§ Data generation

In choosing a strategy to simulate gene expression values, we followed the following principles: (a)  the selection strategy in G.V. should lead to a number of eGene and eVariant discoveries similar to those detected in the real data in <cit.>; (b) there should be some genes that are not true eGenes, so that it is sensible to consider “false discoveries"; (c) there should be eGenes regulated by multiple eVariants; (d) the model should be as simple as possible, to make the interpretation of the results straightforward. After experimenting with a few models, that satisfied criteria (b)-(d), we found the following to be the one that gave results closer to those in <cit.> and adopted it.


For each gene g for which gene expression is available in the real data, and for which V_g cis-variants have been genotyped, we generate a vector Y^(g)∈^n of  synthetic gene expression as follows


  * We randomly set the number of causal variants |𝒮^(g)|  from {0,1,2, ⋯, 9} according to the distribution in <ref>, such that approximately  a third of the  genes  contain at least one true signal. 




  * A set  𝒮^(g) of causal variants of the selected size is  drown randomly from the variants in the cis region.  To assure, however, that the desired number of “independent" signals is present, we want to make sure that 𝒮^(g) does not contain variants that are closely correlated and whose contribution to the gene expression value would be indistinguishable. 
To achieve this goal, we subject the genotyped variants to 
  hierarchical clustering with minimax linkage (an unsupervised pruning technique detailed in Appendix <ref>): we randomly select |𝒮^(g)| clusters,  and randomly assign an element in the cluster to be a causal variant. 
  
We also rely on this same clustering to identify a “pruned” set of variants that will  constitute the input of the penalized regression, which needs to work on non-collinear variables. From each of these clusters, we choose a representative, correlated by at least ρ_0=0.5 with all members of that cluster. We remark that the cluster representative and the possible causal variant residing in the cluster do not necessarily coincide.  
  
  
  * The  expression values Y^(g) are generated according to the following model

    Y^(g) = ∑_k∈𝒮^(g)X^(g)_kβ_k + ϵ^(g), ϵ^(g)∼𝒩(0, I),

where the noise vector ϵ^(g)∈ ^n is independently sampled for each gene, each X^(g)_k is standardized with mean zero and unit variance, and the true effect sizes are  
β_k =3. 

  



  §.§.§ Evaluation metrics

Because of the computational costs associated to the analysis methods, we do not repeat the data generating process multiple times for a given gene, but rather we interpret our results by aggregating across genes with the same number of causal variants. 


The focus of this work is to provide methods that allow correct inference for 
the adaptive target defined in (<ref>), dependent on the selection of variants associated with each eGene.
To evaluate the performance of the three approaches with this respects, we rely on three quantities: the coverage of the confidence intervals, their lengths, and the average empirical risk of the point estimate computed with respect to a quadratic loss function. Note that the parameters to be estimated in the loss metric for risk evaluation are the adaptive targets. 


At the same time, given that the three procedures differ in their selection steps, to  interpret appropriately the results, it is useful to also compare them at the level of selection. In keeping with the hierarchical structure of the selection, there are two levels at which it makes sense to talk about FDR and power: the eGene level and the eVariant level, for selected eGenes. FDR and power are easily defined for eGenes. To  explore the performance at the level of eVariants, we are going to focus on the eVariants for selected eGenes (therefore, eVariants for erroneously missed eGenes are not going to be considered in our power evaluation). Our proposed pipeline and  B.L.L. receive as input only cluster representative SNPs and therefore can only identify these as eVariants: we consider their discoveries correct if they correspond to a cluster that contains a true causal variants;  likewise we consider a causal variant discovered if the cluster that contains it has been selected. We further note that  the selection adjusted confidence intervals can be used to refine the selection of eVariants: while regularized regression might estimate a coefficient as different from zero, if the corresponding adjusted confidence interval covers zero, it make sense to discard the eVariant from the discovery set. Therefore, it seems appropriate to evaluate FDR and power on the basis of this final post-inference set.




  §.§.§ Results

Figures <ref> and <ref> summarize the results of the simulations, the first focusing on the inference on effect sizes, which is the object of the present paper, and the second anchoring these results in the context of FDR and Power.

Figure <ref> clearly illustrates that our pipeline succeeds in producing confidence intervals with the correct 
 empirical coverages, and that the lengths of these confidence intervals are not unduly large.  The failure of G.V. to produce confidence intervals with the right coverage is to be expected: the interest of our results is in showing the extent of the departure. While adjusting for eVariant selection as in  B.L.L. improves coverage, the approach still falls short of the target for variances in erroneously selected eGenes and for eGenes with a small number of causal variants. Moreover, the lengths of these adjusted confidence intervals are substantially longer than the ones we propose, which are just 1.5 times longer than the naive—this is the advantage of randomization.

 


 


Figure <ref> puts these result in context of the selection outputs of the three different methods.
Looking at eGene selections, we can see that  B.L.L. is unduly conservative at the eGene level as its FDR  is much lower than the that of our procedure. It is reassuring to note that Proposed has slightly higher power than the benchmark G.V. at the eGene selection level, while controlling for FDR at the level of 0.1.


At the eVariant selection stage, we report the FDR of the screening procedure standard (used in B.L.L.) and randomized regularized regression (used in Proposed) and the forward-backward selection (used in G.V.). The average number of eVariants per eGenes reported by the methods are shown in Table <ref>. Because the number of declared eVariants (post-inference) are determined by whether their intervals cover 0 or not among the screened evariants (pre-inference) for both B.L.L. and Proposed, we see the decrease in the average number of reported variants with the exception of forward-backward selection. In G.V., the eVariants (reported variants) are selected by an adaptive forward-backward selection and this report is not impacted by inference; thereby, the FDP curve is the same both pre and post inference for G.V.. 

Both the standard and randomized regularized regression do not appear to control the FDR  (Figure <ref> panel (d)):  correction based on the adjusted confidence intervals improves (Figure <ref> panel (e)) the performance, however still is quite inferior to the one of G.V.. This is not surprising as the adjusted approach aims at providing valid p-values, but does not apply any correction for multiplicity encountered at the eVariant level. It is rather reassuring that the selections of G.V., which practically coincide with those reported in <cit.>,  have low FDR: this documents the efforts of a large community of scientists whose focus was precisely the selection of eVariants. 
The better selection performance of G.V. underscores the opportunity of developing alternative  adjusted approaches that condition on selection strategies as the one adopted in GTEx. 
 
 
 





 §.§ Effect sizes in GTEx liver data 




We now turn to the analysis of the real expression data available for liver in <cit.>. We start by noting that the limited sample size (97 observations) makes it impossible to rely on data splitting: a selection based on a 50% of the data results in the detection of only 56 eGenes, a mere 4% of the eGenes reported in <cit.>. We analyze the data using the three procedures we compared via simulation, reporting  the  eVariants whose adjusted-interval estimates does not cover 0. 

The Venn-diagram <ref> compares the  eGenes selected by the three procedures: our pipeline identified 2216 eGenes, with 1683 in common with the G.V., which detects 1831 eGenes in all. The B.L.L. selects 1395 eGenes,  1341 of which are also identified by our proposal.
 


Figure <ref> illustrates the eVariants results: the number of detected eVariants and distribution of the  lengths of their confidence intervals. The pipeline we have constructed detects an average of  4-5 eVariants  per gene, while B.L.L. reports 2-3 eVariants on an average and G.V. 1 eVariant on average per eGene. Our simulations indicate that both our proposed procedure and B.L.L.  tend to have higher FDR than G.V., so that we cannot assume all the additional discoveries to be valid ones. Nevertheless, even assuming a false discovery rate as high as 50%, the  number of extra discoveries is such that we can expect that B.L.L. and our proposed method do allow the identification of a considerable number of additional true discoveries, or, in other words, that  G.V underestimates the number of eVariants. 
  

Table 2 compares the eVariant findings between G.V. and the proposed pipeline in the current work. The table in the left panel gives the number of common and exclusive eVariants reported by each procedure across the common eGenes; a total of 1683 eGenes are reported by both. The table in the right panel compares the eVariant findings across all eGenes reported by any of these proposals. These reported findings show that even, after a deduction of 50% of the discoveries from our proposal, we report potentially new discoveries that are not discovered by the benchmark analysis in G.V..






§ DISCUSSION

In many genomic studies, the first step of data analysis identifies interesting associations between variables. Inference of the parameters governing these associations is attempted only in a second stage.  Naive estimates, that rely on standard methods using the same data that suggested the importance of these associations, do not enjoy satisfactory statistical properties. A common strategy to avoid misleading conclusions is to resort to data splitting, using one portion of the data to identify parameters of interest and another to estimate them reliably.
When the studies, however, rely on scarce biological samples, data splitting (or deferring to a new sample for inference) is not a viable option. The studies of genetic regulation of gene expression in hard-to-access human tissues are a perfect exemplification of these challenges. On the one hand, to identify the parameters of interest one sifts through the possible associations between tens of thousands of genes and hundreds of thousands of DNA polymorphisms. On the other hand, the number of available samples is in the hundreds at the best.

In the present work we have explored to which extent the conceptual framework of conditional inference after selection can be brought to bear to provide researchers with reliable tools to estimate the effect size associated to DNA variants that are discovered to be associated with variation in gene expression. Following what is the standard practice in the eQTL community, we have described a two stage process for the selection of relevant cis variants: first genes under cis regulations are identified, and then the subset of important variants in their vicinity. Recognizing both the appeal that regularized regression has in the statistical community, and the ease with which we can condition on selection events based on the results of lasso, our selection of important variants for a gene of interest is carried out with a regularized regression. Acknowledging that a complete conditioning on the selection event leaves little information in the data for inference, we have exploited randomization techniques.

Deriving appropriate conditional inference in this setting presented a number of challenges. In addition to having to deal with a two-stages selection process, calculation of interval and point estimates is more involved post randomization unlike the easy computation of intervals based on a truncated Gaussian law in <cit.>.
To bypass the fact that the exact selection-adjusted law lacks a closed form expression, we introduce an approximation strategy that is likely to be useful in other contexts. Further, we highlight in this work a clear pay-off in terms of inferential power associated with randomized mining strategies in comparison to <cit.> despite incurring some cost involved in making inference in the randomized framework tractable.

A simulation study carried out starting from real data underscored the dangers of ignoring selection at the inferential stage: naive confidence intervals constructed for the effect sizes of the variants identified with a published strategy missed the target coverage by a wide margin. The pipeline we developed leads to confidence intervals with correct coverage and with lengths that are appreciably shorter than those obtained with out-of-the-box tools for conditional inference. While this results are encouraging, we also noted that the selection procedure  in our pipeline has appreciably worse performance that state-of-the-art variant selection strategies in terms of FDR. Given that we resorted to randomization to reserve information for inference, we expect a loss of precision at the selection stage. However, the size of the difference we observe is likely not to be ascribable only to randomization, but to the property of the selection procedure itself. Specifically, we observe that the regularized regression we adopted does not enjoy FDR control. In fact, we realize that the techniques in this work are amenable to a more general framework of convex learning programs, with the proposed mining pipeline as a specific example in this broader class of selection schemes. It is certainly an interesting direction for future work to study if computation for conditional inference can be carried out for other selection strategies with better FDR control, as SLOPE, the knockoffs etc. 

While we have provided effect size estimates for single tissue eQTL experiments, borrowing strengths from samples across multiple tissues (see <cit.>) is likely to lead to enhanced power in both discoveries and subsequent inference on their effect sizes. The adjusted likelihood provided by our work post an adequate selection strategy for discoveries based on samples, combined across different tissues can be potentially incorporated into the inferential framework. Exploring the proposed pipeline for inference in other genomic applications including the GWA (genome-wide association) studies to effectively calibrate effect sizes post discoveries by optimally using the information in the available data samples remains a challenge ahead. We conclude by remarking that we hope to take on these challenges in other biological contexts with this first attempt of tractably addressing the issue of selection bias in effect size estimation in eQTL studies.




§ ACKNOWLEDGEMENTS

The authors acknowledge that data for the gene expression analysis, described in this manuscript was acquired using dbGaP accession number phs000424.v6.p1. S.P. would like to thank Jonathan Taylor for several helpful discussions regarding this project. J.Z. acknowledges support from Stanford Graduate Fellowship. C.S. acknowledges support from NSF DMS 1712800 and NIH R01MH101782.  

chicago




§ AN UNSUPERVISED PRUNING OF LOCAL VARIANTS


The distance matrix {d_i,j,  i, j ∈  1,2,⋯, V_g} based on which we implement an unsupervised pruning of variants, measured within a cis-window around a gene g in the genome, is defined as
d_i,j = d(X_i, X_j)= 1- ρ(X_i, X_j); ρ denoting the Pearson correlation between variants X_i and X_j. To perform a clustering of variants and obtain a prototype for each cluster, we implement a hierarchical clustering with a minimax linkage, explored in <cit.>. This is based on a distance between clusters 𝒞_1 and 𝒞_2 defined as 

    d(𝒞_1, 𝒞_2) = d_0 (𝒞_1, 𝒞_2), where  d_0(𝒞) = min_X_i∈𝒞max_X_j∈𝒞d(X_i, X_j).
 
The prototype for each cluster 𝒞 is determined as 
    min_X_i∈𝒞max_X_j∈𝒞d(X_i, X_j).

To determine the number of clusters for a set of variants, we cut the minimax tree at a height 1-ρ_0 that  yields a dataset of prototypical variants in which every variant has correlation of at least ρ_0 with one of the prototype SNPs. For our implementations, we use a ρ_0= 0.5 so that the prototypes E^(g) selected by a LASSO analysis can be interpreted as identifying all the local variants that are correlated by 0.5 with any of these representatives to be the set of promising functional variants. Prior to the secondary analysis in both the simulation study in <ref> and the GTEx data analysis in <ref>, we perform a pruning of cis-variants by the above prescribed scheme. In the simulations, pruning is seen to reduce the size of local variants for 1770 selected eGenes with an average number of SNPs of 4497 to 229 SNPs on an average; in the eQTL experiment, the averaged pruned size of variants is 262 for 2261 eGenes with 5101 variants on an average.













§ DECOUPLED SELECTION EVENT: WITH EXTRA CONDITIONING






















Theorem <ref> facilitates a decoupling of the truncated law across eGenes by considering a selection event 

    {g ∈𝒢̂,   K = K_0,    j_(1)^(g)=  j_0,  T_j_(2)^(g)= T_0, sign(T^(g)_j_0) = s_j_0,
    
            Ê(y^(g),ζ^(g))= E^(g), sign(β̂_E^(g)) = s_E^(g)}

with additional information. The decoupled law is computationally less burdensome to handle. We now give a proof of this Lemma.
 Separability result in <ref>
This selection event can be written as


    [0.95]  {g ∈𝒢̂,   K = K_0,    j_(1)^(g)=  j_0,   T_j_(2)^(g)= T_0, sign(T^(g)_j_0) = s_j_0,
    
       [0.95]     Ê(y^(g),ζ^(g))= E^(g), sign(β̂_E^(g)) = s_E^(g)}
       [0.95]= {p̃^(g)≤K_0/G q ,    j_(1)^(g)=  j_0,   T_j_(2)^(g)= T_0, sign(T^(g)_j_0) = s_j_0,
       [0.95]       Ê(y^(g),ζ^(g))= E^(g) , sign(β̂_E^(g)) = s_E^(g),
       [0.95]        p̃^(g')_(k)≥(K_0 + k)/G q   for  g' ∉𝒢};   where {p̃^(g')_(k)} denote ordered p-values amongst non eGenes g' ∉𝒢.
    
    
    
    
    
       =[0.95]{s_j_0 T_j_0^(g)≥ |T_0| ,  s_j_0 T_j_0^(g)≥√(1+γ^2)·Φ^-1(1- K_0/2V_g Gq) 
       [0.95]         Ê(y^(g),ζ^(g))= E^(g), sign(β̂_E^(g)) = s_E^(g)}⋂{p̃^(g')_(k)≥(K_0 + k)/G q   for  g' ∉𝒢}
       =[0.95]{s_j_0(X_j_0^(g)^T y^(g)  + ω_j_0^(g)) ≥max(√(1+γ^2)·Φ^-1(1- K_0/2V_g Gq), |T_0^(g)|), 
       [0.95]        Ê(y^(g),ζ^(g))= E^(g), sign(β̂_E^(g)) = s_E^(g)}⋂{p̃^(g')_(k)≥(K_0 + k)/G q   for  g' ∉𝒢}.

The second step follows by noting that if K_0 is the number of rejections in the BH, then the rejected p-values satisfy

    p̃^(g)≤K_0/G q;

and the genes g' (ordered in terms of their Bonferroni-adjusted p-values) that are discarded at this stage have corresponding  p-values that satisfy

    p̃^(g')_(k)≥(K_0 + k)/G q.

The last equality describes the selection of eGenes in terms of t-statistics.


Define the event {ℋ̂(y^(g), ω^(g), ζ^(g))= ℋ} as 

    [0.95]{s_j_0(X_j_0^(g)^T y^(g)  + ω_j_0) ≥max(√(1+γ^2)·Φ^-1(1- K_0/2V_g Gq), |T_0^(g)|),
       [0.95]         Ê(y^(g),ζ^(g))= E^(g),sign(β̂_E^(g)) = s_E^(g)},

clearly  a function of (y^(g),ω_j_0^(g),ζ^(g)), exclusive to eGene g.
The second component {𝒥̂(y^(g'), ω^(g')), g'∉𝒢) =𝒥} is defined as

    [0.97]{p̃^(g')_(k)≥(K_0 + k)/G q   for  g' ∉𝒢},

a function of (y^(g'), ω^(g')) associated with genes not selected by BH. This completes the proof that the selection region under consideration takes a separable form.




Conditioning on nuisance statistics and marginalizing over randomizations gives a selection-modified family of exponential laws in Lemma <ref> . We include a proof for the same.
Exponential law in <ref>.
With a slight abuse of notation (noting that we use the same ℋ to denote the selection region), we describe the selection event in terms of data and randomization as

    {(b̂_j; E^(g), 𝒰_j^(g), ω^(g), ζ^(g))∈ℋ}.

The joint law of data and randomizations introduced in selections is proportional to

    exp(-(b̂_j;E^(g) -b_j;E^(g))^22σ_j;E^(g)^2)·exp(-Σ_U^-1/2(𝒰_j^(g)-n_j)_2^2)
       ·exp(-ω^(g)_2^22γ^2)·exp(-ζ^(g)_2^22τ^2)· 1_{(b̂_j; E^(g), 𝒰_j^(g), ω^(g), ζ^(g))∈ℋ}

Conditioning on null-statistics 𝒰^(g)_j= u^(g)_j, the marginal law of b̂_j;E^(g) based on the above joint distribution is proportional to

    exp(-(b̂_j;E^(g) -b_j;E^(g))^2/2σ_j;E^(g)^2) ∫exp(-ω^(g)_2^2/2γ^2)·exp(-ζ^(g)_2^2/2τ^2)·
                                                                        1_{(b̂_j; E^(g),ω^(g), ζ^(g)) ∈ℋ(u^(g)_j)} dω d ζ.
       =exp(-(b̂_j;E^(g) -b_j;E^(g))^2/2σ_j;E^(g)^2) ·ℙ[(b̂_j; E^(g),Ω^(g), Z^(g)) ∈ℋ(u^(g)_j)  |b̂_j; E^(g)=b̂_j; E^(g),
                                                                               𝒰^(g)_j= u^(g)_j].


This gives the adjusted marginal law of b̂_j; E^(g), as stated above.




§ K.K.T. MAPPING FOR SELECTION OPTIMIZATIONS


We make explicit the selection maps in (<ref>) and (<ref>) that characterize the identification of a set of eGenes, denoted by 𝒢 and the selection of active variants corresponding to each such eGene, represented by the set E^(g). 

These maps are obtained from the K.K.T. conditions that characterize the solutions of the underlying optimizations/ selections. The K.K.T. maps are, in both the selection of eGenes and active variants with the additionally conditioned information described in <ref>, affine maps in both the target statistic b̂_j;E^(g) and randomizations, ω^(g) and ζ^(g).




  * Selection map for eGenes: 
The selection of eGene g (with additional information) in the event {ℋ(y^(g), ω^(g), ζ^(g))= ℋ} is given by

    {s_j_0(X_j_0^(g)^T y^(g) + ω_j_0) ≥max(√(1+γ^2)·Φ^-1(1- K_0/2V_g Gq), |T_0^(g)|)}.

This is equivalent to a map


    η^(g) = s_j_0(X_j_0^(g)^T y^(g)  + ω_j_0^(g))

where η^(g) satsifies constraints 

    [0.95]𝒞_1 = {η: η≥max( √(1+γ^2)·Φ^-1(1- K_0/2V_g Gq), |T_0^(g)|) }.

Decomposing the data-vector X_j_0^(g)^T y^(g) into the statistic of interest b̂_j; E^(g) and null-statistics that are orthogonal to b̂_j; E^(g), this map can be re-written as follows.

    ω_j_0^(g)   = -X_j_0^(g)^TX̃_E^(g) (X̃_E^(g)^T X̃_E^(g))^-1e_j b̂_j; E^(g)/σ_j;E^(g)^2 + s_j_0η^(g)
             - X_j_0^(g)^T (y - X̃_E^(g) (X̃_E^(g)^T X̃_E^(g))^-1e_jb̂_j; E^(g)/σ_j;E^(g)^2)
       =  P^(g)_j  b̂_j; E^(g)  + s_j_0η^(g) + ℳ^(g)_j,  where


    [0.95]ℳ^(g)_j=- X_j_0^(g)^T (y - X̃_E^(g) (X̃_E^(g)^T X̃_E^(g))^-1e_jb̂_j; E^(g)/σ_j;E^(g)^2) ;


    [0.95]P^(g)_j  = -X_j_0^(g)^TX̃_E^(g) (X̃_E^(g)^T X̃_E^(g))^-1e_j .


Conditioning on ℳ^(g)_j and denoting it as q^(g)_j, the K.K.T. map associated with the selection of eGene g is represented as

    ω_j_0^(g)   = P^(g)_j  b̂_j; E^(g)  + s_j_0η^(g) + q^(g)_j.



  * Selection map for active variants (potential eVariants): 
Fixing notations, we let o_E^(g) and o_-E^(g) represent the active signs and inactive sub-gradients, obtained upon solving the LASSO objective. We denote the vector [  o_E^(g); o_-E^(g) ] as o^(g). The K.K.T map associated with the solution of the randomized LASSO that selects a set of variants, E^(g) with signs s_E^(g) is given by


    ζ^(g)   =  -1.5[  X̃_E^(g)^T X̃_E^(g)                    0; X̃_-E^(g)^T X̃_E^(g)                    I ][                   b̂_E; X̃_-E^T(y - X̃_E b̂_E) ]+ 1.5[ X̃_E^(g)^T X̃_E^(g) +ϵ I                        0;     X̃_-E^(g)^T X̃_E^(g)                        I ][  o_E^(g); o_-E^(g) ] + [ λ s_E^(g);         0 ]
       =-1.5[  X̃_E^(g)^T X̃_E^(g)                    0; X̃_-E^(g)^T X̃_E^(g)                    I ][ (X̃_E^(g)^T X̃_E^(g))^-1 e_j;                            0 ]b̂_j; E^(g)/σ_j;E^(g)^2 + 1.5[ X̃_E^(g)^T X̃_E^(g) +ϵ I                        0;     X̃_-E^(g)^T X̃_E^(g)                        I ][  o_E^(g); o_-E^(g) ]
             + [ λ s_E^(g);         0 ] + 𝒩^(g)_j  
       = A_j^(g)b̂_j; E^(g) + B_j^(g)o^(g) + c_j^(g)


where o^(g) satisfies the constraints

    𝒞_2 ={o: sign(o_E^(g)) = s_E^(g),  o_-E^(g)_∞≤λ}.

The matrices/ vectors  A_j^(g), B_j^(g), c_j^(g) are equal to

    A_j^(g) = -1/σ_j;E^(g)^21.5[  X̃_E^(g)^T X̃_E^(g)                    0; X̃_-E^(g)^T X̃_E^(g)                    I ][ (X̃_E^(g)^T X̃_E^(g))^-1 e_j;                            0 ], B_j^(g)=   1.5[ X̃_E^(g)^T X̃_E^(g) +ϵ I                        0;     X̃_-E^(g)^T X̃_E^(g)                        I ] and


    c_j^(g)=[ λ s_E^(g);         0 ] + 𝒩^(g)_j.

In the above,

    𝒩^(g)_j= -1.5[  X̃_E^(g)^T X̃_E^(g)                    0; X̃_-E^(g)^T X̃_E^(g)                    I ]{[                   b̂_E; X̃_-E^T(y - X̃_E b̂_E) ] -  [ (X̃_E^(g)^T X̃_E^(g))^-1 e_j;                            0 ]b̂_j; E^(g)/σ_j;E^(g)^2 }

representing null statistics in (<ref>) in a saturated model framework. This map also relies on an affine decomposition into the statistic b̂_j; E^(g) and null statistics 𝒩^(g)_j orthogonal to it.


Finally, a decomposition into active (non-zero coordinates) denoted by E^(g) and the non-active coordinates that are shrunk to 0 by the LASSO yields a K.K.T. map in terms of (ζ^(g), b̂_j; E^(g), o^(g)) as

    1.5[  ζ_E^(g); ζ_-E^(g) ] = 1.5[ A_E^(g),jb̂_j; E^(g)+ B_E^(g),j o_E^(g)  + c _E^(g),j;             A_-E^(g),j + B_-E^(g),j o_E^(g)+ o_-E^(g) ].

Here, A_E^(g),j, B_E^(g),j, c _E^(g),j denote the selected subset of rows corresponding to active coordinates E^(g). Similarly, A_-E^(g),j, B_-E^(g),j, c _-E^(g),j represent the subset of rows that are associated with the inactive components or coefficients shrunk to 0 by the LASSO.




§ CHANGE OF VARIABLES FORMULAE


It is easy to see that the maps in (<ref>) and (<ref>) together with the constraints on (η^(g), o^(g)) give rise to half spaces defined in terms of (b̂_j; E^(g), ω^(g), ζ^(g)) as selection-induced regions. That is, they yield a polytope

    {(b̂_j;E^(g), ω^(g), ζ^(g)): (b̂_j;E^(g), ω^(g), ζ^(g)) ∈ℋ}.

The reference measure in the selection-adjusted law of b̂_j;E^(g) is based on computing
 
    𝒫_ℋ(t) = ℙ((b̂_j; E^(g), Ω^(g), Z^(g)) ∈ℋ|b̂_j;E^(g)=t),

 which is the Gaussian volume of a polytope. In general, a multivariate Gaussian volume of a polyhedral region might be hard to compute. However, we can use the K.K.T. maps associated with identifying eGenes and active variants as change of variables formulae from the space of 
(b̂_j; E^(g), ω^(g)_j_0, ζ^(g)) to (b̂_j; E^(g), η^(g), o^(g)). This reduces the calculation of the reference measure to volumes of orthants and cubes, that are geometrically much simpler in nature. 



Adjusted law induced by change of measure:
Based on change of variables formulae

    ω_j_0^(g) = P^(g)_j  b̂_j; E^(g)  + s_j_0η^(g) + q^(g)_j  and


    ζ^(g)=A_j^(g)b̂_j; E^(g) + B_j^(g)o^(g) + c_j^(g);

the marginal selection-adjusted law of b̂_j; E^(g), conditioned on nuisance statistics (ℳ^(g)_j, 𝒩^(g)_j) is proportional to

    [0.95]exp(-(b̂_j;E^(g) -b_j;E^(g))^22σ_j;E^(g)^2) ·𝒫_𝒞_1; 𝒞_2(b̂_j;E^(g))= exp(-(b̂_j;E^(g) -b_j;E^(g))^22σ_j;E^(g)^2) ·ℙ[η^(g)∈𝒞_1, o^(g)∈𝒞_2  |b̂_j; E^(g)].



Conditioning on null statistics together with a change of measure induced by the described selection maps yields a joint law for (b̂_j; E^(g), η^(g), o^(g)) as

    |J|·exp(-(b̂_j;E^(g) -b_j;E^(g))^2/2σ_j;E^(g)^2)·exp(-(P^(g)_j  b̂_j; E^(g)  + s_j_0η^(g) + q^(g)_j)^2/2γ^2)
       ·exp(-A_j^(g)b̂_j; E^(g) + B_j^(g)o^(g) + c_j^(g)_2^2/2τ^2)· 1_{η^(g)∈𝒞_1, o^(g)∈𝒞_2}

with |J| being the Jacobian, a constant in our case. It is trivial to see from here that a marginalization over (η^(g), o^(g)) yields the adjusted law with a reference measure 

    𝒫_𝒞_1; 𝒞_2(t)= ℙ[η^(g)∈𝒞_1, o^(g)∈𝒞_2  |b̂_j; E^(g)=t].





§ APPROXIMATING OPTIMIZATION FOR AFFINE GAUSSIAN VOLUME


Theorem <ref> gives an expression for the probability ℙ(η^(g)∈𝒞_1,    o^(g)∈𝒞_2|b̂_j; E^(g)  = t) as a function of t where 𝒞_1 and 𝒞_2 represent regions based on the selection constraints in (<ref>) and (<ref>). Theorem <ref> derives an upper bound for the same and hence, bounds from the reference measure 𝒫_𝒞_1; 𝒞_2(·) for convex and compact sets 
𝒞_1 ⊂ and 𝒞_2⊂^p_g.



Reference measure in selection-adjusted law:
For the selection regions 𝒞_1⊂ and 𝒞_2⊂^p_g in (<ref>) and (<ref>),
the probability of selection
ℙ(η^(g)∈𝒞_1,    o^(g)∈𝒞_2|b̂_j; E^(g)  = t) equals

    K · C_1(t) ·∫_o_Eexp(-A_E^(g),j t+ B_E^(g),j o_E  + c _E^(g),j_2^2 /2τ^2)C_2(o_E, t )1_{sign(o_E)= s_E^(g)} do_E

where

    C_1(t)  = Φ̅({L+ s_j_0P^(g)_j  t + s_j_0q^(g)_j}/γ)  and


    C_2(o_E, t ) = ∏_k ∉ E^(g){Φ({A_-E^(g),j^kt + B_-E^(g),j^ko_E+ λ}/τ) - Φ({A_-E^(g),j^k t + B_-E^(g),j^k o_E -  λ} /τ)}

and K is some constant that does not depend on t;  A_-E^(g),j^k represents the k-th row of A_-E^(g),j and B_-E^(g),j^k denotes the k-th row of B_-E^(g),j; L is the lower threshold in the set of constraints 𝒞_1. 



We use notations K_1, K_2, ⋯ to denote constants in our proof; these dissolve as constants and do not impact the computation of the reference measure as a function of b̂_j; E^(g).

    ℙ(η^(g)∈𝒞_1,    o^(g)∈𝒞_2|b̂_j; E^(g)  = t)    = K_1·∫_o∫_ηexp(-(P^(g)_j  t + s_j_0η + q^(g)_j)^2/2γ^2)1_{η∈𝒞_1}dη
                ·exp(-A_j^(g)t + B_j^(g)o^(g) + c_j^(g)_2^2/2τ^2) 1_{o ∈𝒞_2}do
       = K_1·∫_ηexp(-(P^(g)_j  t   + s_j_0η + q^(g)_j)^2/2γ^2)1_{η∈𝒞_1}dη
              ·∫_oexp(-A_j^(g)t + B_j^(g)o + c_j^(g)_2^2/2τ^2) 1_{o ∈𝒞_2}do

Denoting the lower threshold in the set of constraints 𝒞_1 as L, we have

    ∫_ηexp(-(P^(g)_j  t  + s_j_0η + q^(g)_j)^2/2γ^2)1_{η∈𝒞_1}dη
       = K_2·Φ̅({L+ s_j_0P^(g)_j  t+ s_j_0q^(g)_j}/γ)= K_2 · C_1(β̂_j; E^(g))

In the computation of 

    ∫_oexp(-A_j^(g)t + B_j^(g)o^(g) + c_j^(g)_2^2/2τ^2) 1_{o ∈𝒞_2}do,

we marginalize over the inactive sub-gradient variables first. That is, this integral equals

    ∫_o_Eexp(-A_E^(g),j t+ B_E^(g),j o_E  + c _E^(g),j_2^2 /2τ^2)1_{sign(o_E)= s_E^(g)}
         ·∫_o_-Eexp(-A_-E^(g),j t + B_-E^(g),j o_E+ o_-E_2^2/2τ^2) 1_{o_-E_∞≤λ} do_-E d o_E
       = K_3 ·∫_o_Eexp(-A_E^(g),j t + B_E^(g),j o_E  + c _E^(g),j_2^2 /2τ^2)1_{sign(o_E)= s_E^(g)}C_2(o_E, b̂_j; E^(g) ) do_E.

where C_2(o_E, b̂_j; E^(g) ) equals

    ∏_k ∉ E^(g){Φ({A_-E^(g),j^kt + B_-E^(g),j^k o_E+ λ}/τ) - Φ({A_-E^(g),j^kt + B_-E^(g),j ^ko_E -  λ} /τ)}.

The above is a consequence of the separability of the inactive sub-gradient equation and separability of the cube region induced by the same. Modulo a constant K, we can thereby write ℙ(η^(g)∈𝒞_1,    o^(g)∈𝒞_2|b̂_j; E^(g)  = t) as 

    C_1(β̂_j; E^(g)) ·∫_o_Eexp(-A_E^(g),jb̂_j; E^(g)+ B_E^(g),j o_E  + c _E^(g),j_2^2 /2τ^2)C_2(o_E, b̂_j; E^(g) )1_{sign(o_E)= s_E^(g)} do_E.






An upper bound on reference:
For convex and compact sets 𝒞_1 ⊂, 𝒞_2⊂^p_g, an upper bound for logℙ(η^(g)∈𝒞_1,  o^(g)∈𝒞_2|b̂_j; E^(g)  = t) is






    log K + log C_1(t)-inf_sign(o_E) = s_E^(g){ A_E^(g), j  t + B_E^(g), j o_E + c_E^(g), j_2^2/2τ^2- log C_2(o_E, t).}

where K is the same constant in the conditional probability in Theorem <ref>.



Based on Lemma <ref>, we have 

    logℙ(η^(g)∈𝒞_1,    o^(g)∈𝒞_2|b̂_j; E^(g)  = t) 
       = log K + log C_1(t) + log𝔼[ C_2(o_E^(g), b̂_j; E^(g))1_{sign(o_E^(g)) = s_E}|b̂_j; E^(g)  = t].

Ignoring the constant above,

    log C_1(t) + log𝔼[ C_2(o_E^(g), b̂_j; E^(g) )1_{sign(o_E^(g)) = s_E}|b̂_j; E^(g)  = t]
       ≤log C_1(t)+ log𝔼[exp{sup_sign(o_E) = s_Elog C_2(o_E, t)-β^T o_E}exp(β^T o_E^(g))|b̂_j; E^(g)  = t]
    
    
       ≤log C_1(t)-sup_β∈^E^(g)inf_sign(o_E) = s_E{β^T o_E-log C_2(o_E, t) }- log𝔼[exp(β^T o_E^(g))|b̂_j; E^(g)  = t]
       = log C_1(t)-inf_sign(o_E) = s_Esup_β∈^E^(g){β^T o_E- log𝔼[exp(β^T o_E^(g))|b̂_j; E^(g)  = t]}-log C_2(o_E, t).

The last step follows by a minimax equality for convex and compact sets 𝒞_1, 𝒞_2.

Finally, using the fact that 
    sup_β∈^E^(g){β^T o_E- log𝔼[exp(β o_E)|b̂_j; E^(g)  = t]}
 is the conjugate of logarithm of the moment generating functions of Gaussian random variable o_E^(g) and equals 

    A_E^(g), j  t + B_E^(g), j o_E + c_E^(g), j_2^2/2τ^2,
 we have the stated bound.



We use a smooth version of the above objective in <ref> through a barrier penalty. The final optimization that we solve to approximate the reference is given by

    log C_1(t) -inf_o_E{ A_E^(g), j  t + B_E^(g), j o_E + c_E^(g), j_2^2/2τ^2- log C_2(o_E, t) + ℬ(o_E)}

with ℬ(·) being a barrier penalty that reflects the sign constraints.

Let the approximate reference measure be 𝒫̂_ℋ(b̂_j;E^(g)). The resulting pseudo selection-adjusted law for b_j;E^(g)  is proportional to

    exp(-(b̂_j;E^(g) -b_j;E^(g))^2/2σ_j;E^(g)^2) ·𝒫̂_𝒞_1;𝒞_2 (b̂_j;E^(g)).

In fact, a pivot based on the approximate law is calculated as

    p̂(b̂_j;E^(g); b, σ_j; E^(g)^2) =∑_t∈χ: ≥b̂_j;E^(g)exp(-(t -b_j;E^(g))^2/2σ_j;E^(g)^2) ·𝒫̂_𝒞_1;𝒞_2 (t)∑_t'∈χexp(-(t' -b_j;E^(g))^2/2σ_j;E^(g)^2) ·𝒫̂_𝒞_1;𝒞_2 (t')

for χ being an appropriate grid on the real line. 
As pointed out in the intractability of the MLE problem, solving for the selection-adjusted MLE based on (<ref>) requires the log-partition function as a function of the target parameter b_j;E^(g). 
Based on an approximated partition function on a grid χ⊂,  the MLE is the solution of the optimization

    b_j;E^(g)∈minimize(b̂_j;E^(g) -b_j;E^(g))^2/2σ_j;E^(g)^2  + log∑_t∈χexp(-(t -b_j;E^(g))_2^2/2σ_j;E^(g)^2)·𝒫̂_𝒞_1;𝒞_2 (t).




§ DATA PROCESSING FOR GTEX CIS-EQTL DATA



The variant calls, which had been pre-processed through GTEx-specific quality control and variant imputation pipelines, were downloaded from dbGaP (accession phs000424.v6.p1),. We followed the processing pipeline for cis-eQTL according to <cit.>: (1) the search for variants was limited to ±1 Mb of the transcription start site of each gene and variants were removed if they occurred in less than 10 samples of had minor allele frequency less than 0.01; (2) a total of 20 covariates, including technical confounders, population background, gender and genotyping platform, were regressed out from the genotypes (i.e., variant matrices X^(g))) as well as phenotypes (i.e., gene expressions y^(g)). We modified the original software fastQTL (<cit.>) to implement these filters and corrections and obtain the final gene-specific genotype matrices. The processed gene expression data was also downloaded from dbGaP, and we regressed out the same 20 covariates as we did for the genotype data for the real cis-eQTL data analysis.

