Joint association analysis of bivariate quantitative and qualitative traits

Univariate genome-wide association analysis of quantitative and qualitative traits has been investigated extensively in the literature. In the presence of correlated phenotypes, it is more intuitive to analyze all phenotypes simultaneously. We describe an efficient likelihood-based approach for the joint association analysis of quantitative and qualitative traits in unrelated individuals. We assume a probit model for the qualitative trait, under which an unobserved latent variable and a prespecified threshold determine the value of the qualitative trait. To jointly model the quantitative and qualitative traits, we assume that the quantitative trait and the latent variable follow a bivariate normal distribution. The latent variable is allowed to be correlated with the quantitative phenotype. Simultaneous modeling of the quantitative and qualitative traits allows us to make more precise inference on the pleiotropic genetic effects. We derive likelihood ratio tests for the testing of genetic effects. An application to the Genetic Analysis Workshop 17 data is provided. The new method yields reasonable power and meaningful results for the joint association analysis of the quantitative trait Q1 and the qualitative trait disease status at SNPs with not too small MAF.


Background
Statistical methods for the univariate association analysis of quantitative and qualitative traits have been well developed in the literature. Complex human diseases are often characterized by multiple traits. These traits tend to be correlated with each other because of common environmental and genetic factors. In the genetic analysis of complex diseases, it is natural to account for the correlations among multiple traits and to model them simultaneously. Joint genetic linkage analysis of multiple correlated phenotypes has been studied by Jiang and Zeng [1], Mangin et al. [2], Amos and Laing [3], Almasy et al. [4], Blangero et al. [5], Wijsman and Amos [6], and Williams et al. [7,8], among others. Joint linkage analysis of multiple correlated traits can potentially improve the power to detect linkage signals at genes that jointly influence a complex disease.
Recently, Liu et al. [9] developed an extended generalized estimating equation method for the bivariate association analysis of continuous and binary traits. Their simulation results demonstrated that, compared with univariate analysis, bivariate analysis may substantially improve power while having comparable type I error rates under certain situations.
In this paper we extend the joint linkage analysis of multivariate qualitative and quantitative traits described by Williams et al. [7,8] to association analysis. Specifically, we assume that a latent variable determines the qualitative trait and that the latent variable and the quantitative trait follow a bivariate normal distribution. With such modeling, we develop likelihood-based inference procedures for testing pleiotropic genetic effects. As an illustration, we perform the joint association analysis of the quantitative trait Q1 and the qualitative trait disease status on chromosome 13 from the Genetic Analysis Workshop 17 (GAW17) data.

Methods
Suppose that the data contain n independent individuals. Let Y i1 be the qualitative trait and Y i2 be the quantitative trait for the ith individual. Let X i denote a vector of covariates, including the intercept and environmental variables, and let Z i denote a vector of genotype score(s) at the major single-nucleotide polymorphism (SNP) locus. We may also include gene by environment interaction terms. We assume that Y i1 is determined by a latent continuous variable Y i1 * such that: We consider the following model: where the b are regression coefficients for the environmental effects, the g are regression coefficients for the genetic effects, and the subscripts 1 and 2 denote the qualitative and quantitative traits, respectively, and the e e e i i i = ( , ) 1 2 are independent and identically distributed bivariate normal random variables with mean 0 and variance-covariance matrix Σ: where σ 1 2 is the variance of ε i1 , σ 2 2 is the variance of ε i2 , and ρ is the correlation between ε i1 and ε i2 .
To ensure the identifiability of the model, we fix s 1 = 1 and a = 0.
Given Y i2 , it can be shown that Y i1 * follows a normal distribution with mean and variance: Therefore the likelihood function for the unknown parameters θ ≡ (b 1 , g 1 , b 2 , g 2 , s 2 , r) takes the form: where: and Φ is the cumulative distribution function of the standard normal distribution.
To obtain the maximum-likelihood estimator (MLE) of θ, denotedq , we maximize the log-likelihood by using the quasi-Newton algorithm described by Press et al. [10]. We can perform various association analyses through the likelihood ratio test (LRT) statistic, given by: where  q is the MLE of θ under the null hypothesis. The LRT follows a chi-square distribution asymptotically with the degrees of freedom being the difference of the number of free parameters under the null hypothesis and the number under the alternative hypothesis.

Results
We first conducted a small simulation study to evaluate the type I error rates of the association test based on the univariate analysis at significance levels of 0.01 and 0.05. Figure 1 presents the results based on 100,000 replicates, where the x-axis denotes the correlation between the two Wald test statistics for the two traits. Figure 1 implies that failing to account for the correlation between traits may lead to inflated type I error rate.
We next performed a joint association analysis of the quantitative trait Q1 and the qualitative trait disease status from the GAW17 data. Before conducting the analyses, we had knowledge of the answers. Sex, Age, and Smoking were included as covariates in the model. We evaluated the additive effect of each of the 11 true SNPs on chromosome 13 on each individual trait as well as the pleiotropic effect on both traits. The genotype score was coded as the number of minor alleles at each SNP locus. For comparison, we also performed a univariate association analysis and used the sum of the LRT statistics for testing genetic effect on an individual trait as the test statistic for testing the pleiotropic effect. Table 1 presents the powers of the tests at the 11 true SNPs on chromosome 13 from the joint association analysis and the powers from the univariate analysis based on 200 replicates. Table 2 presents the type I error rates at 6 randomly chosen SNPs on chromosome 13, excluding the 11 true SNPs. None of these 6 SNPs was in linkage disequilibrium with the 11 true SNPs on chromosome 13. Type I error rates and powers of the association tests appear to depend very much on the minor allele frequency (MAF). Type I error rates of both tests are lower than the nominal significance level when the MAF is small. Essentially no tests have power to detect the SNP effect when the MAF is low and the minor allele is detected in only a few of the 697 individuals. On the other hand, both methods had good power to detect the pleiotropic   effects at true SNPs with MAF > 0.0043. The estimates of the additive effects and the correlations for the 11 true SNPs are shown in Table 3. As expected, a strong correlation exists between quantitative trait Q1 and disease status with an estimated correlation of 0.68.

Discussion and conclusions
In this paper, we described a likelihood-based approach for the joint association analysis of quantitative and qualitative traits in unrelated individuals. As expected from the true data generation model, our association analysis of the GAW17 data reveals that there is strong correlation between the quantitative trait Q1 and the qualitative trait disease status. In this case, the joint association test of pleiotropic effects is more appropriate than the test based on univariate analysis, which ignores the trait correlation. Two factors might have contributed to the small difference between the joint analysis and the univariate analysis using the GAW17 data set. First, both methods had low power with rare variants, for example, when MAF < 0.0014. On the other hand, when the MAF was moderate and the effect size was relatively large, both methods yielded power close to or equal to 100%; therefore we are unlikely to see the difference between the two methods under such situations. In fact, it may not be fair to compare the powers between the two methods because the univariate analysis tends to inflate type I error rates in the presence of correlation between the qualitative trait and the quantitative trait.
In its current form, the single-SNP joint analysis is more ideal for common variants and less ideal for rare variants. To improve the power of the association tests with rare variants, we could apply the approaches of Morris and Zeggini [11] based on accumulations of rare variants within the same functional unit. The idea is to treat the proportion of rare variants at which an individual carries a minor allele as a covariate in the joint analysis model and to test for the covariate effects on the qualitative and quantitative traits jointly.