C-JAMP: Copula-based Joint Analysis of Multiple Phenotypes
Overview
The general goal of C-JAMP (Copula-based Joint Analysis of Multiple Phenotypes) is to jointly model two (or more) traits Y1, Y2 of a phenotype conditional on covariates using copula functions, in order to estimate and test the association of the covariates with either trait.
Copulas are functions used to construct a joint distribution by combining the marginal distributions with a dependence structure (Joe 1997, Nelsen 2006). In contrast to multivariate linear regression models, which model the linear dependence of two random variables coming from a bivariate normal distribution with a certain Pearson’s correlation, more general dependencies can be flexibly investigated using copula functions. This can be helpful when the conditional mean of Yi given Yj is not linear in Yj. Additional properties of copula models are that they allow an investigation of the dependence structure between the phenotypes separately from the marginal distributions. Also, the marginal distributions can come from different families and not necessarily from the normal distribution.
Copula models can be used to make inference about the marginal parameters and to identify markers associated with one or multiple phenotypes. Through the joint modeling of multiple traits, the power of the association test with a given trait can be increased, also if the markers are only associated with one trait. This is of special interest for genetic association studies, and in particular for the analysis of rare genetic variants, where the low power of association tests is one of the main challenges in empirical studies. The potential of C-JAMP to increase the power of association tests is supported by the results of empirical applications of C-JAMP in real-data analyses (Konigorski, Yilmaz & Bull, 2014; Konigorski, Yilmaz, Pischon, 2016).
The functions in this package implement C-JAMP for fitting joint models based on the Clayton and 2-parameter copula, of two normally-distributed traits Y1, Y2 conditional on multiple predictors, and allow obtaining coefficient and standard error estimates of the copula parameters (modeling the dependence between Y1, Y2) and marginal parameters (modeling the effect of the predictors on Y1 and on Y2), as well as performing Wald-type hypothesis tests of the marginal parameters. Summary functions provide a user-friendly output. In addition, likelihood-ratio tests for testing nested copula models are available. Furthermore, functions are available to generate genetic data, and phenotypic data from bivariate normal distrutions and bivariate distributions based on copula functions. Finally, for genetic association analyses, functions are available to calculate the amount of phenotypic variance that can be explained by the genetic markers.
Background on copula functions
In more detail, a d-dimensional copula can be defined as a function C : [0, 1]d → [0, 1], which is a joint cumulative distribution function with uniform marginal cumulative distribution functions. Hence, C is the distribution of a multivariate random vector, and can be considered independent of the margins. For Y1, …, Yd and a covariate vector x, the joint distribution F of Y1, …, Yd, conditional on x, can be constructed by combining the marginal distributions of Y1, …, Yd, F1, …, Fd, conditional on x, using a copula function Cψ with dependence parameter(s) ψ:
F(Y1, …, Yd|x) = Cψ(F1(Y1|x), …, Fd(Yd|x)). For continuous marginal distributions, the copula Cψ is unique and the multivariate distribution can be constructed from the margins in a uniquely defined way (Sklar, 1959).
Popular copula functions include the Clayton family
$$ C_\phi (u_1, \dots, u_d, \phi) = \left( \sum_{i=1}^d u_i^{-\phi} - (d-1) \right)^{ -\frac{1}{\phi}} $$
with ϕ > 0, and the Gumbel-Hougaard family
$$ C_\theta (u_1, \dots, u_d, \theta) = exp\left( - \left[ \sum_{i=1}^d(-log(u_i))^\theta \right]^{\frac{1}{\theta}} \right) $$
with θ > 1. A third family which includes both Clayton and Gumbel-Hougaard copula (for θ = 1 and ϕ → 0) is the 2-paramater copula family
$$ C_\psi (u_1, \dots, u_d, \phi, \theta) = \left( \left[ \sum_{i=1}^d \left( u_i^{-\phi}-1 \right)^\theta \right]^{\frac{1}{\theta}} + 1 \right)^{-\frac{1}{\phi}} $$
with 0 ≤ u1, u2 ≤ 1, and the copula parameters ψ = (ϕ, θ), ϕ > 0, θ ≥ 1, which allows a flexible modeling of both the lower- and upper-tail dependence.
For the 2-parameter copula family, Kendall’s tau can be derived as (Joe, 1997)
$$ \tau = 1 - \frac{2}{\theta(\phi+2)}. $$ Additional dependence parameters of interest are the lower and upper tail dependence measures λL, λU, which explain the amount of dependence between extreme values. For the 2-parameter copula family, these dependence measures become (Joe, 1997)
$$ \lambda_L=2^{-\frac{1}{\theta\phi}}, \
\lambda_U= 2-2^{\frac{1}{\theta}}.$$ Data Y1, Y2
can be sampled from the Clayton copula (with standard normal margins)
using the generate_clayton_copula() function. This uses the
following steps to generate standard normal Y1, Y2
from the Clayton copula Cϕ(Φ(Y1), Φ(Y2)) = Cϕ(U1, U2),
where Φ is the standard normal
cumulative distribution function:
- Generate uniform random variables U1, Ũ2 ∼ Unif(0, 1).
- Generate Y1 from the inverse standard normal distribution of U1, Y1 := Φ−1(U1).
- In order to obtain the second variable Y2, use the conditional distribution of U2 given U1, set $\tilde{U}_2 = \frac{\partial C(U_1, U_2)}{\partial U_1}$ and solve for U2. This yields $U_2 = \left( 1 - U_1^{-\phi} \cdot \left( 1 - \tilde{U}_2^{-\frac{\phi}{1 + \phi}} \right) \right)^{-\frac{1}{\phi}}$.
- Generate Y2 from the inverse standard normal distribution of U2, Y2 := Φ−1(U2).
# Generate phenotype data from the Clayton copula:
set.seed(10)
dat1a <- generate_clayton_copula(n = 1000, phi = 0.5)## [1] "Kendall's tau between Y1, Y2 = 0.2"## [1] "Kendall's tau between Y1, Y2 = 0.5"## [1] "Kendall's tau between Y1, Y2 = 0.8"par(mfrow = c(1, 3))
plot(dat1a$Y1, dat1a$Y2, xlab = "Y1", ylab = "Y2",
main = expression(paste("Scatterplot of ", Y[1], ", ", Y[2], " with ", tau, "=0.2")))
plot(dat1b$Y1, dat1b$Y2, xlab = "Y1", ylab = "Y2",
main = expression(paste("Scatterplot of ", Y[1], ", ", Y[2], " with ", tau, "=0.5")))
plot(dat1c$Y1, dat1c$Y2, xlab = "Y1", ylab = "Y2",
main = expression(paste("Scatterplot of ", Y[1], ", ", Y[2], " with ", tau, "=0.8")))
Data generation functions
In order to generate sample data for simulation studies and to
illustrate C-JAMP, different functions are provided. Firstly, the
functions generate_singleton_data(),
generate_doubleton_data() and
generate_genodata() can be used to generate genetic data in
the form of single nucleotide variants (SNVs).
generate_singleton_data() generates singletons (i.e., SNVs
with one observed minor allele); generate_doubleton_data()
generates doubletons (i.e., SNVs with two observed minor alleles), and
the function generate_genodata() generates n observations of k SNVs with random minor allele
frequencies. The minor allele frequencies can be computed using the
function compute_MAF().
Next, generate_phenodata_1_simple(),
generate_phenodata_1(),
generate_phenodata_2_bvn() and
generate_phenodata_2_copula() can be used to generate
phenotype data based on SNV input data, covariates, and a specification
of the covariate effect sizes. Here, the functions
generate_phenodata_1_simple() and
generate_phenodata_1() generate one normally-distributed or
binary phenotype Y conditional
on provided SNVs and two generated covariates X1, X2. The functions
generate_phenodata_2_bvn() and
generate_phenodata_2_copula() generate two phenotypes Y1, Y2 with dependence
Kendall’s τ conditional on the
provided SNVs and covariates X1, X2 from the bivariate
normal distribution or the Clayton copula function with standard normal
marginal distributions. generate_phenodata_2_copula() is
using the function generate_clayton_copula() described
above.
# Generate genetic data:
set.seed(10)
genodata <- generate_genodata(n_SNV = 20, n_ind = 1000)
compute_MAF(genodata)## SNV1 SNV2 SNV3 SNV4 SNV5 SNV6 SNV7 SNV8 SNV9 SNV10 SNV11 SNV12 SNV13
## 0.254 0.022 0.195 0.253 0.066 0.119 0.421 0.474 0.157 0.049 0.175 0.461 0.168
## SNV14 SNV15 SNV16 SNV17 SNV18 SNV19 SNV20
## 0.064 0.114 0.097 0.102 0.095 0.466 0.195# Generate phenotype data from the bivariate normal distribution given covariates:
phenodata_bvn <- generate_phenodata_2_bvn(genodata = genodata,
tau = 0.5, b1 = 1, b2 = 2)
plot(phenodata_bvn$Y1, phenodata_bvn$Y2, xlab = "Y1", ylab = "Y2",
main = expression(paste("Scatterplot of bivariate normal ", Y[1], ", ",
Y[2], " with ", tau, "=0.5")))
# Generate phenotype data from the Clayton copula given covariates:
phenodata <- generate_phenodata_2_copula(genodata = genodata$SNV1,
MAF_cutoff = 1, prop_causal = 1,
tau = 0.5, b1 = 0.3, b2 = 0.3)## [1] "Kendall's tau between Y1, Y2 = 0.5"
## [1] "Kendall's tau between Y1, Y2 = 0.5"plot(phenodata$Y1, phenodata$Y2, xlab = "Y1", ylab = "Y2",
main = expression(paste("Scatterplot of ", Y[1], ", ", Y[2],
" from the Clayton copula with ", tau, "=0.5")))
For the analysis of genetic associations with the phenotypes, the
amount of variability in the phenotypes that can be explained by the
SNVs might be of interest. This can be computed based on four different
approaches using the function compute_expl_var() (Laird
& Lange, 2011), which is described in more detail in the help pages
of the function.
compute_expl_var(genodata = genodata, phenodata = phenodata$Y1,
type = c("Rsquared_unadj", "Rsquared_adj",
"MAF_based", "MAF_based_Y_adjusted"),
causal_idx = rep(TRUE,20), effect_causal =
c(0.3 * abs(log10(compute_MAF(genodata$SNV1))), rep(0,19)))## $Rsquared_unadj
## [1] 0.02794911
##
## $Rsquared_adj
## [1] 0.008091075
##
## $MAF_based
## [1] 0.01208152
##
## $MAF_based_Y_adjusted
## [1] 0.00867618C-JAMP: Copula-based joint analysis of multiple phenotypes
C-JAMP is implemented for the joint analysis of two phenotypes Y1, Y2 conditional on one or multiple predictors x with linear marginal models. Functions are available to fit the joint model
F(Y1, Y2|x) = Cψ(F1(Y1|x), F2(Y2|x)) for the Clayton and 2-parameter copula, with marginal models
Y1 = αTX + ϵ1, ϵ1 ∼ N(0, σ12)
Y2 = βTX + ϵ2, ϵ2 ∼ N(0, σ22).
Maximum likelihood estimates for all parameters (ψ, α, β, σ1, σ2)
can be obtained using the function get_estimates_naive(),
which is based on the ´lm()´ function.
predictors <- data.frame(X1 = phenodata$X1, X2 = phenodata$X2, SNV = genodata$SNV1)
get_estimates_naive(Y1 = phenodata$Y1, Y2 = phenodata$Y2, predictors_Y1 = predictors,
predictors_Y2 = predictors, copula_param = "both")## log_phi log_theta_minus1 Y1_log_sigma Y2_log_sigma
## 1.021100321 0.327953141 0.015511896 -0.006361787
## Y1_(Intercept) Y1_X1 Y1_X2 Y1_SNV
## -0.011581527 0.492133186 0.583058174 0.129050968
## Y2_(Intercept) Y2_X1 Y2_X2 Y2_SNV
## -0.041396358 0.483807892 0.621047800 0.082656354For these or any other parameter estimates, the log-likelihood (or
rather the minus log-likelihood) of the copula model can be computed
with the function minusloglik() for the Clayton and
2-parameter copula:
predictors <- data.frame(X1 = phenodata$X1, X2 = phenodata$X2, genodata[, 1:5])
estimates <- get_estimates_naive(Y1 = phenodata$Y1, Y2 = phenodata$Y2,
predictors_Y1 = predictors,
predictors_Y2 = predictors, copula_param = "both")
minusloglik(Y1 = phenodata$Y1, Y2 = phenodata$Y2, predictors_Y1 = predictors,
predictors_Y2 = predictors, parameters = estimates, copula = "2param")## [1] 3583.498It should be noted that the minusloglik() function (and
the cjamp(), cjamp_loop() functions below)
assume quantitative predictors and use an additive model. That means
that for the analysis of categorical predictors with more than 2 levels,
dummy variables have to be created beforehand. Accordingly, if single
nucleotide variants (SNVs) are included as predictors, the computation
is based on an additive genetic model if SNVs are provided as 0-1-2
genotypes and on a dominant model if SNVs are provided as 0-1
genotypes.
The lrt_copula() can be used to test the model fit of
different nested copula models with the same marginal models. For this,
likelihood ratio tests are performed. In more detail, to test the fit of
the 1-parameter Clayton copula model fit versus the 2-parameter copula,
i.e. H0 : ϕ → 0, the
likelihood ratio test statistic
logLR = 2 ⋅ (L(α̂, β̂, ψ̂) − L(α̂(ψ0), β̂(ψ0), ψ̂)) can be used, where ψ̂ = (ϕ̂, θ̂), $\hat{\psi_0}=\left(\hat{\phi}(\theta=1), 1\right)$, see Yilmaz & Lawless (2011). Here, α̂(ψ0) and β̂(ψ0) are the maximum likelihood estimates of α and β under the null model ψ = ψ0. P-values are obtained by using that logLR is asymptotically distributed as 0.5 + 0.5χ12 under the null hypothesis ψ = ψ0 and under the assumption that other free parameters in ψ0 don’t lie on the boundary of the parameter space (Self & Liang, 1987).
# Example: Test whether 2-parameter copula model has a better
# model fit compared to Clayton copula (no).
predictors <- data.frame(X1 = phenodata$X1, X2 = phenodata$X2,
SNV = genodata$SNV1)
estimates_c <- get_estimates_naive(Y1 = phenodata$Y1, Y2 = phenodata$Y2,
predictors_Y1 = predictors,
predictors_Y2 = predictors,
copula_param = "phi")
minusloglik_Clayton <- minusloglik(Y1 = phenodata$Y1, Y2 = phenodata$Y2,
predictors_Y1 = predictors,
predictors_Y2 = predictors,
parameters = estimates_c, copula = "Clayton")
estimates_2p <- get_estimates_naive(Y1 = phenodata$Y1, Y2 = phenodata$Y2,
predictors_Y1 = predictors,
predictors_Y2 = predictors,
copula_param = "both")
minusloglik_2param <- minusloglik(Y1 = phenodata$Y1, Y2 = phenodata$Y2,
predictors_Y1 = predictors,
predictors_Y2 = predictors,
parameters = estimates_2p, copula = "2param")
lrt_copula(minusloglik_Clayton, minusloglik_2param)## $chisq
## [1] -2318.222
##
## $pval
## [1] 0.5Log-likelihood ratio tests of marginal parameters within the same
copula model can be performed using the lrt_param()
function:
# Example: Test marginal parameters (alternative model has better fit).
predictors_1 <- data.frame(X1 = phenodata$X1, X2 = phenodata$X2)
estimates_1 <- get_estimates_naive(Y1 = phenodata$Y1, Y2 = phenodata$Y2,
predictors_Y1 = predictors_1,
predictors_Y2 = predictors_1,
copula = "phi")
minusloglik_1 <- minusloglik(Y1 = phenodata$Y1, Y2 = phenodata$Y2,
predictors_Y1 = predictors_1,
predictors_Y2 = predictors_1,
parameters = estimates_1, copula = "Clayton")
predictors_2 <- data.frame(X1 = phenodata$X1, X2 = phenodata$X2,
SNV = genodata$SNV1)
estimates_2 <- get_estimates_naive(Y1 = phenodata$Y1, Y2 = phenodata$Y2,
predictors_Y1 = predictors_2,
predictors_Y2 = predictors_2,
copula = "phi")
minusloglik_2 <- minusloglik(Y1 = phenodata$Y1, Y2 = phenodata$Y2,
predictors_Y1 = predictors_2,
predictors_Y2 = predictors_2,
parameters = estimates_2, copula = "Clayton")
lrt_param(minusloglik_1, minusloglik_2, df=2)## $chisq
## [1] 3.632676
##
## $pval
## [1] 0.1626202Finally, to obtain maximum likelihood estimates of the parameters
(ψ, α, β, σ1, σ2)
as well as standard error estimates under the Clayton or 2-parameter
copula models, the function cjamp() can be used.
cjamp() uses the optimx() function in the
optimx package to maximize the log-likelihood function (i.e., minimize
minusloglik()). For this, the BFGS optimization method is
recommended. In order to deal with convergence problems, several checks
are built-in, and different starting values for the optimization
algorithm are automatically tried. Standard error estimates of the
parameter estimates are obtained from the observed inverse information
matrix. Finally, p-values are computed for the marginal parameter
estimates from hypothesis tests of the absence of effects of each
predictor on each phenotype in the marginal models.
# Restrict example to sample size 100 to decrease running time:
predictors <- data.frame(X1 = phenodata$X1, X2 = phenodata$X2, genodata[, 1:3])[1:100,]
cjamp_res <- cjamp(copula = "2param", Y1 = phenodata$Y1[1:100], Y2 = phenodata$Y2[1:100],
predictors_Y1 = predictors, predictors_Y2 = predictors,
scale_var = FALSE, optim_method = "BFGS", trace = 0,
kkt2tol = 1E-16, SE_est = TRUE, pval_est = TRUE,
n_iter_max = 10)
cjamp_res## $`Parameter point estimates`
## Y1_sigma Y2_sigma Y1_(Intercept) Y1_X1 Y1_X2
## 9.343665e-01 9.105235e-01 8.306169e-02 5.062952e-01 5.641714e-01
## Y1_SNV1 Y1_SNV2 Y1_SNV3 Y2_(Intercept) Y2_X1
## 3.707410e-02 3.854541e-01 -9.356935e-02 -5.692639e-02 3.778275e-01
## Y2_X2 Y2_SNV1 Y2_SNV2 Y2_SNV3 phi
## 4.996719e-01 1.457613e-01 3.059386e-01 -8.816050e-02 2.216355e+00
## theta tau lambda_l lambda_u
## 1.000018e+00 5.256651e-01 7.314427e-01 2.473863e-05
##
## $`Parameter standard error estimates`
## Y1_sigma Y2_sigma Y1_(Intercept) Y1_X1 Y1_X2
## 0.060523619 0.059930322 0.149761801 0.090333461 0.170380948
## Y1_SNV1 Y1_SNV2 Y1_SNV3 Y2_(Intercept) Y2_X1
## 0.126321437 0.473962562 0.155398437 0.155211728 0.085813211
## Y2_X2 Y2_SNV1 Y2_SNV2 Y2_SNV3 phi
## 0.166817053 0.125764605 0.438079460 0.149839103 0.401091962
## theta tau lambda_l lambda_u
## 0.001014005 0.045121585 0.041395510 0.001405642
##
## $`Parameter p-values`
## Y1_(Intercept) Y1_X1 Y1_X2 Y1_SNV1 Y1_SNV2
## 5.791510e-01 2.085729e-08 9.288517e-04 7.691475e-01 4.160699e-01
## Y1_SNV3 Y2_(Intercept) Y2_X1 Y2_X2 Y2_SNV1
## 5.470907e-01 7.137935e-01 1.068104e-05 2.741495e-03 2.464558e-01
## Y2_SNV2 Y2_SNV3
## 4.849500e-01 5.562855e-01
##
## $`Convergence code of optimx function`
## convcode
## 0
##
## $`Karush-Kuhn-Tucker conditions 1 and 2`
## KKT1 KKT2
## TRUE TRUE
##
## $`Maximum log-likelihood`
## maxloglik
## 222.3041
##
## attr(,"class")
## [1] "cjamp"If the goal is to test a large number of predictors in the same
marginal models (such as in genetic association studies), the wrapper
function cjamp_loop() provides an easy use of the
cjamp() function and only extracts the coefficient and
standard error estimates as well as the p-values for the predictor(s) of
interest, and not for all other covariates.
# Restrict example to sample size 100 to decrease running time:
covariates <- data.frame(X1 = phenodata$X1, X2 = phenodata$X2)[1:100,]
predictors <- genodata[1:100,1:5]
cjamp_loop_res <- cjamp_loop(copula = "Clayton", Y1 = phenodata$Y1[1:100],
Y2 = phenodata$Y2[1:100], predictors = predictors,
covariates_Y1 = covariates, covariates_Y2 = covariates,
scale_var = FALSE, optim_method = "BFGS", trace = 0,
kkt2tol = 1E-16, SE_est = TRUE, pval_est = TRUE,
n_iter_max = 10)
cjamp_loop_res## $`Parameter point estimates for effects of predictors on Y1`
## SNV1 SNV2 SNV3 SNV4 SNV5
## 0.03839141 0.42923016 -0.11439386 0.14760786 0.16346598
##
## $`Parameter point estimates for effects of predictors on Y2`
## SNV1 SNV2 SNV3 SNV4 SNV5
## 0.15063307 0.30906230 -0.11098892 0.08360284 0.40747355
##
## $`Parameter standard error estimates for effects on Y1`
## SNV1 SNV2 SNV3 SNV4 SNV5
## 0.1267177 0.4757368 0.1561781 0.1316184 0.2786912
##
## $`Parameter standard error estimates for effects on Y2`
## SNV1 SNV2 SNV3 SNV4 SNV5
## 0.1258775 0.4362601 0.1477147 0.1239290 0.2591875
##
## $`Parameter p-values for effects of predictors on Y1`
## SNV1 SNV2 SNV3 SNV4 SNV5
## 0.7619142 0.3669279 0.4638893 0.2620823 0.5575067
##
## $`Parameter p-values for effects of predictors on Y2`
## SNV1 SNV2 SNV3 SNV4 SNV5
## 0.2314374 0.4786746 0.4524279 0.4999282 0.1159230
##
## $`Convergence code of optimx function`
## [1] 0 0 0 0 0
##
## $`Karush-Kuhn-Tucker conditions 1 and 2`
## [1] TRUE TRUE TRUE TRUE TRUE
##
## $`Maximum log-likelihood`
## [1] 222.8908 223.7511 223.8642 223.5241 222.7414
##
## attr(,"class")
## [1] "cjamp"Both cjamp() and cjamp_loop() return
cjamp objects as output, so that the
summary.cjamp() function can be used through the generic
summary() to provide a reader-friendly formatted output of
the results.
## [1] "C-JAMP estimates of marginal parameters."
## point_estimates SE_estimates pvalues
## Y1_(Intercept) 0.08306169 0.14976180 5.791510e-01
## Y1_X1 0.50629518 0.09033346 2.085729e-08
## Y1_X2 0.56417136 0.17038095 9.288517e-04
## Y1_SNV1 0.03707410 0.12632144 7.691475e-01
## Y1_SNV2 0.38545408 0.47396256 4.160699e-01
## Y1_SNV3 -0.09356935 0.15539844 5.470907e-01
## Y2_(Intercept) -0.05692639 0.15521173 7.137935e-01
## Y2_X1 0.37782751 0.08581321 1.068104e-05
## Y2_X2 0.49967186 0.16681705 2.741495e-03
## Y2_SNV1 0.14576130 0.12576460 2.464558e-01
## Y2_SNV2 0.30593863 0.43807946 4.849500e-01
## Y2_SNV3 -0.08816050 0.14983910 5.562855e-01## [1] "C-JAMP estimates of marginal parameters on Y1."
## point_estimates SE_estimates pvalues
## Y1_SNV1 0.03839141 0.1267177 0.7619142
## Y1_SNV2 0.42923016 0.4757368 0.3669279
## Y1_SNV3 -0.11439386 0.1561781 0.4638893
## Y1_SNV4 0.14760786 0.1316184 0.2620823
## Y1_SNV5 0.16346598 0.2786912 0.5575067
## [1] "C-JAMP estimates of marginal parameters on Y2."
## point_estimates SE_estimates pvalues
## Y2_SNV1 0.15063307 0.1258775 0.2314374
## Y2_SNV2 0.30906230 0.4362601 0.4786746
## Y2_SNV3 -0.11098892 0.1477147 0.4524279
## Y2_SNV4 0.08360284 0.1239290 0.4999282
## Y2_SNV5 0.40747355 0.2591875 0.1159230References
Joe H (1997). Multivariate models and multivariate dependence concepts. London: Chapman & Hall.
Konigorski S, Yilmaz YE, Bull SB (2014). Bivariate genetic association analysis of systolic and diastolic blood pressure by copula models. BMC Proc, 8(Suppl 1): S72.
Konigorski S, Yilmaz YE, Pischon T (2016). Genetic association analysis based on a joint model of gene expression and blood pressure. BMC Proc, 10(Suppl 7): 289-294.
Laird NM, Lange C (2011). The fundamentals of modern statistical genetics. New York: Springer.
Nelsen RB (2006). An introduction to copulas. Springer, New York.
Sklar A (1959). Fonctions de répartition à n dimensions et leurs marges. Publications de l’Institut de statistique de l’Université de Paris 8: 229–231.
Self GS, Liang KY (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. J Am Stat Assoc, 82: 605–610.
Yilmaz YE, Lawless JF (2011). Likelihood ratio procedures and tests of fit in parametric and semiparametric copula models with censored data. Lifetime Data Anal, 17(3): 386-408.