Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Genet Epidemiol. Author manuscript; available in PMC 2010 April 1.
Published in final edited form as:
Genet Epidemiol. 2009 April; 33(3): 217–227.
doi:  10.1002/gepi.20372
PMCID: PMC2745071

Bivariate Association Analyses for the Mixture of Continuous and Binary Traits with the Use of Extended Generalized Estimating Equations


Genome-wide association (GWA) study is becoming a powerful tool in deciphering genetic basis of complex human diseases/traits. Currently, the univariate analysis is the most commonly used method to identify genes associated with a certain disease/phenotype under study. A major limitation with the univariate analysis is that it may not make use of the information of multiple correlated phenotypes, which are usually measured and collected in practical studies. The multivariate analysis has proven to be a powerful approach in linkage studies of complex diseases/traits, but it has received little attention in GWA. In this study, we aim to develop a bivariate analytical method for GWAS, which can be used for a complex situation that a continuous trait and a binary trait measured are under study. Based on the modified extended generalized estimating equation (EGEE) method we proposed herein, we assessed the performance of our bivariate analyses through extensive simulations as well as real data analyses. In the study, to develop an EGEE approach for bivariate genetic analyses, we combined two different generalized linear models corresponding to phenotypic variables using a Seemingly Unrelated Regression (SUR) model. The simulation results demonstrated that our EGEE-based bivariate analytical method outperforms univariate analyses in increasing statistical power under a variety of simulation scenarios. Notably, EGEE-based bivariate analyses have consistent advantages over univariate analyses whether or not there exits a phenotypic correlation between the two traits. Our study has practical importance, as one can always use multivariate analyses as a screening tool when multiple phenotypes are available, without extra costs of statistical power and false positive rate. Analyses on empirical GWA data further affirm the advantages of our bivariate analytical method.


Genome-wide association (GWA) studies offer a powerful tool for identifying genes conferring modest disease risks [Risch and Merikangas 1996]. Currently, statistical methods developed for GWA data analyses were largely for single phenotypes. However, in practice, investigators often collect a set of correlated phenotypes that may share common environmental and/or genetic factors. The multivariate analysis is a natural way to take all phenotypes into consideration simultaneously. Multivariate analyses have been widely adopted in linkage studies [Allison, et al. 1998; Huang and Jiang 2003; Jiang and Zeng 1995; Lange and Whittaker 2001; Liu, et al. 2007a; Williams, et al. 1999], and there is a general consensus that multivariate analyses outperform univariate analyses in increasing statistical power and precision of parameter estimation. However, in GWA, multivariate analyses have not received sufficient attention, although there are some sporadic reports on methodology development of multi-trait analyses [Lange, et al. 2003; Lange, et al. 2002; Verzilli, et al. 2005]. In particular, as to the common situation when both binary (e.g., disease status) and continuous data (e.g., quantitative measures that are risk factors for the disease under study) are involved, little work has been done to evaluate the performance of the multivariate analytical method.

For a mixture of correlated phenotypes involving non-normal traits (such as binary or ordinal data), it is not straightforward to specify a full probability model, making it difficult to directly implement likelihood-based approaches [Lange and Whittaker 2001; Lange, et al. 2002]. To surmount this problem, several methods were developed [Catalano 1997; Gueorguieva and Agresti 2001; Gueorguieva and Sanacora 2006; Regan and Catalano 2000] based on the threshold model assumption, which transform a discrete phenotype to a “latent” normal variable. However, the potential concerns with these methods are the arbitrary and untestable distributional assumptions on the “latent” variable, and the difficulty with computational features of model fitting [Prentice and Zhao 1991].

Since the seminal work of Liang and Zeger [1986], a suite of methods [Hall 2001; Hall and Severini 1998; Liang, et al. 1992; Prentice and Zhao 1991; Zhao and Prentice 1990] based on the Generalized Estimating Equation (GEE) have been developed to address the situation of non-normally distributed variables in multivariate data analyses. GEE can be regarded as a multivariate version of generalized linear model (GLM), a method with reasonable statistical efficiency to analyze correlated data with non-normal distributions. A vast of literature on GEE theory and application has arisen over the past two decades [Zeger and Liang 1992; Ziegler, et al. 1998; Zorn 2001]. The most commonly used GEE-based approaches, referred to as GEE1 and GEE2, were introduced by Liang and Zeger [1986] and Zhao and Prentice [1990], respectively.

GEE1 and GEE2 have their respective limitations. In GEE1, regression and association parameters are presumed to be orthogonal which may not necessarily true, and the association parameter is treated as a nuisance [Liang and Zeger 1986; Liang, et al. 1992]. Because the relationship between the first- and second-order moment parameters is ignored, GEE1 has less efficient estimators than GEE2 if the working covariance model is correctly specified [Hall and Severini 1998; Liang, et al. 1992]. In contrast to GEE1, GEE2 removes the assumption of “orthogonality” initially imposed on GEE1. Thus, GEE2 can perform estimation for regression and association parameters simultaneously under a unified framework, leading to a significant increase in efficiency of parameter estimation. However, a limitation of GEE2 is that it may lose the consistency of regression parameter estimators if the covariance structure is misspecified [Hall 2001; Hall and Severini 1998; Liang, et al. 1992; Zhao and Prentice 1990].

To improve the performance of GEE1 and GEE2, an alternative GEE-based approach, named Extended Generalized Estimating Equation (EGEE), was proposed by Hall and Severini [1998]. EGEE was developed using the ideas of extended quasi-likelihood [McCullagh and Nelder 1989; Nelder and Pregibon 1987]. It retains many of the advantages of the original GEE1 and GEE2 methods while avoiding their potential limitations: (1) EGEE can offer consistent estimations of both regression and association parameters even under an incorrect covariance model; (2) EGEE is not required to assume third- and fourth-order moment “working” models which are necessary in GEE2; and (3) the estimation efficiency of EGEE is comparable to GEE2 and often outperforms GEE1 [Hall 2001; Hall and Severini 1998].

In the present study, using EGEE, being an improved GEE-based approach, we aim to develop bivariate gene/phenotype association analyses for population-based data when both continuous and binary traits are involved in the same studies. Since the original EGEE [Hall 2001; Hall and Severini 1998] cannot deal with mixed continuous and binary outcomes in a straightforward fashion, we modified the EGEE method for such situations. Our strategy is motivated by the usage of Seemingly Unrelated Regression (SUR) [Zellner 1962]. SUR procedure is an extension of ordinary least squares analyses for dealing with the system of linear equations with correlated error terms, which can be treated as a special case of generalized least squares analyses. Based on the idea of SUR models, two GLMs with different link functions can be incorporated into a unified equation system. Consequently, we developed a “specific” form of EGEE, which is well suited for the mixture of two phenotypes with different distributions.

We applied our modified EGEE method to assess the performance of bivariate analyses through extensive simulation studies. The results demonstrated that, compared with univariate analyses, bivariate analyses can substantially improve power while having comparable false-positive rates under almost all the scenarios simulated. The advantages of our EGEE-based bivariate analyses method was further validated in empirical data analyses using our genome-wide association scan data for osteoporosis and obesity. Given the emerging widespread GWA, our proposed method may have practical significance in the genetic epidemiology field in general.


We begin with the situation in which N unrelated random subjects are recruited in a population-based genetic association study, each having observations on two different types of phenotypes: one is in normal distribution (named T1) and the other is a binary trait (named T2). A 2×1 random observation vector for individual i (i = 1, …, N) is denoted as yi = (yi1, yi2)’ with the mean vector μi = E(yi) = (μi1, μi2)’ and the 2×2 (co)variance matrix var(yi) = Vi. Let xi be a L×1 vector of explanatory variables including genetic makers as well as other fixed effects (e.g., age, sex) corresponding to individual i. We assume xi affects both phenotypes through corresponding L-dimensional regression parameter vectors β1 and β2. Similarly, we set a 2L×1 vector β = (β1’, β2’)’. Here, the subscripts 1 and 2 relate the corresponding variables to traits T1 and T2 respectively. Within the GLM framework [McCullagh and Nelder 1989], using an identity link (for T1) and a logit link (for T2), we model the relationship between the explanatory variables and the marginal means of the two phenotypes as follows:



Here, for ease of incorporating the two different types of data into an overall EGEE analysis, we employ SUR model [Zellner 1962] to combine Equations [1] and [2] into a unified framework:


where Xi=(xi00xi), and g(μi) is a compound function vector in the form of (g1(μi1), g2(μi2)) ’.

Following the GLM theory, the respective variance functions of two distributions, i.e., the normal and the binary, are given as:


where μi1and μi2 are variables of the two variance functions corresponding to the two phenotypes. Note that we have ν1(•) [equivalent] 1 here because of a normal distribution chosen by the trait T1.

In the framework of GLM, probability distributions of the random variables are usually parameterized in terms of not only the means (e.g., μi1and μi2 herein) but also dispersion parameters. Accordingly, we denote the dispersion parameters for different distributions of T1 and T2 as ψ1 and ψ2, respectively. Since in the binary outcome context there exists no overdispersion [Hall 2001], we will always have ψ2 [equivalent] 1.

For specifying the marginal (co)variance model Vi, we further define the “working” correlation matrix R(ξ)=(1ξξ1), where ζ is the unknown association parameter falling within (-1.0, 1.0). It is noted that the regressors involved in the SUR model are allowed to have possibly different effects on either of the two regressands considered, with residuals from each regression assumed to be correlated [Verzilli, et al. 2005].

Given the variance functions, dispersion parameters and working correlation matrix defined above, the (co)variance matrix Vi between the two phenotypes for individual i has the form similar to those given by Rochon [1996] and Lange et al. [2002].


On the basis of Equation [4] and the equality ψ2 [equivalent] 1, Equation [5] can be fully described as:


From Equation [6], it is obvious that Vi is determined by β, ζ and ψ1. Since ψ1 should be restricted as positive, we make a “squared” transform, ψ1 = [var phi]2, such that the estimate of ψ1 can fall into its space. Accordingly, Equation [6] is rewritten as:


For ease of presentation, we combine the second-order moment parameters in a 2×1 vector fashion, α~=(ξ,φ).

To examine whether the marker locus is associated with either of the two phenotypes within the EGEE framework, we take two steps (named in turn Estimation Step and Testing Step) on the basis of the SUR regression model and the marginal (co)variance model given in Equations [3] and [7] respectively.

Estimation Step

From the EGEE perspective, the unknown parameters, including the regression vector β and the association vector α~ aforementioned, can be estimated in a set of estimating equations. Here, we directly adopt the same form of the EGEEs as given in the original study [Hall 2001]:


where, the 2×2L matrix Di = [partial differential]μi / [partial differential]β’, and the 4×2 matrix Fi={vecVi1}α~. The operator vec(•) here is a function for generating a column vector via stacking the columns (or rows) of the argument matrix. The 4×1 vectors si and σi are defined for estimating the second-order moment parameters, where si = vec{(yiμi )(yiμi)’} and σi = E(si ) = vecVi.

Obviously, Equation [8] regarding the unknowns β and α~ is hard to resolve straightforwardly due to its nonlinear property. Alternatively, a Fisher scoring algorithm may offer a natural way for obtaining the estimates of β and α~ . Starting with initial values β(0) and α~(0), the updated parameter estimates β(l) and α~l in the lth iteration are given by


The above iteration is repeated until convergence. In our study, we set a sufficiently rigid convergence criterion as


Accordingly, the updated values of β and α~ in the latest iteration will be treated as the solution of the EGEEs given by Equation [8], which are denoted as β~ and α~^, respectively. Accordingly, a sandwich-type estimate of the asymptotic (co)variance matrix of β~ and α~^, var(β~,α~^), is given by



Testing Step

We employ a Wald chi-square statistic to test whether the marker considered has an effect on either the continuous phenotype or the binary trait. With the assumption that the distribution of genotypes at the SNP locus is in Hardy-Weinberg equilibrium, the EGEEs for individual i can include the allele-based index variable xi,m with the regression coefficients βm,1 and βm,2 on traits 1 and 2 respectively. Under the null hypothesis that the SNP has no effect on either of the two traits, an asymptotic chi-square statistic with 2 degree of freedom (df) can be expressed as


where β^m,1 and β^m,2 are the estimates of βm,1 and βm,2 respectively, which are obtained from the estimate vector β^ in the Estimation Step. Similarly, var(β^m,1,β^m,2), the (co)variance matrix between β^m,1 and β^m,2, is a sub-matrix of the matrix var(β~,α~^) given by Equation [10], where only those rows and columns regarding β^m,1 and β^m,2 are selected.

Given a significant threshold such as α = 0.05, if the calculated p-value for χSNPEGEE2 is lower than α, we can declare that the SNP locus affects either or both of the two phenotypes. It is worth noting that we mainly focus on identifying the potential SNP effects on the traits through constructing the statistic χSNPEGEE2 based on the estimates of regression parameters (i.e., β^m,1 and β^m,2). A similar statistic can be easily developed in the same way to test the significance of the association parameter ζ if one is substantially interested in further exploring the correlation between two traits.

For comparison with univariate analyses, we separately detect SNP-phenotype associations based on traditional linear regression analysis (for T1) and logistic regression analysis (for T2). The respective p-values are denoted as p1 and p2. Under this situation, we calculated the adjusted significance level (denoted as αc) for multiple testing by two different methods, Bonferroni correction and Meff method [Nyholt 2004]. The originally Meff method is proposed to adjust for multiple testing of SNPs in linkage disequilibrium (LD) based on the spectral decomposition (SpD) of matrices of pairwise LD between SNPs. Here we directly apply Meff method to correct for multiple testing of two correlated phenotypes. Accordingly, we declare a significant association if min(p1, p2) <αc.



To evaluate the performance of bivariate association analyses based on our developed EGEE method, we conduct extensive simulations for joint analyses of a range of mixtures of one continuous trait and one binary trait in different scenarios.

Under an additive model, we assume a biallelic locus under consideration has effects on both traits, i.e., T1 with normal distribution and T2 in a binary pattern, within a random population. Firstly, for ease of describing the genetic effect on the binary phenotype, we adopt the penetrance function fi as used in Liu et al [2007b], which is the penetrance of the genotype with i copies of causal allele (i = 0, 1 or 2). Similarly, we define the ratio of f1/f0 as heterozygote genotype relative risk (GRR). Then we generate the marker genotypes and the two types of phenotypic observations for each individual under a suite of parameter settings, including sample size N (150 vs. 300), frequency of the liability allele q (0.05 vs. 0.15), the residual error correlation r (-0.45, -0.15, 0, 0.15 vs. 0.45) between two traits, QTL heritability h12 (0, 0.01 vs. 0.03) for trait 1, and heterozygote GRR (1.0, 1.5 vs. 2.0) for trait 2. Under the above settings, a total of 180 simulation scenarios are generated by combining these settings at different levels. For a simple illustration, we set the residual error variance σe,11=1.0 for trait 1. As to trait 2, we assume f0, i.e., the population incidence among individuals carrying zero copy of liability allele, being 0.1, and the ratio of cases and controls of the sample is set to 1:1 for each simulation scenario. According to the above parameter settings, two related phenotypes, both affected by the same locus as well as influenced by the environmental factors, can be generated via the following steps:

  1. Two alleles are sampled and paired to form a random individual according to the liability allele frequency q. Following the theory of quantitative genetics [Falconer and Mackay 1996], under an additive model, the simulated genotype effect on T1 can be expressed as
    where j is the number of liability alleles assigned to the simulated subject.
  2. We assume the locus simulated above also underlies a latent normal variable z2 with mean a2, j and variance σe,22=1, and the residual error correlation between T1 and z2 is set to r.
    The random variable z2 can be truncated to generate the binary trait T2 using a threshold of t. Based on our earlier study [Deng, et al. 2000], a2,j and the threshold t can be determined via the aforementioned simulation parameters, including heterozygote GRR and the penetrance function f0. To be concrete, we provide a detailed description in the APPENDIX section on how to bridge the gap between normal variable z2 and binary data T2 through heterozygote GRR and f0 as well as the threshold t.
  3. For the random individual i, we draw the related variables T1 and z2 conditionally on its genotype Gj from the bivariate normal distribution (T1z2)N((a1,ja2,j),(1rr1)). Accordingly, the respective binary phenotype T2 is truncated to a case (z2t) or a control (z2 < t).
  4. Through steps (1) ~ (3), we can obtain a random individual with a specific marker genotype as well as both the normal and binary phenotypic observations with correlation coefficient r. To ensure the ratio of cases and controls involved in the binary data T2, we repeat these processes till 0.5N cases and 0.5N controls are formed.

Under each of the specific simulation scenarios, we generate 20,000 independent samples for obtaining the robust estimates of statistical power and false-positive rates by the EGEE-based bivariate association analyses as well as by traditional univariate association analyses. For both bivariate and univariate analyses, we define the identical null hypothesis H0 as: the SNP marker has no association with either T1 or T2, i.e., h12=0 and heterozygote GRR =1. The corresponding alternative hypothesis Ha is designed as: the SNP maker is associated with either or both of the two traits, i.e., h12>0 or/and heterozygote GRR >1. Accordingly, the power and the false-positive rate can be calculated as the proportions of significant SNP/phenotype associations reported in 20,000 independent replicates with the true Ha and H0, respectively.


According to the combinations of various parameter settings for simulations, we perform extensive simulations involving 160 scenarios under the hypothesis Ha as well as 20 scenarios under the hypothesis H0, defined above. To study the performance of detecting liability alleles based on the bivariate test, we conduct three association tests for each set of simulated phenotypes employing the EGEE-based bivariate analytical method (denoted as B_EGEE) and separate univariate analyses adjusted for multiple testing by the Bonferroni correction (denoted as U_b) and the Meff method [Nyholt 2004] (denoted as U_e) . For all tests, we set the experiment-wise significance threshold α = 0.05.

Tables Tables11 and and22 present the power estimates of gene-phenotype associations with sample sizes of 150 and 300, respectively, for both bivariate and univariate association analytical methods. Under each setting, we highlight the maximal power in bold to emphasize the best performance among the three tests. It is obvious that the largest increases in power are consistently observed with B_EGEE, with the exception of only 3 data sets. This shows that the EGEE-based bivariate association test is highly preferable to the individual univariate tests adjusted either by the Bonferroni correction or by the Meff method.

Table 1
Power estimates of association tests for the mixture of the continuous and binary traits with the use of EGEE-based bivariate analyses and univariate analyses under the significance level α = 0.05 and the sample size N = 150.
Table 2
Power estimates of association tests for the mixture of the continuous and binary traits with the use of EGEE-based bivariate analyses and univariate analyses under the significance level α = 0.05 and the sample size N = 300.

Comparing the power of B_EGEE across different levels of correlation coefficient r between the two traits, we can see that the performance of bivariate analyses varies along with r in the wide range of situations investigated. The effects of residual error correlation on the performance of B_EGEE are twofold. On one hand, for the situation in which only one phenotype, either T1 or T2, is associated with the locus, the estimated power of B_EGEE has a trend of increasing with a rising absolute value of r from 0 to 0.45, regardless of its direction. Note that we set r > 0 or r < 0 herein, corresponding to the situations where it has either the identical or opposite direction, respectively, as that of correlation induced by the simulated pleiotropic QTL. On the other hand, under the scenario in which both phenotypes are associated with the gene considered, the performance of B_GEE is obviously influenced by both degree and direction of correlation r. Specifically, when the residual error correlation is opposite in sign to the correlation induced by QTL effects, B_EGEE has an increasing gain in power as the degree of r increases; otherwise, the power of B_EGEE decreases with the increase of r. The above trend of the power of B_EGEE is repeated over various design combinations of liability allele frequency and effect of QTL as well as sample size.

In Tables Tables11 and and2,2, comparing the performance of univariate tests with the Meff method and the Bonferroni correction, we can clearly see that the former consistently outperforms the latter across all those simulation scenarios with r > 0. Furthermore, the higher degree r has, the more obvious the advantage of Meff. Since Meff may take non-independence between correlated phenotypes into consideration, it reasonably overcomes the conservativeness of the Bonferroni correction while having the ease of computation.

Results from Tables Tables11 and and22 further demonstrate that the highest power levels for both bivariate and univariate association analyses are achieved under the combinational design of the higher liability allele frequency at q = 0.15, a larger sample size at N = 300 and the highest QTL effects on both phenotypes. This trend logically complies with the routine power profiles in genetic association studies.

Table 3 offers the results of false-positive rate calculations for the three tests, i.e., B_EGEE, U_b and U_e, under various design combinations at the significance level of 0.05. Overall, all estimated levels of false positive rates for both bivariate and univariate analyses are close to the threshold value of 0.05. This aspect further suggests that the increased power in B_EGEE is due primarily to the fact that it considers all phenotype data within a unified test framework rather than focusing on information of each single phenotype separately.

Table 3
False-positive rate estimates of association tests for the mixture of the continuous and binary traits with the use of EGEE-based bivariate analyses and univariate analyses across various simulation scenarios under the significance level α = 0.05. ...

Table 3 also shows that B_EGEE has a more stringent false-positive rate control for lower levels of both allele frequency r and sample size N. This result demonstrates that the convergence aspect of B_EGEE estimates toward their asymptotic distribution is affected, to some extent, by the sample characteristic mentioned above. For the two univariate analyses, it is obvious that U_b is slightly more conservative than U_e (when r≠0) since the former ignores the correlation between tests of the two phenotypes.



We perform real data analyses using our available GWA datasets for osteoporosis and obesity. Osteoporosis is a major public health problem characterized by low bone mineral density (BMD) and increased risk of fracture. Obesity is a condition of extra body weight defined as body mass index (BMI, weight/height2) ≥30 kg/m2. The relationship between osteoporosis and obesity has been widely studied, with conflicting results [Guney, et al. 2003; Radak 2004; Zhao, et al. 2008; Zhao, et al. 2007]. We recently found that osteoporosis and obesity share common genetic factors that may have inverse effects on BMD and BMI. To date, we are not aware of any studies investigating the relationship between obesity and low-trauma osteoporotic fractures (OF), the endpoint of osteoporosis.

Our GWA data are obtained using the Affymetrix 500K arrays on 700 Chinese Han subjects (350 with osteoporotic hip fractures and 350 controls) (Y. Guo, et. al., unpublished data). Among these, 598 also have BMI data.

To further validate the proposed method by assessing its empirical performance, we perform bivariate GWA using B_EGEE on the 598 subjects for which information is available for both OF and BMI. Before the association test, a suite of quality control procedures are performed. First, for the initial full-set of 500,568 SNPs, missing genotypes at each locus are imputed using HelixTree 5.3.1 (Golden Helix, Bozeman, MT), based on Expectation-Maximization (EM) algorithm. Second, we discard SNPs with a call rate < 95% (10,467 SNPs), those deviating from Hardy-Weinberg equilibrium (HWE) (P<0.001, 21,274 SNPs), and those having a minor allele frequency (MAF) < 0.05 in the total sample (141,945 SNPs). Eventually, 326,882 SNPs are available for subsequent analyses. At each locus, the allele-based association test is performed, and age, sex and age2 are included as covariates to adjust their respective effects on raw observations in the models for the B_EGEE bivariate analyses as well as respective univariate analyses.


We test whether either of the two phenotypes is associated with SNPs across the whole genome. For ease of presentation, we only select the top 100 SNPs with the strongest association signals, i.e., the lowest p-values, from the E_GEE-based bivariate analyses. The profile of the p-values (in -log scale) is plotted in Fig. 1. For comparison, for these 100 SNPs, the respective raw p-values of univariate analyses and the p-values adjusted by Bonferroni correction for multiple traits are also presented in Figs. 1a and 1b respectively. We also scrutinize the top 100 most significant SNPs based on the adjusted p-values of two univariate analyses. We find that there is a high overlap between these two sets of top 100 SNPs (89 out of 100), suggesting that the comparisons between the two methods are reasonable in terms of the top 100 SNPs chosen by B_EGEE.

Fig. 1
Profiles of the smallest p-values of top 100 SNPs in the EGEE-based bivariate genome-wide association study for the continuous trait BMI and the binary trait OF

Fig. 1a depicts the p-value profiles of three different analyses, i.e., two separate traditional univariate analyses for the continuous trait BMI and the binary trait OF and EGEE-based bivariate analyses for the mixture of the both traits. It can be seen that the majority of the SNPs (66 out of 100) show the strongest association signals in B_EGEE analyses even compared with the raw p-values in univariate analyses. For the remaining SNPs (34 out of 100), the p-values of B_EGEE fall between the two raw p-values of the univariate analyses, but are closer to the smaller and more significant ones. The trend from Fig. 1a suggests that bivariate association analyses have practical significance. Specifically, when two or more traits are involved in GWA, we may experience the situations shown in Fig. 1a, i.e., some causal variants with modest effects may show moderate signals (e.g., p = 1E-3) in univariate analyses while showing much stronger signals (such as p = 1E-7) in bivariate analyses. If only univariate associations are performed, these signals will probably not reach genome-wide significance. However, if a bivariate GWA screen is performed initially, some of these loci may achieve a genome-wide significance level. In the second stage, using separate univariate analyses to further confirm the association identified in bivariate GWA, we only need to adjust raw p-values of each univariate analysis for the number of phenotypes and specific handful markers screened out in the initial bivariate analyses. Therefore, those SNPs with moderate signals in univariate analyses may remain significant as the issue of multiple testing is minimized.

Fig. 1b further provides comparisons between the p-values of E_GEE statistic in bivariate analyses and the corresponding adjusted p-values for multiple testing of two phenotypes in univariate analyses. It shows that, with the exception of three loci, the p-value of B_EGEE statistic at each SNP locus is apparently lower than the corresponding adjusted p-value of univariate analyses. This further supports the conclusion that the performance of B_EGEE proposed herein greatly outperforms univariate analyses in studies in which multiple phenotypes are involved and multiple testing is adjusted for the number of phenotypes involved.

In addition, we selected 10 most significant SNPs based on p-values of B_EGEE statistic for our bivariate GWA data in Table 4. These identified SNPs could be treated as the potential genetic factors which are associated with the two traits investigated, i.e., BMI and OF. Due to the relatively small sample size of our GWA data, all of the association signals in our analysis did not reach the genome-wide significance threshold, such as 4.2×10-7 proposed by Lencz, et al. [2007]. However, our preliminary findings herein may offer evidence for further follow-up genetic studies on osteoporosis and obesity.

Table 4
The association signals for the 10 most significant SNPs in bivariate analyses on the empirical GWA data using B_EGEE method


In this study, we have developed a bivariate association analytical method based on EGEE for GWA studies. To the best of our knowledge, this is the first study to develop specific methods designed for investigating the performance of bivariate genetic association analyses for population-based GWA data. The extensive simulation experiments and GWA real data analyses consistently demonstrated the advantages of our method over univariate analyses.

Here we consider a complex but common situation where both continuous and discrete phenotypic observations are available, for which, however, no existing analytical methods or public software are readily available. In this study, using a SUR model, we incorporated two GLMs, corresponding to different phenotypic distributions, into the same set of estimating equations; this was followed by parameter estimation and statistical testing based on the EGEE method. As pointed out by Hall and colleagues [Hall 2001; Hall and Severini 1998], EGEE generally performs better in estimation of both regression and association parameters compared to either GEE1 or GEE2 approaches. However, the original EGEE was developed for analyzing multiple related phenotypes with the same type of distributions that have the same dispersion parameter, such as longitudinal data obtained in clinical trials. The novelty of our modified EGEE is that it relaxes the constraint of the original EGEE and can handle response variables having different distributions or dispersion parameters. Hence, our modified EGEE method is well suited for bivariate association analyses on the mixture of a normal distributed continuous trait and a binary trait.

The advantages of multivariate analyses over univariate analyses may lie in two general aspects: (1) Multivariate analyses can incorporate information of correlation between traits into a single test statistic, which renders a gain in power. (2) Multivariate analyses can avoid or alleviate the issue of multiple testing. This becomes particularly important in circumstances in which the traits under investigation can be divided into several subgroups [Lange, et al. 2003]. Another similar situation is that statistical significance (especially when modest) of a potential causal locus may diminish if adjusted for the number of loci in a univariate GWA. However, it is likely that this locus will remain significant if genome-wide multivariate association analysis is performed and the genome-wide significance can be achieved, as described in the section of ANALYSES OF EMPIRICAL GWA DATA.

Two important aspects of this study worth further emphasis. First, in the circumstance in which both phenotypes are affected by the locus, the power of bivariate analyses is affected by both degree and direction of the residual correlation r; a larger power increase is always favored by the residual correlation r with opposite direction to the QTL-induced correlation. This is in agreement with studies of Jiang and Zeng [1995] and Allison et al [1998] for bivariate linkage studies. Further theoretical interpretations on our finding herein can trace to a recently report [Evans 2002], which pinpointed why the sign of corrections between variables affect the power of bivariate analyses. The rationale lies in that the opposite sign between QTL-induced and residual correlations can result in a corresponding positive term involved in the expression of non-centrality parameter (NCP), increasing the magnitude of NCP and the power; otherwise a negative term of NCP arises, decreasing the NCP and the power. As a result, the opposite direction of the correlations always contributes to the greater increase in power in contrast to the correlations with the same direction. Although the above interpretations are originally for bivariate linkage analyses, it should work in the similar way in association studies from statistical perspective. Second, even in the circumstance in which two phenotypes do not share genetic determinations at the locus investigated or residual effects, bivariate analyses still have obvious advantages over univariate analyses across various simulated scenarios due to alleviating the issue of multiple testing. This has also been observed in multivariate family-based association studies [Lange, et al. 2003]. These two features make the B_EGEE procedure an ideal screening tool for detecting gene-phenotype association no matter whether the phenotypic correlation between traits exists or not.

In B_EGEE method, we employed a Fisher scoring algorithm to estimate the unknown parameters due to nonlinear property of the equation system. The program implementing B_EGEE method is written in FORTRAN language, which is available on our website A potential concern is that our method may entail prohibitorily high computational demand for conducting the iterative calculation. To illustrate that the proposed method is computationally practical for GWA studies, we assessed the CPU time required by the program in simulation and empirical GWA data analyses. All the analyses were carried out on a computer with Intel® Pentium® 4 3.4GHz dual processors and 2.0GB RAM. It took 55hrs 16mins to complete simulation analyses for all 180 simulated scenarios and 14hrs 22mins for GWA data analyses. This indicates that the computation time required for simulation and empirical data analyses is acceptable, and thus our method is practical for bivariate analyses in the field of candidate gene or GWA studies.

For simplicity, we only conducted bivariate analyses for two correlated traits, with the purpose of demonstrating the advantages of B_EGEE over univariate analyses. It should be noted that our method can be easily extended to the situation of multiple correlated traits (>2) with more complex distributions, e.g., Poisson distribution, Gamma distribution, etc. The basic idea of this extension is rather straightforward. One only needs to set up the respective GLM for each phenotypic variable via a specific link function according to its distribution feature [McCullagh and Nelder 1989]. Similar to B_EGEE, these GLMs then can be combined into a unified framework through the SUR model.

Another potential extension of B_EGEE is to incorporate haplotype information across multiple loci in the models. It has been shown that the haplotype-based strategy may be preferable in some situations [Akey, et al. 2001; Liu, et al. 2007b]. To test phenotype-haplotype association for multiple traits, the substantial modification needed for B_EGEE is to construct a new design matrix Xi which includes the index variables corresponding to the two haplotypes carried by individual i. For uncertainty of haplotypes from phase-unknown data, we may borrow the idea of haplotype trend regression [Zaykin, et al. 2002] in univariate analysis. We may incorporate the posterior probabilities of all uncertainties into the haplotype design matrix rather than most likely haplotype configurations. This will help to avoid a loss of information and potential bias in subsequent analyses [Liu, et al. 2007b].

As verified in simulation studies, Bonferroni correction is more conservative than the Meff method when performing multiple testing corrections. Even so, in univariate analyses of real GWA data, we did not present the Meff-adjusted p_values for multiple testing of the number of traits involved. The main reasons are as follows. On one hand, the Meff method requires the correlation coefficients among phenotypic variables to estimate an effective number of independent variables [Nyholt 2004]. However, it is unfeasible to calculate the correlation between the continuous observation and the binary disease status using the traditional approach for our real data set, although the correlation coefficient is directly given as a parameter setting in the simulation studies; On the other hand, we borrowed the estimate of the association parameter ζ obtained in B_EGEE as an alternative of correlation coefficient between two phenotypic variables. At 100 most significant loci, the estimates of ζ vary from 0.016 through 0.023. Given these estimates of ζ, we implemented the Meff method to estimate the effective number of independent tests, and the results range within the interval [1.9995, 1.9997]. This means the number of independent phenotypes is quite close to the true number of phenotype. This suggests that Meff-adjusted p_values are almost the same as those by Bonferroni correction. Hence, to simplify our presentation, we just presented p-value comparisons between the proposed B_EGEE test and the adjusted univariate analyses using Bonferroni correction rather than using the Meff method in Fig. 1b.

For population-based association studies, population stratification may cause false positive/negative results [Weiss, et al. 2001]. In both simulation and real data analyses conducted herein using B_EGEE, we did not consider or control population stratification. However, we point out that this can be done using the methods developed and commonly used for univariate analysis, including Structure [Pritchard, et al. 2000] and Principal Component methods [Price, et al. 2006]. For example, we can incorporate estimates of individual ancestry [Price, et al. 2006; Pritchard, et al. 2000] as covariates in each GLM. It has been shown that use of individual ancestry as a covariate is a promising way for adjusting the effects due to population stratification at most marker loci except those with very strong differentiation even among closely related populations and ethnic groups [Chen and Abecasis 2007]. In univariate association analysis, another way to control population stratification is to adjust test statistics using the genomic control method [Roeder, et al. 2006]; however, the direct use of this method for multiple dependent variables is not straightforward. We will explore the feasibility of extending the genomic control approach to multivariate analyses in further endeavors.


We are grateful to the anonymous reviewers for their insightful comments and constructive suggestions that greatly improved our manuscript. We thank Dr. Yongjun Liu for his help in revising our earlier manuscript and valuable discussions which led to improvements in the paper. Investigators of this work were partially supported by grants from NIH (R01 AR050496-01, R21 AG027110, R01 AG026564, and P50 AR055081). The study also benefited from grants from National Science Foundation of China, Huo Ying Dong Education Foundation, HuNan Province, Xi’an Jiaotong University, and the Ministry of Education of China.


How a binary observation is determined by a latent normal distribution variable

Assume a SNP locus with two alleles: A1 and A2, where A2 is the risk allele for the binary trait T2. To keep consistent with previous description, we adopt the same symbols a2,j (j = 0, 1, 2) as before to denote the genotype values corresponding to three genotypes A1A1, A1A2 and A2A2. Under the additive model, we have


For the latent normal variable z2, we assume it follows a mixture of normal distribution z2 ~ N (a2,j, 1) (j = 0, 1, 2 ) conditionally on the genotype carried by each individual. Accordingly, we have the following expressions according to the theory of threshold theory:




where, Φ(•) is the cumulative distribution function of a standard normal variable.

Given the value of penetrance function f0 for the genotype A1A1 and heterozygote GRR, the status of the binary variable T2 can be determined conditionally by the genotype of each individual through the respective penetrance function fi based on the above Equations A1~A4.


  • Akey J, Jin L, Xiong M. Haplotypes vs single marker linkage disequilibrium tests: what do we gain? Eur J Hum Genet. 2001;9(4):291–300. [PubMed]
  • Allison DB, Thiel B, Jean P, Elston RC, Infante MC, Schork NJ. Multiple phenotype modeling in gene-mapping studies of quantitative traits: power advantages. Am J Hum Genet. 1998;63(4):1190–201. [PubMed]
  • Catalano PJ. Bivariate modelling of clustered continuous and ordered categorical outcomes. Stat Med. 1997;16(8):883–900. [PubMed]
  • Chen WM, Abecasis GR. Family-based association tests for genomewide association scans. Am J Hum Genet. 2007;81(5):913–26. [PubMed]
  • Deng HW, Chen WM, Recker RR. QTL fine mapping by measuring and testing for Hardy-Weinberg and linkage disequilibrium at a series of linked marker loci in extreme samples of populations. Am J Hum Genet. 2000;66(3):1027–45. [PubMed]
  • Evans DM. The power of multivariate quantitative-trait loci linkage analysis is influenced by the correlation between variables. Am J Hum Genet. 2002;70(6):1599–602. [PubMed]
  • Falconer DS, Mackay TFC. Introduction to Quantitative Genetics. Longmans Green; Harlow, Essex, UK: 1996.
  • Gueorguieva RV, Agresti A. A Correlated Probit Model for Joint Modeling of Clustered Binary and Continuous Responses. Journal of the American Statistical Association. 2001;96(455):1102–1112.
  • Gueorguieva RV, Sanacora G. Joint analysis of repeatedly observed continuous and ordinal measures of disease severity. Stat Med. 2006;25(8):1307–22. [PubMed]
  • Guney E, Kisakol G, Ozgen G, Yilmaz C, Yilmaz R, Kabalak T. Effect of weight loss on bone metabolism: comparison of vertical banded gastroplasty and medical intervention. Obes Surg. 2003;13(3):383–8. [PubMed]
  • Hall DB. On the Application of Extended Quasi-Likelihood to the Clustered Data Case. The Canadian Journal of Statistics. 2001;29(1):77–97.
  • Hall DB, Severini TA. Extended generalized estimating equations for clustered data. Journal of the American Statistical Association. 1998;93(444):1365–1375.
  • Huang J, Jiang Y. Genetic linkage analysis of a dichotomous trait incorporating a tightly linked quantitative trait in affected sib pairs. Am J Hum Genet. 2003;72(4):949–60. [PubMed]
  • Jiang C, Zeng ZB. Multiple trait analysis of genetic mapping for quantitative trait loci. Genetics. 1995;140(3):1111–27. [PubMed]
  • Lange C, Silverman EK, Xu X, Weiss ST, Laird NM. A multivariate family-based association test using generalized estimating equations: FBAT-GEE. Biostatistics. 2003;4(2):195–206. [PubMed]
  • Lange C, Whittaker JC. Mapping quantitative trait Loci using generalized estimating equations. Genetics. 2001;159(3):1325–37. [PubMed]
  • Lange C, Whittaker JC, Macgregor AJ. Generalized estimating equations: A hybrid approach for mean parameters in multivariate regression models. Statistical Modelling. 2002;2(3):163–181.
  • Lencz T, Morgan TV, Athanasiou M, Dain B, Reed CR, Kane JM, Kucherlapati R, Malhotra AK. Converging evidence for a pseudoautosomal cytokine receptor gene locus in schizophrenia. Mol Psychiatry. 2007;12(6):572–80. [PubMed]
  • Liang K-Y, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika. 1986;73(1):13–22.
  • Liang K-Y, Zeger SL, Qaqish B. Multivariate Regression Analyses for Categorical Data. Journal of the Royal Statistical Society. Series B (Methodological) 1992;54(1):3–40.
  • Liu J, Liu Y, Liu X, Deng HW. Bayesian mapping of quantitative trait loci for multiple complex traits with the use of variance components. Am J Hum Genet. 2007a;81(2):304–20. [PubMed]
  • Liu J, Papasian C, Deng HW. Incorporating single-locus tests into haplotype cladistic analysis in case-control studies. PLoS Genet. 2007b;3(3):e46. [PubMed]
  • McCullagh PJ, Nelder JA. Generalized Linear Models. Chapman and Hall; London: 1989.
  • Nelder JA, Pregibon D. An extended quasi-likelihood function. Biometrika. 1987;74:221–231.
  • Nyholt DR. A simple correction for multiple testing for single-nucleotide polymorphisms in linkage disequilibrium with each other. Am J Hum Genet. 2004;74(4):765–9. [PubMed]
  • Prentice RL, Zhao LP. Estimating Equations for Parameters in Means and Covariances of Multivariate Discrete and Continuous Responses. Biometrics. 1991;47(3):825–839. [PubMed]
  • Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, Reich D. Principal components analysis corrects for stratification in genome-wide association studies. Nat Genet. 2006;38(8):904–9. [PubMed]
  • Pritchard JK, Stephens M, Rosenberg NA, Donnelly P. Association mapping in structured populations. Am J Hum Genet. 2000;67(1):170–81. [PubMed]
  • Radak TL. Caloric restriction and calcium’s effect on bone metabolism and body composition in overweight and obese premenopausal women. Nutr Rev. 2004;62(12):468–81. [PubMed]
  • Regan MM, Catalano PJ. Regression models and risk estimation for mixed discrete and continuous outcomes in developmental toxicology. Risk Anal. 2000;20(3):363–76. [PubMed]
  • Risch N, Merikangas K. The future of genetic studies of complex human diseases. Science. 1996;273(5281):1516–7. [PubMed]
  • Rochon J. Analyzing Bivariate Repeated Measures for Discrete and Continuous Outcome Variables. Biometrics. 1996;52(2):740–750. [PubMed]
  • Roeder K, Bacanu SA, Wasserman L, Devlin B. Using linkage genome scans to improve power of association in genome scans. Am J Hum Genet. 2006;78(2):243–52. [PubMed]
  • Verzilli CJ, Stallard N, Whittaker JC. Bayesian modelling of multivariate quantitative traits using seemingly unrelated regressions. Genet Epidemiol. 2005;28(4):313–25. [PubMed]
  • Weiss ST, Silverman EK, Palmer LJ. Case-control association studies in pharmacogenetics. Pharmacogenomics J. 2001;1(3):157–8. [PubMed]
  • Williams JT, Van Eerdewegh P, Almasy L, Blangero J. Joint multipoint linkage analysis of multivariate qualitative and quantitative traits. I. Likelihood formulation and simulation results. Am J Hum Genet. 1999;65(4):1134–47. [PubMed]
  • Zaykin DV, Westfall PH, Young SS, Karnoub MA, Wagner MJ, Ehm MG. Testing association of statistically inferred haplotypes with discrete and continuous traits in samples of unrelated individuals. Hum Hered. 2002;53(2):79–91. [PubMed]
  • Zeger SL, Liang KY. An overview of methods for the analysis of longitudinal data. Stat Med. 1992;11(14-15):1825–39. [PubMed]
  • Zellner A. An Efficient Method of Estimating Seemingly Unrelated Regressions and Tests for Aggregation Bias. Journal of the American Statistical Association. 1962;57(298):348–368.
  • Zhao LJ, Jiang H, Papasian CJ, Maulik D, Drees B, Hamilton J, Deng HW. Correlation of obesity and osteoporosis: effect of fat mass on the determination of osteoporosis. J Bone Miner Res. 2008;23(1):17–29. [PubMed]
  • Zhao LJ, Liu YJ, Liu PY, Hamilton J, Recker RR, Deng HW. Relationship of obesity with osteoporosis. J Clin Endocrinol Metab. 2007;92(5):1640–6. [PMC free article] [PubMed]
  • Zhao LP, Prentice RL. Correlated binary regression using a quadratic exponential model. Biometrika. 1990;77(3):642–648.
  • Ziegler A, Kastner C, Blettner M. The Generalised Estimating Equations: An Annotated Bibliography. Biometrical Journal. 1998;40(2):115–139.
  • Zorn CJW. Generalized Estimating Equation Models for Correlated Data: A Review with Applications. American Journal of Political Science. 2001;45(2):470–490.