Methods for Analyzing Multivariate Phenotypes in Genetic Association Studies

Multivariate phenotypes are frequently encountered in genetic association studies. The purpose of analyzing multivariate phenotypes usually includes discovery of novel genetic variants of pleiotropy effects, that is, affecting multiple phenotypes, and the ultimate goal of uncovering the underlying genetic mechanism. In recent years, there have been new method development and application of existing statistical methods to such phenotypes. In this paper, we provide a review of the available methods for analyzing association between a single marker and a multivariate phenotype consisting of the same type of components e.g., all continuous or all categorical or different types of components e.g., some are continuous and others are categorical . We also reviewed causal inference methods designed to test whether the detected association with the multivariate phenotype is truly pleiotropy or the genetic marker exerts its effects on some phenotypes through affecting the others.


Introduction
Association studies, where the correlation between a genetic marker and a phenotype is assessed, are useful for mapping genes influencing complex diseases. With reduction of genotyping cost, completion of the HapMap Project 1 , and more recently the 1000 Genomes Project 2 , genome-wide association studies GWAS with several hundred thousands to tens of millions genotyped and/or imputed single nucleotide polymorphisms SNPs have become a common approach nowadays to search for genetic determination of complex traits.
In the study of complex diseases, several correlated phenotypes, or a multivariate MV phenotype with several components, may be measured to study a disorder or trait. For example, hypertension is evaluated using systolic and diastolic blood pressures; 2 Journal of Probability and Statistics a person's cognitive ability is usually measured by tests in domains including memory, intelligence, language, executive function, and visual-spatial function. The tests within and between domains are correlated. Most published GWAS only analyzed each individual phenotype separately, although results on related phenotypes may be reported together. Published single phenotype GWAS have successfully identified a large number of novel genetic variants predisposing to a variety of complex traits 3, 4 . However, majority of the identified genetic variants only explain a small fraction of total heritability defined as between individual phenotype variability attributable to genetic factors 4, 5 . It has been hypothesized that current GWAS may be underpowered to detect many genetic variants of moderate-to-small effects. Joint analysis of correlated phenotypes can exploit the correlation among the phenotypes, which may lead to better power to detect additional genetic variants with small effects across multiple traits or pleiotropy effects. Furthermore, joint analysis avoids multiple testing penalty incurred in analyzing each phenotype separately. Therefore, it is important to identify appropriate methods that fully utilize information in multivariate phenotypes to detect novel genetic loci in genetic association studies.
In addition to discovery of novel loci of potential pleiotropy effects, it is also important to detangle the complex relationship between phenotype components and genetic variants One of the frequently asked questions is whether a genetic variant affects multiple phenotypes simultaneously pleiotropy or affects one phenotype through affecting another phenotype. In this paper, we review methods for both purposes.

Methods for Detecting Association Using Multivariate Phenotypes
For all the methods mentioned in this section, the null hypothesis is no association between a single genetic marker and any components of a multivariate MV phenotype; the alternative hypothesis is the genetic maker associated with at least one phenotype component. Here we review methods for an MV phenotype consisting of all continuous, all categorical, or all timeto-event components, and methods for MV phenotypes consisting of a mixture of different types of components.

Regression Models
Regression models for clustered observations such as linear and generalized mixed effects models, generalized estimating equations, and frailty models can be used to analyze the association of a genetic marker with all continuous, categorical, or survival multivariate phenotypes.

Mixed Effects Models
Mixed effects models such as linear mixed effects model LME and generalized linear mixed effects model GLMM involve using fixed effects for the genetic marker effect and random effects to account for correlation among multivariate phenotypes 6, 7 .
Let y jk denote the kth k 1, . . . , K continuous component of the K-dimensional phenotype of the jth j 1, . . . , J individual. Let g j be the genotype of a genetic marker of Journal of Probability and Statistics 3 the jth individual, and X g j a score of the genotype. The linear mixed effects model takes the following form: where β 0 is the intercept or other genetic or environmental fixed effects; β k is the fixed effect size of X g j on the kth phenotype; η jk k 1, . . . , K ∼ N 0, Σ are the random effects correlated within jth person; e jk is the random errors iid. ∼ N 0, σ 2 e . Between any two individuals, η jk , k 1, . . . , K are independent. Within a person, η jk , k 1, . . . , K are correlated. The null hypothesis that the genetic marker is not associated with any phenotype component corresponds to H 0 : β 1 , . . . , β K 0. The estimation of variance parameters and fixed effect parameters can be obtained using restricted maximum likelihood method REML 8,9 . When y jk is categorical, it can be modeled with generalized mixed effects model GLMM as follows: where μ is a link function and μ −1 is its inverse. For Gaussian distributed traits, μ is the identity link, thus 2.2 is identical to the linear mixed effects model 2.1 ; for binary traits, μ is the logit link μ x ln x/1 − x . For links other than identity function, the likelihood for this model contains integrals without a close form solution. All existing algorithms for likelihood maximization are either based on theoretical or numerical approximation 10, 11 .
The null hypothesis under the LME or GLMM can be tested using the likelihood ratio test or Wald chi-squared test. They can be implemented using SAS PROC Mixed or R lme4 package function lmer(). The Wald chi-squared test statistic takes the form β T cov β −1 β ∼ χ 2 K , where β β 1 , . . . , β K is estimated using 2.1 or 2.2 . For example, Kraja et al. 12 have employed a model similar to 2.1 to the analyses of bivariate continuous metabolic traits. We can also fit a model assuming β 1 · · · β K β, that is, E y jk | η k μ −1 β 0 βX g j η jk , where a single degree-of-freedom df test β/se β can be used to test the null hypothesis. This test can be more powerful than the multi-df Wald chi-squared test if the effect sizes are in the same direction and not very different. It, however, may lack power if the β 1 , . . . , β K are very different, especially have different signs and cancel each other out.

Frailty Models
When the phenotypes are correlated survival times, frailty models can be used to fit the association model. Suppose the survival or censoring times are t kj for the kth k 1, . . . , K phenotype of the jth j 1, . . . , J individual. Let g j be the genotype of a genetic marker of the jth individual, and X g j a score of the genotype as follows: where η kj j 1, . . . , J are subject specific random effects following N 0, Σ , and Σ is a Kdimensional correlation matrix. This is the Gaussian frailty model. There is another class of frailty models where exp η kj follows a gamma distribution. A Gaussian or gamma frailty model assuming an exchangeable correlation within a person can be fitted using coxph() in the survival package of R by including a frailty() term in the regressor. In addition, including a cluster() term in coxph() fits generalized estimating equations GEE type of model that assumes an independent working correlation matrix 13 . Frailty models with an arbitrary prespecified Σ can be fitted with the coxme() in R coxme package for Gaussian random effects model.
Fitting a mixed effects frailty model requires predetermining the correlation matrix Σ of random effects η jk within jth person. The correlation between the phenotypes y jk within a person is attributable to the random effects η jk and the fixed effects of the genetic marker. However, since the fixed effects are unknown, it is impossible to directly infer the correlation among the random effects. Misspecifying the correlation among random effects may result in bias in the inference on fixed effects. But the bias seems to be small for genetic association studies 14, 15 .

Generalized Estimating Equations
Different from mixed effects model is a class of models called marginal models. Instead of having random effects as regressors in addition to random errors to model correlation in multivariable response, marginal models collapse the random effects and random residual errors in the model. Generalized estimating equations GEE 16 solve the quasi-likelihood score function as follows: where V j A 1/2 j R α A 1/2 j , and R α is the working correlation matrix for the residual correlation. The variance and covariance of β is estimated with the so-called robust variance estimator 16 . Similar to the LME, single-or multi-df Wald test statistic can be usually used to test that the genetic marker is not associated with any of the phenotypes.
In our experience, GEE results are inflated with low minor frequency SNPs and not as powerful as LME in general 15, 17 . However, GEE is robust to misspecification of response distribution or association model and thus can be used when the LME shows bias or inflation due to these reasons.

Variable Reduction Method
Variable reduction approaches are in general only applicable to MV phenotype consisting of all continuous phenotypes that are approximately normal distributed. It derives a single or a few new phenotypes that are linear combinations of the original phenotypes, for example, Existing methods include principal components analysis PCA where for the first component, a i , i 1, . . . , K are coefficients that maximize the variance of Y ; principal component of heritability PCH with coefficients maximizing the total heritability of Y 18 Journal of Probability and Statistics 5 and penalized PCH applicable to high-dimensional data 19, 20 ; and principal components of heritability with coefficients maximizing the quantitative trait locus QTL heritability PCQH of Y 21-24 , that is, the variance explained by the genetic marker. The PCQH approaches are designed to maximize the individual phenotype variation explained by the genetic marker and thus may be more powerful than PCA and PCH in genetic association studies.

PCQH Approaches
The approaches proposed by Lange et al. 21,25 and Klei et al. 23 involve using a subset of the sample to estimate the coefficients in 2.5 that maximize the correlation between Y and the genetic marker. Specifically, in the estimation sample, the total phenotype variance is partitioned into QTL variance and residual variance as follows: A that maximizes h 2 A can be obtained by solving the following generalized eigen system 18 : . . , |β K |σ x , σ x is the sample standard deviation of the score of genotype X g across all individuals, β i is estimated using the least squared estimator of Y i α β i X g ε, and 1 sign β 1 , . . . , sign β K .
Lange et al. 21, 25 approaches are only applicable to family-based association design. They suggest using the noninformative families or parental genotypes to estimate A because these data will not contribute directly to the family-based association tests FBAT . Then perform FBAT of Y on X g . However, FBAT has low power in the absence of population stratification 26 compared to population based approaches. Klei et al's. 23 is a populationbased association approach where they randomly split the sample into two subsets: one used to estimate A, the other used to test the association of Y with X g via a linear regression model: Y α βX g ε. This ensures valid P value in the association test.

Canonical Correlation Analysis
Canonical correlation analysis seeks coefficients so that the squared correlation between Y in 2.5 , and the score of genetic marker, X g , is maximized. Here ρ corr Y , X is called where Σ Y Y is the K × K matrix of the variance-covariance matrix of Y , Σ Y X and its transpose Σ XY are K × 1 and 1 × Kmatrix of the covariance matrix between Y and X, Σ XX is the variance of X, a scalar. All these submatrices can be estimated using the respective sample co-variance matrix. The canonical correlation, These tests are implemented in SAS PROC GLM and R function summary.manova(). As part of the PLINK package specifically developed for genetic analysis, Ferreira et al. 24 implemented the Wilks lambda, and its P value is obtained from F-approximation Canonical correlation analysis shares similarity with PCQH 23 in that both estimate a linear combination of original phenotypes, so that the genotype score explains most of the variation in terms of percent of total variance and squared correlation, resp. of the new phenotype. The difference between the two approaches is that the canonical correlation analysis evaluates squared correlation using whole sample, while PCQH estimates the loadings using a subset of the sample and test the association in the rest of the sample. Extensive simulation studies performed in 28 . The author of 28 showed that MANOVA via Wilk's lambda was substantially more powerful than PCQH 23 with K 5 phenotypes.

Combining Test Statistics from Univariate Analysis
An alternative way to analyze multivariate phenotypes is to perform univariate phenotypegenotype association test for each phenotype individually and then combine the test statistics from the univariate analysis. The advantage of such approach is the simplicity, that is, the methods to deal with univariate phenotypes are generally simpler than methods for MV phenotypes. It is especially useful for analyzing multivariate phenotype consisting of components of different types of distributions such as continuous, dichotomous, and survival. Regression methods for analyzing such multivariate phenotype are generally complicated and not trivial to implement for MV phenotype with dimension > 2, see for example, 29, 30 . In recent years, researchers have generated large amount of univariate GWAS results for a variety of complex traits. Methods that combine the univariate results of multiple traits to detect genetic markers associated with multiple phenotypes are appealing.

Methods for Heterogeneous Genetic Effects across Phenotypes
The limitation of O'Brien statistic is that it is not powerful for heterogeneous effects across multiple phenotypes, especially if some effects are of opposite directions. Another class of statistics that takes a quadratic form of the vector of the individual association statistic may overcome the limitation. For example, the following Wald chi-squared type test statistic was mentioned in Xu et al. 32 .

Journal of Probability and Statistics
The difference between 2.10 and 2.11 is that the vector e 1, 1, . . . , 1 T is replaced by the T in 2.11 . S follows a chi-squared distribution with degree of freedom equal to the number of the phenotypes K or rank of Σ if it is not full rank. Due to the "curse of dimensionality," power of 2.11 is diminishing with the increased number of phenotypes. Similar problem has been extensively studied and discussed in high-dimensional data analysis field and most recently in the analyses of multiple rare variants. Borrowing ideas from these fields, we propose the following test statistic that may be more powerful than 2.10 and 2.11 with heterogeneous effects.
The difference between 2.12 and 2.11 is that there is no variance-covariance matrix in 2.12 . This statistic was first proposed by Pan 35 to analyze multiple rare or common variants against a single phenotype, where the t i is the beta coefficient for the ith genetic variant. Different from Pan 35 , here t i is the association statistic for the ith phenotype with a single marker. Based on the groundwork of Zhang 36 , Pan 35 pointed out that the distribution of 2.11 is a mixture of single degree-of-freedom chi-squared variates, where c i s are the eigen values of Σ, that is, the variance-covariate matrix of t i . The distribution of 2.12 can be well approximated by aχ 2 The P value is calculated as p χ 2 d > S Sq − b /a . The degree of freedom of the S sq may be less than K with highly correlated phenotypes. In addition, 2.12 does not have the problem of instability observed for 2.10 and 2.11 when some of the components are highly correlated in one of our applications, a correlation ∼0.7 has resulted in inflated results for 2.10 and 2.11 . We have developed an R package CUMP combining univariate results of multivariate phenotypes that have implemented all the aforementioned combining statistics approaches. The software can be downloaded at http://people.bu.edu/qyang/ , and a short report of this software is submitted 37 .

Identifying Pleiotropy
All the aforementioned methods can be used to detect association that is potentially due to pleiotropy. But they do not answer the question if the detected association is truly pleiotropy, that is, the marker locus affects all components of the MV phenotype directly. The detect association can affect some of the phenotypes and/or mediate through these phenotypes to affect the other phenotypes. Vansteelandt et al. 38 illustrated potential confounding mechanism between the genotype of a genetic marker and a phenotype using a causal diagram Figure 1 : the association between the genotype, denoted as G, and the response phenotype Y can occur through the paths connecting the two variables along all unbroken sequences of edges regardless of the direction of the arrows, given that there are no colliders i.e., variables in which two arrows converge, e.g., variables K and L in Figure 1 in the sequence 39 .
The genotype G may be associated with Y due to 1 direct causal effect, that is, G → Y ; 2 through intermediate phenotype or risk factors, that is, The authors showed that two commonly used approaches to detangle the complex relationship between phenotypes, genotype, and traditional risk factors are flawed. The first commonly used approach derives the residuals of Y regressing on K, say Y Y −βK, and then the association between G and Y is tested. The disadvantage of this approach is that not only the direct causal effect of K on Y is removed but also any indirect effect of K on Y through G e.g., K ← G → Y and Y ← U ← L ← G → K and other factors e.g., K ← L → U → Y . Therefore β may be biased in the presence of confounding factors which leads to biased test of G with Y . The second commonly used approach tests the direct effect of G on Y in a regression model including K and L as covariates. Adjustment of K removes the relationship between G and Y through G → K → Y ; however, because K is a collider Figure 1 , the adjustment of K induces a spurious association 39, 40 along the path G → K → L ← U → Y . Additionally, adjusting for L induces spurious associated through the path G → L ← U → Y .
To overcome the limitation of the two commonly adapted approaches, Vansteelandt et al. 38 proposed a least squared regression model to estimate the direct effect size of K on Y . This regression model includes the suspected intermediate phenotype, the score of the genetic marker genotype, X G , and other common risk factors between the two phenotypes as regressors: The estimated effect size of the phenotype represents the direct effect of the K on Y , that is, not confounded by the effect of X mediated through any of the covariates. Then, a new phenotype is created as the residual of the response subtract the effect of K only Then, whether the G only exerts its effect on Y through K can be tested using any standard association test statistic between the residual and the X. A negative result indicates that G only exerts its effect on Y through K while a positive result indicates that the G has a direct effect on Y and/or a spurious effect through other confounders. Extensions of the method to dichotomous and time-to-event outcomes have been proposed 41, 42 .

Discussion
In this paper, we reviewed methods available for joint analyzing correlated phenotypes in genetic association studies. Some of these methods are designed to detect potential association with multiple phenotypes pleiotropy , while the others are designed to test whether the detected association with the MV phenotype is truly pleiotropy or the genetic marker exerts its effects on some phenotypes through affecting the others. For methods designed to detect association, each method has its own pros and cons. Random effects model requires knowledge of residual correlation, and misspecifying the correlation may incur inflation or power loss. Generalized estimating equations are robust to misspecification of residual correlation, but it is inflated for low-frequency variants and less powerful than random effects model in our experience. Variable reduction approaches are appealing because correlated outcomes are reduced to a single or fewer number of uncorrelated outcomes. However, in the presence of missing data in the outcomes, individuals with missing data do not contribute to the analysis, which may result in power loss. The approaches combining univariate association results are more flexible than the other methods especially when MV phenotypes consist of a mixture of continuous, discrete, and/or time-to-events data. Regression approaches have been developed to deal with such phenotypes. But they are generally complicated and few available software implements these methods. Since univariate association results are used, individuals with incomplete observations still contribute to the analysis of available phenotypes. Simulations on all continuous phenotypes indicated that the power of O'Brien's method, one of the approaches combining univariate association results is similar to regression and variable reduction methods when the effects size are similar across multiple phenotypes 34 .
All the approaches introduced here for population based approaches assume unrelated individuals. When there are related individuals in the data, not accounting for family structure can result in inflation or power loss. Extension of introduced methods to account for family data are possible. For example, one may add a random effect in mixed effects model to account for family structure. For approaches combining univariate association results, a model that account for family structure need be used in the univariate analyses.
In terms of computational cost, mixed effects models may be most time consuming since maximization of likelihood is required.
Finally, it has been shown that traditional causal inference is useful in distinguishing true pleiotropy from other mechanisms that also result in genetic association with multiple phenotypes. A related causal inference in recent genetic literature is Mendelian randomization test 43-45 . This approach can be used to infer whether an intermediate phenotype has a causal effect on an outcome phenotype, using genetic marker s in association with the intermediate phenotype. Unlike a phenotype that is subject to the influence of uncontrolled environmental factors and/or reverse causation of another phenotype, genotype s of genetic marker s is are free of influence of environmental factors and reverse causation. For this approach, marker genotype s is are used as an instrument variable. This test requires that there is no pleiotropy effect of the genetic marker on outcome phenotype. Association of the genotype and outcome phenotype indicates that the intermediate phenotype may causally affect the outcome phenotype.