PaperlayerView exact revision on arXiv

Source-linked AI summary

Efficient Algorithms for Multivariate Linear Mixed Models in Genome-wide Association Studies

Xiang Zhou, Matthew Stephens

arXiv:1305.4366v2q-bio.QMstat.APstat.CO

TL;DR

Computing likelihood ratio tests for multivariate linear mixed models in GWAS is difficult, and existing approaches are limited to approximate tests for two phenotypes. The paper develops efficient algorithms and GEMMA software for fitting mvLMMs and testing associations, making such analyses practical for large samples and several phenotypes. The resulting LRT provides correct null type-I-error control and is more powerful than MTMM's approximate LRT.

  • Problem

    Existing mvLMM methods were not computationally practical for GWAS LRTs, while MTMM provided only an approximate test for two phenotypes.

  • Method

    The paper develops computationally efficient algorithms for fitting mvLMMs and computing GWAS association LRTs, with an implementation in GEMMA.

  • Results

    The LRT provides correct control of type I error and is more powerful than MTMM's approximate LRT.

  • Takeaways & Limitations

    The method provides a computationally practical approach to mvLMM GWAS LRTs for reasonably large samples and more than two phenotypes.

  • Takeaways & Limitations

    The algorithms remain computationally intractable in practice for large n, such as n > 50,000, motivating low-rank approximations.

Abstract

from arXiv · show

Multivariate linear mixed models (mvLMMs) have been widely used in many areas of genetics, and have attracted considerable recent interest in genome-wide association studies (GWASs). However, fitting mvLMMs is computationally non-trivial, and no existing method is computationally practical for performing the likelihood ratio test (LRT) for mvLMMs in GWAS settings with moderate sample size n. The existing software MTMM perform an approximate LRT for two phenotypes, and as we find, its p values can substantially understate the significance of associations. Here, we present novel computationally-efficient algorithms for fitting mvLMMs, and computing the LRT in GWAS settings. After a single initial eigen-decomposition (with complexity O(n^3)) the algorithms i) reduce computational complexity (per iteration of the optimizer) from cubic to linear in n; and ii) in GWAS analyses, reduces per-marker complexity from cubic to quadratic in n. These innovations make it practical to compute the LRT for mvLMMs in GWASs for tens of thousands of samples and a moderate number of phenotypes (~2-10). With simulations, we show that the LRT provides correct control for type I error. With both simulations and real data we find that the LRT is more powerful than the approximate LRT from MTMM, and illustrate the benefits of analyzing more than two phenotypes. The method is implemented in the GEMMA software package, freely available at http://stephenslab.uchicago.edu/software.html

Introduction

mvLMMs are increasingly relevant to GWAS because they model correlated phenotypes and can improve detection of genetic variants, but their LRTs are computationally difficult. The paper introduces efficient algorithms and software that make multivariate GWAS testing practical for large samples and several phenotypes.

  • mvLMMs can increase power to detect variants affecting multiple correlated phenotypes or only one phenotype in the collection.
  • Fitting mvLMMs requires multidimensional optimization of a potentially non-convex likelihood or restricted likelihood.
  • GWAS applications require these optimizations to be repeated for hundreds of thousands of SNPs, making existing mvLMM LRT methods impractical.
  • MTMM provides only an approximate LRT for two phenotypes, and its p values can substantially understate association significance.
  • The paper presents computationally efficient mvLMM-fitting and GWAS-LRT algorithms, implemented in GEMMA, that are practical for about 50,000 individuals and a modest number of phenotypes.

Results

The paper develops computationally efficient mvLMM algorithms and implements them in GEMMA, making GWAS likelihood-ratio testing practical at larger scales. The algorithms reduce computational burden, produce well-calibrated LRT p values, and can improve power relative to MTMM while supporting analyses with more than two phenotypes.

  • Computational contribution: The algorithms provide a computationally practical approach for computing mvLMM likelihood-ratio tests in GWAS, including analyses with more than two phenotypes.The software is described as the first available implementation of an mvLMM test for GWAS with more than two phenotypes.
  • Computational contribution: Existing methods require O(s(t1n^3d^3 + t2n^3d^7)) complexity when repeated across SNPs, which is impractical for large s and moderate n.The expensive estimation procedure must otherwise be repeated for all SNPs.
  • Computational contribution: After one O(n^3) eigen-decomposition of the relatedness matrix, per-SNP complexity increases only quadratically in n.The overall complexity is O(n^3 + n^2d + s(n^2 + t1nd^2 + t2nd^6)).
  • Computational performance: For NFBC1966 with d = 4, GEMMA took about 7 minutes compared with 8 hours for WOMBAT when fitting a single mvLMM.The reported comparison was performed for parameter estimation, and the gains were expected to increase for larger d.
  • Statistical performance: GEMMA LRT p values were well-calibrated under the null, whereas MTMM p values were systematically larger than expected, with extreme values almost an order of magnitude larger.The simulations also reported that GEMMA had greater power than MTMM under alternatives.
  • Comparison with MTMM: GEMMA ran 2-12 times faster than MTMM in the reported GWAS analyses, while MTMM p values were up to 6 fold less significant for HMDP.For NFBC1966, GEMMA completed a two-phenotype analysis in about four hours versus almost two and a half days for MTMM.
  • Analyzing more than two phenotypes: Four-phenotype analyses were consistently at least as powerful as two-phenotype analyses, including settings where only one or two phenotypes were associated.Power could increase when unassociated phenotypes were correlated with associated phenotypes.
  • Analyzing more than two phenotypes: In real data, 16 SNPs were significant only in the four-phenotype analysis, compared with 3 significant only in the univariate analysis.Across the analyses, 45 SNPs from 14 genetic regions passed the Bonferroni-adjusted significance level of 0.05.

Discussion

The paper presents GEMMA algorithms and software that make mvLMM likelihood-ratio testing practical for GWAS, while improving calibration over approximate methods. It also identifies computational, statistical, and data-related boundaries for applying the approach.

  • Contribution: The algorithms reduce repeated computation after an initial eigen-decomposition, enabling multivariate GWAS analyses with moderate phenotype counts.The paper describes reduced computational burden after the eigen-decomposition and applicability to analyses with several phenotypes.
  • Contribution: GEMMA provides a computationally practical mvLMM likelihood-ratio test for GWAS and the first software implementation supporting more than two phenotypes.The implementation is provided in GEMMA, and the method is positioned as the first practical GWAS approach for this test.
  • Statistical performance: GEMMA LRT p values are well-calibrated under the null, unlike p values from existing approximate-LRT methods.In the reported simulations, MTMM p values were systematically larger than expected, with the most significant values almost an order of magnitude larger than expected.
  • Limitations: The method is limited to mvLMMs with one variance component besides residual error, and missing phenotypes require dropping individuals or imputation.The paper also notes that larger phenotype counts can create computational and statistical barriers, including quadratic growth in variance-component parameters.
  • Limitations: For large n, the initial eigen-decomposition can become practically intractable; low-rank relatedness approximations may reduce computation and memory demands.The paper gives n > 50,000 as an example boundary and identifies low-rank approximations as a possible remedy.

Methods

The study analyzes cleaned and transformed genotype–phenotype data from HMDP and NFBC1966, then uses simulations to evaluate calibration, power, missing-phenotype handling, and multivariate phenotype counts.

  • Data preparation: The analyses use four metabolic phenotypes in each dataset and estimate relatedness from centered genotype matrices.HMDP contains 98 strains, 656 individuals, and 108,562 SNPs after filtering; NFBC1966 contains 5,255 individuals and 319,111 SNPs.
  • Data preparation: Each phenotype is quantile-transformed to a standard normal distribution, with covariate adjustment applied in NFBC1966.The transformations were used to guard against model misspecification; NFBC1966 additionally regressed out sex, oral contraceptives, and pregnancy status.
  • Evaluation: Statistical power is measured as the proportion of p values exceeding Bonferroni-corrected genome-wide significance thresholds.The thresholds are p = 4.6 × 10^-7 for HMDP and p = 1.6 × 10^-7 for NFBC1966.
  • Evaluation: Missing-phenotype simulations randomly remove one phenotype for 2.5%, 5%, or 10% of individuals to compare dropping with imputation.Both approaches are evaluated using GEMMA on the simulated phenotypes.
  • Evaluation: Four-phenotype analyses are compared with the minimum p value from six pairwise analyses when simulated SNPs affect one to four phenotypes.The comparison includes effects with concordant or discordant directions across affected traits and uses corrected significance levels.

1.1 Multivariate Linear Mixed Model

The paper formulates mvLMMs for multiple phenotypes and individuals using genetic random effects, residual errors, covariates, and a relatedness matrix. Eigen-decomposition transforms the model so individuals become independent but non-identically distributed multivariate normal observations.

  • Model specification: The mvLMM models d phenotypes for n individuals with covariates, marker effects, genetic random effects, and residual errors.The genetic effects use covariance Vg and relatedness matrix K, while residual errors use covariance Ve and the identity matrix.
  • Model specification: The genetic variance component Vg and environmental variance component Ve are d by d symmetric matrices.The model uses a matrix normal distribution with row covariance across phenotypes and column covariance across individuals.
  • Model transformation: Eigen-decomposition of the relatedness matrix K transforms phenotypes, covariates, and random effects into a computationally convenient representation.The transformation uses eigenvectors Uk and eigenvalues in diagonal matrix Dk.
  • Model transformation: After transformation, the model is Y = BX + G + E, with G and E following matrix normal distributions involving Dk and the identity matrix.Vectorization gives an equivalent multivariate normal formulation with Kronecker-product covariance matrices.
  • Individual-level representation: For each individual l, transformed phenotypes follow independent but non-identical multivariate normal distributions with variance Vl = δlVg + Ve.The individual-specific variance depends on the corresponding eigenvalue δl.

1.2 Optimization Method Overview

Existing mvLMM optimization methods become impractical as sample size and phenotype number grow, especially when optimization is repeated for every SNP. The paper introduces transformations and matrix-structure exploits that reduce fitting and per-SNP costs while combining stable PX-EM with faster Newton-Raphson refinement.

  • Computational challenge: Derivative-free methods become impractical as phenotype number increases because their time complexity grows exponentially with the number of parameters.The cited original derivative-free mvLMM work showed examples only for two phenotypes.
  • Algorithmic strategy: The algorithms combine univariate LMM tricks, simultaneous diagonalization, block-diagonal matrices, and sparse-matrix properties.They provide a multivariate analogue of EMMA and FaSTLMM/GEMMA/CM.
  • Complexity reductions: O(nd^2) replaces O(n^3d^3) per EM iteration, while O(nd^6) replaces O(n^3d^7) per Newton-Raphson iteration.These reductions apply to the multivariate algorithms after the transformation-based computations.
  • Complexity reductions: O(n^2) replaces O(n^3d^3) per SNP, or O(n) when K has low rank.The per-marker reduction makes the multivariate GWAS computation substantially more scalable.
  • Optimization strategy: PX-EM supplies stable optimization but can be slow, whereas Newton-Raphson converges quickly from a good starting point but can fail from a poor one.The implementation therefore combines PX-EM initialization with Newton-Raphson refinement.

1.3 PX-EM Algorithms

The PX-EM framework estimates mvLMM parameters by treating random effects, and in one variant fixed effects, as missing data. Conditional expectation and maximization updates are combined with parameter expansion and transformations that make each iteration computationally efficient.

  • Algorithm overview: The paper describes ECM for maximum likelihood estimation and EM for restricted maximum likelihood estimation in mvLMMs.Parameter-expanded versions of both algorithms are also developed with efficient computations.
  • EM updates: The EM procedure treats G as missing values and uses its conditional distribution given Y and current parameter values.The expected log likelihood is then optimized through conditional maximization updates.
  • ECM updates: ECM updates B conditional on the current variance components and then updates Vg and Ve using conditional expectations.The variance updates incorporate conditional estimates and covariance terms for the random effects and residual errors.
  • Parameter expansion: The PX versions introduce a latent d by d matrix Va so that Vg = VaVg* in a reparameterized model.The expanded parameter is updated using expectations under the conditional distribution of the latent random effects.
  • Efficient computation: The resulting PX-ECM/EM iteration complexity is reduced through efficient evaluation of the transformed individual-level quantities.The method exploits the transformed model and avoids repeatedly handling dense covariance structures directly.
  • Efficient computation: A second transformation converts correlated traits into uncorrelated ones, complementing the transformation that converts correlated individuals into uncorrelated ones.This canonical transformation avoids cubic computation in the number of traits.

1.4 Newton-Raphson’s Algorithms

The Newton-Raphson algorithms estimate mvLMM likelihood or restricted likelihood parameters using analytically derived derivatives and efficient matrix operations. Block-diagonal and sparse structures reduce the cost of evaluating these derivatives and associated target-function quantities.

  • Algorithm overview: Newton-Raphson algorithms are presented for both maximum likelihood and restricted maximum likelihood estimation in mvLMMs.The paper reports that the average-information algorithm is unnecessary when these efficient algorithms are used.
  • Derivative calculations: The log-likelihood and log-restricted likelihood are expressed together with first- and second-order partial derivatives for variance components.Derivatives are specified with respect to elements of Vg and Ve.
  • Derivative calculations: The derivative notation defines variance-component elements, indicator vectors, and matrices used to represent individual covariance parameters.These definitions support the derivative expressions for Vg and Ve.
  • Efficient computation: Efficient calculations evaluate trace terms and vector-matrix-vector products using quantities such as Q, q, and q derived from H^-1 and the data.The target functions and first-order derivatives are reduced to a set of efficiently evaluated quantities.
  • Efficient computation: The key computational trick uses Iij = IiIj^T so that many quantities require only scalar multiplications or rank-one matrix updates.The calculations also use block-diagonal and sparse-matrix properties.
  • Efficient computation: The most time-consuming Newton-Raphson calculation involves Qgg terms, with each requiring O(nc^2d^2) computation time.This identifies the dominant cost in the derivative evaluation described by the paper.

1.5 Test Statistics and p Values

The section defines likelihood ratio, Wald, and score tests for mvLMMs, including their asymptotic calibration and corrections. It reports that the likelihood ratio test provides correct type I error control relative to the other tests.

  • Likelihood ratio test: The likelihood ratio test compares maximized alternative and null log likelihoods using zLR = 2(l1 − l0).The alternative and null models use their respective estimated variance components.
  • Wald test: The Wald test estimates marker effect sizes using variance-component estimates from the alternative model.The transformed genotype and covariance quantities enter the Wald statistic.
  • Score test: The score test computes a d(d + 1)-dimensional score vector using derivatives and the observed information matrix evaluated under the null.These quantities are used to construct the score statistic.
  • Calibration: Under the null hypothesis, all three test statistics asymptotically follow a χ2(d) distribution.When sample size is small or relatedness is strong, the asymptotic distribution may not calibrate p values exactly.
  • Calibration: Marker-specific corrections require alternative-model estimates and partial derivatives; empirically, score tests are often conservative and Wald tests often anti-conservative.The likelihood ratio test behaves between them and has correct type I error control.

1.6 Phenotype Imputation

The method imputes missing phenotypes before association testing so partially observed individuals are not discarded. It estimates null-model parameters from complete cases and uses conditional means for imputation.

  • Motivation: The algorithms require complete or imputed phenotypes, motivating imputation when individuals have partially missing measurements.This avoids dropping individuals with incomplete phenotypes.
  • Method: GEMMA imputes missing phenotypes before association tests.The procedure is designed to retain information from individuals with partially missing phenotypes.
  • Method: The method first estimates null-model fixed effects and covariance components using individuals with completely observed phenotypes.The estimated parameters are then used for imputation.
  • Method: Missing values are estimated by the conditional mean given observed phenotypes and the estimated parameters.Observed and missing values are represented separately in the covariance formulation.

2 Supplementary Results

The supplementary materials describe comparisons of testing algorithms, phenotype counts, missing-phenotype handling, and significant SNPs under multiple-phenotype analyses. They specify the thresholds and encodings used for these comparisons.

  • Figure S1: Figure S1 compares -log10 p values from GEMMA using PX-EM plus NR with MTMM or GEMMA using PX-EM alone across paired HMDP traits.Black denotes the MTMM comparison and red denotes the PX-EM-only comparison.
  • Figure S2: Figure S2 compares power between four-phenotype and two-phenotype GEMMA analyses across PVE, affected-phenotype count, datasets, and multiple-correction thresholds.The four-phenotype thresholds are 4.6 × 10−7 for HMDP-based simulations and 1.6 × 10−7 for NFBC1966-based simulations.
  • Figure S3: Figure S3 compares dropping incomplete individuals with imputing missing phenotypes as the percentage of individuals with one missing phenotype increases.Symbols and line types encode effect direction and PVE for the second phenotype.
  • Table S1: Table S1 lists SNPs significant after Bonferroni correction in four-phenotype mvLMM or univariate LMM analyses of NFBC1966 data.It distinguishes SNPs significant only in one analysis and those more significant in the four-phenotype analysis after test correction.
  • Table S2: Table S2 lists SNPs significant after Bonferroni correction in four-phenotype or two-phenotype mvLMM analyses of NFBC1966 data.No SNP is significant only in the two-phenotype analysis.
Loading 1305.4366v2…