# A copula-based set-variant association test for bivariate continuous or
mixed phenotypes

JULIEN ST-PIERRE∗
Department of Epidemiology, Biostatistics and Occupational Health, McGill
University, Montréal, Québec, Canada
<EMAIL_ADDRESS>
KARIM OUALKACHA
Département de Mathématiques, Université du Québec à Montréal, Montréal,
Québec, Canada

###### Abstract

In genome wide association studies (GWAS), researchers are often dealing with
non-normally distributed traits or a mixture of discrete-continuous traits.
However, most of the current region-based methods rely on multivariate linear
mixed models (mvLMMs) and assume a multivariate normal distribution for the
phenotypes of interest. Hence, these methods are not applicable to disease or
non-normally distributed traits. Therefore, there is a need to develop unified
and flexible methods to study association between a set of (possibly rare)
genetic variants and non-normal multivariate phenotypes. Copulas are
multivariate distribution functions with uniform margins on the $[0,1]$
interval and they provide suitable models to deal with non-normality of errors
in multivariate association studies. We propose a novel unified and flexible
Copula-Based Multivariate Association Test (CBMAT) for discovering association
between a genetic region and a bivariate continuous or mixed phenotype. We
also derive a data-driven analytic p-value procedure of the proposed region-
based score-type test. Through simulation studies, we demonstrate that CBMAT
has well controlled type I error rates and higher power to detect associations
compared with other existing methods, for discrete and non-normally
distributed traits. At last, we apply CBMAT to detect the association between
two genes located on chromosome 11 and several lipid levels measured on 1,477
subjects from the ASLPAC study. copulas; multivariate phenotype; mixed binary-
continuous phenotypes; gene-based tests; rare-variant association; variance
component score test.

00footnotetext: To whom correspondence should be addressed.

##  1 Introduction

Among the thousands of genetic variants that have been identified by genome
wide association studies (GWAS), $4.6\%$ of single-nucleotide polymorphisms
(SNPs) and $16.9\%$ of genes have proved to be associated with multiple
correlated phenotypes (Solovieff _and others_ , 2013). This has led to a
growing interest in discovering novel biology of pleiotropy by jointly
analyzing multiple traits simultaneously. Pleiotropy occurs when several
phenotypes are co-regulated by a common gene (Stearns, 2010; Williams, 1957).
Some phenotypes are intrinsically multivariate, such as blood pressure, which
is measured by both systolic and diastolic pressure. Moreover, analysis of
multiple neurobiological and/or psychological genetically influenced
quantifiable traits, i.e. endophenotypes (Gottesman and Gould, 2003), has been
proposed to facilitate the search of susceptibility genes related to
heterogeneous/complex traits, such as schizophrenia (SZ) and major depressive
disorder (MDD) (Hall, 2010; Iacono _and others_ , 2017). Indeed, multiple
endophenotypes often measure the same underlying trait and might bear a more
direct relationship with the biological etiology of a clinical disorder.
Ignoring the fact that multiple phenotypes are related and might share a
common genetic basis of complex traits, the single-trait study paradigm is
likely to suffer from a potential loss of power and may provide less accurate
prediction of some phenotypes compared to simultaneous analysis of multiple
traits (Yang _and others_ , 2016).

Region-based association methods have become increasingly popular given large-
scale genomic data provided by next generation sequencing technology and are
often of interest in genome wide association studies (GWAS). Such methods
aggregate variants into groups based on pre-defined criteria, such as gene
boundaries or linkage-disequilibrium (LD) blocks, and then test for phenotypic
effects within each group. Region-based methods reduce the multiple testing
burden of single-variant analyses. Furthermore, they can increase power by
combining information across biologically meaningful units (Petersen _and
others_ , 2013). Several approaches have been developed to handle multivariate
association in both single-variant and region-based frameworks. For
association methods between a single genetic variant and a multivariate
phenotype, a nice overview can be found in Yang and Wang (2012) (see also
Schaid _and others_ , 2016; Shriner, 2012; Solovieff _and others_ , 2013;
Zhang _and others_ , 2014). In the region-based framework, many strategies
have been explored for testing multivariate association and can be classified
into three main categories: regression analysis, variable reduction analysis,
and combining approach.

The combining approach consists in aggregating test statistics from either
multiple single-variant-based tests (Van der Sluis _and others_ , 2015) or
test statistics obtained from single-variant-based multivariate tests (Guo
_and others_ , 2013) to obtain one overall, multivariate gene-based p-value.
This approach can lead to low power because it discards the additional
information obtained by combining many outcomes into a unified analysis.
Indeed, modeling multivariate phenotypes may increase the power over analyzing
individual phenotypes separately (Klei _and others_ , 2008). The second
category includes variable reduction methods such as principal component
analysis (PCA) and canonical correlation analysis (CCA). CCA, which extracts
the linear combination of traits that explain the largest possible amount of
covariance between two sets of variables, was initially proposed for
multivariate association analysis by Ferreira and Purcell (2009) and later
extended to perform multivariate gene-based association tests (Tang and
Ferreira, 2012; Seoane _and others_ , 2014). Important limitations of CCA is
that these methods cannot accomodate covariates and are not appropriate for
the analysis of non-normally distributed traits. To this date, work is still
needed to extend the PCA framework for joint analysis of multiple phenotypes
in GWAS to multiple phenotype analysis in region-based association studies
(Liu and Lin, 2019).

Finally, the last category regroups methods based on regression analysis. The
use of functional linear models for the analysis of multiple variants and
multiple traits has been proposed by Wang _and others_ (2015). However, most
of the regression methods are based on multivariate linear mixed models
(mvLMMs) which assume multivariate normal distribution for multiple phenotypes
and a random effect to account for the cross-phenotype correlations (Sun _and
others_ , 2016; Dutta _and others_ , 2019; Lee _and others_ , 2017). A
limitation of all these methods is that they are not applicable to disease or
non-normally distributed traits. Lee _and others_ (2017) extended their
Multivariate Association Analysis using Score Statistics (MAAUSS) method so
that it could also handle correlated binary and continuous phenotypes using
the generalized estimating equation (GEE) framework. In practice, researchers
are often dealing with both non-normally distributed traits and a mixture of
discrete-continuous traits, hence there is a need to develop unified and
flexible methods to study association between a set of rare/common variants
and non-normal multivariate phenotypes.

Copulas are multivariate distribution functions with uniform margins on the
$[0,1]$ interval (Nelsen, 2010). They provide flexible models for joint
distributions by separating specification of the marginal distributions and
the dependence. In other words, they allow for modelling a dependency
structure between several random variables regardless of their marginal
distributions. Hence, they provide suitable models to deal with non normality
of errors in multivariate association studies. Copulas are extensively used in
regression models due to their great flexibility, as they can be incorporated
in models with binary outcomes (Genest _and others_ , 2013) as well as with
mixed discrete and continuous outcomes (de Leon and Wu, 2011). More recently,
they have been incorporated into a single-marker association test to increase
the power of genetic association studies of rare variants (Konigorski _and
others_ , 2020) and into a gene-based association study for bivariate survival
traits based on functional regression (Wei _and others_ , 2019). Of note, the
latter model is only valid for continous survival traits and the authors have
assumed the same effects for the covariates and the genetic region on both
marginal survival times, which is unrealistic in practice when dealing with
different diseases for instance.

In this paper, we propose a novel unified and flexible copula-based
multivariate test, CBMAT (Copula Based Multivariate Association Test), for
discovering association between a genetic region and a bivariate continuous or
mixed phenotype. By using copulas to model the joint distribution of
phenotypes, CBMAT allows a wide range of continuous and discrete marginal
distributions for the phenotypes. Many of the existing multivariate methods
for region-based tests can be viewed as special cases of CBMAT with a Gaussian
copula and particular choices of kernels to relate genetic variants to
multiple phenotypes. Finally, we also derive a data-driven analytic p-value
procedure of the proposed region-based score-type test.

The article is organized as follows: in Section 2, we introduce our proposed
method CBMAT and derive a data-driven analytic p-value procedure of the
proposed region-based score-type test. In Section 3 and Section 4
respectively, we describe our simulation studies and demonstrate that CBMAT
has well controlled type I error rates and higher power to detect associations
compared with other existing methods, for discrete and non-normally
distributed traits. At last, we apply CBMAT to detect the association between
two genes located on chromosome 11 and several lipid levels measured on 1,477
subjects from the ASLPAC study.

##  2 Material and Methods

Our goal in this section is to develop a unified and flexible multivariate
mixed-effects-model for a wide range of discrete and continuous distributions.
First, we present the copula-based model for the joint density of
$(Y_{1},Y_{2})$ when both traits are continuous and when one trait is discrete
binary. Secondly, we briefly review some existing multivariate linear mixed
models for region-based association. Third, we present our proposed CBMAT
score statistics for testing the variance component of the region effect.
Finally, we derive the null distribution of our proposed score statistics and
its data-driven analytic p-value procedure.

###  2.1 Copula-based model for joint density

Assume we observe 2 traits for $N$ independent subjects,
$\bm{Y}_{i}=(y_{i1},y_{i2})^{T}$, $i=1,...,N$, $m$ covariates,
$\bm{X_{i}}=(x_{i1},x_{i2},...,x_{im})^{T}$, and a set of SNPs,
$\bm{G_{i}}=(g_{i1},g_{i2},...,g_{ir})^{T}$ that contains $r$ biallelics
variants coded $0,1,2$, representing the number of minor alleles. The proposed
model assumes that the (conditional) marginal distributions of $Y_{ik}$,
$k=1,2$, are belonging to the exponential family, and it links $\mu_{ik}$, the
conditional mean of $Y_{ik}$, to $\bm{X}_{i}$ and $\bm{G}_{i}$ via a
(conditional) marginal generalized model, as follows

$\displaystyle\text{E}[Y_{ik}|\bm{X}_{i},\bm{G}_{i}]=\mu_{ik}=g_{k}^{-1}\left(\bm{X}_{i}^{T}\bm{\gamma}_{k}+\bm{G}_{i}^{T}\bm{\beta}_{k}\right),$
(1)
$\displaystyle\text{Var}[Y_{ik}|\bm{X}_{i},\bm{G}_{i}]=\phi_{k}\cdot\nu_{k}(\mu_{ik}),$

where, for the $k^{th}$ trait,
$\bm{\gamma}k=(\gamma_{k1},...,\gamma_{km})^{T}$ and
$\bm{\beta}k=(\beta_{k1},...,\beta_{kr})^{T}$ are vectors of regression
coefficients for $\bm{X}i$ and $\bm{G}i$ respectively, $g_{k}(.)$ is the link
function, $\phi_{k}$ a dispersion parameter and $\nu_{k}$ the variance
function. Thus, the traits (conditional) marginal cumulative distribution
functions, $F_{k}(y_{ik}|\bm{X}_{i},\bm{G}_{i})$, are completely specified. We
then specify the joint CDF of $Y_{i1}$ and $Y_{i2}$ as follows

$\bm{F}\left(y_{i1},y_{i2}\mid\bm{X}_{i},\bm{G}_{i}\right)=C_{\theta}\left\\{F_{1}\left(y_{i1}\mid\bm{X}_{i},\bm{G}_{i}\right),F_{2}\left(y_{i2}\mid\bm{X}_{i},\bm{G}_{i}\right)\right\\},$
(2)

where $C_{{\theta}}$ is a copula function, which accounts for dependence
between marginal CDFs $F_{1},F_{2}$ through the parameter ${\theta}$; i.e.,
$C_{{\theta}}$ is a joint distribution function in the unit square

$\displaystyle C_{\bm{\theta}}:[0,1]^{2}$ $\displaystyle\longmapsto$
$\displaystyle[0,1]$ $\displaystyle(u_{1},u_{2})$ $\displaystyle\longmapsto$
$\displaystyle C_{\bm{\theta}}(u_{1},u_{2}).$

In the sequel, we drop the conditional notation with respect to $\bm{X}_{i}$
and $\bm{G}_{i}$ for sake of clarity. For two continuous phenotypes, the
(conditional) joint density function is given as

$\displaystyle
f(y_{i1},y_{i2})=f_{1}(y_{i1})f_{2}(y_{i2})c_{\theta}(F_{1}(y_{i1}),F_{2}(y_{i2})),$
(3)

where $c_{\theta}(u,v)=\partial^{2}C_{\theta}(u,v)/\partial u\partial v$
stands for the copula density. In the mixed phenotype case, that is for
$Y_{i1}$ binary and $Y_{i2}$ continuous, a liability threshold model which
transforms the binary trait into a probit latent marginal variable is
proposed. We define $Y_{i1}^{*}$ such that it has a normal CDF
$F_{1}^{*}(y_{i1}^{*})=\Phi\left(\frac{y_{i1}^{*}-g_{1}(\mu_{i1})}{\sigma_{1}}\right)$
(i.e., $Y_{i1}^{*}\sim N(g_{1}(\mu_{i1}),\sigma_{1}^{2})$), with $Y_{i1}$ and
$Y_{i1}^{*}$ linked in a way such that $Y_{i1}=I(Y_{i1}^{*}>y^{c})$, where
$I(\cdot)$ is the indicator function and $y^{c}$ is a cutpoint. For
convenience and identifiability reasons, we take $y^{c}=0$ and
$\sigma_{1}^{2}=1$. Letting $\bm{F}^{*}$ be the joint CDF of $Y_{i1}^{*}$ and
$Y_{i2}$, we define the following probability of event for $Y_{i1}$ and
$Y_{i2}$ as follows

$\displaystyle
P(Y_{i1}=y_{1},Y_{i2}<y_{2})=\begin{cases}\bm{F}^{*}(0,y_{2})&\text{if
}y_{1}=0,\\\ F_{2}(y_{2})-\bm{F}^{*}(0,y_{2})&\text{if }y_{1}=1.\end{cases}$

Let now $\bm{F}^{*}$ be determined by a copula model, as in (2). This leads to

$\bm{F}^{*}(y_{i1}^{*},y_{i2})=C_{\theta}\big{\\{}F_{1}^{*}(y_{i1}^{*}),F_{2}(y_{i2})\\}.$

Thus, one can write the joint (conditional) density
$f\big{(}y_{i1},y_{i2}\big{)}=\partial P(Y_{i1}=y_{i1},Y_{i2}<y_{i2})/\partial
y_{2}$ as

$f\big{(}y_{i1},y_{i2}\big{)}=f_{2}(y_{i2})\times\Big{[}C_{\theta}^{01}\\{1-\mu_{i1},F_{2}(y_{i2})\\}\Big{]}^{I(y_{i1}=0)}\times\Big{[}1-C_{\theta}^{01}\\{1-\mu_{i1},F_{2}(y_{i2})\\}\Big{]}^{I(y_{i1}=1)},$
(4)

where $\mu_{i1}$ is given in (LABEL:marginalmean), $f_{2}$ is the density
function of $Y_{2}$, and $C_{\theta}^{01}$ is the derivative of the copula
function $C_{\theta}$ over its second argument, $C_{\theta}^{01}(u,v)=\partial
C_{\theta}(u,v)/\partial v$.

Recall that the regression parameters linking both covariates and genotypes to
the phenotypes are kept in the marginal CDFs, and so, they are marginally
meaningful. Moreover, both between- and within-trait polygenic heritability
(i.e. polygenic dependence, (Bauman _and others_ , 2005)) are “margin-free”,
in a way that they are characterized by the copula alone; this is another
advantage of the proposed model.

###  2.2 Existing mvLMM approaches

In the multivariate set-based association framework, the overall effect of the
genomic region on the two phenotypes,
$\bm{\beta}=(\bm{\beta}_{1}^{\top},\bm{\beta}_{2}^{\top})^{\top}$, is assumed
to be a vector of (normally-distributed) random effects with mean zero and
variance-covariance matrix $\eta\Sigma_{P}\otimes\Sigma_{G}$, where $\otimes$
is the Kronecker product, $\eta$ is an unknown variance component,
$\Sigma_{P}$ is a $2\times 2$ matrix that captures the relationship among the
effect sizes of each variant on the $2$ phenotypes, and $\Sigma_{G}$ is a
$r\times r$ matrix that models the relationship among the effect sizes of the
variants on each phenotype.

In our proposed model, we also adopt the random effect assumption of the
region overall-effect on the two phenotypes. Consequently, if one assumes
marginal normal distributions for the two phenotypes and set $C_{{\theta}}$ to
be a Gaussian copula, then the proposed model (2) will be equivalent to the
multivariate mixed-effects-model of MURAT (Sun _and others_ , 2016) and Multi-
SKAT (Dutta _and others_ , 2019), which assume multivariate normality of the
error terms, and thereby Gaussian distributions for the univariate margins.
Thus, in comparison with existing methods, the proposed model allows for
marginal CDFs to be any distribution belonging to the exponential family and
can handle both continuous and discrete phenotypes simultaneously.

Of note, the main difference between most existing mvLMM approaches results
from different specifications of the two kernel matrices, $\Sigma_{G}$ and
$\Sigma_{P}$. For instance, MURAT assumes a common correlation, $\rho$, for
the effects of the same variant on different phenotypes and also that the
effects of different variants are uncorrelated; i.e. MURAT assumes that
${\Sigma_{P}=(1-\rho)I_{K}+\rho\bm{1_{K}}\bm{1_{K}}^{T}}$, where $\bm{1_{K}}$
is a $K$-dimensional vector of ones. The Multi-SKAT method (Dutta _and others_
, 2019) suggests different choices for $\Sigma_{P}$. For instance, if we
assume that the effect sizes of a variant $j$ are the same for all the
phenotypes, in which case $\beta_{1j}=\ldots=\beta_{Kj}$, the homogeneous
kernel is $\Sigma_{P}=\bm{1_{K}}\bm{1_{K}}^{T}$ (i.e. $\rho=1$ in
$\Sigma_{P}$). On the other hand, if we assume that the effect sizes
$(\beta_{1j},...,\beta_{Kj})$ are uncorrelated among themselves, the
heterogeneous kernel is $\Sigma_{P}=I_{K}$ (i.e. $\rho=0$ in $\Sigma_{P}$).
Other suggestions for $\Sigma_{P}$ are given in detail in Dutta _and others_
(2019).

The choice of $\Sigma_{G}$ has been extensively studied in the literature (Wu
_and others_ , 2011; Lee _and others_ , 2012) and a common choice is to assume
$\Sigma_{G}=\bm{W}$, where $\bm{W}=\text{diag}(w_{1},...,w_{r})$ is an
$r\times r$ diagonal matrix of the weights to be used for the r variants.
Typically, the weights are chosen such that they are inversely proportional to
the minor allele frequency (MAF) to up-weight rare variants, for instance by
using $w_{j}=beta(MAF_{j},1,25)$.

###  2.3 Hypothesis testing and parameter estimation

Under the copula-based model (2), the association between a region that
containts $r$ variants and the bivariate phenotype can be tested by evaluating
the null hypothesis $H_{0}:\bm{\beta}=\bm{0}$. As mentionned in the previous
section, we assume $\bm{\beta}$ is a random effects vector that follows an
arbitrary probability distribution

$\displaystyle\bm{\beta}=\left(\bm{\beta}_{1}^{T},\bm{\beta}_{2}^{T}\right)^{T}\sim\bm{H}\left(\bm{0}_{2r},\bm{\Sigma}_{\beta}=\eta\Sigma_{P}\otimes\Sigma_{G}\right).$
(5)

In this context, testing $\bm{\beta}=\bm{0}$ corresponds to testing $\eta=0$.
In what follows, we derive a score-type test statistic for the null hypothesis
$H_{0}:\eta=0$. First, one can see that the conditional likelihood function
given $\bm{\beta}$ is given as

$\displaystyle L(\bm{\xi}|\bm{\beta})$
$\displaystyle=\prod_{i=1}^{n}f(y_{i1},y_{i2}),$ (6)

where $f(y_{i1},y_{i2})$ is given by (3) when both traits are continuous, and
by (4) for the mixed trait case,
$\bm{\xi}=(\theta,\bm{\gamma}1,\bm{\gamma}2,\phi_{1},\phi_{2})^{T}$, where
$\theta$ is the copula parameter, while $\bm{\gamma}1,\bm{\gamma}2,\phi_{1}$,
and $\phi_{2}$ are the parameters of the marginal GLMs given in
(LABEL:marginalmean). Of note, the product factorization of the conditional
likelihood in (6) results from the fact that, conditioning on $\bm{\beta}$,
the subjects are independent. The full likelihood is given by the
$2r$-dimensional integral

$\displaystyle L(\eta,\bm{\xi})=\int_{\bm{\beta}}L(\bm{\xi}|\bm{\beta})\
\bm{H}(\bm{\beta})d\bm{\beta}.$ (7)

Finally, the score test for $H_{0}$: $\eta=0$ is based on the score statistic

$\displaystyle U(\eta)=\frac{\partial}{\partial\eta}\text{log
}L(\eta,\bm{\xi}).$

The problem with the direct computation of this score is the evaluation of the
integral in dimension $\mathbbm{R}^{2r}$ in (7). To solve this computational
problem, we rely on the work proposed by (Lin, 1997) and approximate the
integral using Taylor’s expansion techniques of $L(\bm{\xi}|\bm{\beta})$ in
the neighborhood of $\bm{\beta}=\bm{0}_{2r}$. In Appendix A, we show that

$\displaystyle
U(\eta)\approx\frac{1}{2}\Big{(}\bm{L}^{T}(\Sigma_{P}\otimes\bm{G}\Sigma_{G}\bm{G}^{T})\bm{L}+\text{tr}\left\\{(\Sigma_{P}\otimes\bm{G}\Sigma_{G}\bm{G}^{T})\bm{D}\right\\}\Big{)},$
(8)

where $\bm{L}=(\bm{L}_{1}^{T},\bm{L}_{2}^{T})^{T}$ is a $2n\times 1$ vector,
with $\bm{L}_{1}=(\bm{L}_{11},\ldots,\bm{L}_{1n})^{T}$ and
$\bm{L}_{2}=(\bm{L}_{21},\ldots,\bm{L}_{2n})^{T}$, such that

$\frac{\partial}{\partial\bm{\beta}}l(\bm{\xi}|\bm{\beta})=\sum_{i=1}^{n}\begin{bmatrix}\bm{L}_{1i}\\\\[10.00002pt]
\bm{L}_{2i}\end{bmatrix}\otimes\bm{G}_{i}=\left(I_{2}\otimes\bm{G}^{T}\right)\bm{L},$

and $\bm{G}$ is the $n\times r$ matrix of genotypes, $\bm{D}$ is a block
matrix defined in Appendix A, whose entries are corresponding to the Hessian
of the (conditional) log-likelihood with respect to $g(\bm{\mu})$, such that

$\frac{\partial^{2}}{\partial\bm{\beta}\partial\bm{\beta}^{T}}l(\bm{\xi|\beta})=\left(I_{2}\otimes\bm{G}^{T}\right)\bm{D}\left(I_{2}\otimes\bm{G}\right).$
(9)

The score statistic in (8) depends on the data and the unknown vector
parameter $\bm{\xi}$ through $\bm{L}$ and $\bm{D}$. It depends also on the two
kernel matrices $\Sigma_{P}$ and $\Sigma_{G}$, which need to be specified
based on prior knowledge beforehand. Thus, our plug-in score test statistic is
obtained by substituting the unknown entries of $\bm{\xi}$ in (8) by their
maximum likelihood estimators under the null hypothesis $H_{0}:\eta=0$.

###  2.4 Null distribution correction of the proposed plug-in score statistic

Under the null model, when the unknown parameters are fixed at their true
values, standard asymptotic theory gives that
$\left(I_{2}\otimes\bm{G}^{T}\right)\bm{L}$ follows an $2r$-variate normal
distribution with mean zero and variance-covariance matrix
$-E\left(\left(I_{2}\otimes\bm{G}^{T}\right)\bm{D}\left(I_{2}\otimes\bm{G}\right)\right)$.
Thus, the asymptotic distribution of $U(\eta)$ is equivalent to that of
$\sum_{i=1}^{2r}\lambda_{i}(\chi^{2}_{1}-1)/2$, where $\lambda_{i}$,
$i=1,...,2r$, are eigenvalues of
$\bm{B}^{1/2}(\Sigma_{P}\otimes\Sigma_{G})\bm{B}^{1/2}$, with
$\bm{B}=-\left(I_{2}\otimes\bm{G}^{T}\right)\bm{D}\left(I_{2}\otimes\bm{G}\right)$
is the negative of the Hessian matrix defined in (9). However, ignoring
variability due to the fact that we are using estimates of $\bm{\xi}$, in the
asymptotic variance-covariance matrix, $\bm{B}$, can lead to severe type I
error inflation of the proposed score test. In our work, we rely on (Lu and
Shiou, 2002) and suggest a corrected variance-covariance matrix, say
$\tilde{\bm{B}}$, of the score vector, which accounts for the variability
induced by the plug-in estimates. This corrected matrix satisfies
$\tilde{\bm{B}}^{-1}=\bm{B}^{-1}+\bm{C}_{\hat{\xi}}$, where
$\bm{C}_{\hat{\xi}}$ is a $2r\times 2r$ positive definite symmetric matrix,
which depends on $\hat{\bm{\xi}}$ and the data. Consequently, the corrected
null distribution of the score statistic becomes
$\sum_{i=1}^{2r}\tilde{\lambda}_{i}\chi^{2}_{1}/2-\sum_{i=1}^{2r}\lambda_{i}/2$,
where $\tilde{\lambda}_{i}$, $i=1,...,2r$, are eigenvalues of
$\tilde{\bm{B}}^{1/2}(\Sigma_{P}\otimes\Sigma_{G})\tilde{\bm{B}}^{1/2}$.
Finally, one can approximate the score test distribution under the null using
the Davies method (Davies (1980)), and obtain the corresponding p-value.
Appendix B gives more details about the derivation of the asymptotic
distribution of the proposed score statistic, and outlines steps to obtain the
corrected asymptotic variance-covariance matrix $\tilde{\bm{B}}$.

In this work, we set $\Sigma_{G}=\bm{W}$, an $r\times r$ diagonal matrix of
weights. Several choices of $\Sigma_{P}$ and their effect from the modelling
perspective have been investigated in Dutta _and others_ (2019). If a selected
$\Sigma_{P}$ does not reflect biology, the test may have substantially reduced
power. To avoid such a drawback and achieve robust power, we follow the same
procedure as in MURAT, and set
$\Sigma_{\rho}:=\Sigma_{P}=(1-\rho)I_{K}+\rho\bm{1_{K}}\bm{1_{K}}^{T}$, then
we aggregate results across different $\Sigma_{\rho}$ by selecting, from a
grid of values over $[0,1]$, the optimal value of $\rho$ which maximizes the
power of the test.

###  2.5 Data-driven analytic p-value procedure

We outline here our procedure for obtaining an analytic p-value of the
proposed data-driven score test. To do so, let
$(\rho_{1},\rho_{2},\ldots,\rho_{b})$ be a grid of values of $\rho$ in
$[0,1]$, and let $U_{\rho_{j}}(\eta)$ the corresponding plug-in estimate score
statistic. Deriving p-value for the minimum p-value statistic, based on
$U_{\rho_{j}}(\eta)$’s, can be achieved using permutation or perturbation to
calculate the Monte‐Carlo p-value, which is computationally expensive.
Resampling-based techniques, which require less resampling steps, as in Multi-
SKAT, can also be used to approximate the joint null distribution of
$\left(U_{\rho_{1}}(\eta),\ldots,U_{\rho_{b}}(\eta)\right)$, and then one can
approximate p-value for the minimum p-value statistic. Here, again, we rely on
copulas to develop an analytic p-value procedure, which does not require
permutation nor resampling techniques. In fact, noticing that the second term
in the right-hand side of (8) has small variability compared to the first
term, we consider a series of statistics

$\displaystyle Q_{\rho_{j}}$
$\displaystyle=\bm{L}^{T}(\Sigma_{\rho_{j}}\otimes\bm{G}\Sigma_{G}\bm{G}^{T})\bm{L}$
$\displaystyle=\bm{Z}^{T}\tilde{\bm{K}}_{\rho_{j}}\bm{Z},$ (10)

where
$\bm{Z}=\tilde{\bm{B}}^{-1/2}\left(I_{2}\otimes\bm{G}^{T}\right)\bm{L}\sim
N(\bm{0},\bm{I}_{2r})$, and
$\tilde{\bm{K}}_{\rho_{j}}=\tilde{\bm{B}}^{1/2}(\Sigma_{\rho_{j}}\otimes\Sigma_{G})\tilde{\bm{B}}^{1/2}$
for $j=1,\ldots,b$. Again, one can verify that the null distribution of
$Q_{\rho_{j}}$ is $\sum_{i=1}^{2r}\tilde{\lambda}_{i}^{(j)}\chi^{2}_{1}$,
where $\tilde{\lambda}_{1}^{(j)},\ldots,\tilde{\lambda}_{2r}^{(j)}$ are
eigenvalues of $\tilde{\bm{K}}_{\rho_{j}}$. The p-value of such a test
statistic is $p_{j}=P(Q_{\rho_{j}}>q_{\rho_{j}})=S_{j}(q_{\rho_{j}})$, where
$S_{j}(.)$ and $q_{\rho_{j}}$ are the survival and the observed value of
$Q_{\rho_{j}},j=1,\ldots,b$. Now, assume
$\mathcal{P}_{min}=min(P_{1},\ldots,P_{b})$, with $P_{j}=S_{j}(Q_{\rho_{j}})$.
Our data-driven procedure rejects the null hypothesis if $\mathcal{P}_{min}$
is large, and its corresponding p-value is $P(\mathcal{P}_{min}<p_{min})$,
where $p_{min}$ is the observed value of $\mathcal{P}_{min}$, which is given
as follows

$\displaystyle P\left(\mathcal{P}_{\min}<p_{\min}\right)$
$\displaystyle=1-P\left\\{P_{1}>p_{\min},\ldots,P_{b}>p_{\min}\right\\}$
$\displaystyle=1-P\left\\{P_{1}<1-p_{\min},\ldots,P_{b}<1-p_{\min}\right\\}$
$\displaystyle=1-P\left\\{Q_{\rho_{1}}>S_{1}^{-1}\left(1-p_{\min}\right),\ldots,Q_{\rho_{b}}>S_{b}^{-1}\left(1-p_{\min}\right)\right\\}.$
(11)

In order to derive a closed-form expression for the p-value in (11), we
approximate the unknown joint null distribution of
$(Q_{\rho_{1}},\ldots,Q_{\rho_{b}})$ by copulas. More precisely, we assume

$P\left(Q_{\rho_{1}}>q_{\rho_{1}},\ldots,Q_{\rho_{b}}>q_{\rho_{b}}\right)=C_{\Gamma}\left\\{S_{1}\left(q_{\rho_{1}}\right),\ldots,S_{b}\left(q_{\rho_{b}}\right)\right\\},$

where $C_{\Gamma}$ is a Gaussian copula indexed by a correlation matrix
$\Gamma$. Under such an assumption, we have the following closed-form
expression for the p-value

$\displaystyle
P\left(\mathcal{P}_{\min}<p_{\min}\right)=1-C_{\Gamma}\left(1-p_{\min},\ldots,1-p_{\min}\right).$

The p-value calculation procedure is completed by specifying the correlation
matrix $\Gamma$. In this work, we take advantage of the fact that, for $1\leq
j,l\leq b$, we have (Magnus, 1978)

$Cov(Q_{\rho_{j}},Q_{\rho_{l}})=2\text{tr}(\tilde{\bm{K}}_{\rho_{j}}\tilde{\bm{K}}_{\rho_{l}}),$
(12)

and we set $\Gamma$ to be the Pearson correlation matrix of
$(Q_{\rho_{1}},\ldots,Q_{\rho_{b}})$.

Of note, Kendall’s tau correlations might be more suitable as entries of
$\Gamma$, however analytical derivation of such correlations from (12)
requires solving $b(b-1)/2$ double-integrals equations (Sun _and others_ ,
2019). In our simulation studies in Section 3, we show that the use of Pearson
correlation matrix to model dependence between the test statistics has no
impact on type I error. Moreover, in Appendix D of the Supplementary
Materials, we compare CBMAT p-values obtained using these two strategies, that
is, either an approximation of the matrix $\Gamma$ using the rank-based
dependence between the test statistics (i.e. Kendall’s tau coefficients)
calculated based on a resampling procedure, or using the Pearson correlation
matrix, as described above. We show that these two procedures yield similar
power for both continuous and mixed bivariate phenotypes.

##  3 Simulation Studies

We carried out simulation studies to evaluate and compare the proposed
methodology for region-based association testing of bivariate phenotypes.
Comparisons were made with MURAT and Multi-SKAT, using available packages in
R. For Multi-SKAT, we assumed that effect sizes of a variant on different
phenotypes were either homogeneous (Multi-SKAT Hom) or uncorrelated (Multi-
SKAT Het).

Data generation: Genotype data of subjects with European ancestry from the
1000 Genomes Project (The 1000 Genomes Project Consortium, 2015) were used to
mimic real data situations. From this public database, we extracted genetic
variants within 500kbs of the BRCA1 gene, for which several known mutations
are associated with a higher risk of developing breast, ovarian and prostate
cancers (Foulkes and Shuen, 2013). To avoid multicollinearity issues, we
pruned variants in strong LD ($r^{2}>0.7$). Simulated data of $N=503$ subjects
were generated according to scenarios defined by varying combinations of
values for parameters (Kendall’s $\tau$, $\rho$, $\eta$), for a total of 21
scenarios: 3 × 1 × 1 = 3 scenarios under $H_{0}$, and 3 × 3 × 2 = 18 scenarios
under the alternative hypothesis, $H_{a}$. For a given scenario and for each
individual, $i$, independently, we simulated the data according to the
following steps:
1\. $\mathbf{G}_{i}$ is the genotype vector of a randomly selected genetic
region of $r=30$ consecutive common/rare variants from the extracted large
region around BRCA1 gene (the selected region is the same for all subjects);
2\. $\mathbf{X}_{i}=(1,X_{i1},X_{i2})^{T}$, with
$X_{i1}\sim\textit{Bernoulli}(0.5)$ and $X_{i2}\sim\textit{N}(0,1)$;
3\. $\bm{\beta}=\begin{pmatrix}\bm{\beta}1\\\ \bm{\beta}2\end{pmatrix}\sim
N_{2r}\left(\bm{0}_{2r},\eta\begin{bmatrix}\bm{W}&\rho\bm{W}\\\
\rho\bm{W}&\bm{W}\end{bmatrix}\right)$;
4\. $(U_{i},V_{i})\sim C_{\theta}$
5\. $Y_{i1}=F^{-1}_{1}(U_{i}|\mu_{i1},\phi_{1},\tau)$, and
$Y_{i2}=F^{-1}_{2}(V_{i}|\mu_{i2},\phi_{2},\tau),$
with
$\mu_{ik}=g_{k}^{-1}(\bm{X}_{i}^{T}\bm{\gamma}_{k}+\bm{G}_{i}^{T}\bm{\beta}_{k}),$
where $F_{k}^{-1}(\cdot)$ and $g_{k}^{-1}(\cdot)$ are respectively the inverse
of the CDF function and the inverse link function of $Y_{k}$, $k=1,2$.

We fixed the covariates’ effects at $\bm{\gamma}_{1}=(-0.20,0.33,0.78)^{T}$
and $\bm{\gamma}_{2}=(1.52,1.25,1.86)^{T}$, and we set the dispersion
parameters to $\phi_{1}=\phi_{2}=1$. We fixed the weights matrix
$\bm{W}=\bm{I}_{r}$. The dependency structure between $Y_{1}$ and $Y_{2}$ was
induced using a selected Kendall’s $\tau$ value and either a Gaussian, a
Frank, or a Clayton copula. Our selected $\tau$ values ($0.05,0.2,0.4$)
correspond to weak, moderate, and strong phenotype dependency, respectively.
Under $H_{0}$, we set $\eta=0$. Under $H_{a}$, the value of $\eta$ depends on
the fraction of causal variants $v$ and on the traits heritability $h^{2}$,
which is defined as the fraction of the total phenotypic variability that is
attributable to genetic variability. We set $v$ equals to $10\%$ and $20\%$
respectively, and $h^{2}$ equals to 2$\%$. Of note, as the binary trait was
generated based on a latent probit model, heritability was calculated for its
underlying latent continuous trait. As there exists no closed form
relationship between $\eta$, $v$ and $h^{2}$ for non-linear models,
calculations were made based on the parameter relationships derived for the
normal traits. More details are given in Appendix C of the Supplementary
Materials. Under $H_{a}$, we varied the region-specific pleiotropy parameter
as $\rho=(0,0.4,0.8)$, which correspond to no, moderate, and strong pleiotropy
effect, respectively. We refer the reader to Table 1 for the specific values
of $v$, $\rho$ and Kendall’s $\tau$ considered.

We assessed the performance of the CBMAT and its competitors with respect to
the parameter of interest $\eta$ (i.e. $\bm{\beta}$) in terms of the type I
error and power. Under the null hypothesis, 10 000 random samples were
generated as described above in all simulation scenarios. Under the
alternative hypothesis, 5000 random samples were generated in all the
simulation scenarios, to evaluate the power of the methods. All type I error
and power results were calculated with significance threshold $\alpha=1\%$.
For CBMAT, we used the package optim in R to estimate the copula parameter as
well as the parameters of the marginal distributions, under the null
hypothesis of no genetic association.

Analysis of simulated data: Three simulation settings were investigated to
evaluate the methods performance. Setting 1 aims to evaluate the robustness of
all methods against the non-normality assumption of the marginal
distributions, Setting 2 illustrates the performance of CBMAT in the case of
the mixed-trait case, and Setting 3 focuses on the impact of misspecifying the
true copula model for the CBMAT approach.

Setting 1: In this setting, both traits were assumed to be continuous, and we
considered two cases for the marginal distributions: 1) $(F_{1},F_{2})$ are
set to be (Exponential, Exponential) CDFs with the log-link function for the
covariates, 2) $F_{1}$ is the Student-t distribution with $\nu=3$ degrees of
freedom and we used the identity-link function for the covariates, and $F_{2}$
is set to be an Exponential distribution (log-link). The dependence between
the traits is generated from either a Normal, a Frank, or a Clayton copula.
For CBMAT, we evaluated its performance with the correct copula, and with both
the correct marginal distributions of the phenotypes (true), and when the
margins are assumed to be unknown. In the latter case, selection of the
marginal models was based on the best fitted GLM models in
(LABEL:marginalmean) with the lowest AIC, among the Gaussian, Gamma and
Student-t distributions. For the other methods, as these approaches
theoretically require the traits normality assumption, we compared their
performance when an inverse normal transformation (INT) was applied over the
residuals of each trait after adjusting for covariates (Beasley _and others_ ,
2009), and without normalizing the traits.

Setting 2: For the mixed binary-continuous bivariate phenotype, $Y_{1}$ is
assumed to follow a probit latent distribution and $Y_{2}$ follows an
Exponential distribution (log-link). In this setting, we assessed the
performance of CBMAT when the fitted marginal models are the true ones, and
when the marginal distribution of $Y_{2}$ is considered unknown, but it is
selected based on the AIC criterion, in a similar way as in Setting 1.
Although both MURAT and Multi-SKAT are designed to handle only continuous
phenotypes, we fitted both models in order to evaluate their performance in
the mixed-trait case. For both methods, we analysed the traits with and
without the INT transformation over the residuals.

Of note, in both Settings 1 and 2, CBMAT was fitted with the true copula that
generates the dependence structure between the traits. The impact of
misspecifying the traits’ dependence structure on the performance of CBMAT is
given next in the Setting 3.

Setting 3: We aimed to evaluate the robustness of CBMAT when we misspecify the
dependence structure (i.e. true copula) of the two phenotypes. More
specifically, we used a Normal copula to simulate either a bivariate
continuous or a mixed binary-continuous phenotype and then evaluated the type
I error of CBMAT when it is fitted with a Normal, a Clayton or a Frank copula
to model the traits dependence. Also, to mimic analysis of real data where the
phenotypes’ true copula model is unknown, we evaluated the performance of
CBMAT where the fitted copula was chosen following the AIC criterion. Thus,
AIC criterion is derived based on the joint model (3) when both traits are
continuous, or based on the joint model (4) for the mixed-trait case. The
traits’ marginal distributions are set to be Exponential (log-link) or
Student-t ($\nu=3$, identity link) for continuous traits and probit latent
model for the binary traits.

##  4 Results

###  4.1 Type I Error Rates

Setting 1: Empirical type I error rates for continuous bivariate traits are
presented in Table 2, when both traits are distributed following Exponential
distributions, and in Table 3 when $Y_{1}$ follows an Exponential distribution
and $Y_{2}$ follows a Student-t distribution, under the null hypothesis
($\eta=0$). From these tables, one can see that both MURAT and Multi-SKAT
methods have inflated type I error rates when the traits are not normally
distributed (no INT), except for Multi-SKAT with the homogeneous effects
assumption (Multi-SKAT Hom) in the case where $Y_{2}$ follows the Student-t
distribution (Table 3), which shows a conservative type I error rate. Although
the INT transformation seems to be effective for correcting the type I error
rate of both MURAT and Multi-SKAT, the methods become relatively conservative
after the traits’ transformation, except for Multi-SKAT with the heterogeneous
effects assumption, when both traits are distributed following Expenential
distributions (Table 2). The dependence structure of the phenotypes (Normal,
Frank, or Clayton) seem to have similar impact on the results of the two
methods. In contrast, the type I error rate of CBMAT, in all scenarios of
Setting 1, is well controlled when both correct marginal distributions (true)
are used and also when the traits’ margins are chosen based on the AIC
criterion.

Setting 2: Empirical type I error rates for the mixed-trait case are outlined
in Table 4 where $Y_{1}$ follows a Probit marginal model and $Y_{2}$ follows a
marginal Exponential distribution. From this table, one can notice that Multi-
SKAT, in general, shows conservative type I error rates, for both cases: with
and without INT transformation of the traits. MURAT shows a type I error rate
inflation when the traits are not transformed, and it becomes relatively
conservative after the traits’ INT. Again, in this setting, CBMAT has well
controlled type I error rates, in all scenarios.

Setting 3: Empirical type I error rates for potential misspecifications of the
copula model are presented in Table 5. One can notice that using the AIC
criterion for selecting a copula to model the two traits dependence
effectively controls type I error rates of CBMAT. On the other hand,
mispecifiying a copula in the case of weakly correlated traits (Kendall’s
$\tau=0.05$), has no impact on the type error I rate. For the scenarios with
moderate and high residual correlation, the inflation becomes quite important
when the true copula is misspecified, especially for the case of using a
Clayton copula instead of the (correct) Gaussian copula. This can be explained
by the fact that the Clayton copula, in contrary to the Normal copula, models
asymmetric dependence with greater dependence in the negative tail than in the
positive (Cuvelier and Fraiture, 2005).

###  4.2 Power

Settings 1 and 2: We show the empirical power results for the continuous-
traits case in Figures 1 and 2, and for the mixed-traits case in Figure 3. The
traits dependence structure is generated from the Gaussian copula. Empirical
results, when the dependence structure is generated from a Clayton copula are
presented in Appendix E of the Supplementary Materials. For the CBMAT
approach, results were reported when the marginal distributions are selected
based on the best fitted marginal models using the AIC criterion. For MURAT
and Multi-SKAT, we reported their empirical power after INT transformation of
the residuals since both methods have shown inflated type I error rates when
the traits are not normally distributed.

From these figures, one can notice that CBMAT has greater power compared to
MURAT and Multi-SKAT for almost all scenarios. Moreover, CBMAT power to detect
association is not affected when the margins are chosen based on the AIC
criterion compared to when CBMAT is fitted with the true/correct marginal
distributions, as shown in Supplementary Tables of Appendix E.

###  4.3 ALSPAC data analysis

We illustrate our method on the Avon Longitudinal Study of Parents and
Children (ALSPAC) data. This transgenerational prospective study is concerned
with the genetic, epigenetic, biological, psychological and social factors
influencing various indicators of health and social development over a
lifetime (Boyd _and others_ , 2013). The cohort is composed of children born
between 1990 and 1992 in the Bristol area in the United Kingdom. For analysis,
we focused on 1477 subjects for whom the entire genome has been sequenced. In
addition, several clinical phenotypes have been measured for these subjects,
including high-density lipoprotein (HDL) and low-density lipoprotein (LDL),
responsible for cholesterol transport to the liver, and triglycerides (Trigl),
apolipoproteins B (ApoB) and apolipoproteins A1 (ApoA1), which are
respectively the main constituents of low and high density lipoproteins. We
are interested in the association between these phenotypes and two sets of
SNPs falling within the APOC3 and APOA1 genes, located on chromosome 11.
Previous studies have shown that mutations in the APOA1 gene cause a decrease
in HDL levels in the blood, potentially increasing the risk of cardiovascular
disease (Eriksson _and others_ , 2009) while mutations in the APOC3 gene are
associated with low plasma triglyceride levels and reduced risk of ischemic
cardiovascular disease (TG and HDL Working Group of the Exome Sequencing
Project, National Heart, Lung, and Blood Institute _and others_ , 2014.

Pearson correlation coefficient values for each pair of phenotypes are
presented in Appendix E of the Supplementary Materials. As expected, the
correlation between HDL and ApoA1 (0.828) is very high. In addition, HDL is
also strongly correlated with Trigl. Thus, the pairs of phenotypes (HDL,
ApoA1) and (HDL, Trigl) are considered for association analysis. Box plots for
HDL, Trigl, and ApoA1 measurements before and after logarithmic transformation
are presented in Appendix F of the Supplementary Materials. Thus, assuming
ApoA1 and Trigl traits follow log-normal distributions, we normalized them
using a logarithmic transformation. In order to mimic the mixed-trait case of
the analysis of this data, we dichotomized HDL using the sample median as a
cutoff value to create a binary trait. For MURAT and Multi-SKAT, we first
fitted a linear regression model of the HDL binary trait and adjusted for Sex
as a covariate, before applying an INT transformation of the residuals.

The CBMAT approach is fitted directly to the traits without transformations as
the method handles both dichotomous and non-normally-distributed traits. More
precisely, we selected marginal models of the continuous traits, ApoA1 and
Trig, among Gaussian, Gamma and Student-t distributions, based on AIC
criterion. This procedure suggests Student-t distributions with 12 and 3
degrees of freedom as best marginal models for ApoA1 and Trigl, respectively.
The choice of the copula to model the traits dependence structure is also
conducted based on AIC criterion under the null hypothesis. The Frank and
Gaussian copulas are selected as the best joint models for (HDL, ApoA1) and
(HDL, Trigl), respectively. The estimated copula parameters, the corresponding
Kendall’s $\tau$, and pleiotropy coefficient for each bivariate phenotype are
presented in Table 6.

The two sets of SNPs that we analysed consist of 27 SNPs located in the APOA1
gene (MAF ranging from 0.0003 to 0.45), and 59 SNPs located in the APOC3 gene
(MAF ranging from 0.0003 to 0.3). The p-values for the Multi-SKAT, MURAT and
CBMAT approaches are presented in Table 7. As observed in the simulations, the
power of Multi-SKAT decreases more importantly when the pleiotropy is low and
homogeneous effect sizes on different phenotypes are assumed. This may explain
why the association p-values for (HDL,Trigl) are significantly higher than for
(HDL,ApoA1). For all four association tests, CBMAT is more powerful than both
both MURAT and Multi-SKAT, except for the analysis of (HDL, ApoA1), where
Multi-SKAT has smaller p-values. However, Multi-SKAT results need to be
interpreted with caution as the method has shown inflation and/or power loss
in our simulation studies. Finally, recall that, as opposed to the other
methods, CBMAT allows to test for association and provides covariates’ effects
estimation directly on the original phenotypes, without any transformation.

##  5 Discussion

In this study, we have introduced a new copula based method for rare‐variant
association tests for both continuous and mixed bivariate phenotypes, where
binary traits are incorporated through the use of a marginal latent probit
model. As demonstrated, CBMAT is more flexible with regard to modeling the
dependence structure between the traits through the use of copulas. Moreover,
our proposed test is applicable for a wide range of distributions as it does
not assume normality of the traits’ margins. The simulation studies show that
the type I error of our proposed data-driven score-type test is well
controlled for both continuous and mixed bivariate phenotypes, even when both
the true copula and marginal distributions were not known. The simulation
studies investigation also shows that methods which assume multivariate
normality, such as MURAT and Multi-SKAT, are prone to type I error inflation
and a loss of power to detect association when applied to non-normally
distributed outcomes.

Our method is a generalization of MURAT when bivariate normality is not
assumed. Hence, when one uses normal marginal distributions and a Gaussian
copula to model the joint distribution of the traits, CBMAT reduces to MURAT.
One of the main advantages of our method, though, is that we model the
dependence structure between traits regardless of their marginal
distributions. Moreover, the choice among different copula families to model
the traits dependence allows for more modeling flexibility. Furthermore, our
method circumvent inherent complications in terms of analysis and
interpretation resulting from transformation of variables prior to analysis.
Indeed, transformation of the original data may not always result in accurate
interpretation of the estimated coefficients for the covariates and/or the
genetic effect, especially since several clinical measurements can only be
modeled by distributions taking positive values. Our method is also applicable
for studies where there is clinical utility in dichotomizing a trait or when
one outcome consists of disease status, and we aim to analyse such a trait
with other related continuous phenotypes. Most existing region-based genetic
association tests do not have current implementations for the mixture of
binary-continuous traits.

We assumed a bivariate phenotype throughout the paper; however, in many
applications, researchers have access to more than two traits. Although the
mathematical derivation of the CBMAT approach in the case where $K>2$ is
straightforward, work still needs to be made for the method implementation and
its computational efficiency in order to scale to GWAS data sets. Also,
Archimedean copulas in this case, such as Clayton and Frank copulas, might not
be suitable anymore to model dependence between the traits as these two copula
families only allow for exchangeable dependence structure between the $K$
phenotypes. Finally, another useful extension of our methodology would be to
take advantage of the great flexibility provided by copulas to allow for
several forms of clustering, such as in family-based or longitudinal studies.
For instance, in famliy-based designs, one needs to capture both between- and
within-subjects dependencies. Thus, again marginal GLM models could link the
traits to the covariates and the genetic region, and Vine copulas (Joe, 2014)
could be exploited to model the additional within-family source of dependence
and develop a copula-based extension of CBMAT for related subjects.

##  6 Software

A software R package, together with a toy data example for illustration and
complete documentation, are available at https://github.com/julstpierre/CBMAT.

## Acknowledgments

This research supported by both the Fonds de recherche Québec-Santé through
individual grant $\\#$ 267074 and by the Natural Sciences and Engineering
Research Council of Canada through individual discovery research grants to
Karim Oualkacha. This study makes use of data generated by the UK10 K
Consortium, derived from samples from ALSPAC, under data access agreement
ID2250. A full list of the investigators who contributed to the generation of
the data is available from www.UK10K.org. Funding for UK10 K was provided by
the Wellcome Trust under award WT091310.

## Data availability statement

Data is available under request from UK10K consortium or from ALSPAC
consortium. We have access to this data through UK10K consortium since Dr
Oualkacha is a member of the statistics group of the UK10K consortium.

Conflict of Interest: None declared.

## References

  * Bauman _and others_ (2005) Bauman, L. E., Almasy, L., Blangero, J., Duggirala, R., Sinsheimer, J. S. and Lange, K. (2005, September).  Fishing for pleiotropic QTLs in a polygenic sea.  Annals of Human Genetics 69(5), 590–611.
  * Beasley _and others_ (2009) Beasley, T. M., Erickson, S. and Allison, D. B. (2009, September).  Rank-Based Inverse Normal Transformations are Increasingly Used, But are They Merited?  Behavior Genetics 39(5), 580–595.
  * Boyd _and others_ (2013) Boyd, A., Golding, J., Macleod, J. _and others_. (2013, February).  Cohort Profile: The ’children of the 90s’–the index offspring of the Avon Longitudinal Study of Parents and Children.  International Journal of Epidemiology 42(1), 111–127.
  * Cuvelier and Fraiture (2005) Cuvelier, Etienne and Fraiture, Monique Noirhomme. (2005).  Clayton copula and mixture decomposition.  In: Janssen, Jacques and Lenca, Philippe (editors), Applied Stochastic Models and Data Analysis (ASMDA 2005), Brest, 17-20 May 2005.  Publication editors : Jacques Janssen and Philippe Lenca.
  * Davies (1980) Davies, Robert B. (1980).  Algorithm AS 155: The Distribution of a Linear Combination of $\chi$ 2 Random Variables.  Applied Statistics 29(3), 323.
  * de Leon and Wu (2011) de Leon, A. R. and Wu, B. (2011, January).  Copula-based regression models for a bivariate mixed discrete and continuous outcome.  Statistics in Medicine 30(2), 175–185.
  * Dutta _and others_ (2019) Dutta, D., Scott, L., Boehnke, M. and Lee, S. (2019, February).  Multi-SKAT: General framework to test for rare-variant association with multiple phenotypes: DUTTA et al. Genetic Epidemiology 43(1), 4–23.
  * Eriksson _and others_ (2009) Eriksson, M., Schönland, S., Yumlu, S. _and others_. (2009, May).  Hereditary Apolipoprotein AI-Associated Amyloidosis in Surgical Pathology Specimens.  The Journal of Molecular Diagnostics 11(3), 257–262.
  * Ferreira and Purcell (2009) Ferreira, M. A. R. and Purcell, S. M. (2009, January).  A multivariate test of association.  Bioinformatics 25(1), 132–133.
  * Foulkes and Shuen (2013) Foulkes, W. D. and Shuen, A. Y. (2013, August).  In Brief: _BRCA1_ and _BRCA2_.  The Journal of Pathology 230(4), 347–349.
  * Genest _and others_ (2013) Genest, C., Nikoloulopoulos, A. K., Rivest, L. and Fortin, M. (2013, August).  Predicting dependent binary outcomes through logistic regressions and meta-elliptical copulas.  Brazilian Journal of Probability and Statistics 27(3), 265–284.
  * Gottesman and Gould (2003) Gottesman, I. I. and Gould, T. D. (2003, April).  The Endophenotype Concept in Psychiatry: Etymology and Strategic Intentions.  American Journal of Psychiatry 160(4), 636–645.
  * Guo _and others_ (2013) Guo, X., Liu, Z., Wang, X. and Zhang, H. (2013, January).  Genetic association test for multiple traits at gene level.  Genetic Epidemiology 37(1), 122–129.
  * Hall (2010) Hall, M. (2010, January).  A New Role for Endophenotypes in the GWAS Era: Functional Characterization of Risk Variants.  Harvard Review of Psychiatry 18(1), 67–74.
  * Iacono _and others_ (2017) Iacono, W. G., Malone, S. M. and Vrieze, S. I. (2017, January).  Endophenotype best practices.  International Journal of Psychophysiology 111, 115–144.
  * Joe (2014) Joe, Harry. (2014).  Dependence modeling with copulas. Chapman and Hall/CRC.
  * Klei _and others_ (2008) Klei, L., Luca, D., Devlin, B. and Roeder, K. (2008, January).  Pleiotropy and principal components of heritability combine to increase power for association analysis.  Genetic Epidemiology 32(1), 9–19.
  * Konigorski _and others_ (2020) Konigorski, S., Yilmaz, Y. E., Janke, J., Bergmann, M. M., Boeing, H. and Pischon, T. (2020, January).  Powerful rare variant association testing in a copula-based joint analysis of multiple phenotypes.  Genetic Epidemiology 44(1), 26–40.
  * Lee _and others_ (2012) Lee, S., Emond, M. J., Bamshad, M. J., Barnes, K. C., Rieder, M. J., Nickerson, D. A., NHLBI GO Exome Sequencing Project—ESP Lung Project Team, Christiani, D. C., Wurfel, M. M. and Lin, X. (2012, August).  Optimal unified approach for rare-variant association testing with application to small-sample case-control whole-exome sequencing studies.  American Journal of Human Genetics 91(2), 224–237.
  * Lee _and others_ (2017) Lee, S., Won, S., Kim, Y. J., Kim, Y., T2D-Genes Consortium, Kim, B. and Park, T. (2017, April).  Rare variant association test with multiple phenotypes.  Genetic Epidemiology 41(3), 198–209.
  * Lin (1997) Lin, X. (1997, June).  Variance component testing in generalised linear models with random effects.  Biometrika 84(2), 309–326.
  * Liu and Lin (2019) Liu, Z. and Lin, X. (2019, July).  A Geometric Perspective on the Power of Principal Component Association Tests in Multiple Phenotype Studies.  Journal of the American Statistical Association 114(527), 975–990.
  * Lu and Shiou (2002) Lu, T. and Shiou, S. (2002, January).  Inverses of 2 × 2 block matrices.  Computers & Mathematics with Applications 43(1-2), 119–129.
  * Magnus (1978) Magnus, J.R. (1978).  The moments of products of quadratic forms in normal variables.  Statistica Neerlandica 32, 201–210.
  * Nelsen (2010) Nelsen, R. B. (2010).  An Introduction to Copulas, 2. ed. 2006. Corr. 2. pr. Softcover version of original hardcover edition 2006 edition., Springer Series in Statistics. New York, NY: Springer New York.
  * Petersen _and others_ (2013) Petersen, A., Alvarez, C., DeClaire, S. and Tintle, N. L. (2013).  Assessing methods for assigning SNPs to genes in gene-based tests of association using common variants.  PloS One 8(5), e62161.
  * Schaid _and others_ (2016) Schaid, D. J., Tong, X., Larrabee, B. _and others_. (2016, October).  Statistical Methods for Testing Genetic Pleiotropy.  Genetics 204(2), 483–497.
  * Seoane _and others_ (2014) Seoane, J. A., Campbell, C., Day, Ian N. M. _and others_. (2014, October).  Canonical Correlation Analysis for Gene-Based Pleiotropy Discovery.  PLoS Computational Biology 10(10), e1003876.
  * Shriner (2012) Shriner, D. (2012).  Moving toward System Genetics through Multiple Trait Analysis in Genome-Wide Association Studies.  Frontiers in Genetics 3.
  * Solovieff _and others_ (2013) Solovieff, N., Cotsapas, C., Lee, P. H. _and others_. (2013, July).  Pleiotropy in complex traits: Challenges and strategies.  Nature Reviews Genetics 14(7), 483–495.
  * Stearns (2010) Stearns, Frank W. (2010, November).  One Hundred Years of Pleiotropy: A Retrospective.  Genetics 186(3), 767–773.
  * Sun _and others_ (2016) Sun, J., Oualkacha, K., Forgetta, V. _and others_. (2016, September).  A method for analyzing multiple continuous phenotypes in rare variant association studies allowing for flexible correlations in variant effects.  European Journal of Human Genetics 24(9), 1344–1351.
  * Sun _and others_ (2019) Sun, Jianping, Oualkacha, Karim, Greenwood, Celia M.T. and Lakhal-Chaieb, Lajmi. (2019, March).  Multivariate association test for rare variants controlling for cryptic and family relatedness.  Canadian Journal of Statistics 47(1), 90–107.
  * Tang and Ferreira (2012) Tang, C. S. and Ferreira, M. A. R. (2012, March).  A gene-based test of association using canonical correlation analysis.  Bioinformatics 28(6), 845–850.
  * TG and HDL Working Group of the Exome Sequencing Project, National Heart, Lung, and Blood Institute _and others_ (2014) TG and HDL Working Group of the Exome Sequencing Project, National Heart, Lung, and Blood Institute _and others_. (2014, July).  Loss-of-function mutations in APOC3, triglycerides, and coronary disease.  The New England Journal of Medicine 371(1), 22–31.
  * The 1000 Genomes Project Consortium (2015) The 1000 Genomes Project Consortium. (2015, October).  A global reference for human genetic variation.  Nature 526(7571), 68–74.
  * Van der Sluis _and others_ (2015) Van der Sluis, S., Dolan, . V., Li, J. _and others_. (2015, April).  MGAS: A powerful tool for multivariate gene-based genome-wide association analysis.  Bioinformatics (Oxford, England) 31(7), 1007–1015.
  * Wang _and others_ (2015) Wang, Y., Liu, A., Mills, J. L. _and others_. (2015, May).  Pleiotropy analysis of quantitative traits at gene level by multivariate functional linear models.  Genetic Epidemiology 39(4), 259–275.
  * Wei _and others_ (2019) Wei, Y., Liu, Y., Chen, W. and Ding, Y. (2019, April).  Gene-based Association Analysis for Bivariate Time-to-event Data through Functional Regression with Copula Models.  arXiv:1904.01116 [stat].
  * Williams (1957) Williams, G. C. (1957, December).  Pleiotropy, Natural Selection, and the Evolution of Senescence.  Evolution 11(4), 398.
  * Wu _and others_ (2011) Wu, M. C., Lee, S., Cai, T. _and others_. (2011, July).  Rare-variant association testing for sequencing data with the sequence kernel association test.  American Journal of Human Genetics 89(1), 82–93.
  * Yang _and others_ (2016) Yang, J. J., Li, J., Williams, L. K. and Buu, A. (2016, December).  An efficient genome-wide association test for multivariate phenotypes based on the Fisher combination function.  BMC Bioinformatics 17(1), 19.
  * Yang and Wang (2012) Yang, Q. and Wang, Y. (2012).  Methods for Analyzing Multivariate Phenotypes in Genetic Association Studies.  Journal of Probability and Statistics 2012, 1–13.
  * Zhang _and others_ (2014) Zhang, Y., Xu, Z., Shen, X. _and others_. (2014, August).  Testing for association with multiple traits in generalized estimation equations, with application to neuroimaging data.  NeuroImage 96, 309–325.

Table 1: Values of $v$, $\rho$, $\tau$ in the simulations (leading to 21 distinct scenarios). $H_{0}$ and $H_{a}$ stand for the null and the alternative hypotheses of the association between the genetic region and the two phenotypes. a fraction of variants that are causal and induce the association between the genetic region and the two traits; b region-specific pleiotropy effect; c Kendall’s $\tau$ in the CBMAT copula model. | $v^{a}$ | $\rho^{b}$ | Kendall’s $\tau^{c}$
---|---|---|---
$H_{0}$ | 0 | 0 | {0.05, 0.20, 0.40}
$H_{a}$ | {10%, 20%} | {0, 0.4, 0.8} | {0.05, 0.20, 0.40}
Figure 1: Power $(\%)$ to detect SNPs/phenotype association. Results are computed from 5 000 data sets generated under a Gaussian copula model when either $10\%$ or $20\%$ of the variants are causal and when both traits are continuous and both follow Exponential marginal distributions. The variance component parameter $\eta$ was fixed such that the traits heritabilty under a linear model is approximately $2\%$. The compared methods are CBMAT fitted using marginal models chosen based on AIC, MURAT and Multi-SKAT fitted with INT transformation after adjusting for the covariates; Multi-SKAT is fitted with homogeneous (Multi-SKAT Hom) and uncorrelated (Multi-SKAT Het) effect sizes. Figure 2: Power $(\%)$ to detect SNPs/phenotype association. Results are computed from 5 000 data sets generated under a Gaussian copula model when either $10\%$ or $20\%$ of the variants are causal and when both traits are continuous and follow Exponential/Student-t marginal distributions. The variance component parameter $\eta$ was fixed such that the traits heritabilty under a linear model is approximately $2\%$. The compared methods are CBMAT fitted using marginal models chosen based on AIC, MURAT and Multi-SKAT fitted with INT transformation after adjusting for the covariates; Multi-SKAT is fitted with homogeneous (Multi-SKAT Hom) and uncorrelated (Multi-SKAT Het) effect sizes. Figure 3: Power $(\%)$ to detect SNPs/phenotype association. Results are computed from 5 000 data sets generated under a Gaussian copula model when either $10\%$ or $20\%$ of the variants are causal and when traits are mixed binary-continuous and follow Probit/Exponential marginal distributions. The variance component parameter $\eta$ was fixed such that the traits heritabilty under a linear model is approximately $2\%$. The compared methods are CBMAT fitted using marginal models chosen based on AIC, MURAT and Multi-SKAT fitted with INT transformation after adjusting for the covariates; Multi-SKAT is fitted with homogeneous (Multi-SKAT Hom) and uncorrelated (Multi-SKAT Het) effect sizes. Table 2:  Empirical type 1 error rate $(\%)$ under the null hypothesis of no SNPs/phenotype association $(\eta=0)$ where the data are generated under either Gaussian, Frank or Clayton copula models. Results are computed from 10 000 data sets generated under Setting 1 where both traits are continuous and both follow Exponential marginal distributions. The compared methods are CBMAT fitted with true marginal models (true), and fitted using marginal models chosen based on AIC (AIC); MURAT and Multi-SKAT are fitted without INT transformation (no INT) and with INT transformation after adjusting for the covariates (INT); Multi-SKAT is fitted with homogeneous (Multi-SKAT Hom) and uncorrelated (Multi-SKAT Het) effect sizes. |  | CBMAT | MURAT | Multi-SKAT Het | Multi-SKAT Hom
---|---|---|---|---|---
copula | Kendall’s $\tau$ | true | AIC | no INT | INT | no INT | INT | no INT | INT
Normal | 0.05 | 1.11 | 1.11 | 2.34 | 0.84 | 4.10 | 1.08 | 4.08 | 0.60
| 0.20 | 1.19 | 1.18 | 2.27 | 0.84 | 4.13 | 1.18 | 4.10 | 0.59
| 0.40 | 1.19 | 1.19 | 2.27 | 0.90 | 4.12 | 1.16 | 4.10 | 0.45
Franck | 0.05 | 1.14 | 1.14 | 2.26 | 0.80 | 4.21 | 1.10 | 4.19 | 0.58
| 0.20 | 1.22 | 1.22 | 2.31 | 0.80 | 4.24 | 1.14 | 4.18 | 0.58
| 0.40 | 1.14 | 1.14 | 2.31 | 0.89 | 4.31 | 1.18 | 4.26 | 0.48
Clayton | 0.05 | 1.06 | 1.06 | 2.33 | 0.76 | 4.15 | 1.09 | 4.16 | 0.59
| 0.20 | 1.18 | 1.18 | 2.30 | 0.88 | 4.12 | 1.10 | 4.12 | 0.59
| 0.40 | 1.06 | 1.06 | 2.26 | 0.82 | 4.09 | 1.05 | 4.07 | 0.56
Table 3: Empirical type 1 error rate $(\%)$ under the null hypothesis of no SNPs/phenotype association $(\eta=0)$ where the data are generated under either Gaussian, Frank or Clayton copula models. Results are computed from 10 000 data sets generated under Setting 1 where both traits are continuous and follow Exponential/Student-t marginal distributions. The compared methods are CBMAT fitted with true marginal models (true), and fitted using marginal models chosen based on AIC (AIC); MURAT and Multi-SKAT are fitted without INT transformation (no INT) and with INT transformation after adjusting for the covariates (INT); Multi-SKAT is fitted with homogeneous (Multi-SKAT Hom) and uncorrelated (Multi-SKAT Het) effect sizes. |  | CBMAT | MURAT | Multi-SKAT Het | Multi-SKAT Hom
---|---|---|---|---|---
copula | Kendall’s $\tau$ | true | AIC | no INT | INT | no INT | INT | no INT | INT
Normal | 0.05 | 0.93 | 0.96 | 1.57 | 0.93 | 3.50 | 0.76 | 0.02 | 0.77
| 0.20 | 1.00 | 1.02 | 1.35 | 0.91 | 3.95 | 0.67 | 0.01 | 0.89
| 0.40 | 1.07 | 1.10 | 1.37 | 0.89 | 4.48 | 0.86 | 0.02 | 0.85
Franck | 0.05 | 0.92 | 0.95 | 1.61 | 0.90 | 3.46 | 0.75 | 0.02 | 0.79
| 0.20 | 0.91 | 0.93 | 1.48 | 0.91 | 3.94 | 0.73 | 0.02 | 0.88
| 0.40 | 0.98 | 1.04 | 1.53 | 0.95 | 4.50 | 0.88 | 0.02 | 0.88
Clayton | 0.05 | 0.86 | 0.90 | 1.43 | 0.89 | 3.54 | 0.80 | 0.02 | 0.81
| 0.20 | 0.95 | 0.99 | 1.36 | 0.96 | 3.90 | 0.69 | 0.02 | 0.80
| 0.40 | 1.01 | 1.03 | 1.37 | 0.83 | 4.36 | 0.83 | 0.03 | 0.85
Table 4: Empirical type 1 error rate $(\%)$ under the null hypothesis of no SNPs/phenotype association $(\eta=0)$ where the data are generated under either Gaussian, Frank or Clayton copula models. Results are computed from 10 000 data sets generated under Setting 2 for mixed-trait case where $Y_{1}$ follows a Probit marginal model and $Y_{2}$ follows a marginal Exponential distribution. The compared methods are CBMAT fitted with true marginal models (true), and fitted using marginal models chosen based on AIC (AIC); MURAT and MUlti-SKAT are fitted without INT transformation (no INT) and with INT transformation after adjusting for the covariates (INT); Multi-SKAT is fitted with homogeneous (Multi-SKAT Hom) and uncorrelated (Multi-SKAT Het) effect sizes. |  | CBMAT | MURAT | Multi-SKAT Het | Multi-SKAT Hom
---|---|---|---|---|---
copula | Kendall’s $\tau$ | true | AIC | no INT | INT | no INT | INT | no INT | INT
Normal | 0.05 | 1.14 | 1.14 | 2.17 | 0.78 | 0.03 | 0.40 | 0.03 | 0.64
| 0.20 | 1.15 | 1.15 | 2.19 | 0.81 | 0.03 | 0.40 | 0.02 | 0.55
| 0.40 | 1.19 | 1.19 | 2.13 | 0.82 | 0.02 | 0.44 | 0.02 | 0.62
Franck | 0.05 | 1.10 | 1.10 | 2.16 | 0.78 | 0.03 | 0.40 | 0.03 | 0.63
| 0.20 | 1.11 | 1.11 | 2.20 | 0.83 | 0.02 | 0.37 | 0.02 | 0.55
| 0.40 | 1.05 | 1.05 | 2.11 | 0.79 | 0.02 | 0.43 | 0.02 | 0.61
Clayton | 0.05 | 1.12 | 1.12 | 2.18 | 0.82 | 0.03 | 0.43 | 0.03 | 0.61
| 0.20 | 1.08 | 1.08 | 2.10 | 0.78 | 0.02 | 0.42 | 0.02 | 0.60
| 0.40 | 1.06 | 1.06 | 2.18 | 0.80 | 0.01 | 0.37 | 0.01 | 0.69
Table 5: Empirical type 1 error rate $(\%)$ under the null hypothesis of no SNPs/phenotype association $(\eta=0)$. Results are computed from 10 000 data sets generated under Setting 3 where we used a Normal copula to simulate either a bivariate continuous or a mixed binary-continuous phenotype and then fitted CBMAT with a Normal, Clayton or Frank copula, or a model where the copula was chosen based on AIC. |  |  | Copula
---|---|---|---
$F_{1}$ | $F_{2}$ | Kendall’s $\tau$ | Normal | Clayton | Frank | AIC
Binomial | Exponential | 0.05 | 1.00 | 1.03 | 1.08 | 1.05
|  | 0.20 | 1.07 | 1.24 | 1.14 | 1.08
|  | 0.40 | 1.12 | 1.44 | 1.31 | 1.09
Exponential | Exponential | 0.05 | 1.31 | 1.32 | 1.29 | 1.27
|  | 0.20 | 1.34 | 1.53 | 1.42 | 1.34
|  | 0.40 | 1.26 | 2.10 | 1.77 | 1.26
| Student | 0.05 | 1.07 | 1.19 | 1.11 | 1.11
|  | 0.20 | 1.13 | 1.47 | 1.25 | 1.24
|  | 0.40 | 1.12 | 1.66 | 1.52 | 1.51
Table 6: Estimated copula parameter $\hat{\theta}$, corresponding Kendall’s $\hat{\tau}$ and pleitropic coefficient $\hat{\rho}$ for each mixed bivariate phenotype. Phenotype | Copula | $\hat{\theta}$ | Kendall’s $\hat{\tau}$ | $\hat{\rho}$
---|---|---|---|---
(HDL,ApoA1) | Frank | 7.05 | 0.56 | 0.4
(HDL,Trigl) | Gaussian | -0.47 | -0.31 | 0
Table 7: P-values for Multi-SKAT, MURAT and CBMAT under the null hypothesis of no gene/phenotype association for APOA1 and APOC3 genes. HDL trait has been dichotomized, while ApoA1 and Trigl traits have been log-transformed prior to analysis. All association tests have been adjusted for sex as a potential confounder. Gene |
---|---
Phenotype | Multi-SKAT Het | Multi-SKAT Hom | MURAT | CBMAT
APOA1 |  |  |  |
\- (HDL,ApoA1) | $2.44\times 10^{-3}$ | $9.40\times 10^{-4}$ | $5.38\times 10^{-2}$ | $5.99\times 10^{-3}$
\- (HDL,Trigl) | $1.10\times 10^{-2}$ | $1.77\times 10^{-2}$ | $4.49\times 10^{-4}$ | $6.82\times 10^{-5}$
APOC3 |  |  |  |
\- (HDL,ApoA1) | $1.99\times 10^{-3}$ | $6.15\times 10^{-4}$ | $9.36\times 10^{-3}$ | $1.86\times 10^{-4}$
\- (HDL,Trigl) | $2.61\times 10^{-2}$ | $2.50\times 10^{-2}$ | $6.00\times 10^{-4}$ | $2.29\times 10^{-4}$