#### Abstract

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.

#### 1. 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; 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.

#### 2. 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 time-to-event components, and methods for MV phenotypes consisting of a mixture of different types of components.

##### 2.1. 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.

###### 2.1.1. 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 denote the *k*th continuous component of the -dimensional phenotype of the *j*th individual. Let be the genotype of a genetic marker of the *j*th individual, and a score of the genotype. The linear mixed effects model takes the following form:
where is the intercept or other genetic or environmental fixed effects; is the fixed effect size of on the *k*th phenotype; are the random effects correlated within *j*th person; is the random errors iid. ~. Between any two individuals, are independent. Within a person, are correlated. The null hypothesis that the genetic marker is not associated with any phenotype component corresponds to . The estimation of variance parameters and fixed effect parameters can be obtained using restricted maximum likelihood method (REML) [8, 9].

When is categorical, it can be modeled with generalized mixed effects model (GLMM) as follows: where is a link function and 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 . 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 , where 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 , that is, , where a single degree-of-freedom (df) test 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 are very different, especially have different signs and cancel each other out.

###### 2.1.2. 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 for the *k*th phenotype of the *j*th individual. Let be the genotype of a genetic marker of the *j*th individual, and a score of the genotype as follows:
where are subject specific random effects following , and is a *-*dimensional correlation matrix. This is the Gaussian frailty model. There is another class of frailty models where 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 within *j*th person. The correlation between the phenotypes within a person is attributable to the random effects 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].

###### 2.1.3. 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 , and 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.

##### 2.2. 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, are coefficients that maximize the variance of ; principal component of heritability (PCH) with coefficients maximizing the total heritability of [18] 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 [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.

###### 2.2.1. 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 and the genetic marker. Specifically, in the estimation sample, the total phenotype variance is partitioned into QTL variance and residual variance as follows: where is the total phenotype variance-covariance matrix, the QTL variance matrix, and the residual variance matrix. Let , then the variance of explained by the genetic marker is A that maximizes can be obtained by solving the following generalized eigen system [18]: can be approximated by , where , is the sample standard deviation of the score of genotype across all individuals, is estimated using the least squared estimator of , and .

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 because these data will not contribute directly to the family-based association tests (FBAT). Then perform FBAT of on . However, FBAT has low power in the absence of population stratification [26] compared to population based approaches. Klei et alβs. [23] is a population-based association approach where they randomly split the sample into two subsets: one used to estimate , the other used to test the association of with via a linear regression model: . This ensures valid value in the association test.

###### 2.2.2. Canonical Correlation Analysis

Canonical correlation analysis seeks coefficients so that the squared correlation between in (2.5), and the score of genetic marker, , is maximized. Here is called estimated canonical correlation. To obtain , the covariance matrix of and is partitioned as follows: where is the matrix of the variance-covariance matrix of , and its transpose are and matrix of the covariance matrix between and , is the variance of , a scalar. All these submatrices can be estimated using the respective sample co-variance matrix. The canonical correlation, is solved as the squared root of the largest eigenvalue of , and the corresponding eigenvector contains the coefficients for constructing . multivariate analysis of variance (MANOVA) tests correspond to evaluating canonical correlation. Table 1 details the relationship between and commonly reported test statistics in MANOVA of a multivariate phenotype on [27].

These tests are implemented in SAS PROC GLM and function *summary.manova()*. As part of the PLINK package specifically developed for genetic analysis, Ferreira et al. [24] implemented the Wilks lambda, and its value is obtained from -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 phenotypes.

##### 2.3. Combining Test Statistics from Univariate Analysis

An alternative way to analyze multivariate phenotypes is to perform univariate phenotype-genotype 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.

###### 2.3.1. Methods for Homogeneous Genetic Effects across Phenotypes

Assume that is a vector of test statistics obtained from association analyses of each individual component phenotype against the genetic marker. Assume that follows a multivariate normal distribution with mean and a nonsingular covariance matrix . For example, ββcan be the coefficients from least squared regression model for individual components or the *-*test statistics from the regression models. The null hypothesis of no association to any phenotypes is H_{0}: . OβBrien [31β33] suggested the following linear combination of , with weight of length :
when (2.10) is the most powerful test among a class of tests statistics that are linear combinations of . Under the null hypothesis, follows the normal distribution with mean 0 and variance . To estimate with GWAS results, Yang et al. [34] suggested using the sample covariance matrix of the statistics on a large number of SNPs genomewide with little or no linkage disequilibrium among them (say HapMap ).

The power of OβBrienβs method depends on the assumption . When the means are very different or with opposite signs, OβBrienβs method may not be efficient. Yang et al. proposed a sample splitting approach that replaces the uniform weight by weights ββestimated using a portion of the sample and only used the remaining sample to estimate ββin (2.10), that is, . βTo overcome the variability introduced by a random sample splitting, Yang et al. also evaluated a cross-validation approach that averages the test statistics of 10 random splitting samples. The results showed that when are of different magnitude or in opposite directions, OβBrienβs method is less powerful than Yang et al., which indicates room for improvements for OβBrienβs method. However, the sample splitting and cross-validation methods are less powerful than OβBrienβs method with homogeneous effect sizes.

###### 2.3.2. 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].
The difference between (2.10) and (2.11) is that the vector is replaced by the in (2.11). follows a chi-squared distribution with degree of freedom equal to the number of the phenotypes 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 is the beta coefficient for the *i*th genetic variant. Different from Pan [35], here is the association statistic for the *i*th 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 s are the eigen values of Ξ£, that is, the variance-covariate matrix of . The distribution of (2.12) can be well approximated by with
The value is calculated as . The degree of freedom of the may be less than 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].

#### 3. 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 , and the response phenotype 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 and in Figure 1) in the sequence [39].

The genotype may be associated with *Y *due to direct causal effect, that is, ; through intermediate phenotype or risk factors, that is, ; because of confounding factors, that is, or .

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 , say , and then the association between and is tested. The disadvantage of this approach is that not only the direct causal effect of on is removed but also any indirect effect of on through (e.g., and ) and other factors (e.g., ). Therefore may be biased in the presence of confounding factors which leads to biased test of with .

The second commonly used approach tests the direct effect of on in a regression model including and as covariates. Adjustment of *K* removes the relationship between and through ; however, because is a collider (Figure 1), the adjustment of induces a spurious association [39, 40] along the path . Additionally, adjusting for induces spurious associated through the path .

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 on . This regression model includes the suspected intermediate phenotype, the score of the genetic marker genotype, , and other common risk factors between the two phenotypes as regressors: The estimated effect size of the phenotype represents the direct effect of the on , that is, not confounded by the effect of mediated through any of the covariates. Then, a new phenotype is created as the residual of the response subtract the effect of only . Then, whether the only exerts its effect on through can be tested using any standard association test statistic between the residual and the . A negative result indicates that only exerts its effect on through while a positive result indicates that the has a direct effect on and/or a spurious effect through other confounders. Extensions of the method to dichotomous and time-to-event outcomes have been proposed [41, 42].

#### 4. 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.

#### 5. URLs to Software Mentioned in This Paper

βSAS: http://www.sas.com/,βR: http://www.r-project.org/,βCUMP: http://cran.r-project.org/web/packages/CUMP/index.html,βcoxme: http://cran.r-project.org/web/packages/coxme/index.html,βgee: http://cran.r-project.org/web/packages/gee/index.html,βsurvival: http://cran.r-project.org/web/packages/survival/index.html,βlme4: http://cran.r-project.org/web/packages/lme4/index.html,βPLINK: http://pngu.mgh.harvard.edu/~purcell/plink/.

#### Acknowledgments

Q. Yangβs work is supported by the National Heart, Lung, and Blood Instituteβs Framingham Heart Study (Contract no. N01-HC-25195) and Grant no. R01HL093328 and R01HL093029. Y. Wangβs work is supported by NIH Grants nos. R03AG031113-01A2 and 1R01NS073671-01 1.