  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
