# Associating Multivariate Traits with Genetic Variants Using Collapsing and Kernel Methods with Pedigree- or Population-Based Studies

Authors: Li-Chu Chien
Journal: Computational and Mathematical Methods in Medicine
 (2021)

---

## Abstract

In genetic association analysis, several relevant phenotypes or multivariate traits with different types of components are usually collected to study complex or multifactorial diseases. Over the past few years, jointly testing for association between multivariate traits and multiple genetic variants has become more popular because it can increase statistical power to identify causal genes in pedigree- or population-based studies. However, most of the existing methods mainly focus on testing genetic variants associated with multiple continuous phenotypes. In this investigation, we develop a framework for identifying the pleiotropic effects of genetic variants on multivariate traits by using collapsing and kernel methods with pedigree- or population-structured data. The proposed framework is applicable to the burden test, the kernel test, and the omnibus test for autosomes and the X chromosome. The proposed multivariate trait association methods can accommodate continuous phenotypes or binary phenotypes and further can adjust for covariates. Simulation studies show that the performance of our methods is satisfactory with respect to the empirical type I error rates and power rates in comparison with the existing methods.


---

## Body

## 1. Introduction

Genome-wide association studies (GWAS) intend to find genetic variants such as single nucleotide polymorphisms (SNPs) associated with common traits or with complex diseases [1, 2]. Association studies, where the correlation relationship between a genetic variant and a trait is evaluated, are helpful for mapping genes influencing complex diseases [3]. In the study of complex diseases, data on several correlated phenotypes or a multivariate phenotype with several components are often collected to get a better understanding of the disease [1, 3, 4]. Multivariate correlated traits are influenced through multiple variants simultaneously. Therefore, by a suitable joint or multivariate analysis framework of multivariate traits, we can not only gain more statistical power to identify pleiotropic effects of genetic variants on multivariate traits [3, 5–12] but also can further understand the genetic architecture of the disease of interest [5, 13]. Thus, recently, the joint analysis of multivariate traits has become popular because it can increase statistical power over analyzing only one trait at a time [1, 4].Several statistical methods have been developed to identify the association between multivariate traits and a genetic variant [1, 5]. Current multivariate methods can be classified into three groups [1, 2, 5]: regression methods [14–16], variable reduction methods [11, 13, 17, 18], and combining analysis [9, 19–23]. However, many of the existing methods for multivariate association analysis cannot be straightaway extended to rare variant analyses, due to their enormous numbers causing the problems of multiple comparison or multiple testing and their low minor allele frequencies [2, 5, 24]. Moreover, sparsity of data could lead to problems on estimating regression parameters and fitting regression models [2]. Hence, it is necessary for proposing statistical methods for identifying the association between multivariate traits and multiple genetic variants (common and/or rare variants) [5]. In recent years, various statistical techniques have been proposed for this purpose in GWAS [8, 17, 25–27]. Furthermore, several approaches have been extendedly developed for the investigation of rare variants associated with multivariate traits [2, 28–38].Although these new developments keep many benefits, existing methods have some potential limitations [39]. Most current methods are constructed under some specific assumptions about the effects of genetic variants on multivariate traits [39]. These current approaches suffer a severe loss in power once the model assumptions are violated [26, 39].In this investigation, we develop the statistical methods for identifying pleiotropic effects of genetic variants on multivariate traits using collapsing and kernel methods with pedigree- or population-structured data. The proposed multivariate trait association method is able to handle binary phenotypes or continuous phenotypes and further can adjust for covariates. Moreover, the proposed multivariate trait association method not only can leverage the dependence on the phenotypes but also can account for the sample relatedness in the pedigree-based or population-based structured data.The rest of the article is organized as follows. In the materials and methods section, we construct the multivariate effect model using the joint GEE model formulation (JGEE) [40]. We apply the JGEE to pedigree- or population-structured data and introduce a retrospective framework for analyzing multivariate traits in genetic association studies. The proposed framework is applicable to the burden test, the kernel test, and the omnibus test for autosomes and the X chromosome. In the simulation studies, we examine the finite sample size performance of the proposed multivariate association methods and evaluate the comparison results with the existing method, Multi-SKAT [39]. Concluding remarks and future possibilities for continuity are given in the conclusion section and the limitation section.


## 2. Materials and Methods

### 2.1. Notations

To describe the proposed multivariate trait association method based on the pedigree- or population-based structured data, we suppose that there areN independent pedigrees and each pedigree has ni subjects. We assume that the ni subjects have been sequenced in a genetic region of interest (e.g., a gene) that contains p variants. Let yik=yi1k,yi2k,⋯,yinikT be the ni×1 phenotype vector for the kth phenotype of the ith pedigree. Let yi=yi1,yi2,⋯,yiK be the ni×K×1 response vector for the K phenotypes that we are interested in. Let xim=xi1m,xi2m,⋯,xinimT be the ni×1 vector for the mth covariate of the ith pedigree. Let xi=xi0,xi1,⋯,xiq be the ni×q+1 covariate matrix for the q+1 nongenetic covariates that we want to adjust for. Let αk=α0k,α1k,⋯,αqkT be the q+1×1 vector of regression coefficients of the q+1 nongenetic covariates with the element αmk being the effect of the mth covariate on the kth trait. Let gi=gi1,gi2,⋯,gip be the ni×p genetic matrix for p genetic variants in a target region of interest where gil=gi1l,gi2l,⋯,ginilT is the ni×1 vector for a genetic variant l (gijl=0, 1, or 2 for 0, 1, or 2 copies of the minor allele, respectively). Let βk=β1k,β2k,⋯,βpkT be the p×1 vector of regression coefficients of the p genetic variants with the element βlk being the effect of the lth genetic variant on the kth trait.


### 2.2. Multitrait Regression-Based Tests for Pedigree Data

We letXi=IK⊗xi be the ni×K×q+1×K covariate matrix and Gi=IK⊗gi be the ni×K×p×K genotype matrix for the ith pedigree where IK is an identity matrix of dimension K×K and ⊗ stands for the Kronecker product. According to the generalized linear model [41], we assume that the marginal density of yijk is fyijk=expyijkθijk−aθijk+byijk/ϕ with two moments, μijk=Eyijk=∂aθijk/∂θijk and Varyijk=∂μijk/∂θijkϕ, where ϕ is a scale parameter. Let θi=θi1T,θi2T,⋯,θiKTT be the ni×K×1 vector with the components θik=θi1k,θi2k,⋯,θinikT,k=1,2,⋯,K and ηi=ηi1T,ηi2T,⋯,ηiKTT be the ni×K×1 vector with the components ηik=ηi1k,ηi2k,⋯,ηinikT=xiαk+giβk,k=1,2,⋯,K for the kth trait of the ith pedigree.Based on the joint GEE model formulation [40], we construct the multivariate linear model for describing the association relationship between K correlated traits and genetic variants, which is given as follows:
(1)μi=gηiandηi=Xiα+Giβ,where g−1• is the inverse function of g• and is a response-specific link function [40], μi=μi1T,μi2T,⋯,μiKTT is the ni×K×1 vector of the expected mean of the multivariate traits yi=yi1T,yi2T,⋯,yiKTT, α=α1T,α2T,⋯,αKTT is the q+1×K×1 vector of regression coefficients of the q+1 nongenetic covariates for the K correlated traits, and β=β1T,β2T,⋯,βKTT is the p×K×1 vector of regression coefficients of the p genetic variants for the K correlated traits.LetRniφ and RKγ be the ni×ni within-in cluster correlation matrix and the K×K multivariate-response cluster correlation matrix, which depend on a vector of parameters φ and γ, respectively. The ni×K×ni×K working (or approximate) covariance matrix of yi is given by [40].
(2)Vi=Ai1/2RKγ⊗RniφAi1/2ϕ,where Ai=diagAi1,Ai2,⋯,AiK is a ni×K×ni×K block diagonal matrix with the components Aik=diag∂μi1k/∂θi1k,∂μi2k/∂θi2k,⋯,∂μinik/∂θinik,k=1,2,⋯,K being the ni×ni diagonal matrices. According to equation (1), under the null hypothesis of no association between genotypes and phenotypes, we propose the multivariate association methods including the homogeneous kernel statistic (HoK), the heterogeneous kernel statistic (HeK), and burden test (BT). Moreover, we propose the homogeneous omnibus test (HoO) and heterogeneous omnibus test (HeO) by combining the HoK with the BT and by combining the HeK with the BT, respectively.


#### 2.2.1. Kernel Statistic

We letH be a p×p correlation matrix of genotype scores with element Hll′ for markers l and l′. Let ml denote the minor allele frequency (MAF) of the lth marker. Let Si=V^i−1yi−μ^i=Si1T,Si2T,⋯,SiKTT be the ni×K×1 vector of the standard residuals with components Sik=Si1k,Si2k,⋯,SinikT,k=1,2,⋯,K, where V^i−1 is the inverse matrix of V^i. Here, V^i and μ^i are the estimates of Vi and μi. Here and henceforth, all estimates are calculated based on the null hypothesis of the genetic effects β equal to zero. All unknown parameters and the working within-in and multivariate-response cluster correlation matrices are estimated by the R package JGEE [42].(1) Homogeneous Kernel Statistic. We suppose that wl is a marker-specific weight of the lth variant and assume that the genetic effects on the K different phenotypes are homogeneous (i.e., β1=β2=⋯=βK. Based on the JGEE model with the genotype as random variables considered, we propose the homogeneous quadratic (kernel) association statistic (HoK) as follows:
(3)κHo=∑l=1pwl∑k=1K∑i=1NgilTΔ∧ikA∧ikSik2=∑l=1pwl∑k=1KZlk2=∑l=1pwlZl2=∑l=1pZ~l2,where Z~l=wlZl, Zl=∑k=1KZlk, Zlk=∑i=1NgilTΔ^ikA^ikSik,A^ik is the estimate of Aik, and Δ^ik is the estimate of Δik=diag∂θi1k/∂ηi1k,∂θi2k/∂ηi2k,⋯,∂θinik/∂ηinik that is a ni×ni diagonal matrix for the kth phenotype of the ith pedigree. The null distribution of κHo asymptotically follows a mixture chi-square distribution ∑l=1pλlχl,12, where χl,12s are independent random variables following a chi-square distribution with one degree of freedom and λ1,λ2,⋯,λp are nonzero eigenvalues of the null covariate matrix of Cov0Z~l,Z~l′=2CHowlwl′Hll′ml1−mlml′1−ml′ where CHo=∑i=1N∑k=1KSikTA^ikΔ^ikΩi∑k=1KΔ^ikA^ikSik and Ωi is a ni×ni matrix of genetic correlations for all ni individuals in the ith pedigree, which has the same definition given by Schaid et al. [43] and can be calculated by the R package kinship2 [44]. When the genetic relationship between subjects j and j′ in the ith pedigree is unknown, the elements of the genetic correlation Ωi can be estimated through genomic data [43, 45], and its estimate is given by [43]
(4)Ω^i=1p∑l=1pgijl−2mlgij′l−2ml2ml1−ml.(2) Heterogeneous Kernel Statistic. We assume that the genetic effects on the K different phenotypes are heterogeneous (i.e., β1≠β2≠⋯≠βK. The heterogeneous quadratic (kernel) association statistic (HeK) is defined by
(5)κHe=∑l=1p∑k=1Kwlk∑i=1NgilTΔ∧ikA∧ikSik2=∑l=1p∑k=1KwlkZlk2=∑l=1p∑k=1KZ~lk2,where Z~lk=wlkZlk and wlk is a marker-specific weight of the lth variant of the kth trait. The null distribution of κHe asymptotically follows a mixture chi-square distribution ∑l=1p×Kλlχl,12, where χl,12s are independent random variables following a chi-square distribution with one degree of freedom, and λ1,λ2,⋯,λp×K are nonzero eigenvalues of the null covariate matrix of Cov0Z~lk,Z~l′k′=2CHewlkwl′k′Hll′ml1−mlml′1−ml′, where CHe=∑i=1NSikTA^ikΔ^ikΩiΔ^ik′A^ik′Sik′.Theoreticalp values of κHo and κHe are approximately calculated by Kuonen’s saddlepoint method [46] and obtained by the R package pchisqsum. A theory for the derivation of the HoK test κHo and the HeK test κHe is shown in Appendix S1.


#### 2.2.2. Burden Test

We letg~iT=∑l=1pwlgilT be a weighted average of genotype scores for the ith pedigree. On the basis of the HoK test κHo and the HeK test κHe in equations (3) and (5) with the same marker-specific weight of the lth variant for each trait k (i.e., wl=wlk,k=1,2,⋯,K), we propose the burden test (BT) as follows:
(6)BT=∑i=1N∑k=1KSikTA∧ikΔ∧ikg~i2∑i=1N∑k=1KSikTA^ikΔ^ikCov0g~i∑k=1KΔ^ikA^ikSik,where the null covariance matrix of g~i is given by
(7)Cov0g~i=Cov0∑l=1pwlgil=∑l=1pwl2Cov0gil,gil+2∑l=1p∑l′=l+1pwlwl′Cov0gil,gil′=Ωi∑l=1p∑l′=1p2wlwl′Hll′ml1−mlml′1−ml′.Then,(8)BT=∑i=1N∑k=1KSikTA∧ikΔ∧ikg~i22∑l=1p∑l′=1pwlwl′Hll′ml1−mlml′1−ml′CHo.The null distribution ofBT asymptotically follows a chi-square distribution with one degree of freedom.


#### 2.2.3. Omnibus Test

LetpHo, pHe, and pBT denote the p values obtained by the HoK, HeK, and BT statistics. Based on the idea of the p value combination method through the Cauchy distribution [47–49], we propose the homogeneous omnibus test (HoO) and heterogeneous omnibus test (HeO).(1) Homogeneous Omnibus Test. Combining the pHo with the pBT, we construct the homogeneous omnibus test (HoO) as follows:
(9)OHo=−12FC−1pHo+FC−1pBT,where FC−1 stands for the inverse cumulative distribution function of the standard Cauchy distribution.(2) Heterogeneous Omnibus Test. Combining the pHe with the pBT, we construct the heterogeneous omnibus test (HeO) as follows:
(10)OHe=−12FC−1pHe+FC−1pBT.The null distributions of theOHo test and the OHe test asymptotically follow a standard Cauchy distribution [47–49]. The p values of the OHo test and the OHe test are calculated by the R package RNOmni [50].The kernel statistic, the burden test, and the omnibus test are also applicable to the X chromosome. Additional technical information for extensions to the X chromosome is shown in AppendixS2.


## 2.1. Notations

To describe the proposed multivariate trait association method based on the pedigree- or population-based structured data, we suppose that there areN independent pedigrees and each pedigree has ni subjects. We assume that the ni subjects have been sequenced in a genetic region of interest (e.g., a gene) that contains p variants. Let yik=yi1k,yi2k,⋯,yinikT be the ni×1 phenotype vector for the kth phenotype of the ith pedigree. Let yi=yi1,yi2,⋯,yiK be the ni×K×1 response vector for the K phenotypes that we are interested in. Let xim=xi1m,xi2m,⋯,xinimT be the ni×1 vector for the mth covariate of the ith pedigree. Let xi=xi0,xi1,⋯,xiq be the ni×q+1 covariate matrix for the q+1 nongenetic covariates that we want to adjust for. Let αk=α0k,α1k,⋯,αqkT be the q+1×1 vector of regression coefficients of the q+1 nongenetic covariates with the element αmk being the effect of the mth covariate on the kth trait. Let gi=gi1,gi2,⋯,gip be the ni×p genetic matrix for p genetic variants in a target region of interest where gil=gi1l,gi2l,⋯,ginilT is the ni×1 vector for a genetic variant l (gijl=0, 1, or 2 for 0, 1, or 2 copies of the minor allele, respectively). Let βk=β1k,β2k,⋯,βpkT be the p×1 vector of regression coefficients of the p genetic variants with the element βlk being the effect of the lth genetic variant on the kth trait.


## 2.2. Multitrait Regression-Based Tests for Pedigree Data

We letXi=IK⊗xi be the ni×K×q+1×K covariate matrix and Gi=IK⊗gi be the ni×K×p×K genotype matrix for the ith pedigree where IK is an identity matrix of dimension K×K and ⊗ stands for the Kronecker product. According to the generalized linear model [41], we assume that the marginal density of yijk is fyijk=expyijkθijk−aθijk+byijk/ϕ with two moments, μijk=Eyijk=∂aθijk/∂θijk and Varyijk=∂μijk/∂θijkϕ, where ϕ is a scale parameter. Let θi=θi1T,θi2T,⋯,θiKTT be the ni×K×1 vector with the components θik=θi1k,θi2k,⋯,θinikT,k=1,2,⋯,K and ηi=ηi1T,ηi2T,⋯,ηiKTT be the ni×K×1 vector with the components ηik=ηi1k,ηi2k,⋯,ηinikT=xiαk+giβk,k=1,2,⋯,K for the kth trait of the ith pedigree.Based on the joint GEE model formulation [40], we construct the multivariate linear model for describing the association relationship between K correlated traits and genetic variants, which is given as follows:
(1)μi=gηiandηi=Xiα+Giβ,where g−1• is the inverse function of g• and is a response-specific link function [40], μi=μi1T,μi2T,⋯,μiKTT is the ni×K×1 vector of the expected mean of the multivariate traits yi=yi1T,yi2T,⋯,yiKTT, α=α1T,α2T,⋯,αKTT is the q+1×K×1 vector of regression coefficients of the q+1 nongenetic covariates for the K correlated traits, and β=β1T,β2T,⋯,βKTT is the p×K×1 vector of regression coefficients of the p genetic variants for the K correlated traits.LetRniφ and RKγ be the ni×ni within-in cluster correlation matrix and the K×K multivariate-response cluster correlation matrix, which depend on a vector of parameters φ and γ, respectively. The ni×K×ni×K working (or approximate) covariance matrix of yi is given by [40].
(2)Vi=Ai1/2RKγ⊗RniφAi1/2ϕ,where Ai=diagAi1,Ai2,⋯,AiK is a ni×K×ni×K block diagonal matrix with the components Aik=diag∂μi1k/∂θi1k,∂μi2k/∂θi2k,⋯,∂μinik/∂θinik,k=1,2,⋯,K being the ni×ni diagonal matrices. According to equation (1), under the null hypothesis of no association between genotypes and phenotypes, we propose the multivariate association methods including the homogeneous kernel statistic (HoK), the heterogeneous kernel statistic (HeK), and burden test (BT). Moreover, we propose the homogeneous omnibus test (HoO) and heterogeneous omnibus test (HeO) by combining the HoK with the BT and by combining the HeK with the BT, respectively.


### 2.2.1. Kernel Statistic

We letH be a p×p correlation matrix of genotype scores with element Hll′ for markers l and l′. Let ml denote the minor allele frequency (MAF) of the lth marker. Let Si=V^i−1yi−μ^i=Si1T,Si2T,⋯,SiKTT be the ni×K×1 vector of the standard residuals with components Sik=Si1k,Si2k,⋯,SinikT,k=1,2,⋯,K, where V^i−1 is the inverse matrix of V^i. Here, V^i and μ^i are the estimates of Vi and μi. Here and henceforth, all estimates are calculated based on the null hypothesis of the genetic effects β equal to zero. All unknown parameters and the working within-in and multivariate-response cluster correlation matrices are estimated by the R package JGEE [42].(1) Homogeneous Kernel Statistic. We suppose that wl is a marker-specific weight of the lth variant and assume that the genetic effects on the K different phenotypes are homogeneous (i.e., β1=β2=⋯=βK. Based on the JGEE model with the genotype as random variables considered, we propose the homogeneous quadratic (kernel) association statistic (HoK) as follows:
(3)κHo=∑l=1pwl∑k=1K∑i=1NgilTΔ∧ikA∧ikSik2=∑l=1pwl∑k=1KZlk2=∑l=1pwlZl2=∑l=1pZ~l2,where Z~l=wlZl, Zl=∑k=1KZlk, Zlk=∑i=1NgilTΔ^ikA^ikSik,A^ik is the estimate of Aik, and Δ^ik is the estimate of Δik=diag∂θi1k/∂ηi1k,∂θi2k/∂ηi2k,⋯,∂θinik/∂ηinik that is a ni×ni diagonal matrix for the kth phenotype of the ith pedigree. The null distribution of κHo asymptotically follows a mixture chi-square distribution ∑l=1pλlχl,12, where χl,12s are independent random variables following a chi-square distribution with one degree of freedom and λ1,λ2,⋯,λp are nonzero eigenvalues of the null covariate matrix of Cov0Z~l,Z~l′=2CHowlwl′Hll′ml1−mlml′1−ml′ where CHo=∑i=1N∑k=1KSikTA^ikΔ^ikΩi∑k=1KΔ^ikA^ikSik and Ωi is a ni×ni matrix of genetic correlations for all ni individuals in the ith pedigree, which has the same definition given by Schaid et al. [43] and can be calculated by the R package kinship2 [44]. When the genetic relationship between subjects j and j′ in the ith pedigree is unknown, the elements of the genetic correlation Ωi can be estimated through genomic data [43, 45], and its estimate is given by [43]
(4)Ω^i=1p∑l=1pgijl−2mlgij′l−2ml2ml1−ml.(2) Heterogeneous Kernel Statistic. We assume that the genetic effects on the K different phenotypes are heterogeneous (i.e., β1≠β2≠⋯≠βK. The heterogeneous quadratic (kernel) association statistic (HeK) is defined by
(5)κHe=∑l=1p∑k=1Kwlk∑i=1NgilTΔ∧ikA∧ikSik2=∑l=1p∑k=1KwlkZlk2=∑l=1p∑k=1KZ~lk2,where Z~lk=wlkZlk and wlk is a marker-specific weight of the lth variant of the kth trait. The null distribution of κHe asymptotically follows a mixture chi-square distribution ∑l=1p×Kλlχl,12, where χl,12s are independent random variables following a chi-square distribution with one degree of freedom, and λ1,λ2,⋯,λp×K are nonzero eigenvalues of the null covariate matrix of Cov0Z~lk,Z~l′k′=2CHewlkwl′k′Hll′ml1−mlml′1−ml′, where CHe=∑i=1NSikTA^ikΔ^ikΩiΔ^ik′A^ik′Sik′.Theoreticalp values of κHo and κHe are approximately calculated by Kuonen’s saddlepoint method [46] and obtained by the R package pchisqsum. A theory for the derivation of the HoK test κHo and the HeK test κHe is shown in Appendix S1.


### 2.2.2. Burden Test

We letg~iT=∑l=1pwlgilT be a weighted average of genotype scores for the ith pedigree. On the basis of the HoK test κHo and the HeK test κHe in equations (3) and (5) with the same marker-specific weight of the lth variant for each trait k (i.e., wl=wlk,k=1,2,⋯,K), we propose the burden test (BT) as follows:
(6)BT=∑i=1N∑k=1KSikTA∧ikΔ∧ikg~i2∑i=1N∑k=1KSikTA^ikΔ^ikCov0g~i∑k=1KΔ^ikA^ikSik,where the null covariance matrix of g~i is given by
(7)Cov0g~i=Cov0∑l=1pwlgil=∑l=1pwl2Cov0gil,gil+2∑l=1p∑l′=l+1pwlwl′Cov0gil,gil′=Ωi∑l=1p∑l′=1p2wlwl′Hll′ml1−mlml′1−ml′.Then,(8)BT=∑i=1N∑k=1KSikTA∧ikΔ∧ikg~i22∑l=1p∑l′=1pwlwl′Hll′ml1−mlml′1−ml′CHo.The null distribution ofBT asymptotically follows a chi-square distribution with one degree of freedom.


### 2.2.3. Omnibus Test

LetpHo, pHe, and pBT denote the p values obtained by the HoK, HeK, and BT statistics. Based on the idea of the p value combination method through the Cauchy distribution [47–49], we propose the homogeneous omnibus test (HoO) and heterogeneous omnibus test (HeO).(1) Homogeneous Omnibus Test. Combining the pHo with the pBT, we construct the homogeneous omnibus test (HoO) as follows:
(9)OHo=−12FC−1pHo+FC−1pBT,where FC−1 stands for the inverse cumulative distribution function of the standard Cauchy distribution.(2) Heterogeneous Omnibus Test. Combining the pHe with the pBT, we construct the heterogeneous omnibus test (HeO) as follows:
(10)OHe=−12FC−1pHe+FC−1pBT.The null distributions of theOHo test and the OHe test asymptotically follow a standard Cauchy distribution [47–49]. The p values of the OHo test and the OHe test are calculated by the R package RNOmni [50].The kernel statistic, the burden test, and the omnibus test are also applicable to the X chromosome. Additional technical information for extensions to the X chromosome is shown in AppendixS2.


## 2.2.1. Kernel Statistic

We letH be a p×p correlation matrix of genotype scores with element Hll′ for markers l and l′. Let ml denote the minor allele frequency (MAF) of the lth marker. Let Si=V^i−1yi−μ^i=Si1T,Si2T,⋯,SiKTT be the ni×K×1 vector of the standard residuals with components Sik=Si1k,Si2k,⋯,SinikT,k=1,2,⋯,K, where V^i−1 is the inverse matrix of V^i. Here, V^i and μ^i are the estimates of Vi and μi. Here and henceforth, all estimates are calculated based on the null hypothesis of the genetic effects β equal to zero. All unknown parameters and the working within-in and multivariate-response cluster correlation matrices are estimated by the R package JGEE [42].(1) Homogeneous Kernel Statistic. We suppose that wl is a marker-specific weight of the lth variant and assume that the genetic effects on the K different phenotypes are homogeneous (i.e., β1=β2=⋯=βK. Based on the JGEE model with the genotype as random variables considered, we propose the homogeneous quadratic (kernel) association statistic (HoK) as follows:
(3)κHo=∑l=1pwl∑k=1K∑i=1NgilTΔ∧ikA∧ikSik2=∑l=1pwl∑k=1KZlk2=∑l=1pwlZl2=∑l=1pZ~l2,where Z~l=wlZl, Zl=∑k=1KZlk, Zlk=∑i=1NgilTΔ^ikA^ikSik,A^ik is the estimate of Aik, and Δ^ik is the estimate of Δik=diag∂θi1k/∂ηi1k,∂θi2k/∂ηi2k,⋯,∂θinik/∂ηinik that is a ni×ni diagonal matrix for the kth phenotype of the ith pedigree. The null distribution of κHo asymptotically follows a mixture chi-square distribution ∑l=1pλlχl,12, where χl,12s are independent random variables following a chi-square distribution with one degree of freedom and λ1,λ2,⋯,λp are nonzero eigenvalues of the null covariate matrix of Cov0Z~l,Z~l′=2CHowlwl′Hll′ml1−mlml′1−ml′ where CHo=∑i=1N∑k=1KSikTA^ikΔ^ikΩi∑k=1KΔ^ikA^ikSik and Ωi is a ni×ni matrix of genetic correlations for all ni individuals in the ith pedigree, which has the same definition given by Schaid et al. [43] and can be calculated by the R package kinship2 [44]. When the genetic relationship between subjects j and j′ in the ith pedigree is unknown, the elements of the genetic correlation Ωi can be estimated through genomic data [43, 45], and its estimate is given by [43]
(4)Ω^i=1p∑l=1pgijl−2mlgij′l−2ml2ml1−ml.(2) Heterogeneous Kernel Statistic. We assume that the genetic effects on the K different phenotypes are heterogeneous (i.e., β1≠β2≠⋯≠βK. The heterogeneous quadratic (kernel) association statistic (HeK) is defined by
(5)κHe=∑l=1p∑k=1Kwlk∑i=1NgilTΔ∧ikA∧ikSik2=∑l=1p∑k=1KwlkZlk2=∑l=1p∑k=1KZ~lk2,where Z~lk=wlkZlk and wlk is a marker-specific weight of the lth variant of the kth trait. The null distribution of κHe asymptotically follows a mixture chi-square distribution ∑l=1p×Kλlχl,12, where χl,12s are independent random variables following a chi-square distribution with one degree of freedom, and λ1,λ2,⋯,λp×K are nonzero eigenvalues of the null covariate matrix of Cov0Z~lk,Z~l′k′=2CHewlkwl′k′Hll′ml1−mlml′1−ml′, where CHe=∑i=1NSikTA^ikΔ^ikΩiΔ^ik′A^ik′Sik′.Theoreticalp values of κHo and κHe are approximately calculated by Kuonen’s saddlepoint method [46] and obtained by the R package pchisqsum. A theory for the derivation of the HoK test κHo and the HeK test κHe is shown in Appendix S1.


## 2.2.2. Burden Test

We letg~iT=∑l=1pwlgilT be a weighted average of genotype scores for the ith pedigree. On the basis of the HoK test κHo and the HeK test κHe in equations (3) and (5) with the same marker-specific weight of the lth variant for each trait k (i.e., wl=wlk,k=1,2,⋯,K), we propose the burden test (BT) as follows:
(6)BT=∑i=1N∑k=1KSikTA∧ikΔ∧ikg~i2∑i=1N∑k=1KSikTA^ikΔ^ikCov0g~i∑k=1KΔ^ikA^ikSik,where the null covariance matrix of g~i is given by
(7)Cov0g~i=Cov0∑l=1pwlgil=∑l=1pwl2Cov0gil,gil+2∑l=1p∑l′=l+1pwlwl′Cov0gil,gil′=Ωi∑l=1p∑l′=1p2wlwl′Hll′ml1−mlml′1−ml′.Then,(8)BT=∑i=1N∑k=1KSikTA∧ikΔ∧ikg~i22∑l=1p∑l′=1pwlwl′Hll′ml1−mlml′1−ml′CHo.The null distribution ofBT asymptotically follows a chi-square distribution with one degree of freedom.


## 2.2.3. Omnibus Test

LetpHo, pHe, and pBT denote the p values obtained by the HoK, HeK, and BT statistics. Based on the idea of the p value combination method through the Cauchy distribution [47–49], we propose the homogeneous omnibus test (HoO) and heterogeneous omnibus test (HeO).(1) Homogeneous Omnibus Test. Combining the pHo with the pBT, we construct the homogeneous omnibus test (HoO) as follows:
(9)OHo=−12FC−1pHo+FC−1pBT,where FC−1 stands for the inverse cumulative distribution function of the standard Cauchy distribution.(2) Heterogeneous Omnibus Test. Combining the pHe with the pBT, we construct the heterogeneous omnibus test (HeO) as follows:
(10)OHe=−12FC−1pHe+FC−1pBT.The null distributions of theOHo test and the OHe test asymptotically follow a standard Cauchy distribution [47–49]. The p values of the OHo test and the OHe test are calculated by the R package RNOmni [50].The kernel statistic, the burden test, and the omnibus test are also applicable to the X chromosome. Additional technical information for extensions to the X chromosome is shown in AppendixS2.


## 3. Simulation Studies

We conduct the numerical simulation studies to assess the finite sample performance of the proposed methods and evaluate the comparison results with two existing methods, the minimump value SKAT statistic (mPK), and the minimum p value burden statistic (mPB) [39]. The two existing methods are implemented by the R package Multi-SKAT [39]. Based on the similar simulation set-up as those usually considered from existing genetic association tests [39, 43, 51], we investigate the effect of the proposed methods, HoK, HeK, BT, HoO, and HeO, for identifying genetic variants that are associated with multiple traits. We simultaneously generate 10,000 European-like (EUR) and 10,000 admixed African American-like (AA) haplotypes of length 200 kb using a calibrated human demographic model through the COSI software [51, 52]. A 3 kb region is randomly selected in our numerical simulations. We generate a total of 10,000 databases for each simulation scenario in our studies.


### 3.1. Type I Error Rate and Power Simulations

In the heterogeneous population with nuclear family data considered, continuous and binary phenotypes for traitk for individual j in the ith family are generated from the multivariate linear model in equation (1) with K=2 and ni=3. More precisely, continuous and binary phenotypes are generated by the following linear and logit models, respectively:
(11)yi=Xiα+Giβ+εi,(12)logitPyi=1=Xiα+Giβ,where yi=yi1T,yi2TT, Xi=I2⊗xi, Gi=I2⊗gi, α=α1T,α2TT, β=β1T,β2TT, and εi=εi11,εi21,εi31,εi12,εi22,εi32T. Here, the elements xi0=xi10,xi20,xi30T of the covariance matrix xi=xi0,xi1,xi2 is a 3×1 vector of all ones. The elements xi1=xi11,xi21,xi31T of xi are independently generated with an equal probability of being 0 or 1. The elements xi2=xi12,xi22,xi32T of xi are generated from a multivariate normal distribution with a mean of 0.5 and a covariance matrix with diagonal entries of 1 and all off-diagonal entries of 0.1. The regression coefficients of the covariate matrix xi for the kth correlated trait are given by αk=α0k,α1k,α2kT=0.01,0.1,0.1T and αk=α0k,α1k,α2kT=−1.4,0.1,0.1T, respectively, for continuous traits and binary traits for k=1,2.For continuous traits, the error termsεi=εi11,εi21,εi31,εi12,εi22,εi32T in equation (11) follow a multivariate normal distribution having a mean of zero, a within-in cluster correlation matrix (i.e., Corεijk,εij′k) with diagonal entries of 1 and all off-diagonal entries of 0.2 and a subject-across-response correlation matrix (i.e., Corεijk,εij′k′) with diagonal entries of 0.3 and all off-diagonal entries of 0.1. Similarly, binary traits yi in equation (12) are generated with the same within-in cluster correlation matrix (i.e., Coryijk,yij′k) and the same subject-across-response correlation matrix (i.e., Coryijk,yij′k′) as the continuous traits yi in equation (11). These correlated phenotypes are generated by the R package BinNor [53].For type I error simulations, the regression coefficients of genetic variants,β=β1T,β2TT, in equations (11) and (12) are equal to zero under the null hypothesis. For power simulations, under the alternative hypothesis, we simulate that 35% of low variants with the MAF<0.03 are causal. For each setting, either all causal SNPs have a positive effect, or 80% of causal SNPs are positive, and 20% of causal SNPs are negative. The regression coefficients of genetic variants, β=β1T,β2TT, are set by 0.095×∣log10ml∣ or 0.095×log10ml corresponding to the risk or protective variant l,l=1,2,⋯,p [51]. Under the assumption that the genetic effects on the two different phenotypes are heterogeneous (i.e., β1≠β2, the genetic effects β1 for the first traits yi1 are set as described above, while the genetic effects β2 for the second traits yi2 are set by zero. On the other hand, under the assumption that the genetic effects on the two different phenotypes are homogeneous (i.e., β1=β2, the genetic effects β1 and β2 for the first and second traits have the same settings as described above.We simulate 1,400 nuclear families with 800 nuclear families from the European samples and 600 nuclear families from African-American samples. The marker-specific weightwl for variant l is considered as the beta density function wl=Betaml,λ1,λ2 with shape parameters λ1>0 and λ2>0 [51]. To study the effect of the marker-specific weight wl of variant l on the phenotypes, we consider the unweighted marker-specific weight with wl=Betaml,1,1=1 and the weighted marker-specific weight with wl=Betaml,1,25 [51]. The empirical type I error rates based on fifty thousand replicates and the empirical power rates based on two thousand replicates are reported for all simulation results. The “exchangeable” and “unstructured” structures are considered for the working within-cluster and multivariate-response correlation matrices for the proposed methods, HoK, HeK, and BT, respectively.


## 3.1. Type I Error Rate and Power Simulations

In the heterogeneous population with nuclear family data considered, continuous and binary phenotypes for traitk for individual j in the ith family are generated from the multivariate linear model in equation (1) with K=2 and ni=3. More precisely, continuous and binary phenotypes are generated by the following linear and logit models, respectively:
(11)yi=Xiα+Giβ+εi,(12)logitPyi=1=Xiα+Giβ,where yi=yi1T,yi2TT, Xi=I2⊗xi, Gi=I2⊗gi, α=α1T,α2TT, β=β1T,β2TT, and εi=εi11,εi21,εi31,εi12,εi22,εi32T. Here, the elements xi0=xi10,xi20,xi30T of the covariance matrix xi=xi0,xi1,xi2 is a 3×1 vector of all ones. The elements xi1=xi11,xi21,xi31T of xi are independently generated with an equal probability of being 0 or 1. The elements xi2=xi12,xi22,xi32T of xi are generated from a multivariate normal distribution with a mean of 0.5 and a covariance matrix with diagonal entries of 1 and all off-diagonal entries of 0.1. The regression coefficients of the covariate matrix xi for the kth correlated trait are given by αk=α0k,α1k,α2kT=0.01,0.1,0.1T and αk=α0k,α1k,α2kT=−1.4,0.1,0.1T, respectively, for continuous traits and binary traits for k=1,2.For continuous traits, the error termsεi=εi11,εi21,εi31,εi12,εi22,εi32T in equation (11) follow a multivariate normal distribution having a mean of zero, a within-in cluster correlation matrix (i.e., Corεijk,εij′k) with diagonal entries of 1 and all off-diagonal entries of 0.2 and a subject-across-response correlation matrix (i.e., Corεijk,εij′k′) with diagonal entries of 0.3 and all off-diagonal entries of 0.1. Similarly, binary traits yi in equation (12) are generated with the same within-in cluster correlation matrix (i.e., Coryijk,yij′k) and the same subject-across-response correlation matrix (i.e., Coryijk,yij′k′) as the continuous traits yi in equation (11). These correlated phenotypes are generated by the R package BinNor [53].For type I error simulations, the regression coefficients of genetic variants,β=β1T,β2TT, in equations (11) and (12) are equal to zero under the null hypothesis. For power simulations, under the alternative hypothesis, we simulate that 35% of low variants with the MAF<0.03 are causal. For each setting, either all causal SNPs have a positive effect, or 80% of causal SNPs are positive, and 20% of causal SNPs are negative. The regression coefficients of genetic variants, β=β1T,β2TT, are set by 0.095×∣log10ml∣ or 0.095×log10ml corresponding to the risk or protective variant l,l=1,2,⋯,p [51]. Under the assumption that the genetic effects on the two different phenotypes are heterogeneous (i.e., β1≠β2, the genetic effects β1 for the first traits yi1 are set as described above, while the genetic effects β2 for the second traits yi2 are set by zero. On the other hand, under the assumption that the genetic effects on the two different phenotypes are homogeneous (i.e., β1=β2, the genetic effects β1 and β2 for the first and second traits have the same settings as described above.We simulate 1,400 nuclear families with 800 nuclear families from the European samples and 600 nuclear families from African-American samples. The marker-specific weightwl for variant l is considered as the beta density function wl=Betaml,λ1,λ2 with shape parameters λ1>0 and λ2>0 [51]. To study the effect of the marker-specific weight wl of variant l on the phenotypes, we consider the unweighted marker-specific weight with wl=Betaml,1,1=1 and the weighted marker-specific weight with wl=Betaml,1,25 [51]. The empirical type I error rates based on fifty thousand replicates and the empirical power rates based on two thousand replicates are reported for all simulation results. The “exchangeable” and “unstructured” structures are considered for the working within-cluster and multivariate-response correlation matrices for the proposed methods, HoK, HeK, and BT, respectively.


## 4. Results

### 4.1. Empirical Type I Error Rates

Table1 reports the results of a simulation comparison on empirical type I error rates when the phenotypes are considered to be continuous. Table 1 displays that the proposed methods, HoK, HoO, HeK, HeO, and BT, well control the empirical type I error rates regardless of the weight of the marker-specific weight. Similarly, the existing methods, mPK and mPB, have good performance on controlling the empirical type I error rates. Our simulation results show that the seven competing methods, HoK, HoO, HeK, HeO, BT, mPK, and mPB, reasonably control the empirical type I error rates for autosome analyses with continuous traits. The seven competing approaches display similar performance in terms of the empirical type I error rates for the X chromosome analyses with continuous traits (Appendix S3: Table S1).Table 1
Empirical type I errors of the seven competing methods with continuous traits.
Marker-specific weightwlNominal levelWorking correlationMethodHoK3HoOHeKHeOBTmPK4mPBUnweighted marker-specific weight10.05U/U20.048760.049600.050360.052280.049140.043520.04692E/E0.048660.049940.050160.052160.049140.01U/U0.009180.010120.010160.010300.010340.008540.01036E/E0.009240.009940.010080.010220.010280.001U/U0.000780.000820.000860.000700.000840.000840.00088E/E0.000800.000780.000840.000700.000820.0001U/U0.000080.000020.000060.000080.000080.000060.00014E/E0.000060.000020.000060.000080.00008Weighted marker-specific weight0.05U/U0.050300.049980.051580.051340.046960.046040.04536E/E0.050540.050100.051760.051220.047140.01U/U0.009920.009420.010800.009720.008880.009780.01008E/E0.009920.009440.010880.009780.008860.001U/U0.000780.000860.001260.000980.000820.001240.00134E/E0.000760.000880.001220.001020.000800.0001U/U0.000060.000080.000060.000060.000100.000020.00010E/E0.000060.000080.000080.000060.000101The unweighted marker-specific weight is given by wl=Betaml,1,1=1; the weighted marker-specific weight is given by wl=Betaml,1,25. 2U/U represents the structures of the working within-cluster and multivariate-response correlation matrices considered by the unstructured structures; E/E represents the structures of the working within-cluster and multivariate-response correlation matrices considered by the exchangeable structures. 3HoK, HoO, HeK, HeO, and BT are our proposed methods. 4mPK and mPB are executed by the R package Multi-SKAT [39].Table2 reports the empirical type I error rates based on the proposed methods, HoK, HeK, BT, HoO, and HeO, for the binary data. The two existing methods, mPK and mPB, aren’t included for comparison. This reason is that implementing the two existing methods, mPK and mPB, via the R package Multi-SKAT [39], the MPMM (multiple phenotype mixed model) function in the R package PHENIX [54–56] is a necessary tool for this process. However, the MPMM function is suitable for the continuous phenotypes [56] or is suitable for the binary phenotypes with the condition that the number of cases is sufficiently large [39]. In other words, in some sense, the two existing methods, mPK and mPB, are limited to continuous phenotypes [39].Table 2
Empirical type I errors of the five competing methods with binary traits.
Marker-specific weightwlNominal levelWorking correlationMethodHoK3HoOHeKHeOBTUnweighted marker-specific weight10.05U/U20.049440.051540.050860.052800.04952E/E0.049300.051440.050680.053180.049460.01U/U0.009740.009940.009820.010260.01000E/E0.009740.009980.009840.010280.009980.001U/U0.000680.000840.001000.000980.00106E/E0.000660.000840.001020.000940.001040.0001U/U0.000080.000020.000120.000100.00000E/E0.000080.000020.000120.000100.00002Weighted marker-specific weight0.05U/U0.051700.049000.052560.049220.04576E/E0.051680.049200.052320.049300.045560.01U/U0.010280.009760.009960.009720.00886E/E0.010240.009820.009860.009760.008840.001U/U0.001100.000800.000960.000900.00088E/E0.001120.000760.000960.000880.000900.0001U/U0.000040.000080.000100.000120.00006E/E0.000060.000080.000100.000120.000081The unweighted marker-specific weight is given by wl=Betaml,1,1=1; the weighted marker-specific weight is given by wl=Betaml,1,25. 2U/U represents the structures of the working within-cluster and multivariate-response correlation matrices considered by the unstructured structures; E/E represents the structures of the working within-cluster and multivariate-response correlation matrices considered by the exchangeable structures. 3HoK, HoO, HeK, HeO, and BT are our proposed methods.Table2 shows that the proposed methods appropriately control the type I error rates when the marker-specific weight is considered for wl=Betaml,1,1 or wl=Betaml,1,25 for variant l for binary traits. On the other hand, the empirical type I error rates of the proposed methods for X chromosome analyses with binary traits are depicted in Table S2 in Appendix S3. These empirical type I error rates show similar results as that for autosome analyses.In summary, our simulation results show that the proposed multivariate trait association methods, HoK, HoO, HeK, HeO, and BT, have reasonable control of type I error rates for continuous traits or binary traits whether the marker is X chromosomal or autosomal. On the other hand, the existing methods, mPK and mPB, yield well-controlled type I error rates for the autosome analyses or the X chromosome analyses with continuous traits (Table1 or Table S1), regardless of the weight of the marker-specific weight.


### 4.2. Empirical Power

Figure1 exhibits the comparison results of the empirical power rates for the autosome analyses with continuous traits, when the working within-cluster and multivariate-response correlation matrices of the proposed methods, HoK, HeK, and BT, are considered to be exchangeable. As expected, the empirical power rates of the seven competing methods with a weighted marker-specific weight of wl=Betaml,1,25 are higher than that with an unweighted marker-specific weight of wl=Betaml,1,1=1. The heterogeneous kernel statistic (HeK) has slightly greater empirical power rates than other methods, when the genetic effects on the different phenotypes are heterogeneous (i.e., β1≠β2), and causal SNPs have positive effects or negative effects on phenotypes. On the other hand, the existing method, mPB, has bigger empirical power rates, when the genetic effects on the different phenotypes are heterogeneous (i.e., β1≠β2), and all causal SNPs have a positive association on phenotypes. Moreover, the empirical power rates of the homogeneous omnibus test (HoO) are larger than that of the other six competing methods, when the genetic effects on the different phenotypes are homogeneous (i.e., β1=β2). Evidently, the seven competing methods have their respective advantages in identifying the association between genetic effects and multiple continuous traits for autosome analyses.Figure 1
Power comparisons of the seven competing methods with continuous traits for each scenario at the nominal level of 0.001. (a) Unweighted marker-specific weight:wl=Betaml,1,1=1. (b) Weighted marker-specific weight: wl=Betaml,1,25.Similar empirical power rates are obtained from the working within-cluster and multivariate-response correlation matrices of the proposed methods, HoK, HeK, and BT, considered to be unstructured. Hence, these empirical power rates are not shown in order to save space. On the other hand, the seven competing approaches display a similar performance in testing for the X chromosome analyses with continuous traits (AppendixS3: Figure S1).Figure2 exhibits the comparison results of empirical power rates for the autosome analyses with binary traits when the working within-cluster and multivariate-response correlation matrices of the proposed methods, HoK, HeK, and BT, are considered to be exchangeable. As a similar reason for investigating the empirical type I error rates with binary traits, the two existing methods, mPK and mPB, aren’t included for power comparison.Figure 2
Power comparisons of the five competing methods with binary traits for each scenario at the nominal level of 0.001. (a) Unweighted marker-specific weight:wl=Betaml,1,1=1. (b) Weighted marker-specific weight: wl=Betaml,1,25.Figure2 shows that the heterogeneous kernel statistic (HeK) and the heterogeneous omnibus test (HeO) outperform over other methods in terms of the empirical power rates, when the genetic effects on the different phenotypes are heterogeneous (i.e., β1≠β2). On the other hand, the empirical power rates of the homogeneous omnibus test (HoO) are bigger than that of the other competing methods, when the genetic effects on the different phenotypes are homogeneous (i.e., β1=β2). As expected, in general, the heterogeneous kernel statistic (HeK) is more powerful than the homogeneous kernel statistic (HoK), when the genetic effects on the different phenotypes are heterogeneous (i.e., β1≠β2). On the other hand, the homogeneous kernel statistic (HoK) is more powerful than the heterogeneous kernel statistic (HeK), when the genetic effects on the different phenotypes are homogeneous (i.e., β1=β2). In a word, the proposed methods, HoK, HoO, HeK, HeO, and BT, have their respective merits in examining the association between genetic effects and multiple binary traits for autosome analyses.Similarly, when the working within-cluster and multivariate-response correlation matrices of the proposed methods, HoK, HeK, and BT, are considered to be unstructured, the empirical power rates have similar results and thus they are omitted. On the other hand, the empirical power rates of the proposed methods for X chromosome analyses with binary traits are presented in FigureS2 in Appendix S3. These empirical power rates show similar results as that discussed in Figure 2.In summary, the seven competing methods, HoK, HoO, HeK, HeO, BT, mPK, and mPB, have their respective merits in diagnosing whether genetic effects are associated with multiple continuous traits for autosome analyses or the X chromosome analyses. Similarly, the proposed methods, HoK, HoO, HeK, HeO, and BT, have their respective advantages in examining whether there are associations between genetic effects and multiple binary traits for autosome analyses or the X chromosome analyses.To furthermore examine the performance of the proposed methods, additional simulation studies for continuous traits and binary traits are presented in AppendixS4 and Appendix S5 with higher correlations of phenotypes and higher dimensions of phenotypes considered, respectively. In general, these competing methods based on higher correlations of phenotypes or higher dimensions of phenotypes can provide a bigger empirical power rate for the analysis of continuous traits or binary traits. However, we note that these competing methods based on higher correlations of phenotypes or higher dimensions of phenotypes more easily have empirical type I error rate inflation at a smaller nominal level, especially for binary data analysis (Appendix S5: Tables S5-S6 and Appendix S6: Table S7), in comparison with these methods based on lower correlations of phenotypes or lower dimensions of phenotypes. A detailed discussion of these additional simulation results is given in Appendixes S4 and S5.However, we note that the proposed methods have a high computational cost, especially for binary data. Under our simulation setting and framework, we carry out a single simulated data set by using a computer based on one CPU core at 2.1 GHz. The average computational times of the homogeneous and heterogeneous tests with a weighted marker-specific weightwl=Betaml,1,25 under the alternative hypothesis for continuous data are 0.83 and 0.91 minutes, respectively, while that for binary data are 4.77 and 4.80 minutes, respectively. Therefore, in the current version, such a framework algorithm implementation is unsatisfactory for analyzing a large-scale high-dimensional data set in practice.


## 4.1. Empirical Type I Error Rates

Table1 reports the results of a simulation comparison on empirical type I error rates when the phenotypes are considered to be continuous. Table 1 displays that the proposed methods, HoK, HoO, HeK, HeO, and BT, well control the empirical type I error rates regardless of the weight of the marker-specific weight. Similarly, the existing methods, mPK and mPB, have good performance on controlling the empirical type I error rates. Our simulation results show that the seven competing methods, HoK, HoO, HeK, HeO, BT, mPK, and mPB, reasonably control the empirical type I error rates for autosome analyses with continuous traits. The seven competing approaches display similar performance in terms of the empirical type I error rates for the X chromosome analyses with continuous traits (Appendix S3: Table S1).Table 1
Empirical type I errors of the seven competing methods with continuous traits.
Marker-specific weightwlNominal levelWorking correlationMethodHoK3HoOHeKHeOBTmPK4mPBUnweighted marker-specific weight10.05U/U20.048760.049600.050360.052280.049140.043520.04692E/E0.048660.049940.050160.052160.049140.01U/U0.009180.010120.010160.010300.010340.008540.01036E/E0.009240.009940.010080.010220.010280.001U/U0.000780.000820.000860.000700.000840.000840.00088E/E0.000800.000780.000840.000700.000820.0001U/U0.000080.000020.000060.000080.000080.000060.00014E/E0.000060.000020.000060.000080.00008Weighted marker-specific weight0.05U/U0.050300.049980.051580.051340.046960.046040.04536E/E0.050540.050100.051760.051220.047140.01U/U0.009920.009420.010800.009720.008880.009780.01008E/E0.009920.009440.010880.009780.008860.001U/U0.000780.000860.001260.000980.000820.001240.00134E/E0.000760.000880.001220.001020.000800.0001U/U0.000060.000080.000060.000060.000100.000020.00010E/E0.000060.000080.000080.000060.000101The unweighted marker-specific weight is given by wl=Betaml,1,1=1; the weighted marker-specific weight is given by wl=Betaml,1,25. 2U/U represents the structures of the working within-cluster and multivariate-response correlation matrices considered by the unstructured structures; E/E represents the structures of the working within-cluster and multivariate-response correlation matrices considered by the exchangeable structures. 3HoK, HoO, HeK, HeO, and BT are our proposed methods. 4mPK and mPB are executed by the R package Multi-SKAT [39].Table2 reports the empirical type I error rates based on the proposed methods, HoK, HeK, BT, HoO, and HeO, for the binary data. The two existing methods, mPK and mPB, aren’t included for comparison. This reason is that implementing the two existing methods, mPK and mPB, via the R package Multi-SKAT [39], the MPMM (multiple phenotype mixed model) function in the R package PHENIX [54–56] is a necessary tool for this process. However, the MPMM function is suitable for the continuous phenotypes [56] or is suitable for the binary phenotypes with the condition that the number of cases is sufficiently large [39]. In other words, in some sense, the two existing methods, mPK and mPB, are limited to continuous phenotypes [39].Table 2
Empirical type I errors of the five competing methods with binary traits.
Marker-specific weightwlNominal levelWorking correlationMethodHoK3HoOHeKHeOBTUnweighted marker-specific weight10.05U/U20.049440.051540.050860.052800.04952E/E0.049300.051440.050680.053180.049460.01U/U0.009740.009940.009820.010260.01000E/E0.009740.009980.009840.010280.009980.001U/U0.000680.000840.001000.000980.00106E/E0.000660.000840.001020.000940.001040.0001U/U0.000080.000020.000120.000100.00000E/E0.000080.000020.000120.000100.00002Weighted marker-specific weight0.05U/U0.051700.049000.052560.049220.04576E/E0.051680.049200.052320.049300.045560.01U/U0.010280.009760.009960.009720.00886E/E0.010240.009820.009860.009760.008840.001U/U0.001100.000800.000960.000900.00088E/E0.001120.000760.000960.000880.000900.0001U/U0.000040.000080.000100.000120.00006E/E0.000060.000080.000100.000120.000081The unweighted marker-specific weight is given by wl=Betaml,1,1=1; the weighted marker-specific weight is given by wl=Betaml,1,25. 2U/U represents the structures of the working within-cluster and multivariate-response correlation matrices considered by the unstructured structures; E/E represents the structures of the working within-cluster and multivariate-response correlation matrices considered by the exchangeable structures. 3HoK, HoO, HeK, HeO, and BT are our proposed methods.Table2 shows that the proposed methods appropriately control the type I error rates when the marker-specific weight is considered for wl=Betaml,1,1 or wl=Betaml,1,25 for variant l for binary traits. On the other hand, the empirical type I error rates of the proposed methods for X chromosome analyses with binary traits are depicted in Table S2 in Appendix S3. These empirical type I error rates show similar results as that for autosome analyses.In summary, our simulation results show that the proposed multivariate trait association methods, HoK, HoO, HeK, HeO, and BT, have reasonable control of type I error rates for continuous traits or binary traits whether the marker is X chromosomal or autosomal. On the other hand, the existing methods, mPK and mPB, yield well-controlled type I error rates for the autosome analyses or the X chromosome analyses with continuous traits (Table1 or Table S1), regardless of the weight of the marker-specific weight.


## 4.2. Empirical Power

Figure1 exhibits the comparison results of the empirical power rates for the autosome analyses with continuous traits, when the working within-cluster and multivariate-response correlation matrices of the proposed methods, HoK, HeK, and BT, are considered to be exchangeable. As expected, the empirical power rates of the seven competing methods with a weighted marker-specific weight of wl=Betaml,1,25 are higher than that with an unweighted marker-specific weight of wl=Betaml,1,1=1. The heterogeneous kernel statistic (HeK) has slightly greater empirical power rates than other methods, when the genetic effects on the different phenotypes are heterogeneous (i.e., β1≠β2), and causal SNPs have positive effects or negative effects on phenotypes. On the other hand, the existing method, mPB, has bigger empirical power rates, when the genetic effects on the different phenotypes are heterogeneous (i.e., β1≠β2), and all causal SNPs have a positive association on phenotypes. Moreover, the empirical power rates of the homogeneous omnibus test (HoO) are larger than that of the other six competing methods, when the genetic effects on the different phenotypes are homogeneous (i.e., β1=β2). Evidently, the seven competing methods have their respective advantages in identifying the association between genetic effects and multiple continuous traits for autosome analyses.Figure 1
Power comparisons of the seven competing methods with continuous traits for each scenario at the nominal level of 0.001. (a) Unweighted marker-specific weight:wl=Betaml,1,1=1. (b) Weighted marker-specific weight: wl=Betaml,1,25.Similar empirical power rates are obtained from the working within-cluster and multivariate-response correlation matrices of the proposed methods, HoK, HeK, and BT, considered to be unstructured. Hence, these empirical power rates are not shown in order to save space. On the other hand, the seven competing approaches display a similar performance in testing for the X chromosome analyses with continuous traits (AppendixS3: Figure S1).Figure2 exhibits the comparison results of empirical power rates for the autosome analyses with binary traits when the working within-cluster and multivariate-response correlation matrices of the proposed methods, HoK, HeK, and BT, are considered to be exchangeable. As a similar reason for investigating the empirical type I error rates with binary traits, the two existing methods, mPK and mPB, aren’t included for power comparison.Figure 2
Power comparisons of the five competing methods with binary traits for each scenario at the nominal level of 0.001. (a) Unweighted marker-specific weight:wl=Betaml,1,1=1. (b) Weighted marker-specific weight: wl=Betaml,1,25.Figure2 shows that the heterogeneous kernel statistic (HeK) and the heterogeneous omnibus test (HeO) outperform over other methods in terms of the empirical power rates, when the genetic effects on the different phenotypes are heterogeneous (i.e., β1≠β2). On the other hand, the empirical power rates of the homogeneous omnibus test (HoO) are bigger than that of the other competing methods, when the genetic effects on the different phenotypes are homogeneous (i.e., β1=β2). As expected, in general, the heterogeneous kernel statistic (HeK) is more powerful than the homogeneous kernel statistic (HoK), when the genetic effects on the different phenotypes are heterogeneous (i.e., β1≠β2). On the other hand, the homogeneous kernel statistic (HoK) is more powerful than the heterogeneous kernel statistic (HeK), when the genetic effects on the different phenotypes are homogeneous (i.e., β1=β2). In a word, the proposed methods, HoK, HoO, HeK, HeO, and BT, have their respective merits in examining the association between genetic effects and multiple binary traits for autosome analyses.Similarly, when the working within-cluster and multivariate-response correlation matrices of the proposed methods, HoK, HeK, and BT, are considered to be unstructured, the empirical power rates have similar results and thus they are omitted. On the other hand, the empirical power rates of the proposed methods for X chromosome analyses with binary traits are presented in FigureS2 in Appendix S3. These empirical power rates show similar results as that discussed in Figure 2.In summary, the seven competing methods, HoK, HoO, HeK, HeO, BT, mPK, and mPB, have their respective merits in diagnosing whether genetic effects are associated with multiple continuous traits for autosome analyses or the X chromosome analyses. Similarly, the proposed methods, HoK, HoO, HeK, HeO, and BT, have their respective advantages in examining whether there are associations between genetic effects and multiple binary traits for autosome analyses or the X chromosome analyses.To furthermore examine the performance of the proposed methods, additional simulation studies for continuous traits and binary traits are presented in AppendixS4 and Appendix S5 with higher correlations of phenotypes and higher dimensions of phenotypes considered, respectively. In general, these competing methods based on higher correlations of phenotypes or higher dimensions of phenotypes can provide a bigger empirical power rate for the analysis of continuous traits or binary traits. However, we note that these competing methods based on higher correlations of phenotypes or higher dimensions of phenotypes more easily have empirical type I error rate inflation at a smaller nominal level, especially for binary data analysis (Appendix S5: Tables S5-S6 and Appendix S6: Table S7), in comparison with these methods based on lower correlations of phenotypes or lower dimensions of phenotypes. A detailed discussion of these additional simulation results is given in Appendixes S4 and S5.However, we note that the proposed methods have a high computational cost, especially for binary data. Under our simulation setting and framework, we carry out a single simulated data set by using a computer based on one CPU core at 2.1 GHz. The average computational times of the homogeneous and heterogeneous tests with a weighted marker-specific weightwl=Betaml,1,25 under the alternative hypothesis for continuous data are 0.83 and 0.91 minutes, respectively, while that for binary data are 4.77 and 4.80 minutes, respectively. Therefore, in the current version, such a framework algorithm implementation is unsatisfactory for analyzing a large-scale high-dimensional data set in practice.


## 5. Conclusion

In this investigation, we develop a retrospective framework for identifying the pleiotropic effects of genetic variants on multivariate traits by using collapsing and kernel methods with pedigree- or population-structured data. The proposed framework, corresponding to the burden test, the kernel test, and the omnibus test, provides a sound basis for genetic association analyses for autosomes and the X chromosome. The proposed multivariate trait association methods based on the JGEE model can flexibly accommodate continuous phenotypes or binary phenotypes and further can adjust for covariates.One critical advantage of the proposed methods is that the homogeneous kernel statistic (HoK), the heterogeneous kernel statistic (HeK), and the burden test (BT) retain all of the benefits of the retrospective tests proposed by Schaid et al. [43] who treated the genotype data as random variables by conditioning the phenotypes as constants. On the other hand, the homogeneous omnibus test (HoO) and the heterogeneous kernel statistic (HeO) keep the advantages of the Cauchy combination tests proposed by Liu and Xie [48] who showed that the Cauchy combination tests are robust to model misspecification and robustly protect the type I error rates [49].Another important benefit of the proposed method is that the HoK test, the HeK test, and the BT test keep the benefits of the JGEE model that validly account for complex correlations between subjects within the cluster (within-cluster correlations) and between different phenotypes from the same subjects (multivariate-response correlations). Moreover, the proposed test statistics, HoK, HeK, and BT, based on the JGEE model can efficaciously account for covariate adjustment whether the phenotypes are continuous or binary.Our simulation studies show that an unweighted marker-specific weightwl=Betaml,1,1=1 and an exchangeable structure of the working within-cluster and multivariate-response correlations are recommended for the practical data analysis if the data cannot sufficiently provide valid information for estimating the structures of the working within-cluster and multivariate-response correlations before the start of the data analysis. Moreover, the homogeneous kernel statistic (HoK) is more robust than the heterogeneous statistic (HeK) in controlling the empirical type I errors, because the null distribution of the HeK statistic asymptotically follows a mixture chi-square distribution with a larger degree of freedom, in comparison with the null distribution of the HoK statistic. However, the HeK statistic is more powerful than the HoK statistic when the genetic effects on the different phenotypes are heterogeneous.On the other hand, our simulation results show that for the autosome analyses or the X chromosome analyses with continuous traits, the seven competing methods, HoK, HoO, HeK, HeO, BT, mPK, and mPB, show good performance with well-controlled type I errors, while the seven competing methods have their respective merits for identifying the association between the genetic effects and multiple continuous traits. In addition, our simulation results show that for the autosome analyses or the X chromosome analyses with binary traits, the proposed methods, HoK, HoO, HeK, HeO, and BT, can control empirical type I errors with lower correlations of phenotypes or with lower dimensions of phenotypes (Table2 and Table S2), while these proposed methods have their respective advantages for identifying the genetic variants associated with multiple binary traits. However, we observe that the proposed methods, HoK, HoO, HeK, HeO, and BT, with higher correlations of phenotypes or with higher dimensions of phenotypes, more easily have the infection of empirical type I errors at a smaller nominal level (Appendix S5: Tables S5-S6 and Appendix S6: Table S7), although these method under such situations have higher empirical power rates.


## 6. Limitation

The proposed multivariate trait association methods have their limitations. First, these proposed methods cannot simultaneously include the continuous traits and binary traits in analysis. Thus, future studies are needed to extend the idea of the proposed multivariate trait association methods for simultaneously considering continuous traits and binary traits in analysis. Second, the multivariate trait association methods, based on higher correlations of phenotypes or higher dimensions of phenotypes, easily suffer from the problem of the inflated type I errors, especially when the binary traits are considered (AppendixS5: Tables S5-S6 and Appendix S6: Table S7). Although the JGEE model provides an efficient algorithm for estimating the structure of the working within-cluster and multivariate-response correlations, a large-scale pedigree study always suffers from a more complex and high-dimensional structure of the within-cluster and multivariate-response correlations in pedigree database analysis. Hence, in the future, a more effective algorithm for estimating the complicated and high-dimensional (or higher correlational) structure of the working within-cluster and multivariate-response correlations is necessary to be proposed, especially when the analysis focuses on the binary traits. Third, in comparison with the null distribution of the homogeneous kernel statistic, the null distribution of the heterogeneous kernel statistic follows a larger degree of freedom test, which easily causes such a heterogeneous test to suffer from the problem of the type I error inflation. Therefore, overcoming the problem of the type I error inflation from the heterogeneous test is an essential part of the future work. Fourth, the proposed methods, which have a high computational cost especially for binary data, are inappropriate for analyzing large-scale high-dimensional data in practice. Thus, a more effective algorithm for reducing computational cost is needed to be proposed in further research. Moreover, the software of the proposed methods is computationally inconvenient and particularly inadequate for the mass GWAS data in practice. Therefore, the software of the proposed methods, which is convenient to be used, is a further work in the future. Fifth, our current work focuses mainly on the low- and common-frequency variants. Extension of the proposed methods to the rare variants deserves further works.



---

*Source: 8812282-2021-02-10.xml*