Computational cost of BOLT-LMM versus existing methods

To analyze the computational performance of BOLT-LMM, we simulated data sets of sizes ranging from N=3,750 to 480,000 individuals and M=300,000 SNPs. We used genotypes from the WTCCC2 data set³⁰ analyzed in ref.¹², which contains 15,633 individuals of European ancestry, to form mosaic chromosomes, and we used a phenotype model in which 5,000 SNPs explained 20% of phenotypic variance (Supplementary Note).

We benchmarked BOLT-LMM against existing mixed model association methods, running each method for up to 10 days on machines with 96GB of memory. BOLT-LMM completed all analyses through N=480,000 individuals within these constraints, whereas previous methods could only analyze a maximum of N=7,500–30,000 individuals (Fig. 1 and Supplementary Table 1). All previous methods require O(MN²) running time (for M>N), whereas the running time of BOLT-LMM scales roughly with MN^1.5 (Fig. 1a and Supplementary Fig. 1a). We also observed substantial savings in memory use with BOLT-LMM (Fig. 1b and Supplementary Fig. 1b), which requires little more than the MN/4 bytes of memory needed to store raw genotypes (as in GenABEL software³¹).

Open in a separate window

Figure 1

Computational performance of mixed model association methods

Log-log plots of (a) run time and (b) memory as a function of sample size (N). Slopes of the curves correspond to exponents of power-law scaling with N. Benchmarking was performed on simulated data sets in which each sample was generated as a mosaic of genotype data from 2 random “parents” from the WTCCC2 data set (N=15,633, M=360K) and phenotypes were simulated with M_causal=5,000 SNPs explaining h²_causal=0.2 of phenotypic variance. Reported run times are medians of five identical runs using one core of a 2.27 GHz Intel Xeon L5640 processor. We caution that running time comparisons may vary by a small constant factor as a function of computing environment. FaST-LMM-Select (resp. GCTA-LOCO, EMMAX) memory usage exceeded the 96GB available at N=15K (resp. 30K, 60K). GEMMA encountered a runtime error (segmentation fault) at N=30K. Software versions: FaST-LMM-Select, v2.07; GCTA-LOCO, v1.24; EMMAX, v20120210; GEMMA, v0.94. Numerical data are provided in Supplementary Table 1.

The running time of BOLT-LMM depends not only on the cost of matrix arithmetic, which scales linearly with M and N, but also the number of O(MN)-time iterations required for convergence, which empirically scales roughly as N^0.5 (Supplementary Fig. 1) and also varies with heritability, relatedness, and population structure (Supplementary Note and Supplementary Fig. 2). These observations apply both to the full Gaussian mixture modeling performed by BOLT-LMM and to the subset of the computation (Steps 1a and 1b) needed to compute BOLT-LMM-inf infinitesimal mixed model association statistics, which in our benchmarks required ≈40% of the full BOLT-LMM run time (Fig. 1a and Supplementary Fig. 1a). Our results show that even on very large data sets, BOLT-LMM is efficient enough to enable mixed model analysis using a Gaussian mixture prior, which we recommend because of its potential to increase power.

Power and false positive control of BOLT-LMM in simulations

To assess the power of BOLT-LMM to detect associated loci, we performed additional simulations using real genotypes from the WTCCC2 data set, which is an ancestry-stratified sample containing both Northern and Southern European samples. We simulated phenotypes with 1,250–10,000 causal SNPs^13,14 explaining 50% of phenotypic variance and an additional 60 standardized effect SNPs explaining 2% of variance. We included the latter category of SNPs to allow direct power comparisons across different simulation setups, as the 60 standardized effect SNPs always explain the same total amount of variance regardless of other simulation parameters. We further introduced environmental differences in ancestry by including a phenotypic component aligned with the top principal component that explained an additional 1% of variance. (We note that principal component analysis is not part of BOLT-LMM; it is unnecessary to perform PCA when running mixed model association methods¹².) We chose causal SNPs randomly from the first halves of chromosomes, leaving the second halves of chromosomes to contain only non-causal SNPs (Supplementary Note).

We computed χ² association statistics using linear regression with 10 principal components (PCA)³², GCTA-LOCO¹², BOLT-LMM-inf, and BOLT-LMM. We were unable to test FaST-LMM-Select¹⁵ on this data set because of its memory requirements (Fig. 1). For each method, we computed means of its χ² statistics over standardized effect SNPs and compared these means across simulations involving different numbers of causal SNPs (Fig. 2a and Supplementary Table 2). We observed that BOLT-LMM achieved power gains by modeling non-infinitesimal architectures. For the sparsest genetic architecture (1,250 causal SNPs plus 60 standardized effect SNPs), we observed a 25% increase in mean BOLT-LMM χ² statistics at standardized effect SNPs compared to GCTA-LOCO and BOLT-LMM-inf infinitesimal mixed model χ² statistics. This metric is readily interpretable as corresponding to a 25% increase in effective sample size; for completeness, we also computed traditional power curves at two significance thresholds (Supplementary Fig. 3). The power gain of the Gaussian mixture model decreased with increasing numbers of causal SNPs (Fig. 2a). This behavior is expected because the advantage of the Gaussian mixture lies in its ability to more accurately model a small fraction of SNPs with larger effects amid a majority of SNPs with near-zero effects. Larger numbers of causal SNPs explaining a fixed proportion of variance result in smaller effect sizes per causal SNP, giving BOLT-LMM less opportunity for power gain. In contrast, all methods other than BOLT-LMM had performance independent of the number of causal SNPs, consistent with the fact that none of these methods model non-infinitesimal genetic architectures. GCTA-LOCO and BOLT-LMM-inf mean χ² statistics at standardized effect SNPs were essentially identical and slightly exceeded PCA, consistent with theory¹². We also tested EMMAX³ and GEMMA⁶, which are vulnerable to proximal contamination^5,9,12; these methods suffered loss of power relative to PCA (Supplementary Fig. 4a), consistent with theory¹².

Open in a separate window

Figure 2

BOLT-LMM increases power to detect associations in simulations

Mean χ² at standardized effect SNPs as a function of (a) number of causal SNPs, (b) proportion of variance explained by causal SNPs, (c) number of samples. Simulations used real genotypes from the WTCCC2 data set (N=15,633, M=360K) and simulated phenotypes with the specified number of causal SNPs explaining the specified proportion of phenotypic variance and 60 more standardized effect SNPs explaining an additional 2% of the variance. Error bars, s.e.m., 100 simulations. We verified on the first 5 simulations that the BOLT-LMM-inf and GCTA-LOCO statistics are nearly identical (Supplementary Table 7). Numerical data are provided in Supplementary Table 2.

To further explore the relationship between the magnitude of Gaussian mixture model power gain and other parameters of the data set, we also varied the proportion of variance explained by causal SNPs (Fig. 2b) and the number of individuals (Fig. 2c). We observed that the power boost of BOLT-LMM over infinitesimal mixed model analysis (GCTA-LOCO, BOLT-LMM-inf) increased with each of these parameters. In further simulations using data sets of size N=30,000 and N=60,000 (Supplementary Note) and simulated phenotypes with M_causal=250–15,000 causal SNPs explaining 15–35% of the variance, we observed that the effectiveness of the Gaussian mixture model is closely tied to h_g²N/M_causal (where h_g² is the heritability parameter estimated by BOLT-LMM; see Online Methods for interpretation); intuitively, this quantity measures the effective number of samples per causal SNP (Supplementary Fig. 5). These results are consistent with theory (Supplementary Note and Supplementary Table 2 of ref.¹²), which explains that even in the absence of confounding, mixed model analysis provides a power gain over marginal regression by conditioning on the estimated effects of other SNPs when testing a candidate SNP^9,12. As sample size increases, the power gain of both methods approaches an asymptote corresponding to an increase in effective sample size of 1/(1−h_g²), but for sparse genetic architectures, the Gaussian mixture model approaches this asymptote much faster.

To verify that BOLT-LMM is correctly calibrated and robust to confounding, we also computed mean χ² statistics across SNPs on the second halves of chromosomes, simulated to all have zero effect (“null SNPs”). Because our simulated phenotypes included an ancestry effect, linear regression without correcting for population stratification suffered 35% inflation. In contrast, the BOLT-LMM and BOLT-LMM-inf statistics were both well-calibrated (Supplementary Fig. 4b, Supplementary Table 3, and Supplementary Table 4). We further verified that Type I error was properly controlled (Online Methods and Supplementary Table 5) and that the distribution of statistics at null SNPs did not deviate noticeably from a 1 d.o.f. chi-squared distribution (Supplementary Fig. 6a,b). Genomic inflation factors³³ for BOLT-LMM and BOLT-LMM-inf exceeded 1 in these simulations (Supplementary Fig. 6c,d), consistent with polygenicity of the simulated phenotype and use of a mixed model statistic that successfully avoids proximal contamination^12,13. In contrast, EMMAX and GEMMA had deflated test statistics (Supplementary Fig. 4b).

To examine the tightness of the variational approximation used by BOLT-LMM for Bayesian model fitting and to enable comparison with FaST-LMM-Select, we ran a small-scale simulation using the same setup as above but only one-third of the samples (N=5,211). We simulated genetic architectures with 1,250 causal SNPs explaining 70% of phenotypic variance (and 60 additional standardized effect SNPs explaining 2% of variance and ancestry explaining 1%, as before). We ran PCA, BOLT-LMM-inf, BOLT-LMM, FaST-LMM-Select, and a modified version of BOLT-LMM in which we replaced the variational iteration of Step 2b with a Markov chain Monte Carlo (MCMC) Gibbs sampler. In the limit of infinite sampling iterations, MCMC would produce exact versions of the posterior approximations computed by BOLT-LMM. In these simulations, the variational iteration (i.e., standard BOLT-LMM) achieved statistically identical results to MCMC (Supplementary Table 6a), supporting the choice of variational Bayes for BOLT-LMM. We also observed that while BOLT-LMM-inf achieved a power gain over PCA and BOLT-LMM achieved a further power gain over BOLT-LMM-inf (consistent with previous simulations), FaST-LMM-Select achieved lower power than BOLT-LMM-inf and BOLT-LMM (Supplementary Table 6a). Upon repeating this experiment with the number of causal SNPs reduced to 500, we observed that FaST-LMM-Select achieved a power gain in between BOLT-LMM-inf and BOLT-LMM (Supplementary Table 6b). Finally, we observed that the LD Score calibration approach used by BOLT-LMM also worked well when applied to FaST-LMM-Select, validating this calibration approach (Supplementary Table 6).

Lastly, we investigated the similarity between the BOLT-LMM-inf mixed model statistic and existing methods at the individual SNP level. Despite its use of an infinitesimal model, the BOLT-LMM-inf statistic is not identical to any existing mixed model statistic because it is an approximate test statistic and avoids proximal contamination (Online Methods and Table 1). Nonetheless, we observed that BOLT-LMM-inf statistics very nearly match GCTA-LOCO statistics (which use the standard prospective model), with R²>0.999 (Supplementary Table 7 and Supplementary Fig. 7).

BOLT-LMM-inf mixed model statistic

The BOLT-LMM-inf infinitesimal mixed model statistic is slightly different:

where c_inf is a constant calibration factor estimated as

so that

In practice, for computational efficiency, we take means over 30 pseudorandom SNPs not significantly associated with the phenotype (χ²<5 estimated with the GRAMMAR statistic⁴³). We have observed empirically that 30 random SNPs are enough to estimate the calibration factor to within 1% (Supplementary Table 14).

We can view the BOLT-LMM-inf statistic either as an approximation of the standard prospective statistic (which treats phenotypes as random) or as a retrospective statistic (which treats genotypes as random and builds a null model on SNPs). The first perspective is motivated by the observation that in human genetics applications, the denominator of the prospective statistic in equation (5), x_test'V⁻¹x_test, is nearly independent of the SNP x_test being tested¹⁰. From this perspective, BOLT-LMM-inf is similar to GRAMMAR-Gamma¹⁰, with two key differences: (1) BOLT-LMM-inf is computed via much faster algorithms (described below) for performing initial variance parameter estimation and estimating the calibration constant, and (2) BOLT-LMM-inf avoids proximal contamination via LOCO analysis. Alternatively, we can also view BOLT-LMM-inf as a retrospective quasi-likelihood score test similar to T^SCORE–R (ref.⁴⁴) and MASTOR²³ (Supplementary Note).

BOLT-LMM Gaussian mixture model association statistic

We now generalize BOLT-LMM-inf by observing that the vector V_LOCO⁻¹y appearing in equation (8) is a scalar multiple of the residual phenotype vector σ_e²V_LOCO⁻¹y from best linear unbiased prediction (BLUP). Thus, the χ²_BOLT-LMM-inf statistic is equivalent to computing (and then calibrating) squared correlations between SNPs x_test and BLUP residuals. The power of mixed model association is driven by the fact that SNPs x_test are tested against these “de-noised” residual phenotypes from which other SNP effects estimated by the mixed model have been conditioned out^9,12.

We may generalize this approach by defining

where y_resid-LOCO denotes a generalized residual phenotype vector obtained after fitting a Gaussian mixture extension of the standard LMM (using SNPs not on the same chromosome as x_test) and c denotes a calibration factor, estimated so that the LD Score regression intercept²⁴ of χ²_BOLT-LMM matches that of the (properly calibrated) χ²_BOLT-LMM-inf statistic. Under the infinitesimal model, y_resid-LOCO is proportional to V_LOCO⁻¹y, so χ²_BOLT-LMM reduces to χ²_BOLT-LMM-inf. The general χ²_BOLT-LMM statistic can still be interpreted as a retrospective quasi-likelihood score test and is thus asymptotically chi-squared distributed.

To define the Gaussian mixture LMM extension, it is helpful to first frame the standard LMM in a Bayesian formulation. The null model of BOLT-LMM-inf is

where SNP effects β_m (m indexing SNPs not on the left-out chromosome) are independently drawn from the Gaussian prior distribution

and environmental effects e_n (n indexing samples) are independently drawn from e_n ~ N(0, σ_e²). Performing best linear unbiased prediction amounts to computing the posterior mean of the genetic effect X_LOCOβ_LOCO.

To generalize this model to non-infinitesimal genetic architectures, we replace the Gaussian prior on SNP effect sizes with a more general distribution; this approach has been extensively applied by the “Bayesian alphabet” of genomic prediction methods in the animal breeding literature^17–19. In BOLT-LMM, we use a spike-and-slab mixture of two Gaussians¹⁹ as the prior:

This mixture more flexibly models the heavier-tailed distributions of genetic effects of typical (non-infinitesimal) phenotypes. Explicitly, if p1 and σ_β,1² σ_β,2², the first component of the mixture is a “slab” that models the existence of a small number of relatively large-effect loci, while the second component is a “spike” that models the assumption that most SNPs have near-zero—but not exactly zero—effect on the phenotype. (Note, however, that all SNPs are assigned the same mixture prior; i.e., SNPs are not individually allocated to one or the other component.) It is important that the spike component have nonzero variance so as to capture genome-wide effects on phenotype such as ancestry or relatedness; then, when testing SNPs for association, these genome-wide effects are conditioned out from residual phenotypes, protecting against confounding. The prior could in principle be further generalized; we chose to use a mixture of two Gaussians to keep the model fairly simple and because Gaussian distributions produce convenient analytical formulas during model-fitting.

Under this generalized model, posterior means no longer correspond to BLUP, but we can still approximately fit the Bayesian model (once per left-out chromosome) and obtain residuals

where β_LOCO are estimated posterior mean effect sizes. Plugging these residuals into equation (11) gives the BOLT-LMM Gaussian mixture model association test statistic.

Fast iterative algorithm

The BOLT-LMM software performs a four-step computation for mixed model association analysis, stopping after the first two steps when specialized to the infinitesimal model. We outline the algorithm here and provide full details and pseudocode in the Supplementary Note.

Step 1a: Estimate variance parameters

A key feature of BOLT-LMM is estimation of variance parameters σ_g² and σ_e² using only linear-time iterations without building or decomposing any covariance matrices. We use a Monte Carlo REML approach^26,27 that eliminates all O(MN²) and O(N³)-time matrix computations, requiring only the solution of linear systems of mixed model equations. We solve the mixed model equations using conjugate gradient iteration, which requires only O(MN)-time matrix-vector products^28,29 (Supplementary Note).

Step 1b: Compute and calibrate BOLT-LMM-inf statistics

Having variance parameter estimates from Step 1a, it is straightforward to compute (for each LOCO rep) the quantity V_LOCO⁻¹y in the numerator of the BOLT-LMM-inf statistic, equation (8), using conjugate gradient iteration as above. Completing the computation of the numerator of χ²_BOLT-LMM-inf then just amounts to calculating one dot product per SNP x_test, which requires only O(MN) additional cost across all SNPs. Moreover, this computation can easily be performed for additional SNPs not included in the mixed model but at which association statistics are desired; BOLT-LMM handles imputed “dosage” data in this way. To compute the calibration constant c_inf in equation (9), BOLT-LMM rapidly computes the prospective statistic χ²_LMM-LOCO from equation (7) at 30 random SNPs by applying conjugate gradient iteration to compute V_LOCO⁻¹x_test for each of the 30 selected SNPs x_test. Finally, in addition to computing χ² association statistics, BOLT-LMM also computes effect size estimates for all SNPs tested (Supplementary Note).

There is a slight mismatch between the variance parameters estimated in Step 1a, which BOLT-LMM computes once using all SNPs—not leaving any chromosomes out—and the theoretically optimal parameter estimates that would be obtained by refitting once per left-out chromosome. However, we have observed in simulations that slight mis-specification of the variance parameters has a negligible impact (<0.5%) on the calibration of the BOLT-LMM-inf and BOLT-LMM statistics (Supplementary Table 4). Because very slight miscalibration is not a concern for confounding from population stratification at highly differentiated markers (Supplementary Table 12) and has little impact on Type I error (Supplementary Table 5), the BOLT-LMM software does not by default refit variance parameters for each LOCO rep. If extremely precise calibration is desired, we provide a runtime option to refit variance parameters for each LOCO rep, at the cost of a factor of 2–3 in running time. We believe that LOCO strikes a good balance in terms of achieving ≈95% of the potential power gain (by jointly fitting ≈95% of markers that are not in LD with the candidate marker) while keeping run time down¹², but we also provide a runtime option to partition the genome more finely (e.g., into 100 segments rather than 22), again at the cost of a factor of 2–3 in running time.

Step 2a: Estimate Gaussian mixture prior parameters

The first step of BOLT-LMM Gaussian mixture model association analyis is to estimate parameters of the generalized prior on SNP effect sizes. As written in equation (14), this mixture has three parameters: σ_β,1² and σ_β,2², the variances of the two Gaussians, and p, the probability of drawing from the first Gaussian. To reduce the complexity of parameter estimation, we constrain the total variance of the mixture to equal the variance σ_g²/M estimated under the infinitesimal model in Step 1a:

We reparameterize the remaining two degrees of freedom using the parameters p and f₂, where f₂ denotes the proportion of the total mixture variance within the second Gaussian (the “spike” component that models small genome-wide effects):

Because the model fit is insensitive to the precise values of the mixture parameters, we test a discrete set of model parameter combinations: f₂{0.5,0.3,0.1}, p{0.5,0.2.0.1,0.05,0.02,0.01}. Note that f₂=0.5, p=0.5 corresponds to the infinitesimal model: when f₂=1−p, the two Gaussians are identical and the mixture is degenerate. We bound f₂ from below to ensure that at least a small amount (10%) of the mixture variance is assigned to the spike component, protecting against confounding from genome-wide effects. We bound p from below to prevent the model from trying to fit too strongly to a few SNPs, which makes model-fitting computationally difficult and also increases susceptibility to confounding. BOLT-LMM performs model selection among the 18 possible parameter pairs (f₂, p) by performing cross-validation to optimize mean-squared prediction R².

BOLT-LMM uses a variational approximation to fit Bayesian linear regressions with Gaussian mixture priors. Approximation methods are necessary for Bayesian inference in this setting because exact posterior means involve intractable integrals. We apply a fully factored variational approximation^21,22,38 that repeatedly loops through the SNPs, updating the estimated effect size of each SNP with its posterior mean conditional on current estimates of all other SNP effects. This iteration has also previously been termed “iterative conditional expectation (ICE)”²⁰. The variational Bayes framework puts this iteration on a sound theoretical footing as an optimization of an approximate log likelihood function; the iteration monotonically increases this function and is guaranteed to converge⁴⁵. In fact, we show that the optimization can be reformulated as cyclic coordinate descent applied to a penalized regression problem arising from Bayesian linear regression using a transformed prior (Supplementary Note). The approximate log likelihood also serves as a convenient convergence criterion: BOLT-LMM stops the iteration when the increase in approximate log likelihood over one full update cycle drops below 0.01.

While the core variational iteration that BOLT-LMM uses is identical to previous methods^20–22,38 up to the choice of SNP effect size prior, BOLT-LMM uses cross-validation to estimate hyperparameters¹⁵ rather than doing so within the variational iteration^22,38 or based on variational approximate log likelihoods²¹. We found this approach to be more robust to slackness of the variational approximation caused by linkage disequilibrium.

Step 2b: Compute and calibrate BOLT-LMM Gaussian mixture model statistics

After inferring parameters of the mixture prior in Step 2a, BOLT-LMM uses the same variational iteration to estimate posterior mean residuals y_resid-LOCO (independently for each left-out chromosome). The numerators of the BOLT-LMM Gaussian mixture model statistic from equation (11) are then easily obtained as dot products with test SNPs, leaving only the constant calibration factor c in the denominator to be calculated. Unlike the case of the infinitesimal model, here we do not have a prospective statistic to calibrate against, so we instead apply LD Score regression²⁴ (Supplementary Note). In practice, the calibration factor is usually quite close to 1 (e.g., 1.00 to two decimal places for all WGHS traits; see Supplementary Table 15).

WGHS data set

The Women’s Genome Health Study (WGHS) is a prospective cohort of initially healthy, female North American health care professionals. We analyzed 23,294 individuals with self-reported European ancestry with genotyping at 324,488 SNPs after QC (Supplementary Note).

We are grateful to M. Lipson, S. Simmons, A. Gusev, K. Galinsky, J. Yang, P. Visscher, Z. Zhu, and D. Gudbjartsson for helpful discussions. This research was supported by NIH grant R01 HG006399 and NIH fellowship F32 HG007805. H. K. F. was supported by the Fannie and John Hertz Foundation. The WGHS is supported by HL043851 and HL080467 from the National Heart, Lung, and Blood Institute and CA047988 from the National Cancer Institute, the Donald W. Reynolds Foundation and the Fondation Leducq, with collaborative scientific support and funding for genotyping provided by Amgen.