{'fraction_non_alphanumeric': 0.035509352410760864, 'fraction_numerical': 0.015463149265966168, 'mean_word_length': 4.8085705010280275, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient inference for genetic association studies with multiple outcomes

Combined inference for heterogeneous high-dimensional data is critical in modern biology, where clinical and various kinds of molecular data may be available from a single study. Classical genetic association studies regress a single clinical outcome on many genetic variants one by one, but there is an increasing demand for joint analysis of many molecular outcomes and genetic variants in order to unravel functional interactions. Unfortunately, most existing approaches to joint modelling are either too simplistic to be powerful or are impracticable for computational reasons. Inspired by Richardson et al. (2010, Bayesian Statistics 9), we consider a sparse multivariate regression model that allows simultaneous selection of predictors and associated responses. As Markov chain Monte Carlo (MCMC) inference on such models can be prohibitively slow when the number of genetic variants exceeds a few thousand, we propose a variational inference approach which produces posterior information very close to that of MCMC inference, at a much reduced computational cost. Extensive numerical experiments show that our approach outperforms popular variable selection methods and tailored Bayesian procedures, dealing within hours with problems involving hundreds of thousands of genetic variants and tens to hundreds of clinical or molecular outcomes.


Introduction
Much current research in genetics focuses on combining heterogeneous data for the same samples. This is prompted by the increasing availability of diverse molecular data types from a single study and should lead to a more complete understanding of biological systems, as combined inference for data from different molecular layers might unravel regulatory interactions within and across layers (Civelek and Lusis, 2014). An example is protein quantitative trait locus (pQTL) analyses to detect associations between hundreds of thousands of genetic variants and hundreds of proteomic expression levels. When associated with disease genetic variants, proteins are often regarded as intermediate phenotypes or molecular proxies for the disease of interest, as they may provide direct insights on biological processes underlying the clinical condition, and likewise for other molecules such as genes or metabolites involved in so-called eQTL or mQTL analyses.
An important goal of molecular QTL studies is to detect genetic variants with pleiotropic effects, i.e., variants that regulate the expression levels of several molecules, as the genome location where they lie may initiate essential functional mechanisms (Breitling et al., 2008). As pleiotropy has been acknowledged as a central property of genetic variants causing phenotypic variation, efforts are being made to uncover its activity patterns and gain understanding of the shared biological processes it induces (Sivakumaran et al., 2011;Solovieff et al., 2013). It is also of interest to identify those cases where the same molecule is simultaneously influenced by several genetic variants. This dual task requires a model that allows flexible selection, ideally performed jointly on the genetic variants and the molecular outcomes.
Although in practice univariate analyses still dominate, several proposals for joint modelling of multiple outcomes in genetic association problems have recently been made. Flutre et al. (2013) and Zhou and Stephens (2014) model the outcomes as multivariate responses having a matrixvariate distribution. The former rely on a Bayes factor framework to uncover the associations between a given genetic variant and any subset of outcomes, whereas the latter propose a linear mixed model whose random effect accounts for relatedness among individuals. While these methods show improvements over the fully marginal regression approach, they are restricted to problems with a few outcomes, because their models involve unstructured covariance matrices. O'Reilly et al.
(2012)'s MultiPhen method reverses the classical regression setup and fits a succession of models where each genetic variant is regressed on several outcomes. This eliminates the need to model large covariance matrices, but also entails a marginal treatment of the genetic variants. Moreover, MultiPhen does not penalize model complexity, which may cause instabilities when many outcomes are modelled. Molecular QTL problems are particularly complex, because in addition to the socalled p n paradigm, whereby the number of covariates (genetic variants) p greatly exceeds the number of samples n, there is also a "large d" characteristic, as the number of responses (expression levels) d is usually large. Methodologies to accommodate this are needed, especially when joint response modelling is sought.
Two-stage procedures are natural approaches to association problems with both p and d large.
2HiGWAS (Jiang et al., 2015) is essentially an implementation of the screening strategy of Fan and Lv (2008) in the context of longitudinal outcomes. It consists of a dimension reduction step, where each genetic variant is tested against each outcome, followed by a penalized regression recast into functional mapping in which only the genetic variants remaining after screening are involved.
While the second stage of 2HiGWAS is an interesting approach to joint covariate and response modelling, the method is not tailored to typical molecular QTL analysis, as it is designed for outcomes measured over time. Also, the fully marginal first stage screening, if too stringent, may cause important predictors to drop out of the analysis and thus lead to false negatives. Instead of pruning the covariate set, Wang et al. (2016) summarize information at outcome level in the context of eQTL analyses. They precluster the expression levels using a block-mixture model and test for association between each genetic variant and the resulting clustering. This approach requires that the stability of the clustering and its functional relevance are carefully checked, as the group pattern chosen is critical to subsequent analysis. More generally, two-stage procedures often rely on ad-hoc thresholding decisions at the first stage which influence the conclusions of the second stage.
The approaches sketched above use very diverse strategies to model predictor and outcome variables in high-dimensional settings. Trade-offs between realistic modelling and computational efficiency are inevitable, but two components seem critical to practical and powerful analyses: involving all variables in a single multivariate model, and maintaining interpretable inference.
Unified selection of genetic variants and associated outcomes is then possible, unlike for most existing methods, whose focus is on selecting either predictors or outcomes. The Bayesian framework seems particularly suitable, as it offers flexible modelling possibilities in which biological beliefs may be naturally incorporated. The hierarchical regression approach of Richardson et al. (2010), HESS (hierarchical evolutionary stochastic search), is an appealing example of such approaches; it can identify associations between hundreds of covariates and up to a few thousand responses from a single model. Jia and Xu (2007) and Scott-Boyer et al. (2012) propose methods called BAYES and iBMQ (integrated Bayesian hierarchical model for eQTL mapping) based on models similar to that of HESS, but a major drawback of all three approaches is the lack of scalability of their Markov chain Monte Carlo inference procedures. Problems whose size corresponds to actual genome-wide association studies with molecular outcomes (with hundreds of thousands of genetic variants and hundreds to thousands of outcomes and individuals) are out of reach even for HESS, which is based on adaptive parallel tempering/evolutionary Monte Carlo techniques. We are unaware of any fully multivariate approach that can deal with such data within a reasonable time.
In this paper, we propose a Bayesian inference strategy that avoids sampling, via an algorithm that is fast and whose convergence is easy to monitor, while having performance comparable with Markov chain Monte Carlo approaches. We describe a variational inference procedure for a model similar to that of Richardson et al. (2010), comprising a series of parallel linear regressions, one for each response, combined in a hierarchical manner to leverage shared information. A spike-and-slab prior (Ishwaran and Rao, 2005) is used to induce sparsity of the regression coefficients, and the probability that a given covariate affects any response is modelled through a parameter that is shared across responses. Interpretable posterior quantities, such as the probability of association of each covariate-response pair, are produced. An efficient algorithm for this model is crucial, as its number of parameters can be very large. Our variational approach can update the parameters jointly in a tractable manner. Carbonetto and Stephens (2012) provide a good discussion of variational inference in the context of genetic association and propose a variational regression method called "varbvs", which can be seen as a single-response counterpart of our approach. In our multiple response setting, we show by simulation that our procedure is accurate and reliable, and is better than existing methods at selecting variables for very large problems. Hence, it offers clear added-value in practice: it enables complex and flexible Bayesian inference based on large batches of genetic data, without having to prune them beforehand.
The paper is organized as follows. Section 2 describes the model, discusses its relation to earlier proposals, and presents a procedure to allow for sparsity control at both covariate and response levels. Section 3 gives an overview of variational Bayes approaches and describes our inference strategy. Section 4 compares variational and Markov chain Monte Carlo inferences on the same model, also using direct approximations of posterior quantities. Section 5 describes numerical experiments for larger problems, comparing our method with several predictor selection methods, including the varbvs approach of Carbonetto and Stephens (2012), and with methods performing combined covariate and response selection, namely HESS (Richardson et al., 2010) and iBMQ (Scott-Boyer et al., 2012). The section also presents a permutation-based comparison of our method with varbvs on a real mQTL problem. Section 6 summarizes the discussion and highlights further possible developments.
Although the applicability of our method is not restricted to any particular context, all the numerical experiments presented in this paper use settings tailored to genome-wide association studies. The data-generation schemes are designed to embody common biological assumptions and are described in Section 5.1. The method is implemented in the publicly available R package locus.

Model and earlier proposals
Let y = (y 1 , . . . , y d ) be a n × d matrix of d centered responses and let X be a n × p matrix of p covariates, for each of n samples. The covariates are p genetic variants, more precisely single nucleotide polymorphisms, SNPs, and the responses might represent d gene, protein or metabolite expression levels for n individuals, depending on whether an eQTL, pQTL or mQTL problem is considered. Our model is intended to accommodate all the constraints entailed by molecular QTL analyses; it is adapted from that of HESS (Richardson et al., 2010). Suppose that where δ 0 is the Dirac distribution. Each response, y t , is related linearly to the covariates and has a specific precision, τ t . The responses are conditionally independent across the regressions, but dependence among responses associated with the same covariates is captured through the prior specification of parameters ω s and σ −2 , which are shared across the responses. This formulation circumvents modelling the covariance between the responses, which is infeasible when d is large.
Each covariate-response pair has its own regression parameter, β st , for which sparsity is induced using a spike-and-slab prior. The binary parameter γ st acts as a "pair selection" indicator; covariate X s is associated with response y t if and only if γ st = 1. The parameter σ represents the typical size of nonzero effects and is modulated by the residual variance, τ −1 t , of the response concerned by the effect. The parameters γ s1 , . . . , γ sd specify the response(s) associated with X s and are identically distributed as Bernoulli with common parameter ω s . Thus, ω s controls the proportion of responses associated with covariate X s . The goal of inference is variable selection. Selection of predictors can be performed by ranking the posterior means of {ω s } and selection of covariate-response pairs can be performed by ranking the posterior probabilities of inclusion, i.e., the posterior means of {γ st }.
Our model differs from that of Richardson et al. (2010) in two respects. One concerns the treatment of the regression coefficient parameters β st : we use independent priors, whereas Richardson et al. rely on g-priors (Zellner, 1986). The main motivation for our choice is that the effects of genetic variants on a given outcome can be understood as causal, since no retroactive process can affect the variants, and they can take place at locations of the genome that are far apart, so their correlation structure need not reflect the spatial correlation of the SNPs; see Guan and Stephens (2011). Jia and Xu (2007) also rely on independent priors for the regression coefficients of BAYES, but they model the latter with a mixture of two normal distributions rather than a spike-and-slab prior and impose a residual variance parameter that is common to all responses. This stringent assumption may represent a weakness of their proposal.
The second difference concerns the third level of the model. Richardson et al. (2010) opt for a quite complex specification, in which In their case, the inclusion probability of X s for response y t is modelled through ω t ; it is specific to that response but can be regulated using the parameter ρ s , common to all responses. Jia and Xu (2007) and Scott-Boyer et al. (2012) propose other variants for this prior. The former choose a treatment similar to ours, with ω st ≡ ω s ∼ Dirichlet(1, 1), and the latter consider an additional level of hierarchy, with a s ∼ Exp(λ a ) and b s ∼ Exp(λ b ). Our choice ω st ≡ ω s ∼ Beta(a s , b s ) is partly driven by our wish to design a simpler model and partly by practical considerations, since it ensures a closed form for our variational algorithm, unlike with (1). While such a formulation was mentioned by Richardson et al. (2010) and by Scott-Boyer et al. (2012), they did not pursue it because of concerns regarding its ability to control for multiplicity. Indeed, our model inherently enforces sparse associations as the number of responses, d, increases, but no control is achieved when the number of covariates, p, grows. We address this below by providing a procedure to induce a correction through the prior of ω s .
Part of the flexibility of our model comes from the fact that the hyperparameters, a, b ∈ R p + (for ω's Beta prior), λ, ν ∈ R + (for σ −2 's Gamma prior) and η, κ ∈ R d + (for τ 's Gamma prior) are readily interpreted. One option is to set them based on external information regarding the likelihood of given associations, if available. For instance, to favour associations with covariate X s , one can set a s and b s so that the prior proportion of responses affected by X s , E(ω s ), is large.
The use of such assumptions may be very efficient, but it may also skew the inference towards existing knowledge. In the simulations presented in this paper, we assume that the regression and variance parameters are exchangeable, i.e., that all covariates, and responses have the same prior propensity to be involved in associations, by selecting a single value for all components of a, b, η and κ. Without favouring any covariate or response, however, we can control signal sparsity at the level of covariates by specifying (possibly through cross-validation) a prior average number of covariates, p * , expected to be included in the model. Setting the prior probability that X s is associated with at least one response is and simpler models are favoured as p increases. To see this, one can consider the prior odds ratio representing the support for a model to have an additional response associated with X s , i.e., Figure 1: Prior odds ratios, POR(qs − 1 : qs), for qs = 1, . . . , 5, as and bs as in (3), d = 100, p * = 2, and for a total number of covariates ranging from p = 5 to 5, 000; see Scott and Berger (2010) for a similar visualization of prior odds ratios in a single response context. Clearly, penalties arise and increase with the total number of responses in the model, d. Figure 1 displays (4) for q s = 1, . . . , 5 as a function of p and indicates that, when a s and b s are specified as in (3), the penalties also increase with the total number of covariates, p, therefore naturally adjusting for multiplicity. Moreover, the penalties are not uniform when moving from one to two responses associated with X s , or from four to five, for instance.
The experiment reported in Table 1 confirms that adjustment takes place in practice. It considers problems with p 0 = 20 "active" covariates, i.e., those associated with at least one response, and an increasing number of "noise" covariates and it compares the regime with a s and b s set according to (3) to an "uncorrected" regime with a s ≡ 1, b s ≡ 2d − 1, so that the prior mean number of responses associated with X s is 0.5, i.e., E(ω s ) ≡ (2d) −1 . The number of false positives, based on a posterior probability of inclusion greater than 0.5, grows linearly with p when the uncorrected model is used but remains roughly constant close to zero with correction (3), giving a clear multiplicity adjustment. Other experiments confirmed strong sparsity control; the reported findings on real data should therefore be plausible when (3) is used.
Section 2 described some differences between our model and those of Jia and Xu (2007, BAYES), Richardson et al. (2010, HESS) and Scott and Berger (2010, iBMQ), but a more fundamental distinction concerns the inference procedure. The three earlier methods rely on Markov chain Monte Carlo (MCMC) techniques and require massive computing resources when the dimensionality of the problem is large. We instead employ a variational inference procedure, which is deterministic and hence can be much cheaper (Ormerod and Wand, 2010).
Instead of sampling from the joint posterior probability p(θ | y) of the parameter vector of interest, θ, variational approaches proceed by replacing it by a tractable analytical approximation, q(θ). We focus on so-called mean-field variational formulations (Xing et al., 2002;Attias, 2000) to construct such a class, i.e., we assume that q(θ) factorizes over some partition of θ, {θ j } j=1,...,J , no further assumption is made about the distribution, and in particular no constraint is imposed on the functional forms of the q j (θ j ). Here, we consider the factorization and turn the inference into an optimization problem where q(θ) is obtained by minimizing its Kullback-Leibler divergence KL (q p) from the target distribution, p(θ | y). Because the marginal log-likelihood may be written as log p(y) = L(q) + KL (q p) , where minimizing the Kullback-Leibler divergence amounts to maximizing L(q), which represents a lower bound for log p(y). To this end, we observe that where cst is constant with respect to θ j and where we introduced the distribution with E −j {·} denoting the expectation with respect to the distributions q k over all variables θ k , k = j. The right-hand side of (7) corresponds to the negative Kullback-Leibler divergence between q j (θ j ) and p −j (θ j ; y), plus a constant. Hence, assuming that the q k (θ k ), k = j, are fixed, the The relations (8) give rise to cyclic dependencies among the densities q j (θ j ). This suggests an iterative algorithm whose convergence can easily be monitored by evaluating changes in the lower bound L(q). Our choice (5) ensures that the coordinate updates can be derived in closed form; in particular, the semi-conjugacy of our model implies that the prior densities of all parameters are preserved by the variational densities. For instance, a spike-and-slab distribution with modified parameters is recovered at posterior level, st are to be updated iteratively. Convergence is ensured by the convexity of L(q) in each of the q j (θ j ) (Boyd and Vandenberghe, 2004, § § 3.1.5, 3.2.4, 3.2.5). The algorithm and its derivation are given in Appendix B.
4 Empirical quality assessment of the variational approximation 4.1 Tightness of the marginal log-likelihood lower bound. In this section we evaluate the closeness of the variational density q to the target posterior distribution by approximating the Kullback-Leibler divergence KL (q p). Because of relation (6), this amounts to assessing the tightness of the variational lower bound for the marginal log-likelihood, L(q). For small problems, the likelihood p(y) may be accurately approximated using simple Monte Carlo sums. We have see Appendix C for details. As no closed form is available for the remaining integrals, we use where we independently generate  covariates, d = 6 responses and increasing sample sizes, n. In the left panel, the covariates are independent of each other, and so are the responses. In the right panel, the covariates are equicorrelated with correlation coefficient ρ = 0.75, and so are the responses. In both cases, the mean relative difference is below 1% with n = 50 and seems to decrease as n grows. Although we are not aware of any such study with which to benchmark our results, these values seem very small, suggesting that our variational distribution q adequately reflects the target distribution p, at least for small problems. Likewise, the variational lower bound L(q) may be used as a proxy for the marginal log-likelihood when performing model selection; this use will be illustrated in Section 5.3.
The fact that the variational lower bound remains tight in the correlated data case is reassuring, as it suggests that the independence assumptions underlying the mean-field factorization of q may only weakly impact the quality of the approximation.

Comparison with Markov Chain Monte Carlo.
We complement our quality assessment by comparing several variational posterior quantities with those for MCMC inference on problems of moderate size. A fair comparison is not straightforward, as these two types of inference rely on stopping rules and convergence diagnostics of very different natures. While the convergence criterion for variational inference comes down to a tolerance to be prescribed, the ability of MCMC sampling to adequately explore the model space for a given chain length can be difficult to evaluate, and usually varies greatly with the problem size. To alleviate the risk of inaccurate MCMC inference, we run 10 5 iterations and discard the first half. We also support our comparison with selected quantities approximated by simple Monte Carlo sums, namely, the posterior probability of inclusion of a covariate X s for a response y t , and the posterior mean of ω s , controlling the proportion of responses associated with covariate X s , with the samples σ −2 (i) and {ω (i) s } generated as in (9), with I = 2 × 10 5 draws. Table 2 reports the variational, MCMC and simple Monte Carlo estimates of β and ω for a problem with p = 8 covariates, d = 5 responses for n = 250 samples, and with each nonzero association explaining on average 13.5% of response variance. The estimates all agree closely. Those of the five active regression coefficients, β 1,2 , β 2,1 , β 3,2 , β 4,1 and β 4,2 , are significantly different from zero, unlike the average estimate of the inactive coefficients. A plot of the MCMC and variational posterior densities (the latter obtained in closed form), given in Figure  Both experiments indicate equivalent performance for MCMC and variational inferences.
5 Statistical performance 5.1 Predictor selection. The problems considered in Section 4.2 were small enough to allow accurate and tractable MCMC inference. In this section, we assess the performance of our approach on larger problems by comparing it to popular variable selection methods; i.e., with joint modelling of outcomes and covariates (elastic net for multivariate Gaussian responses), with joint modelling of covariates only (Bayesian multiple regression based on MCMC inference, "BAS", or variational inference, "varbvs"), or with fully marginal modelling (univariate ordinary least squares and "lmBF" Bayesian regressions). Complete descriptions and references are in Appendix D.2. The methods are compared by measuring their ability to detect the active covariates, i.e., to determine which covariates are associated with at least one response. For our variational approach, this task is achieved by ranking the posterior means of the ω s , which control the proportion of responses associated with a given covariate.
Our data-generation design is based on generally accepted principles of population genetics. We simulate SNPs under Hardy-Weinberg equilibrium from a binomial distribution with probabilities corresponding to minor allele frequencies of common variants, chosen in the interval (0.05, 0.5) uniformly at random, and we generate outcomes from Gaussian distributions with specific error variances. The dependence structure of these variables is either enforced block-wise with preselected auto-or equicorrelation coefficients or chosen to be that of real data. The labels of the active SNPs and outcomes are picked randomly, and each active SNP is associated to one (randomly selected) active outcome and to each of the remaining active outcomes with a prescribed probability; some outcomes are therefore under pleiotropic control. The proportions of outcome variance explained per SNP are simulated for all associations from a positively skewed Beta distribution to favour the generation of smaller effects, and they are then rescaled to match a given average proportion. To mimic the result of natural selection, the effect sizes are inversely related to the SNP minor allele frequencies. For more details, see Appendix D.1.
We perform 48 replications for each of three simulation configurations. The first configuration has moderate numbers of covariates (p = 5, 000) and outcomes (d = 50), and allows time-consuming methods to run within hours. The second has many outcomes (d = 20, 000) and the third has many covariates (p = 150, 000); these numbers approach those encountered in molecular QTL studies.
The remaining settings (numbers of active outcomes and covariates, of observations, effect sizes, etc) are detailed in the caption to Figure 3.
The ROC curves in Figure 3 indicate that our approach outperforms the other methods. It is appreciably more powerful for low false positive rates, which are of particular interest for the highly sparse scenarios typically expected for genome-wide association studies. Despite the correlation among the covariates and the outcomes, our method does not seem to suffer from the independence assumptions implied by the mean-field approximation, as suggested by the results of Section 4.1. The marginal ordinary least squares and marginal lmBF regressions appear to miss many associations because of their univariate modelling of covariates, but jointly accounting for the covariates may not suffice, as suggested by the rather poor performances of the Bayesian multiple regression approaches, BAS and varbvs, which apply separate multiple linear regressions for each outcome.
It appears that the ability of our approach to exploit the similarity across outcomes yields more power to detect their shared associations. Finally, even though the multivariate elastic net models jointly the covariate and outcome variables, its inference suffers from the assumption that to each covariate corresponds a single regression coefficient, shared for all responses. As a consequence, regression estimates of covariates with weak or few associations with the responses may be shrunk to zero.

Combined selection of predictors and outcomes.
Unlike the classical variable selection methods used as comparators in Section 5.1, our approach and those of Richardson et al. (2010), HESS, and Scott-Boyer et al. (2012), iBMQ, are tailored to molecular QTL problems: they quantify the associations between each covariate-response pair in a single model, and thus provide flexible and unified frameworks for detecting pairs of associated SNP-molecules, as well as pleiotropic SNPs associated with many molecular outcomes. In this section, we compare the three approaches in terms of the posterior quantities used to perform such selection. As both HESS and iBMQ rely on MCMC sampling, we consider smaller problems than in Section 5.1 in order to ensure convergence within a reasonable time. The simulated datasets have p = 250 covariates, of which p 0 = 50 are active, and d = 100 outcomes, of which d 0 = 50 are active, the probability of association being     Table 3, which gathers sensitivity and specificity measures based on median probability models (Barbieri and Berger, 2004) (consisting of those covariate-response pairs whose posterior inclusion probability is higher than 0.5). As discussed in Section 2, control of signal sparsity can be induced through the prior for ω s , still, rather than median probability models, one may prefer to use a data-driven false discovery threshold in order to prescribe a desired level of false discoveries.   interest, weight maintenance. We also use this illustration on real data to further highlight the benefits of modelling the outcomes jointly via an extensive permutation-based comparison with the single-response variational method varbvs (Carbonetto and Stephens, 2012). In order to adjust for multiplicity, we specify the hyperparameters for ω according to the discussion of Section 2 and choose the prior average number of active SNPs, p * , by grid search within a 3-fold cross-validation procedure that maximizes the variational lower bound. After hyperparameter selection, the algorithm converged in 83 iterations, taking about 10 hours on an Intel Xeon CPU at 2.60 GHz with 512 GB RAM. The posterior means E q (ω s | y) suggest the presence of several active SNPs, spread across the chromosomes (Figure 6), but we use the marginal posterior inclusion probabilities, E q (γ st | y), to declare pairwise associations and active SNPs.
We compare varbvs and our method on real data based on the number of associations declared by each method at specific false discovery rates estimated by permutations. We apply Efron's Bayesian interpretation of the false discovery rate (Efron, 2008) to posterior probabilities of inclusion, and use an empirical null distribution based on B = 400 permutations to compute the estimate FDR(τ ) = median b=1,...,B #{PPI for a grid of thresholds τ ; we then fit a cubic spline to the resulting false discovery rates to find thresholds for specific rates. The analysis suggests that our method is more powerful, with 89 associations declared at an estimated FDR of 25%, against 47 for varbvs; the superiority of our method is further highlighted by Table 4. Figure 6 displays the associations declared by our method and the overlap with those declared by varbvs at estimated FDR of 25%; the associations detected also largely agree with those obtained with marginal screening at Benjamini-Hochberg FDR of 25%.
Database searches on the functional relevance of the detected associations give hints of promising biological functions related to metabolic activities for 12 of the 25 SNPs declared as active by our procedure. For instance, the most outstanding SNP Figure 6, rs4316911, shows many associations with triglyceride levels, and turns out to be located less than 150kb from the protein coding gene ITGA6 known to be linked to diabetic kidney disease (Iyengar et al., 2015). The second most prominent pleiotropic SNP, rs174535, is declared by our approach to be associated with phospholipids, more precisely with 14 different phosphatidylcholine levels, of which four are ether-linked/plasmalogen (PC-O). Interestingly, this latter SNP has been recently reported to be related to metabolite levels; among others, it was found to be associated with trans fatty acid levels and plasma phospholipid levels (Mozaffarian et al., 2015), in line with our findings. Moreover, it was found to be an eQTL for the fatty acid desaturase genes FADS1 and FADS2. The SNP rs3903703 too has been identified as associated with very long-chain fatty acid levels . This seems to agree with our findings, in which rs3903703 exhibits associations with sphingomyelin, a type of lipid containing fatty acids of different chain lengths. The complete subset of SNPs with metabolism-related links found by our procedure is given in Appendix E. Additional details on this real data study, as well as on its replication using simulated data, are also provided there.

Conclusions
We have described a scalable and efficient approach to joint variable selection from large numbers of candidate predictor and outcome variables. As it exploits the similarity across outcomes through a flexible hierarchical structure, our procedure outperforms the most popular predictor selection approaches in high-dimensional set-ups. The variational approximation on which our approach relies provides accurate posterior quantities, with reduced computational effort relative to MCMC procedures; in particular, it yields inferences comparable to those of the MCMC procedure HESS (Richardson et al., 2010). Convergence control is automatic, whereas convergence assess-ment for MCMC algorithms can be difficult, especially in high dimensions. Our simulations also show that variable selection remains powerful when the predictors and outcomes are correlated, notwithstanding the independence assumptions underlying the mean-field factorization.
The key added-value of our approach is its applicability to molecular QTL datasets without the need for prior dimension reduction. In an application, our approach recovered several previously reported SNP-metabolite associations, and declared more associations than the single-outcome method "varbvs" (Carbonetto and Stephens, 2012) at prescribed false discovery rates, thus highlighting the benefits of jointly modelling the outcomes. To the best of our knowledge, no competing Bayesian approach for joint inference on two high-dimensional sets of variables can deal with the problem sizes typically encountered in molecular QTL analyses.
Bioinformatics is moving towards whole-genome analyses, for which several million genetic variants need to be considered, so it seems worthwhile to consider further speed-up strategies for our approach. One option is to use new optimization procedures. So-called natural gradient methods, which rely on the Riemannian structure of variational approximate distributions, seem particularly attractive, as they can be orders of magnitude faster than conventional gradient algorithms (Honkela et al., 2008). Another possibility is a "split-and-merge" strategy, i.e., first partitioning the

Software
The algorithm and the data-generation functions used in this paper are implemented in the publicly available R package locus.

A Multiplicity control at covariate level
Sparsity control at covariate level can be induced through the prior distribution of ω, by carefully selecting its hyperparameters. The prior probability that X s is "active" (i.e., associated with at least one response) is and this after a little algebra equals , so assuming exchangeability and setting implies that where p * is interpreted as a prior average number of active covariates. The choice (11) yields a multiplicity adjustment as suggested by a plot (Figure 1) of the prior odds ratios, indicating the penalty induced by the prior when moving from q s − 1 to q s responses associated with X s ,

B Derivation of the variational algorithm
B.1 Variational distributions. We provide the detailed derivation of our variational algorithm, which is given in Appendix B.3. We have where Let θ = β, γ, τ, σ −2 , ω and consider the following mean-field form for the variational approximation, We obtain each component of this factorization using the formula log q j (θ j ) = E −j {log p(y, θ)} + cst, j = 1, . . . , J, with p(y, θ) given in (12), where E −j is the expectation with respect to the distribution q k over all variables θ k (k = j) and where cst is a constant with respect to θ j . Writing θ (r) j for the r th moment with respect to the approximate posterior distribution q j of θ j , we have where cst is a constant with respect to β st and γ st . Completing the square yields .
We therefore observe that We now compute the variational approximate distribution for the error variance of each y t : Therefore we have where Since τ t has a Gamma distribution, the expectation E (log τ t ) appearing in (13) can be rewritten in terms of η * t and κ * t using the digamma function (Abramowitz and Stegun, 1965), We also find that and, as before, we now have Finally, we have and As ω s has a Beta distribution, we also get (Abramowitz and Stegun, 1965) B.2 Lower bound of the marginal log-likelihood. We provide the computational details for the lower bound, L(q), of the marginal log-likelihood, log p(y). It is evaluated at each iteration of our algorithm, in order to monitor its convergence: where E q (log τ t ), E q log σ −2 , E q (log ω s ) and E q {log(1 − ω s )} are given by (14), (15) and (16).

B.3 Variational algorithm.
Algorithm 1 inputs: y (centered), X (standardized using the usual unbiased estimator of the variance), a, b, η, κ, λ, ν, tol, maxit initialize: for s = 1, . . . , p do for t = 1, . . . , d do 2) based on the current parameter updates until |L(q) − L old (q)| < tol or it = maxit The symbols and are the Hadamard operators standing for element-wise multiplication and division of two matrices of the same dimension.
B.4 Computational details. The convergence characteristics of our procedure can be described in terms of those of any deterministic iterative algorithm. In our experiments, we set the tolerance for the stopping criterion to 10 −6 , as our empirical tests suggest that when using a smaller tolerance, the additional time required until convergence does not yield noticeably better inferences.
While MCMC sampling may require thousands of iterations to converge, our algorithm usually converges in tens of iterations. As suggested by the runtime profiling provided in Appendix D.3, inference for typical genome-wide association problems with multiple outcomes is usually completed in hours. Our algorithm requires the initialization of the variational parameters, but unfortunately comes with no guarantee that it will attain a global minimum for the Kullback-Leibler divergence, KL (q p). This drawback can be alleviated by using several different initializations, at the price of increasing the computational effort. In practice we did not encounter situations where different starting points gave different optima. The source code can be found in the publicly available R package locus.
C Details on the empirical quality assessment of the variational approximation C.1 Marginal likelihood computation. We have and one can obtain a closed form expression for p y t | γ t , σ −2 , after integrating out β t and then τ t . Indeed, proceeding similarly as in George and McCulloch (1997), Hence, if q γt = 0, If q γt = 0, then  each explains on average 13.5% of response variance. We use the software OpenBUGS (Spiegelhalter et al., 2007) and the R package coda (Plummer et al., 2006) for the MCMC inference and convergence diagnostics.   D Details on the simulation studies D.1 Data generation design. We provide here some implementation details on the datageneration settings described in Section 5.1. The same general design is used for all the numerical experiments; it is tailored to genetic association studies with multiple outcomes. The dependence structure of simulated SNPs and molecular outcomes is by blocks. As we assume Hardy-Weinberg equilibrium for the SNPs, their marginal distribution is a binomial B(2, m) with probability m equal to their respective minor allele frequencies, so we model dependence within each block using realisations from a multivariate Gaussian latent variable whose correlation matrix describes a desired dependence structure (autocorrelation with prescribed correlation coefficient or correlation structure corresponding to that of real SNPs at disposal). SNPs are then obtained using a quantile thresholding rule involving their preselected minor allele frequencies; this approach is along the lines of copula-based dependence modelling (without having to resort to copulas). When simulating SNPs that emulate real SNPs data by approximating their minor allele frequencies and correlation structure, the empirical covariance of real SNP blocks may not be positive definite, in which case we approximate it by the closest positive definite covariance matrix (in Frobenius norm) using the algorithm of Higham (2002) implemented in the R package Matrix (Bates and Maechler, 2015). Outcomes are associated with SNPs under an additive dose-effect scheme (uniform and linear increase in risk for each copy of the minor allele) and the proportions of outcome variance explained per active SNP are drawn from a Beta(α, β) distribution, with shape parameters α = 2 and β = 5 chosen to give more weight to smaller effect sizes. These proportions are rescaled to match a preselected average proportion of outcome variance explained per active SNP (assuming that genetic and external environmental effects are uncorrelated). The magnitude of SNP effects derives from these choices and the sign of these effects is switched with probability 0.5. Such a design implies an inverse relationship between minor allele frequencies and effect sizes, which is expected to occur since selection against SNPs with large penetrance is stronger (see, e.g., Park et al., 2011). Additional information regarding the generation of a pleiotropic association pattern is given in the paper, as it may vary with the test case considered.
D.2 Competing predictor selection methods. The performance of our approach in terms of predictor selection is compared with the following regression procedures: • univariate ordinary least squares: each response y t is regressed on each covariate X s and the statistics max t=1,...,d t st , where t st is the t-statistic for the significance of β st , are gathered and ranked; • elastic net regression for multivariate Gaussian responses (α = 0.5) with 10-fold crossvalidation for choosing the tuning parameter λ (glmnet, Friedman et al., 2009). The estimates |β s | (s = 1, . . . , p), where β s is the regression coefficient estimate for X s and common to all responses, are gathered and ranked; • univariate Bayesian regressions, lmBF (Morey and Rouder, 2015): each response y t is regressed on each covariate X s with all computations made analytically. The (average) Bayes factors, d t=1 BF st /d (s = 1, . . . , p), are gathered and ranked; • d Bayesian multiple regressions, BAS (Clyde, 2016), one for each response, using MCMC inference. A g-prior is used for the regression coefficients. The (average) Bayes factors, d t=1 BF st /d (s = 1, . . . , p), are gathered and ranked; and • d Bayesian multiple regressions, varbvs (Carbonetto and Stephens, 2012), one for each response, using variational inference. The posterior probabilities of inclusion are summed across responses and ranked.
D.3 Runtime profiling. We report a runtime profiling of the different methods for a range of problem sizes (n × p × d). Figure 11 displays the run times in minutes, averaged over 24 replications. We do not aim to provide precise and exhaustive comparisons, since the methods all depend on parallelism and convergence characteristics that are not directly comparable. All methods were run serially, in an attempt to treat them on an equal footing. This choice could be challenged, as some approaches are more parallelizable than others, but the number of cores used for the latter represents an additional setting that would come into play with a potentially large impact on the measures. The number of chains for MCMC inferences also matters: HESS (Richardson et al., 2010) is run with three chains (following its authors' choice made in their simulations); the other MCMC inferences are based on a single chain. Finally, the runtime may also greatly vary depending on the chosen chain length: the latter is adaptively selected by the Bayesian multiple regression method BAS (Clyde, 2016), and, based on preliminary convergence diagnostics, it is set to 50, 000 samples for HESS, iBMQ (Scott-Boyer et al., 2012) and the MCMC inference on our model using OpenBUGS (Spiegelhalter et al., 2007). In practice, the number of samples needed until convergence may increase greatly with the problem size, a fact that was not accounted for in this profiling. Hence, by adequately increasing the chain lengths when considering larger dimensionalities, we expect the curves of Figure 11 corresponding to the MCMC approaches to deviate more widely from that of our method. With these serial settings, our approach is faster than all Bayesian and frequentist methods. MCMC inference for our model is the slowest of all tested methods, which underlines the intractability of MCMC sampling for large problems.
Inference for HESS is faster, but still more than 650 times slower than our variational approach.
Our method is also about 10 times faster than d applications (one for each outcome) of the varbvs method (Carbonetto and Stephens, 2012) but the runtime of the varbvs procedure can be reduced using multiple cores. E Details on the real data problem E.1 Clinical study and preprocessing of the mQTL data analysis. Diogenes (Larsen et al., 2010) is a large clinical dietary intervention study, which enrolled roughly 1, 000 overweight subjects from 8 European centers. The subjects were assigned to a 8-week weight-loss program followed by 6 months of supervised ad-libitum diet, meant as a weight-maintenance phase. For the 6-month period, individuals were randomized into five intervention groups, whose diets differed in their macronutrient and glycemic load content. The clinical trial had the broad objective of elucidating whether certain macronutrient compositions are more successful in weight maintenance than others. It gathered data on genetic variants, gene, protein and metabolite expression. In addition, more than 7, 000 clinical, anthropometric and behavioural variables were made available for each individual. When meaningful to do so, the data were collected at three time points, reflecting different stages of the dietary treatment.
In Section 5.3, we performed an mQTL analysis which involves tag SNPs from Illumina Human-Core chips (about 300k SNPs) and metabolites from plasma obtained by liquid chromatographymass spectrometry (LC-MS). We excluded SNPs having missing values, minor allele frequency below 5%, call rate below 95%, or violating the Hardy-Weinberg equilibrium. Moreover, subjects having gender discrepancy (e.g., subject recorded as male but being homozygous for each X chromosome marker), abnormal autosomal heterozygosity or whose genomes were too close to each other (IBS> 95%) were also excluded. Additional subjects were removed after applying the Tukey method for outlier detection (Tukey, 1949), based on the metabolomic data. Metabolite expression levels were log2-transformed, and had no missing values. After these quality checks, the data consisted of p = 215 907 SNPs and d = 125 metabolite expression levels, for n = 317 individuals.
E.2 Permutation-based Bayesian false discovery rate estimation. We provide details on the false discovery rate estimation procedure applied in Section 5.3 to compare our method with the varbvs method of Carbonetto and Stephens (2012) on the real data. We use the two-group mixture approach proposed by Efron (2008) in the context of microarray data analysis: we simultaneously consider N null hypotheses and their corresponding test statistics, which we assume to follow a mixture distribution where π 0 is the prior probability for a null case, and F 0 and F 1 are the null and non-null cumulative distribution functions. The "Bayesian false discovery rate" for some threshold τ is the posterior probability that a rejected null hypothesis (i.e., test statistic exceeding τ ) is a false positive, whereF = 1 − F andF 0 = 1 − F 0 . We derive estimates of (20) based on the posterior probabilities of inclusion quantifying the associations between each covariate-response pair. Specifically, we obtain an empirical null distribution by running our algorithm (with the same hyperparameters as those used for the actual inference) on B datasets with permuted outcome sample labels and compute, for a grid of thresholds 0 < τ 1 < . . . < τ K < 1, FDR(τ k ) = median b=1,...,B #{PPI (b) st > τ k ; s = 1, . . . , p; t = 1, . . . , d} #{PPI st > τ k ; s = 1, . . . , p; t = 1, . . . , d} , k = 1, . . . , K, (21) where we conservatively set π 0 to 1 in (20). As the posterior probabilities of inclusion obtained by applying varbvs d times (one multiple regression for each outcome) are not identically calibrated across all d estimations, we use adaptive thresholds on the columns of the varbvs posterior probabilities of inclusion matrix, τ k (y t ) = median s (PPI st ) median s,t (PPI st ) τ k t = 1, . . . , d, k = 1, . . . , K.
To find thresholdsτ corresponding to preselected false discovery rates, we fit a cubic spline to the estimates (21) previously obtained for τ 1 , . . . , τ K .
E.3 Biological evidence for the mQTL analysis findings. We support the findings obtained by our method for the mQTL data analysis with external association results from the following online databases: • GWAS catalog (Welter et al., 2014); • UCSC genome browser (Karolchik et al., 2003); • GTEx (Lonsdale et al., 2013); and • GeneCards (Rebhan et al., 1998). We find direct or indirect links with metabolic activities for 12 of the 25 SNPs declared as "active' by our method: SNPs rs3820711, rs4316911, rs4909818, rs4744227, rs174535, rs680379, rs8012466, rs4906771, rs573922, rs3903703, rs8114788 and rs6001093 have been identified as associated with BMI, diverse diabetic or obesity diseases, fatty acid, sphingolipid or phospholipid levels, either from direct genome-wide association analyses or through protein coding genes for which they were reported as eQTLs.
E.4 Replication of the mQTL data analysis. In this appendix, we replicate the real data analysis of Section 5.3 on a simulated dataset (with twice as many outcomes and slightly more observations), which can be found on Figshare (https://dx.doi.org/10.6084/m9.figshare. 4509755.v1). The analysis can be reproduced using the code available on GitHub (https:// github.com/hruffieux/mQTL_analysis_example), provided that a machine with adequate RAM memory (about 400G for this large problem) is used.
To best mimic the real mQTL data used in Section 5.3, we simulate p = 215, 907 SNPs based on the sample minor allele frequencies of the real tag SNPs and we reproduce their dependence structure by blocks of 1, 000 consecutive SNPs according to the discussion of Appendix D.1. We also simulate d = 250 normally distributed outcomes with equicorrelation by blocks using correlation coefficient ρ ∈ {0.5, 0.6, 0.7}; this block structure is similar to that of the real metabolite data. To induce a realistic pleiotropic pattern, we simulate associations between 750 SNPs and 175 metabolites (randomly chosen) taking the block-wise dependence structure of the latter into account: the probability that a given active SNP is associated with a given active outcome varies across blocks, so correlated metabolites are either all highly likely or all less likely to be associated with the SNP. The average proportion of metabolic variance explained with each association is 2.5% (but more associations explain less than this, as discussed in Appendix D.1). We generate n = 350 observations.
We apply the permutation analysis described in Section 5.3 to varbvs and our method and again obtain a more powerful selection for our method compared to varbvs at estimated FDR of 20% and 25% (Table 5 and Figure 13). Because simulated data are used, we can further support the overall superiority of our method with ROC curves, see Figure 12.  Figure 13: Comparison of the associations declared by VB (left) and varbvs (right) at FDR of 20% estimated using B = 200 permutations for simulated mQTL data; see also Table 5. The simulated association pattern is overlaid (red crosses). The large number of false positives on the right plot indicates that the permutation-based FDR estimates are somewhat underestimated for varbvs; an improved estimation based on more permutations would further emphasize the improvement of our method over varbvs. Figure 14: Manhattan plots of SNP association and evidence of pleiotropy obtained with VB (left) and varbvs (right) for simulated mQTL data.