Polygenic Modeling with Bayesian Sparse Linear Mixed Models

Both linear mixed models (LMMs) and sparse regression models are widely used in genetics applications, including, recently, polygenic modeling in genome-wide association studies. These two approaches make very different assumptions, so are expected to perform well in different situations. However, in practice, for a given dataset one typically does not know which assumptions will be more accurate. Motivated by this, we consider a hybrid of the two, which we refer to as a “Bayesian sparse linear mixed model” (BSLMM) that includes both these models as special cases. We address several key computational and statistical issues that arise when applying BSLMM, including appropriate prior specification for the hyper-parameters and a novel Markov chain Monte Carlo algorithm for posterior inference. We apply BSLMM and compare it with other methods for two polygenic modeling applications: estimating the proportion of variance in phenotypes explained (PVE) by available genotypes, and phenotype (or breeding value) prediction. For PVE estimation, we demonstrate that BSLMM combines the advantages of both standard LMMs and sparse regression modeling. For phenotype prediction it considerably outperforms either of the other two methods, as well as several other large-scale regression methods previously suggested for this problem. Software implementing our method is freely available from http://stephenslab.uchicago.edu/software.html.

Despite the considerable overlap in their applications, in the context of polygenic modeling, LMMs and sparse regression models are based on almost diametrically opposed assumptions.Precisely, applications of LMMs to polygenic modeling (e.g.[22]) effectively assume that every genetic variant affects the phenotype, with effect sizes normally distributed, whereas sparse regression models, such as Bayesian variable selection regression models (BVSR) [18,19], assume that a relatively small proportion of all variants affect the phenotype.The relative performance of these two models for polygenic modeling applications would therefore be expected to vary depending on the true underlying genetic architecture of the phenotype.However, in practice, one does not know the true genetic architecture, so it is unclear which of the two models to prefer.Motivated by this observation, we consider a hybrid of these two models, which we refer to as the "Bayesian sparse linear mixed model", or BSLMM.This hybrid includes both the LMM and a sparse regression model, BVSR, as special cases, and is to some extent capable of learning the genetic architecture from the data, yielding good performance across a wide range of scenarios.By being "adaptive" to the data in this way, our approach obviates the need to choose one model over the other, and attempts to combine the benefits of both.
The idea of a hybrid between LMM and sparse regression models is, in itself, not new.Indeed, models like these have been used in breeding value prediction to assist genomic selection in animal and plant breeding programs [30][31][32][33][34][35], gene selection in expression analysis while controlling for batch effects [36], phenotype prediction of complex traits in model organisms and dairy cattle [37][38][39][40], and more recently, mapping complex traits by jointly modeling all SNPs in structured populations [41].Compared with these previous papers, our work makes two key contributions.First, we consider in detail the specification of appropriate prior distributions for the hyper-parameters of the model.We particularly emphasize the benefits of estimating the hyper-parameters from the data, rather than fixing them to pre-specified values to achieve the adaptive behavior mentioned above.Second, we provide a novel computational algorithm that exploits a recently described linear algebra trick for LMMs [8,9].The resulting algorithm avoids ad hoc approximations that are sometimes made when fitting these types of model (e.g.[37,41] ), and yields reliable results for data sets containing thousands of individuals and hundreds of thousands of markers.(Most previous applications of this kind of model involved much smaller data sets.) Since BSLMM is a hybrid of two widely used models, it naturally has a wide range of potential uses.Here we focus on its application to polygenic modeling for genome-wide association studies, specifically two applications of particular interest and importance: PVE estimation (or "chip heritability" estimation) and phenotype prediction.Estimating the PVE from large-scale genotyped marker sets (e.g.all the SNPs on a genome-wide genotyping chip) has the potential to shed light on sources of "missing heritability" [42] and the underlying genetic architecture of common diseases [19,[22][23][24]43].And accurate prediction of complex phenotypes from genotypes could ultimately impact many areas of genetics, including applications in animal breeding, medicine and forensics [27][28][29][37][38][39][40].Through simulations and applications to real data, we show that BSLMM successfully combines the advantages of both LMMs and sparse regression, is robust to a variety of settings in PVE estimation, and outperforms both models, and several related models, in phenotype prediction.

Background and Motivation
We begin by considering a simple linear model relating phenotypes y to genotypes X: ǫ ∼ MVN n (0, I n τ −1 ). ( Here y is an n-vector of phenotypes measured on n individuals, X is an n × p matrix of genotypes measured on the same individuals at p genetic markers, β is a p-vector of (unknown) genetic marker effects, 1 n is an n-vector of 1s, µ is a scalar representing the phenotype mean, and ǫ is an n-vector of error terms that have variance τ −1 .MVN n denotes the n-dimensional multivariate normal distribution.Note that there are many ways in which measured genotypes can be encoded in the matrix X.We assume throughout this paper that the genotypes are coded as 0, 1 or 2 copies of a reference allele at each marker, and that the columns of X are centered but not standardized; see Text S1 Detailed Methods.
A key issue is that, in typical current datasets (e.g.GWAS), the number of markers p is much larger than the number of individuals n.As a result, parameters of interest (e.g.β or PVE) cannot be estimated accurately without making some kind of modeling assumptions.Indeed, many existing approaches to polygenic modeling can be derived from (1) by making specific assumptions about the genetic effects β.For example, the LMM approach from [22], which has recently become commonly used for PVE estimation (e.g.[24,[44][45][46]), is equivalent to the assumption that effect sizes are normally distributed, such that [Specifically, exact equivalence holds when the relatedness matrix in the LMM is computed from the genotypes as K = 1 p XX T (e.g.[47]).[22] use a matrix in this form, with X centered and standardized, and with a slight modification of the diagonal elements.]For brevity, in this paper we refer to the regression model that results from this assumption as the LMM (note that this is relatively restrictive compared with the usual definition); it is also commonly referred to as "ridge regression" in statistics [48].The estimated combined effects (Xβ), or equivalently, the estimated random effects, obtained from this model are commonly referred to as Best Linear Unbiased Predictors (BLUP) [49].
An alternative assumption, which has also been widely used in polygenic modeling applications [18,19,34], and more generally in statistics for sparse high-dimensional regression with large numbers of covariates [50,51], is that the effects come from a mixture of a normal distribution and a point mass at 0, also known as a point-normal distribution: where π is the proportion of non-zero β and δ 0 denotes a point mass at zero.[This definition of π follows the convention from statistics [19,50,51], which is opposite to the convention in animal breeding [27,[32][33][34]40].]We refer to the resulting regression model as Bayesian Variable Selection Regression (BVSR), because it is commonly used to select the relevant variables (i.e.those with non-zero effect) for phenotype prediction.Although (4) formally includes (3) as a special case when π = 1, in practice ( 4) is often used together with an assumption that only a small proportion of the variables are likely to be relevant for phenotype prediction, say by specifying a prior distribution for π that puts appreciable mass on small values (e.g.[19]).In this case, BVSR and LMM can be viewed as making almost diametrically opposed assumptions: the LMM assumes every variant has an effect, whereas BVSR assumes that a very small proportion of variants have an effect.(In practice, the estimate of σ b under LMM is often smaller than the estimate of σ a under BVSR, so we can interpret the LMM as assuming a large number of small effects, and BVSR as assuming a small number of larger effects.) A more general assumption, which includes both the above as special cases, is that the effects come from a mixture of two normal distributions: Setting π = 0 yields the LMM (3), and σ b = 0 yields BVSR (4).we can interpret this model as assuming that all variants have at least a small effect, which are normally distributed with variance σ 2 b /(pτ ), and some proportion (π) of variants have an additional effect, normally distributed with variance σ 2 a /(pτ ).The earliest use of a mixture of two normal distributions for the regression coefficients that we are aware of is [52], although in that paper various hyper-parameters were fixed, and so it did not include LMM and BVSR as special cases.
Of the three assumptions on the effect size distributions above, the last ( 5) is clearly the most flexible since it includes the others as special cases.Obviously other assumptions are possible, some still more flexible than the mixture of two normals: for example, a mixture of three or more normals.Indeed, many other assumptions have been proposed, including variants in which a normal distribution is replaced by a t distribution.These variants include the "Bayesian alphabet models" -so-called simply because they have been given names such as BayesA, BayesB, BayesC, etc. -that have been proposed for polygenic modeling, particularly breeding value prediction in genomic selection studies.Table 1 summarizes these models, and some other effect size distributions that have been proposed, together with relevant references (see also [53] and the references there in).Among these, the models most closely related to ours are BayesCπ [34] and BayesR [35].Specifically, BayesCπ without a random effect is BVSR, and with a random effect is BSLMM (which we define below).BayesR models effect sizes using a mixture of three normal components plus a point mass at zero, although the relative variance for each normal distribution is fixed.
Given the wide range of assumptions for effect size distributions that have been proposed, it is natural to wonder which are the most appropriate for general use.However, answering this question is complicated by the fact that even given the effect size distribution there are a number of different ways that these models can be implemented in practice, both in terms of statistical issues, such as treatment of the hyper-parameters, and in terms of computational and algorithmic issues.Both these types of issues can greatly affect practical performance.For example, many approaches fix the hyper-parameters to specific values [27,32,33,40] which makes them less flexible [34,54].Thus, in this paper we focus on a particular effect size distribution (5), which while not the most flexible among all that could be considered, is certainly more flexible than the one that has been most used in practice for estimating PVE (LMM), and admits certain computational methods that could not be applied in all cases.We consider in detail how to apply this model in practice, and the resulting advantages over LMM and BVSR (although we also compare with some other existing approaches).A key contribution of this paper is to provide new approaches to address two important practical issues: the statistical issue of how to deal with the unknown hyper-parameters (π, σ a , σ b ), and the computational issue of how to fit the model.Notably, with the computational tools we use here, fitting the model ( 5) becomes, for a typical data set, less computationally intensive than fitting BVSR, as well as providing gains in performance.
With this background, we now turn to detailed description of the model, its prior specification and its computation algorithm.

The Bayesian Sparse Linear Mixed Model
In this paper we focus on the simple linear model (1) with mixture prior (5) on the effects.However, the computational and statistical methods we use here also apply to a more general model, which we refer to as the Bayesian Sparse Linear Mixed Model (BSLMM), and which includes the model ( 1) with (5) as a special case.
The BSLMM consists of a standard linear mixed model, with one random effect term, and with sparsity inducing priors on the regression coefficients: where u is an n-vector of random effects with known n × n covariance matrix K.In referring to u as the "random effects" we are following standard terminology from LMMs.Standard terminology also refers to the coefficients β as "fixed effects", but this phrase has a range of potential meanings [55] and so we avoid it here.Instead we use the term "sparse effects" for these parameters to emphasize the sparsity-inducing prior.
It is straightforward to show that when K = 1 p XX T , BSLMM is equivalent to the simple linear model (1) with mixture prior (5) on the effects.However, our discussion of prior specification, computational algorithms, and software, all apply for any K.
When we say that ( 6) is equivalent to (1) with ( 5), this equivalence refers to the implied probability model for y given X and the hyper-parameters µ, τ, π, σ a , σ b .However, β and β are not equivalent (explaining our use of two different symbols): in (6) the random effect u captures the combined small effects of all markers, whereas in (1) these small effects are included in β.Since both our applications focus on the relationship between y and X, and not on interpreting estimates of β or β, we do not concern ourselves any further with this issue, although it may need consideration in applications where individual estimated genetic effects are of more direct interest (e.g.genetic association mapping).A related issue is the interpretation of the random effect u in BSLMM: from the way we have presented the material u is most naturally interpreted as representing a polygenic component, specifically the combined effect of a large number of small effects across all measured markers.However, if there are environmental factors that influence the phenotype and are correlated with genotype (e.g.due to population structure), then these would certainly affect estimates of u, and consequently also affect estimates of other quantities, including the PVE.In addition, phenotype predictions from BSLMM will include a component due to unmeasured environmental factors that are correlated with measured genotypes.These issues are, of course, not unique to BSLMM -indeed, they apply equally to the LMM; see [56] and the response from [57] for relevant discussion.
Finally, given the observation that a mixture of two normals is more flexible than a point-normal, it might seem natural to consider modifying (6) by making the assumption that β comes from a mixture of two normal distributions rather than a point-normal.However, if K = 1 p XX T then this modification is simply equivalent to changing the values of σ a , σ b .

Prior Specification
The BSLMM involves (hyper-)parameters, µ, τ, π, σ a , and σ b .Before considering prior specification for these parameters, we summarize their interpretations as follows: • µ and τ −1 control the phenotype mean and residual variance.
• π controls the proportion of β values in (6) that are non-zero.
• σ a controls the expected magnitude of the (non-zero) β.
• σ b controls the expected magnitude of the random effects u.
Appropriate values for these parameters will clearly vary for different data sets, so it seems desirable to estimate them from the data.Here we accomplish this in a Bayesian framework by specifying prior distributions for the parameters, and using Markov chain Monte Carlo (MCMC) to obtain approximate samples from their posterior distribution given the observed data.Although one could imagine instead using maximum likelihood to estimate the parameters, the Bayesian framework has several advantages here: for example, it allows for incorporation of external information (e.g. that most genetic markers will, individually, have small effects), and it takes into account of uncertainty in parameter estimates when making other inferences (e.g.phenotype prediction).
For the mean µ and the inverse of error variance, τ , we use the standard conjugate prior distributions: where κ 1 and κ 2 denote, respectively, shape and rate parameters of a Gamma distribution.Specifically we consider the posterior that arises in the limits κ 1 → 0, κ 2 → 0 and σ 2 µ → ∞.These limits correspond to improper priors, but the resulting posteriors are proper, and scale appropriately with shifting or scaling of the phenotype vector y [58].In particular, these priors have the property that conclusions will be unaffected by changing the units of measurement of the phenotype, which seems desirable for a method intended for general application.
Prior specification for the remaining hyper-parameters π, σ 2 a , σ 2 b is perhaps more important.Our approach is to extend the prior distributions for BVSR described in [19].
Following [19] we place a uniform prior on log π: where p is total number of markers being analyzed.The upper and lower limits of this prior were chosen so that π (the expected proportion of markers with non-zero β) ranges from 1/p to 1.A uniform prior on log π reflects the fact that uncertainty in π in a typical GWAS will span orders of magnitude.A common alternative (see e.g.[18,34]) is a uniform distribution on π, but as noted in [19] this puts large weight on large numbers of markers having non-zero β (e.g. it would correspond to placing 50% prior probability to the event that more than half of the markers have non-zero β, and correspond to placing 90% prior probability to the event that more than 10% of the markers have non-zero β).
To specify priors for σ a and σ b , we exploit the following idea from [19]: prior specification is easier if we first re-parameterize the model in terms of more interpretable quantities.Specifically we extend ideas from [19] to re-parameterize the model in terms of the (expected) proportion of phenotypic variance explained by the sparse effects and by the random effects.
To this end, we define PVE (the total proportion of variance in phenotype explained by the sparse effects and random effects terms together) and PGE (the proportion of genetic variance explained by the sparse effects terms) as functions of β, u and τ : where the function V (x) is defined as These definitions ensure that both PVE and PGE must lie in the interval [0, 1].PVE reflects how well one could predict phenotype y from the available SNPs if one knew the optimal β as well as the random effects u; together with PVE, PGE reflects how well one could predict phenotype y using β alone.
Since PVE and PGE are functions of ( β, u, τ ), whose distributions depend on hyper-parameters π, σ a , σ b , the prior distribution for PVE and PGE depends on the priors assigned to these hyperparameters.In brief, our aim is to choose priors for the two hyper-parameters σ 2 a and σ 2 b so that the induced priors on both PVE and PGE are roughly uniform on 0 and 1. (Other distributions could be chosen if desired, but we consider this uniform distribution one reasonable default.)However, because the relationship between the distribution of PVE, PGE and the hyper-parameters is not simple, we have to make some approximations.
Specifically, we introduce h, ρ as approximations (they are ratios of expectations rather than expectations of ratios) to the expectations of PVE and PGE, respectively: where s a is the average variance of genotypes across markers, and s b is the mean of diagonal elements in K.In other words, ii , where x ij and k ij are the ijth elements of matrices X and K, respectively.See Text S1 Detailed Methods for derivations.Intuitively, the term pπs a σ 2 a captures the expected genetic variance contributed by the sparse effects term (relative to the error variance), because pπ is the expected number of causal markers, σ 2 a is the expected effect size variance of each causal marker (relative to the error variance), and s a is the average variance of marker genotypes.Similarly, the term s b σ 2 b captures the expected genetic variance contributed by the random effects term (relative to the error variance), because σ 2 b is the expected variance of the random effects (relative to the error variance) when the relatedness matrix has unit diagonal elements, while s b properly scales it when not.
The parameter ρ provides a natural bridge between the LMM and BVSR: when ρ = 0 BSLMM becomes the LMM, and when ρ = 1 BSLMM becomes BVSR.In practice, when the data favors the LMM, the posterior distribution of ρ would mass near 0, and when the data favors BVSR, the posterior distribution of ρ would mass near 1.
In summary, the above re-parameterizes the model in terms of (h, ρ, π) instead of (σ a , σ b , π).Now, instead of specifying prior distributions for σ a , σ b , we rather specify prior distributions for h, ρ.Specifically we use uniform prior distributions for h, ρ: independent of one another and of π.Since h and ρ approximate PVE and PGE, these prior distributions should lead to reasonably uniform prior distributions for PVE and PGE, which we view as reasonable defaults for general use.(If one had specific information about PVE and PGE in a given application then this could be incorporated here.)In contrast it seems much harder to directly specify reasonable default priors for σ a , σ b (although these priors on h, ρ, π do of course imply priors for σ a , σ b ; see Text S1 Detailed Methods).Note that we treat h and ρ as approximations to PVE and PGE only for prior specification; when estimating PVE and PGE from data we do so directly using their definitions ( 13) and ( 14) (see below for details).
We use the resulting sample from the posterior distribution (23) to estimate PVE and PGE as follows.For each sampled value of (h, ρ, π, γ), we sample a corresponding value for (τ, β, u) from the conditional distribution P (τ, β, u|y, h, ρ, π, γ).We then use each sampled value of (τ, β, u) to compute a sampled value of PVE and PGE, using equations ( 13) and (14).We estimate the posterior mean and standard deviation of PVE, PGE, from these samples.
The novel part of our algorithm is a new efficient approach to evaluating the likelihood P (y|h, ρ, π, γ) that considerably reduces the overall computational burden of the algorithm.The direct naive approach to evaluating this likelihood involves a matrix inversion and a matrix determinant calculation that scale cubically with the number of individuals n, and this cost is incurred every iteration as hyper parameter values change.Consequently, this approach is impractical for typical association studies with large n, and ad hoc approximations are commonly used to reduce the burden.For example, both [37] and [41] fix σ 2 b to some pre-estimated value.As we show later, this kind of approximation can reduce the accuracy of predicted phenotypes.Here, we avoid such approximations by exploiting recently developed computational tricks for LMMs [8,9].The key idea is to perform a single eigen-decomposition and use the resulting unitary matrix (consisting of all eigen vectors) to transform both phenotypes and genotypes to make the transformed values follow independent normal distributions.By extending these tricks to BSLMM we evaluate the necessary likelihoods much more efficiently.Specifically, after a single n 3 operation at the start of the algorithm, the per iteration computational burden is linear in n (the same as BVSR), allowing large studies to be analyzed.
Full details of the sampling algorithm appear in Text S2.

URLs
Software implementing our methods is included in the GEMMA software package, which is freely available at http://stephenslab.uchicago.edu/software.html.

Simulations: PVE Estimation
Both the LMM and BVSR have been used to estimate the PVE [19,22].Since the LMM assumes that all SNPs have an effect, while BVSR assumes that only a small proportion of SNPs have an effect, we hypothesize that BVSR will perform better when the true underlying genetic structure is sparse and LMM will perform better when the true genetic structure is highly polygenic.Further, because BSLMM includes both as special cases, we hypothesize that BSLMM will perform well in either scenario.
To test these hypotheses, we perform a simulation using real genotypes at about 300,000 SNPs in 3,925 Australian individuals [22], and simulate phenotypes under two different scenarios.In Scenario I we simulate a fixed number S of causal SNPs (with S = 10, 100, 1000, 10000), with effect sizes coming from a standard normal distribution.These simulations span a range of genetic architectures, from very sparse to highly polygenic.In Scenario II we simulate two groups of causal SNPs, the first group containing a small number of SNPs of moderate effect (S = 10 or S = 100), plus a second larger group of 10, 000 SNPs of small effect representing a "polygenic component".This scenario might be considered more realistic, containing a mix of small and larger effects.For both scenarios we added normally-distributed errors to phenotype to create datasets with true PVE =0.6 and 0.2 (equation 13).We simulate 20 replicates in each case, and run the algorithms with all SNPs, including the simulated causal variants, so that the true PVE for typed markers is either 0.6 or 0.2 (if we excluded the causal variants then the true PVE would be unknown).
Figures 1A and 1C, show the root of mean square error (RMSE) of the PVE estimates obtained by each method, and Figure 1B and 1D summarize the corresponding distributions of PVE estimates.In agreement with our original hypotheses, BVSR performs best (lowest RMSE) when the true model is sparse (e.g.Scenario I, S = 10 or S = 100 in Figures 1A, 1C).However, it performs very poorly under all the other, more polygenic, models.This is due to a strong downward bias in its PVE estimates (Figures 1B, 1D).Conversely, under the same scenarios, LMM is the least accurate method.This is because the LMM estimates have much larger variance than the other methods under these scenarios (Figure 1B,1D), although, interestingly, LMM is approximately unbiased even in these settings where its modeling assumptions are badly wrong.As hypothesized, BSLMM is robust across a wider range of settings than the other methods: its performance lies between LMM and BVSR when the true model is sparse, and provides similar accuracy to LMM when not.
Of course, in practice, one does not know in advance the correct genetic architecture.This makes the stable performance of BSLMM across a range of settings very appealing.Due to the poor performance of BVSR under highly polygenic models, we would not now recommend it for estimating PVE in general, despite its good performance when its assumptions are met.

Simulations: Phenotype Prediction
We also compare the three methods on their ability to predict phenotypes from genotypes, using the same simulations.
To measure prediction performance, we use relative prediction gain (RPG; see Text S1 Detailed Methods).In brief, RPG is a standardized version of mean square error: RPG=1 when accuracy is as good as possible given the simulation setup, and RPG=0 when accuracy is the same as simply predicting everyone to have the mean phenotype value.RPG can be negative if accuracy is even worse than that.
Figure 2 compares RPG of different methods for simulations with PVE=0.6 (results for PVE=0.2 are qualitatively similar, not shown).Interestingly, for phenotype prediction, the relative performance of the methods differs from results for PVE estimation.In particular, LMM performs poorly compared with the other two methods in all situations, except for Scenario I with S = 10, 000, the Scenario that comes closest to matching the underlying assumptions of LMM.As we expect, BSLMM performs similarly to BVSR in scenarios involving smaller numbers of causal SNPs (up to S = 1, 000 in Scenario I), and outperforms it in more polygenic scenarios involving large numbers of SNPs of small effect (e.g.Scenario II).This is presumably due to the random effect in BSLMM that captures the polygenic component, or, equivalently, due to the mixture of two normal distributions in BSLMM that better captures the actual distribution of effect sizes.The same qualitative patterns hold when we redo these simulation comparisons excluding the causal SNPs (Figure S1) or use correlation instead of RPG to assess performance (Figure S2 and S3).
For a potential explanation why LMM performs much less well for phenotype prediction than for PVE estimation, we note the difference between these two problems: to accurately estimate PVE it suffices to estimate the "average" effect size reliably, whereas accurate phenotype prediction requires accurate estimates of individual effect sizes.In situations where the normal assumption on effect sizes is poor, LMM tends to considerably underestimate the number of large effects, and overestimate the number of small effects.These factors can cancel one another out in PVE estimation, but both tend to reduce accuracy of phenotype prediction.

Estimating PVE in Complex Human Traits
To obtain further insights into differences between LMM, BVSR and BSLMMM, we apply all three methods to estimate the PVE for five traits in two human GWAS data sets.The first data set contains height measurements of 3,925 Australian individuals with about 300,000 typed SNPs.The second data set contains four blood lipid measurements, including high-density lipoprotein (HDL), low-density lipoprotein (LDL), total cholesterol (TC) and triglycerides (TG) from 1,868 Caucasian individuals with about 550,000 SNPs.The narrow sense heritability for height is estimated to be 0.8 from twin-studies [22,59].The narrow sense heritabilities for the lipid traits have been estimated, in isolated founder populations, to be 0.63 for HDL, 0.50 for LDL, 0.37 for TG in Hutterites [60], and 0.49 for HDL, 0.42 for LDL, 0.42 for TC and 0.32 for TG in Sardinians [61].
Table 2 shows PVE estimates from the three methods for the five traits.PVE estimates from BVSR are consistently much smaller than those obtained by LMM and BSLMM, which are almost identical for two traits and similar for the others.Estimates of PVE from both LMM and BSLMM explain over 50% of the narrow sense heritability of the five traits, suggesting that a sizable proportion of heritability of these traits can be explained, either directly or indirectly, by available SNPs.
These results, with LMM and BSLMM providing similar estimates of PVE, and estimates from BVSR being substantially lower, are consistent with simulation results for a trait with substantial polygenic component.One feature of BSLMM, not possessed by the other two methods, is that it can be used to attempt to quantify the relative contribution of a polygenic component, through estimation of PGE, which is the proportion of total genetic variance explained by "large" effect size SNPs (or more precisely, by the additional effects of those SNPs above a polygenic background).Since the PGE is defined within an inevitably over-simplistic model, specifically that effect sizes come from a mixture of two normal distributions, and also because it will be influenced by unmeasured environmental factors that correlate with genetic factors, we caution against over-interpreting the estimated values.We also note that estimates of PGE for these data (Table 2) are generally not very precise (high posterior standard deviation).Nonetheless, it is interesting that the estimated PGE for height, at 0.12, is lower than for any of the lipid traits (ranging from 0.18 for TG to 0.46 for TC), and that all these estimates suggest a substantial contribution from small polygenic effects in all five traits.

Predicting Disease Risk in the WTCCC Data Set
To assess predictive performance on real data, we turn to the Wellcome trust case control consortium (WTCCC) 1 study [62], which have been previously used for assessing risk prediction [63][64][65].These data include about 14,000 cases from seven common diseases and about 3,000 shared controls, typed at a total of about 450,000 SNPs.The seven common diseases are bipolar disorder (BD), coronary artery disease (CAD), Crohn's disease (CD), hypertension (HT), rheumatoid arthritis (RA), type 1 diabetes (T1D) and type 2 diabetes (T2D).
We compared the prediction performance of LMM, BVSR and BSLMM for all seven diseases.Following [64], we randomly split the data for each disease into a training set (80% of individuals) and a test set (remaining 20%), performing 20 such splits for each disease.We estimated parameters from the training set by treating the binary case control labels as quantitative traits, as in [5,9].[This approach can be justified by recognizing the linear model as a first order Taylor approximation to a generalized linear model; we discuss directly modeling binary phenotypes in the Discussion section.]We assess prediction performance in the test set by area under the curve (AUC) [66].
Figure 3 shows AUC for the three methods on all seven diseases.As in our simulations, we find BSLMM performs as well as or better than either of the other two methods for all seven diseases.Indeed, the performance of BSLMM appears to compare favorably with previous methods applied to the same data set [63][64][65] (a precise comparison with previous results is difficult, as some studies use slightly different splitting strategies [63,65] and some do not perform full cross validation [64]).As might be expected from the simulation results, BVSR performs better than LMM in diseases where a small number of relatively strong associations were identified in the original study [62] (CD, RA and T1D) and worse in the others.We obtained qualitatively similar results when we measured performance using the Brier score instead of AUC (Text S3; Figure S4).
Finally, we caution that, although BSLMM performs well here relative to other methods, at the present time, for these diseases, its prediction accuracy is unlikely to be of practical use in human clinical settings.In particular, in these simulations the number of cases and controls in the test set is roughly equal, which represents a much easier problem than clinical settings where disease prevalence is generally low even for common diseases (see [64] for a relevant discussion).

Predicting Quantitative Phenotypes in Heterogeneous Stock Mice
In addition to the WTCCC data set, we also assess perdition performance using a mouse data set [67], which has been widely used to compare various phenotype prediction methods [37][38][39].The mouse data set is substantially smaller than the human data (n ≈ 1000, p ≈ 10, 000, with exact numbers varying slightly depending on the phenotype and the split).This makes it computationally feasible to compare with a wider range of other methods.Therefore, we include in our comparisons here five other related approaches, some of which have been proposed previously for phenotype prediction.Specifically we compare with: 1. LMM-Bayes, a Bayesian version of the LMM, which we fit by setting ρ = 0 in BSLMM using our software.
3. BayesA-Flex, our own modification of BayesA, which assumes a t distribution for the effect sizes.
Our modification involves estimating the scale parameter associated with the t 4 distribution from the data (Text S1 Detailed Methods).Although the original BayesA performs poorly in this data set [38], this simple modification greatly improves its prediction performance (emphasizing again the general importance of estimating hyper-parameters from the data rather than fixing them to arbitrary values).We modified the R package BLR [39] to obtain posterior samples from this model.

4.
BayesCπ, implemented in the online software GenSel [34].This implementation does not allow random effects, and therefore uses the same model as our BVSR, although with different prior distributions.
5. BSLMM-EB (where EB stands for empirical Bayes1 ), an approximation method to fit BSLMM.The method fixes the variance component σ 2 b to its REML (restricted maximum likelihood) estimate obtained with β = 0, which is one of several strategies used in previous studies to alleviate the computational burden of fitting models similar to BSLMM [37,41].We sample posteriors from this model using our software.
See Text S1 Detailed Methods for further details.
Following previous studies that have used these data for prediction [37][38][39] we focused on three quantitative phenotypes: CD8, MCH and BMI.These phenotypes have very different estimated narrow sense heritabilities: 0.89, 0.48, and 0.13 respectively [69].Table S1 lists estimates of some key quantities for the three traits -including PVE, PGE and log 10 (π) -obtained from LMM, BVSR and BSLMM.All three methods agree well on the PVE estimates, suggesting that the data is informative enough to overwhelm differences in prior specification for PVE estimation.
Following [37,38], we divide the mouse data set roughly half and half into a training set and a test set.As the mice come from 85 families, and individuals within a family are more closely related than individuals from different families, we also follow previous studies and use two different splits of the data: the intra-family split mixes all individuals together and randomly divides them into two sets of roughly equal size; the inter-family split randomly divides the 85 families into two sets, where each set contains roughly half of the individuals.We perform 20 replicates for each split of each phenotype.It is important to note that the intra-family split represents an easier setting for phenotype prediction, not only because individuals in the test set are more related genetically to those in the training set, but also because the individuals in the test set often share a similar environment with those in the training set (specifically, in the intra-family split, many individuals in the test set share a cage with individuals in the training set, but this is not the case in the inter-family split).
We apply each method using genotypes only, without other covariates.We obtain effect size estimates in the training data set, and assess prediction performance using these estimates in the test set by root of mean square error (RMSE), where the mean is across individuals in the test set.We contrast the performance of other methods to BSLMM by calculating the RMSE difference, where a positive number indicates worse performance than BSLMM.We perform 20 inter-family splits and 20 intra-family splits for each phenotype.
Figure 4 summarizes the prediction accuracy, measured by RMSE, of each method compared against BSLMM.Measuring prediction performance by correlation gives similar results (Figure S5).For the lowheritability trait BMI, where no measured SNP has a large effect, all methods perform equally poorly.For both the more heritable traits, CD8 and MCH, BSLMM consistently outperformed all other methods, which seem to split into two groups: LMM, LMM-Bayes and Bayesian Lasso perform least well and similarly to one another on average; BVSR, BayesA-Flex, BayesCπ and BSLMM-EB perform better, and similarly to one another on average.A general trend here is that accuracy tends to increase as model assumptions improve in their ability to capture both larger genetic effects, and the combined "polygenic" contribution of smaller genetic effects (and possibly also confounding environmental effects correlating with genetic background).In particular, the t 4 distribution underlying BayesA-Flex, which has a long tail that could capture large effects, performs noticeably better than either the normal or double-exponential distributions for effect sizes underlying LMM and Bayesian Lasso.
Comparisons of pairs of closely-related methods yield additional insights into factors that do and do not affect prediction accuracy.The fact that BSLMM performs better than BSLMM-EB illustrates how the approximation used in BSLMM-EB can degrade prediction accuracy, and thus demonstrates the practical benefits of our novel computational approach that avoids this approximation.Similarly, the superior performance of BayesA-Flex over BayesA (which performed poorly; not shown) also illustrates the benefits of estimating hyper parameters from the data, rather than fixing them to pre-specified values.
The similar performance between BVSR and BayesCπ, which fit the same model but with different priors, suggests that, for these data, results are relatively robust to prior specification.Presumably, this is because the data are sufficiently informative to overwhelm the differences in prior.

Computational Speed
The average computational time taken for each method on the Mouse data is shown in Table 3.Some differences in computational time among methods may reflect implementational issues, including the language environment in which the methods are implemented, rather than fundamental differences between algorithms.In addition, computing times for many methods will be affected by the number of iterations used, and we did not undertake a comprehensive evaluation of how many iterations suffice for each algorithm.Nonetheless, the results indicate that our implementation of BSLMM is competitive in computing speed with the other (sampling-based) implementations considered here.
In particular, we note that BSLMM is computationally faster than BVSR.This is unexpected, since BSLMM is effectively BVSR plus a random effects term, and the addition of a random effects term usually complicates computation.The explanation for this is that the (per-iteration) computational complexity of both BSLMM and BVSR depends, quadratically, on the number of selected SNPs in the sparse effects term (|γ|), and this number can be substantially larger with BVSR than with BSLMM, because with BVSR additional SNPs are included to mimic the effect of the random effects term in BSLMM.The size of this effect will vary among data sets, but it can be substantial, particularly in cases where there are a large number of causal SNPs with small effects.
To illustrate this, Table 4 compares mean computation time for BSLMM vs BVSR for all data sets used here.For simulated data with a small number of causal SNPs, BSLMM and BVSR have similar computational times.However, in other cases (e.g.PVE=0.6,S=10,000 in Scenario I) BSLMM can be over an order of magnitude faster than BVSR.In a sense, this speed improvement of BSLMM over BVSR is consistent with its hybrid nature: in highly polygenic traits, BSLMM tends to behave similarly to LMM, resulting in a considerable speed gain.

Discussion
We have presented novel statistical and computational methods for BSLMM, a hybrid approach for polygenic modeling for GWAS data that simultaneously allows for both a small number of individually large genetic effects, and combined effects of many small genetic effects, with the balance between the two being inferred from the data in hand.This hybrid approach is both computationally tractable for moderately large data sets (our implementation can handle at least 10,000 individuals with 500,000 SNPs on our moderately-equipped modern desktop computer), and is sufficiently flexible to perform well in a wide range of settings.In particular, depending on the genetic architecture, BSLMM is either as accurate, or more accurate, than the widely-used LMM for estimating PVE of quantitative traits.And for phenotype prediction BSLMM consistently outperformed a range of other approaches on the examples we considered here.By generalizing two widely-used models, and including both as special cases, BSLMM should have many applications beyond polygenic modeling.Indeed, despite its increased computational burden, we believe that BSLMM represents an attractive alternative to the widely-used LASSO [70] for general regression-based prediction problems.
Although it was not our focus here, BSLMM can be easily modified to analyze binary phenotypes, including for example, a typical human case-control GWAS.For PVE estimation, one can directly apply BSLMM, treating the 1/0 case-control status as a quantitative outcome, and then apply a correction factor derived by [24] to transform this estimated PVE on the "observed scale" to an estimated PVE on a latent liability scale.This correction, for which we supply an alternative derivation in Text S3, corrects for both ascertainment and the binary nature of case-control data.For phenotype prediction, one can again directly apply BSLMM, treating the 1/0 case-control status as a quantitative outcome, and interpret the resulting phenotype predictions as the estimated probability of being a case.Although in principle one might hope to improve on this by modifying BSLMM to directly model the binary outcomes, using a probit link for example, we have implemented this probit approach and found that not only is it substantially more computationally expensive (quadratic in n instead of linear in n), but it performed slightly worse than treating the binary outcomes as quantitative, at least in experiments based on the mouse phenotypes considered here (Text S3).This may partly reflect inaccuracies introduced due to the greater computational burden.
The computational innovations we introduce here, building on work by [8,9], make BSLMM considerably more tractable than it would otherwise be.Nonetheless, the computational burden, as with other posterior sampling based methods, remains heavy, both due to memory requirements (e.g. to store all genotypes) and CPU time (e.g. for the large number of sampling iterations required for reasonable convergence).Although more substantial computational resources will somewhat increase the size of data that can be tackled, further methodological innovation will likely be required to apply BSLMM to the very large data sets that are currently being collected.
In addition to providing a specific implementation that allows BSLMM to be fitted to moderately large data sets, we hope that our work also helps highlight some more general principles for improving polygenic modeling methodology.These include: 1.The benefits of characterizing different methods by their effect size distribution assumptions.While this point may seem obvious, and is certainly not new (e.g.[40,71]), polygenic modeling applications often focus on the algorithm used to fit the model, rather than the effect size distribution used.While computational issues are very important, and often interact with modeling assumptions, we believe it is important to distinguish, conceptually, between the two.One benefit of characterizing methods by their modeling assumptions is that it becomes easier to predict which methods will tend to do well in which settings.
flexible distributions such as a single normal (LMM) or a point-normal (BVSR).More generally, in our phenotype prediction experiments, methods with more flexible effect size distributions tended to perform better than those with less flexible distributions.
3. The importance of estimating hyper-parameters from data, rather than fixing them to pre-specified values.Here we are echo-ing and reinforcing similar themes emphasized by [34] and [19].Indeed, our comparison between BSLMM and BSLMM-EB for phenotype prediction illustrates the benefits not only of estimating hyper-parameters from the data, but of doing so in an integrated way, rather than as a two-step procedure.
4. The attractiveness of specifying prior distributions for hyper-parameters by focusing on the proportion of variance in phenotypes explained by different genetic components (e.g.PVE and PGE in our notation).This idea is not limited to BSLMM, and could be helpful even with methods that make use of other effect size distributions.
One question to which we do not know the answer is how often the mixture of two normal distributions underlying BSLMM will be sufficiently flexible to capture the actual effect size distribution, and to what extent more flexible distributional assumptions (e.g. a mixture of more than two normals, or a mixture of t distributions with degrees of freedom estimated from the data) will produce meaningful gains in performance.It seems likely that, at least in some cases, use of a more flexible distribution will improve performance, and would therefore be preferable if it could be accomplished with reasonable computational expense.Unfortunately some of the tricks we use to accomplish computational gains here may be less effective, or difficult to apply, for more flexible distributions.In particular, the tricks we use from [8] and [9] may be difficult to extend to allow for mixtures with more than two components.In addition, for some choices of effect size distribution, one might have to perform MCMC sampling on the effect sizes β directly, rather than sampling γ, integrating β out analytically, as we do here.It is unclear whether this will necessarily result in a loss of computational efficiency: sampling β reduces computational expense per update at the cost of increasing the number of updates necessary (sampling γ by integrating over β analytically ensures faster mixing and convergence [72,73]).Because of these issues, it is difficult to predict which effect size distributions will ultimately provide the best balance between modeling accuracy and computational burden.Nonetheless, compared with currently available alternatives, we believe that BSLMM strikes an appealing balance between flexibility, performance, and computational tractability.

GWAS Datasets
We used four GWAS data sets in the present study.
The first data set contains height measurements for 3925 Australian individuals [22].We used the data set for simulation and PVE estimation.The phenotypes were already regressed out of age and sex effects, and were quantile normalized to a standard normal distribution afterwards [22].A total of 294,831 SNPs were available after stringent quality control [22].We imputed missing SNPs using IMPUTE2 [74], and further excluded five SNPs that have minor allele frequencies below 1% after imputation.
The second data set consists of blood lipid measurements for 1868 individuals [75].We used this data for PVE estimation.The individuals came from two study groups: the Cholesterol and Pharmacogenetics (CAP) group [76] and the Pravastatin Inflammation/CRP Evaluation (PRINCE) group [77].The PRINCE study consists of two cohorts: one contains individuals with history of cardiovascular diseases (CVD) and the other contains individuals with no history of CVD.From both study groups, we selected all 1868 individuals that have complete low-density lipoprotein (LDL) subfraction measurements.We selected four different blood lipid measurements as phenotypes in the present study: LDL, high-density lipoprotein (HDL), total cholesterol (TC) and triglycerides (TG).Phenotypes were quantile-normalized to a standard normal distribution within each group, corrected for covariates including BMI (body mass index), age, sex, and smoking status effects, and quantile-normalized again [75].Individuals were typed on two different SNP arrays (Illumina HumanHap300 and HumanQuad610 bead chips, Illumina, San Diego, CA).We used all SNPs that appeared in either of the arrays and we imputed missing genotypes using IMPUTE2 [74].We obtained a total of 582,962 SNPs and we used 555,601 polymorphic SNPs with minor allele frequency above 1% for analysis.
The third data set is from the Wellcome trust case control consortium (WTCCC) 1 study [62].We used this data set to assess phenotype prediction performance.The data set consists of about 14,000 cases of seven common diseases, including 1868 cases of bipolar disorder (BD), 1926 cases of coronary artery disease (CAD), 1748 cases of Crohn's disease (CD), 1952 cases of hypertension (HT), 1860 cases rheumatoid arthritis (RA), 1963 cases of type 1 diabetes (T1D) and 1924 cases of type 2 diabetes (T2D), as well as 2938 shared controls.We obtained quality controlled genotypes from WTCCC and we further imputed missing genotypes using BIMBAM [78], which resulted in a total of 458,868 shared SNPs.All polymorphic SNPs with minor allele frequency above 1% in the training data were used for prediction (about 400,000 SNPs; depending on the disease and the split).
The fourth data set comes from a genetically heterogeneous stock of mice, consisting of 1904 individuals from 85 families, all descended from eight inbred progenitor strains [67].We used this data set to assess phenotype prediction performance of several methods.Multiple phenotype measurements are available for the data set, and we selected three phenotypes among them: percentage of CD8+ cells (CD8, n = 1410), mean corpuscular hemoglobin (MCH, n = 1580) and body mass index (BMI, n = 1828).We selected these phenotypes because they were previously used for comparing prediction performance of various methods [37][38][39], and they represent a wide range of narrow sense heritability: CD8 has a high heritability, MCH has a median heritability and BMI has a low heritability [69].All phenotypes were already corrected for sex, age, body weight, season and year effects [67], and we further quantile normalized the phenotypes to a standard normal distribution.A total of 12,226 autosomal SNPs were available for all mice.For individuals with missing genotypes, we imputed missing values by the mean genotype of that SNP in their family.All polymorphic SNPs with minor allele frequency above 1% in the training data were used for prediction (about 10,000 SNPs; depending on the phenotype and the split).

Simulations
We used genotypes from the human height data set [22] described above and simulated phenotypes from the simple linear model (1) with different assumptions for the distribution of effect sizes β.We consider two simulation scenarios where the true PVE is known and we simulated 20 independent sets of phenotype data in each case.
Scenario I: the effect sizes of causal SNPs come from a normal distribution.We randomly chose a fixed number of causal SNPs (10, 100, 1000, 10000) and simulated their effect sizes from a N(0, 1) distribution.We drew the errors from a normal distribution with variance chosen to achieve a given PVE (0.2 and 0.6).
Scenario II: the effect sizes of causal SNPs come from a mixture of two normal distributions, such that a small group of causal SNPs have additional effects.We first randomly chose a large number of causal SNPs (10000), and among them, we further selected a small number of medium effect size SNPs (10 or 100) and used what were left as small effect size SNPs (9990 or 9900).We simulated the small effect sizes for all causal SNPs (10000) from a N(0, 1) distribution.Afterwards, we drew additional effect sizes (in addition to the small effects already drawn) for those medium effect SNPs (10 or 100) from a N(0, 1) distribution, and scaled these additional effect sizes further so that together they explained a fixed proportion of genetic variance, or PGE (0.1 and 0.2, for 10 and 100 medium effect size SNPs, respectively).Once we obtained the final effect sizes for all causal SNPs, we drew errors to achieve a given PVE (0.2 and 0.6).

MSPE and RPG
We assess prediction accuracy mainly using mean square prediction error (MSPE), and a rescaled version of MSPE called relative predictive gain (RPG).The MSPE for predicting a future observation in the simple linear model (1) depends on comparing an estimated value of β, β, the true value of β, and the error variance τ , as follows: = E(( where y is the phenotype for a future observation, x is the corresponding p-vector of genotypes, r ij = E(x i x j ) is the covariance between SNP i and j.In practice, we approximate r ij with the sample covariance s ij = 1 n n k=1 x ik x jk , and we only consider s ij for neighboring SNPs that satisfy |i − j| ≤ 20.This is because linkage dis-equilibrium (LD) decays with distance and remote SNPs are approximately independent with each other.The above definition of MSPE extends the definition in [19] to take into account correlations among neighboring SNPs.
We denote MSPE 0 as the MSPE obtained using only the mean of the phenotype (i.e.ȳ) for prediction, and we define RPG as the rescaled version of MSPE following [19]: When applying these formulae for LMM or BSLMM note that we use estimates β for β in (1), and not for β in (6).This ensures that the resulting predictions take account of both the sparse effects β and the random effect u in (6).The way we obtain these estimates β is detailed below.

Correlation
We also assess prediction accuracy using correlation.The correlation between the predicted value and the true value for a future observation in the simple linear model (1) depends on β, β and τ as follows: Again, when applying these formulae for LMM or BSLMM note that we use estimates β for β in (1), and not for β in (6).

Centering X and K
We assume that the genotypes in X have been measured on bi-allelic markers, and that the genotypes at each marker are coded as 0, 1 or 2 copies of some reference allele.(For imputed genotypes we use the posterior mean genotype [78].)It occasionally simplifies the algebra to assume that each column of X is centered to have mean 0; since the results will be the same with or without centering, we make this assumption throughout.It is also common to standardize the columns of X to have unit variance.In contrast to centering, standardizing the columns will affect the results, and we do not standardize the columns in our applications here, although all our methods could be applied with the matrix standardized in this way.(Standardizing the columns of X corresponds to making an assumption that rarer variants tend to have larger effects than common variants, and precisely that marker effect sizes tend to decay with the inverse of the genotype variance; see [79,80] for relevant discussion.)In summary, throughout this paper X ij = (x ij − xj ) where x ij is the number of copies of the reference allele at marker j in individual i, and xj := (1/n) i x ij .
To facilitate prior specification, in addition to centering the genotype matrix X, we also assume that the relatedness matrix K is "centered", in the sense that the random effects have mean zero: n i=1 u i = 0.This holds automatically for K ∝ XX T , with X centered.More generally it can be achieved by multiplying the relatedness matrix with a projection matrix on both sides: MKM, where The resulting transformed relatedness matrix is positive-semidefinite as long as the original relatedness matrix is positive-semidefinite. where x ij and k ij are the ijth elements of matrices X and K, respectively.The above derivation assumes centered genotypes and relatedness matrix.
Induced Priors on σ 2 a and σ 2 b Solving ( 16) and ( 17) for σ 2 a and σ 2 b as functions of h, ρ and π gives: The independent priors ( 18), ( 19), ( 12) for (h, ρ, π) induce a joint prior on (σ a , σ b , π): which is heavy tailed for (σ 2 a , σ 2 b ) marginally (i.e.tail has a polynomial decay), a feature desirable in association studies [80].The priors also have another nice property that the marker effect size variance σ 2 a tends to decrease as the proportion of markers with an effect (π) increases.

PVE Estimation with LMM, BVSR and BSLMM
We could estimate PVE using the posterior mean of the MCMC samples for all three models, and we do so for both BVSR and BSLMM.For LMM, we follow previous studies [22] and use an approximation to the PVE.In particular, we consider the LMM defined by equations ( 6)-( 7) with β = 0, and we estimate PVE by the ratio of expectations where σ2 b is the REML estimate for the variance component.This formula can be viewed as a generalization of the form used in [22], and is valid for any choice of centered relatedness matrix K.
To obtain the standard error of the above estimate, we compute the second derivative of the log restricted likelihood function with respect to σ 2 b [9] and evaluate it at σ2 b : where l r denotes the log restricted likelihood, Despite its complicated form, the second derivative can be easily and efficiently evaluated using recursions in [9].As the variance of σ2 b is asymptotically v(σ 2 b ) = −1/l ′′ (σ 2 b ), using the delta method, we approximate the standard error of PVE by As a check on correctness, we also used the posterior mean for PVE obtained from LMM-Bayes to estimate PVE.This gave almost identical results in all cases considered here and therefore only results from LMM are presented.

Real Data
In the real data, we can perform prediction by estimating both the sparse effects (µ, β) and the random effects (u) in (6) from the training set, and use these to predict phenotypes in the test set.
Let μo , βo and ûo denote estimates for the sparse and random effects obtained from the observed (training) sample.For BSLMM and BVSR these estimates are the posterior means for these parameters, estimated from MCMC samples (for BVSR ûo ≡ 0).For LMM ûo is obtained as the conditional posterior mean of u o given the REML estimate for σ 2 b (i.e.BLUP).We then obtain predictions for a future (test) sample as follows.For a general relatedness matrix K, we assume that the random effects for the observed and future samples follow a multivariate normal distribution where n o and n f are the sample size for observed (training) and future (test) data, respectively.Standard multivariate normal theory gives the conditional distribution We use the conditional mean as an estimate for u f and thus the predicted phenotypes for future observations are ŷf =

Simulated Data
In the simulated data we use RPG and correlation to assess prediction accuracy.To compute RPG and correlation we need to obtain estimates for β in the simple linear model (1).To do this, in the special case when the relatedness matrix K = 1 p XX T , we rewrite the model ( 6) as where we can think of the p-vector α as representing the "small" effect sizes present at every locus.The special case of π = 0 ( β ≡ 0) gives LMM, and the special case of σ 2 b = 0 (α ≡ 0) gives BVSR.Note that α + β = β, so summing estimates of α and β yields an estimate for β in (1).
For LMM we estimate α by its conditional expectation where σ2 b is the REML estimate of the variance component in the observed data.Since β ≡ 0 in LMM, this estimate for α provides the required estimate for β in (1).
For BVSR, we use the posterior mean of β (since α ≡ 0).For BSLMM, we use Rao-Blackwellisation to obtain an approximation to the posterior mean for α (Text S2), and then add this to the (approximate) posterior mean for β obtained from the MCMC sampler to obtain an approximation for the posterior mean of β in (1).

BVSR:
We fit BVSR by fixing ρ = 1 in BSLMM using our software.This gives slightly better results, and is faster than the BVSR software piMASS (version 0.90) [19], in all examples considered here.

LMM-Bayes:
We fit this by fixing ρ = 0 in BSLMM using our software.

5.
BayesA-Flex: This assumes a scaled t-distribution for each coefficient β i in (1): where IG stands for the inverse gamma distribution.Following previous studies, we set the degree of freedom parameter ν to 4 [27,34,40] and set κ 3 = 1/2 and κ 4 = 1/2 [39].We also consider the posterior distribution where κ 1 → 0 and κ 2 → 0. The above model is similar to BayesA [27], but with a key difference in the way the scaling parameter σ 2 is treated: BayesA fixes σ 2 to some pre-specified value, whereas here we specify a prior for σ 2 and allow it to be estimated from the data (and hence the name "BayesA-Flex").Using BayesA in this data set gives poor results [38] (and data not shown); but estimating the scaling parameter greatly improves prediction performance.We modified the R package BLR [39] to obtain posterior samples from this model.The modified code is freely available online.
7. BSLMM-EB: we fit this using our BSLMM software, fixing σ 2 b to the REML estimate σ2 b from the null model (i.e.LMM).This approximation avoids updating σ 2 b in each iteration of the MCMC, and is one of the several approximation strategies used by previous studies to alleviate the computation burden of models similar to BSLMM [37,41].Intuitively, by fixing the variance component to its null estimate, BSLMM-EB discourages the inclusion of large effect SNPs into the sparse effects term and risks underestimating their effect sizes.Therefore, this approximation may reduce the prediction performance of BSLMM, especially when there are large effect SNPs.We confirm this in the real data set.
For all MCMC based methods except for BayesCπ, we run 2.1 million iterations with the first 0.1 million iterations as burn-in steps.For BayesCπ, due to web server restriction, we run 1.1 million iterations with the first 0.1 million iterations as burn-in.
In the above equation, we explored the fact that the parameters β, u and τ can be integrated out analytically to compute the marginal likelihood P (y|h, ρ, π, γ).The marginal likelihood is where a is a function of h, ρ and π, while σ 2 b is a function of h and ρ.To efficiently evaluate the marginal likelihood, we perform an eigen decomposition of the relatedness matrix K = UDU T at the beginning of the MCMC, where U is the matrix of eigen vectors and D is a diagonal matrix of eigen values.We transform both the phenotype vector and the genotype matrix by multiplying the eigen matrix and calculate U T y and U T X. Afterwards, as has been shown previously, the calculations of the determinant and the inverse of matrix H, as well as the vector-matrix-vector form y T Py, in each iteration of the MCMC, are easy to perform [8,9].
We use a standard Metropolis-Hastings algorithm to draw posterior samples of the hyper-parameters (h, ρ, π, γ) based on the above marginal likelihood.Following [19], we use a rank based proposal distribution for γ, and use random walk proposals based on uniform distributions for h, ρ and log(π).In particular, we first obtain single-SNP p values using a standard LMM with GEMMA algorithm [9], and then rank SNPs based on these p values from small to large.Our aim is to use a proposal distribution for γ that puts more weights on SNPs that are ranked higher by the single SNP tests, and to do this we consider a mixture distribution Q p = 0.3U p + 0.7G p , where U p is a uniform distribution on 1, • • • , p and G p is a geometric distribution truncated to 1, • • • , p with its parameter chosen to give a mean of 2000.We denote γ + = {i : γ i = 1} and we propose the new γ randomly choose one of the following steps: • add a covariate with probability 0.4: generate r from Q p until the covariate with rank r is not in γ + , then add this covariate to γ + • remove a covariate with probability 0.4: pick a covariate in γ + uniformly at random and remove it from γ + • switch a pair of covariates with probability 0.2: pick up two covariates by the above two steps, and switch their indicator values For the other hyper-parameters, we update log(π) by adding a random variable from U(-0.05, 0.05) to the current value, and update h and ρ by adding a random variable from U(-0.1, 0.1) to the current values.New values of h and ρ that lie outside the boundary [0,1] are reflected back.
In addition to the above local proposal distributions, we also use a "small world proposal" which improves theoretical MCMC convergence [81].In brief, with probability 0.33 in each iteration, we make a longer-range proposal by compounding many local moves, where the number of compounded local moves is draw uniformly from 1 to 20.
For each sampled values of (h, ρ, π, γ), we further obtain samples of τ and β using the conditional distributions P (τ |y, h, ρ, π, γ) and P ( β|y, h, ρ, π, γ, τ ) listed below: values that are above 1 to be exact 1 and those below 0 to be exact 0. We contrast the performance of the probit BSLMM against BSLMM by calculating the Brier score difference, where a positive value indicates worse performance than BSLMM.
Figure S6 shows Brier score differences for the three traits.Interestingly, for the three traits here, treating binary values as quantitative traits using BSLMM works better than modeling them directly using the probit BSLMM.This may partly reflect numerical inaccuracies due to the greater computational burden of fitting the probit BSLMM.

Correction factor for estimating PVE for case-control data
Here, we provide an alternative way to derive the correction factor, that appeared in [24], for transforming PVE estimate in the observed scale back to that in the liability scale.Our approach is based on Taylor series approximation.To simplify notation, we denote k p = P p (y i = 1) as the case proportion in the population, k a = P a (y i = 1) as the case proportion in the ascertained case-control sample, Φ as the normal CDF (cumulative distribution function), φ as the normal PDF (probability distribution function), µ p satisfies Φ(µ p ) = k p , µ a satisfies Φ(µ a ) = k a , and z p = φ(µ p ).
First, we assume, following [24], a probit model on the population scale: where N is the population sample size.The conditional distribution in the ascertained case-control sample can be derived by Bayes theorem where n is the case-control sample size.We notice that P a (x i , u i |y i = 1, β) = P p (x i , u i |y i = 1, β) holds for ideal case-control studies, as cases in the ascertained sample are selected randomly from all cases in the population.We assume further P p (y i = 1| β) ≈ P p (y i = 1) = k p and P a (y i = 1| β) ≈ P a (y i = 1) = k a , that the probability of being a case does not depend on parameters, an assumption commonly made (see e.g.[83]) and likely hold when parameters are close to their true values.We further denote a normalizing constant Z = Pa(xi,ui| β) Pp(xi,ui| β) , and we have and similarly which give the normalizing constant We expand the above two likelihoods using Taylor series expansion with respect to x i β + u i at 0. If we use the linear term only for approximation, we obtain In other words, the expected value of individual binary label can be approximated by which suggests using a linear mixed model to treat binary values as quantitative traits to infer the parameters.The estimated PVE on the observed scaling using a linear mixed model is where PVE o is the PVE estimate on the observed scale, and PVE l is the true PVE on the liability scale.
Therefore, we can use the correction factor  Tables Table 1.Summary of some effect size distributions that have been proposed for polygenic modeling.The reference list contains only a selection of relevant publications.Abbreviations: DE denotes double exponential distribution, NEG denotes normal exponential gamma distribution, and other abbreviations can be found in the main text.In the scaled t-distribution, ν and σ 2 a are the degree of freedom parameter and scale parameter, respectively.In the DE distribution, θ is the scale parameter.In the NEG distribution, κ and θ are the shape and scale parameters, respectively.Notes: 1.Some applications of these methods combine a particular effect size distribution with a random effects term, with covariance matrix K, to capture sample structure or relatedness.If K ∝ XX T then this is equivalent to adding a normal distribution to the effect size distribution.The listed effect size distributions in this table do not include this additional normal component.2. BayesC has been used to refer to models with different effect size distributions in different papers.3.In some papers, keywords listed here have been used to refer to fitting techniques rather than effect size distributions.Values in parentheses are standard deviations.Means and standard deviations are calculated based on 2.1 million MCMC iterations in 120 replicates: 20 intra-family and 20 inter-family splits for three phenotypes.Computations were performed on a single core of an Intel Xeon L5420 2.50 GHz CPU.Since computing times for many methods will vary with number of iterations used, and we did not undertake a comprehensive evaluation of how many iterations suffice for each algorithm, these results provide only a very rough guide to the relative computational burden of different methods.

Figure 4 .Figure S2 .Figure S6 .
Figure 4. Comparison of prediction performance of several models with BSLMM for three traits in the heterogenous stock mouse data set.Performance is measured by RMSE difference with respect to BSLMM, where a positive value indicates worse performance than BSLMM.The x-axis shows two different ways to split the data into a training set and a test set, each with 20 replicates.The mean RMSE of BSLMM for the six cases are 0.70, 0.80, 0.79, 0.90, 0.98 and 0.99, respectively.

2 p
to transform the PVE estimate on the observed scale back to that on the liability scale.

number of causal SNPs number of causal SNPs with moderate/small effects
Comparison of prediction performance of LMM (blue), BVSR (red) and BSLMM (purple) in two simulation scenarios, where all causal SNPs are included in the data.Performance is measured by Relative Predictive Gain (RPG).True PVE=0.6.Means and standard deviations (error bars) are based on 20 replicates.The x-axis show the number of causal SNPs (Scenario I) or the number of medium/small effect SNPs (Scenario II).

Table 2 .
PVE and PGE estimates for five human traits.

Table 3 .
Mean computation time, in hours, of various methods for the mouse data set.