Stochastic Lanczos estimation of genomic variance components for linear mixed-effects models

Abstract

Background: Linear mixed-effects models (LMM) are a leading method in conducting genome-wide association studies (GWAS) but require residual maximum likelihood (REML) estimation of variance components, which is computationally demanding. Previous work has reduced the computational burden of variance component estimation by replacing direct matrix operations with iterative and stochastic methods and by employing loose tolerances to limit the number of iterations in the REML optimization procedure. Here, we introduce two novel algorithms, stochastic Lanczos derivative-free REML (SLDF_REML) and Lanczos first-order Monte Carlo REML (L_FOMC_REML), that exploit problem structure via the principle of Krylov subspace shift-invariance to speed computation beyond existing methods. Both novel algorithms only require a single round of computation involving iterative matrix operations, after which their respective objectives can be repeatedly evaluated using vector operations. Further, in contrast to existing stochastic methods, SLDF_REML can exploit precomputed genomic relatedness matrices (GRMs), when available, to further speed computation.Results: Results of numerical experiments are congruent with theory and demonstrate that interpreted-language implementations of both algorithms match or exceed existing compiled-language software packages in speed, accuracy, and flexibility. Conclusions:Both the SLDF_REML and L_FOMC_REML algorithms outperform existing methods for REML estimation of variance components for LMM and are suitable for incorporation into existing GWAS LMM software implementations.

Assortative Mating Biases Marker-based Heritability Estimators

“…To our knowledge, existing software is incapable of efficient REML analysis using adjusted genotypes (analogous to dosages) in large samples; e.g., BOLT-REML requires hard-calls as input, whereas GCTA and LDAK have cubic complexity in the number of individuals and markers and would require multiple weeks to run on a high thread-count server. We therefore utilized a modified Python implementation of the REML algorithm presented in [31] (code available by request). We used LDAK 5.0 to obtain adjusted HE regression estimates [32].…”

“…We used GCTA v1.91.3b [14] to construct genomic related matrices and perform HE regression. We obtained REML heritability estimates using BOLT-LMM v2.3.4 [30] for computational efficiency; though BOLT-LMM uses a randomized algorithm, its numerical accuracy is comparable to that of the exact algorithm implemented GCTA [31].…”

“…We run on a high thread-count server. We therefore utilized a modified Python implementation of the REML algorithm presented in [31] (code available by request). We used LDAK 5.0 to obtain adjusted HE regression estimates [32].…”