Robust Joint Analysis with Data Fusion in Two-Stage Quantitative Trait Genome-Wide Association Studies

Abstract

Genome-wide association studies (GWASs) in identifying the disease-associated genetic variants have been proved to be a great pioneering work. Two-stage design and analysis are often adopted in GWASs. Considering the genetic model uncertainty, many robust procedures have been proposed and applied in GWASs. However, the existing approaches mostly focused on binary traits, and few work has been done on continuous (quantitative) traits, since the statistical significance of these robust tests is difficult to calculate. In this paper, we develop a powerful 𝐹-statistic-based robust joint analysis method for quantitative traits using the combined raw data from both stages in the framework of two-staged GWASs. Explicit expressions are obtained to calculate the statistical significance and power. We show using simulations that the proposed method is substantially more robust than the 𝐹-test based on the additive model when the underlying genetic model is unknown. An example for rheumatic arthritis (RA) is used for illustration.

1. Introduction

Genome-wide association studies (GWASs) have identified a large number of genomic regions (especially single-nucleotide polymorphisms (SNPs)) with a wide variety of complex traits/diseases. In a GWAS, two most common types of data, qualitative (or binary) and quantitative (or continuous) traits, are analyzed and two contentious points are often faced; one is how to construct the test statistic considering the genetic model uncertainty and the other is how to evaluate the statistical significance for controlling the false positive rates efficiently (e.g., [1, 2]). Considering these issues, a lot of work has been done on the binary trait in the past 10 years (e.g., [3–7]). Computer algorithms have also been developed to calculated the significance level of robust tests in GWASs, taking into account the genetic model uncertainty [8]. However, few work has been done on continuous traits, only recently So and Sham [9] proposed a MAX3 based on score test statistics, and Li et al. [10] gave a MAX3 based on 𝐹-test statistics. Note that these tests just focus on single-marker analysis in one-stage analysis.

Although the costs of whole-genome genotyping are decreasing with the high-throughput biological technology, the total costs for a GWAS are still very expensive due to the thousands of sampling units and huge amounts of single-nucleotide polymorphisms. In order to save the costs, the two-stage design and the corresponding statistical analysis where all the SNPs are genotyped in Stage 1 on a portion of the samples and the promising SNPs with small 𝑃-values (e.g., <0.001) based on some efficient tests are further screened on the remaining subjects, are often adopted in practice (e.g., [11–15]).

In genetic association studies, especially GWASs, genetic markers are routinely tested under the assumption of additive effects. Although convenient to use, those tests are optimal only when the true underlying genetic model is additive so that they are not robust against the genetic model misspecification. To our best knowledge, few work has been done on the two-stage joint analysis for quantitative trait GWASs allowing for genetic model uncertainty. Here, we attempt to develop a joint analysis method with data fusion in the two-stage design using 𝐹-statistic, since 𝐹-test is commonly employed from the linear regression model for quantitative trait, and Li et al. [10] show that MAX3 based on 𝐹-statistics is more powerful than So and Sham’s method by extensively numerical simulation.

The content of this paper is organized as follows. In Section 2, we give some notations and the proposed robust joint test statistics. Further, we derive the asymptotic distribution of the test statistics under the null and the alternative hypotheses. In Section 3, we show that the proposed joint analysis method is substantially more robust than the additive-model-based 𝐹-test from the numerical results of power comparison when the real genetic model is unknown. After that, an illustrative example for rheumatic arthritis (RA) is presented. Finally, we give some discussion of this paper in Section 4.

2. Methods

2.1. Notations

Assume that 𝑛 individuals are randomly selected to be genotyped in a two-staged GWAS for a certain quantitative trait and that 𝜋 is the sampling proportion in Stage 1. Let 𝑛₁=𝑛𝜋 and 𝑛₂=𝑛(1−𝜋) be the sample sizes for Stages 1 and 2, respectively. Consider a biallelic marker with two alleles G and g. Without loss of generality, we assume that G is the minor or high-risk allele. We suppose that the total 𝑚 SNPs are genotyped on the samples of Stage 1, and SNPs with 𝑃-values less than 𝛾 in Stage 1 will be further genotyped and tested in Stage 2. Let the significance level be 𝛼, and then the genome-wide significance level per SNP is 𝛼/𝑚 with the Bonferroni adjustments. Let 𝐘₁=(𝑦₁,𝑦₂,…,𝑦_𝑛1)^′ and 𝐘₂=(𝑦_𝑛1+1,𝑦_𝑛1+2,…,𝑦_𝑛)^′ be the observed quantitative outcome vectors for Stage 1 and Stage 2, respectively. Without loss of generality, we assume that the first 𝑛₁₀ individuals in Stage 1 have the genotype gg, the second 𝑛₁₁ individuals in Stage 1 have the genotype Gg, and the last 𝑛₁₂ subjects in Stage 1 possess the genotype GG. Similarly, the first 𝑛₂₀ subjects in Stage 2 have the genotype gg, the second 𝑛₂₁ individuals in Stage 2 have the genotype Gg, and the last 𝑛₂₂ subjects in Stage 2 possess the genotype GG. Let 𝟎_𝑘=(0,0,…,0)^′_𝑘×1 and 𝟏_𝑘=(1,1,…,1)^′_𝑘×1, and let 𝐎_𝑘×𝑗 be the 𝑘×𝑗 matrix with all its entries being zero and 𝐈_𝑛 be the 𝑛×𝑛 identity matrix.

2.2. 𝐹-Statistic-Based Robust Joint Analysis

We firstly briefly introduce 𝐹-statistic-based MAX3 by Li et al. [10] just using the data from Stage 1. Consider the following linear regression model: 𝑦_𝑖=𝛽₀+𝑔_𝑖𝛽₁+𝜀_𝑖,𝜀_𝑖∼𝑁0,𝜎²,𝑖=1,2,…,𝑛₁, (1) where 𝛽₀ is the nuisance parameter for the intercept, 𝛽₁ is the parameter of interest for genetic effect, and 𝑔_𝑖 is the genotype value, which takes 0, 1, or 2 corresponding to the count of G at a marker locus for the 𝑖th subject, 𝑖=1,2,…,𝑛₁. The hypotheses of interest are 𝐻₀∶𝛽₁=0⟷𝐻₁∶𝛽₁≠0. (2) The variable 𝑔_𝑖 in the previously stated equation is coded differently for the three common genetic models. Let 𝐗_1𝑅=(𝟏_𝑛1,𝐆_1𝑅), 𝐗_1𝐴=(𝟏_𝑛1,𝐆_1𝐴), and 𝐗_1𝐷=(𝟏_𝑛1,𝐆_1𝐷) be the design matrices under three commonly used genetic models, where 𝐆_1𝑅=(𝟎^′_{𝑛10+𝑛11},𝟏^′_𝑛12)^′ corresponds to the recessive model, 𝐆_1𝐴=(𝑔₁,𝑔₂,…,𝑔_𝑛1)^′ corresponds to the additive model, and 𝐆_1𝐷=(𝟎^′_𝑛10,𝟏^′_{𝑛11+𝑛12})^′ is for the dominant model. Denote 𝐗₁=(𝟏_𝑛1,𝐱₁₁,𝐱₁₂), where 𝐱₁₁=(𝟎^′_𝑛10,𝟏^′_𝑛11,𝟎^′_𝑛12)^′ and 𝐱₁₂=(𝟎^′_𝑛10,𝟎^′_𝑛11,𝟏^′_𝑛12)^′. The modified 𝐹-test statistics under the recessive, additive, and dominant models for Stage 1 are given by 𝐹^𝑅₁=𝐘^′₁𝐗_1𝑅𝐗^′_1𝑅𝐗_1𝑅⁻¹𝐗^′_1𝑅−𝟏_𝑛1𝟏^′_𝑛1𝟏_𝑛1⁻¹𝟏^′_𝑛1𝐘₁𝐘^′₁𝐈_𝑛1−𝐗₁𝐗^′₁𝐗₁⁻¹𝐗^′₁𝐘₁/𝑛₁=𝑍−3_𝑅12RSS₁/𝑛₁,𝐹−3_𝐴1=𝐘^′₁𝐗_1𝐴𝐗^′_1𝐴𝐗_1𝐴⁻¹𝐗^′_1𝐴−𝟏_𝑛1𝟏^′_𝑛1𝟏_𝑛1⁻¹𝟏^′_𝑛1𝐘₁𝐘^′₁𝐈_𝑛1−𝐗₁𝐗^′₁𝐗₁⁻¹𝐗^′₁𝐘₁/𝑛₁=𝑍−3_𝐴12RSS₁/𝑛₁,𝐹−3_𝐷1=𝐘^′₁𝐗_1𝐷𝐗^′_1𝐷𝐗_1𝐷⁻¹𝐗^′_1𝐷−𝟏_𝑛1𝟏^′_𝑛1𝟏_𝑛1⁻¹𝟏^′_𝑛1𝐘₁𝐘^′₁𝐈_𝑛1−𝐗₁𝐗^′₁𝐗₁⁻¹𝐗^′₁𝐘₁/𝑛₁=𝑍−3_𝐷12RSS₁/𝑛₁,−3 (3) where 𝑍^𝑅₁=𝑛₁₀+𝑛₁₁𝑛₁₂𝑛₁𝑦_{𝑛10+𝑛11}−𝑦_𝑛12,𝑍^𝐷₁=𝑛₁₀𝑛₁₁+𝑛₁₂𝑛₁𝑦_𝑛10−𝑦_{𝑛11+𝑛12},𝑍^𝐴₁=𝑛₁₀𝑛₁₁+2𝑛₁₂𝑦_𝑛10−𝑛₁₁𝑛₁₀−𝑛₁₂𝑦_𝑛11−𝑛₁₂2𝑛₁₀+𝑛₁₁𝑦_𝑛12×𝑛₁𝑛₁₀𝑛₁₁+4𝑛₁₂+𝑛₁₁𝑛₁₂⁻¹,𝑦_𝑛10=1𝑛^𝑛₁₀¹⁰_𝑗=1𝑦_𝑗,𝑦_𝑛11=1𝑛^𝑛₁₁^10+𝑛11_{𝑗=𝑛10+1}𝑦_𝑗,𝑦_𝑛12=1𝑛^𝑛₁₂¹_{𝑗=𝑛10+𝑛11+1}𝑦_𝑗,𝑦_{𝑛10+𝑛11}=1𝑛₁₀+𝑛^𝑛₁₁^10+𝑛11_𝑗=1𝑦_𝑗,𝑦_{𝑛11+𝑛12}=1𝑛₁₁+𝑛^𝑛₁₂¹_{𝑗=𝑛10+1}𝑦_𝑗. (4) The robust test statistic in Stage 1 is 𝐹^MAX₁=max𝐹^𝑅₁,𝐹^𝐴₁,𝐹^𝐷₁. (5)

We now give the proposed robust joint analysis. In the framework of two-stage design GWAS of quantitative traits, the SNPs with 𝑃-values less than 𝛾 will be genotyped on the remaining 𝑛₂ subjects in Stage 2. Following the previous notation for Stage 1, corresponding to the recessive, additive, and dominant models, the genotype data in Stage 2 are denoted by 𝐆_2𝑅=(𝟎^′_{𝑛20+𝑛21},𝟏^′_𝑛22)^′, 𝐆_2𝐴=(𝑔_𝑛1+1,𝑔_𝑛1+2,…,𝑔_𝑛)^′, and 𝐆_2𝐷=(𝟎^′_𝑛20,𝟏^′_{𝑛21+𝑛22})^′, respectively, and the design matrices are 𝐗_2𝑅=(𝟏_𝑛2,𝐆_2𝑅), 𝐗_2𝐴=(𝟏_𝑛2,𝐆_2𝐴), and 𝐗_2𝐷=(𝟏_𝑛2,𝐆_2𝐷), respectively. Denote 𝐗₂=(𝟏_𝑛2,𝐱₂₁,𝐱₂₂), where 𝐱₂₁=(𝟎^′_𝑛20,𝟏^′_𝑛21,𝟎^′_𝑛22)^′ and 𝐱₂₂=(𝟎^′_𝑛20,𝟎^′_𝑛21,𝟏^′_𝑛22)^′. Then, we can obtain three modified 𝐹-test statistics under the recessive, additive, and dominant models for Stage 2 similarly, and denote them by 𝐹^𝑅₂, 𝐹^𝐴₂, and 𝐹^𝐷₂. Let 𝐘=(𝐘^′₁,𝐘^′₂)^′, 𝐆_𝑅=(𝐆^′_1𝑅,𝐆^′_2𝑅)^′, 𝐆_𝐴=(𝐆^′_1𝐴,𝐆^′_2𝐴)^′, and 𝐆_𝐷=(𝐆^′_1𝐷,𝐆^′_2𝐷)^′. Denote 𝑁₀=𝑛₁₀+𝑛₂₀, 𝑁₁=𝑛₁₁+𝑛₂₁, and 𝑁₂=𝑛₁₂+𝑛₂₂ for the combined sample sizes from two stages, corresponding to three genotypes. Then the proposed 𝐹-test statistics under three genetic models on the basis of the combined data are as follows: 𝐹^𝑅_𝐽=𝐘^′𝐗_𝑅𝐗^′_𝑅𝐗_𝑅⁻¹𝐗^′_𝑅−𝟏_𝑛𝟏^′_𝑛𝟏_𝑛⁻¹𝟏^′_𝑛𝐘𝐘^′𝐈_𝑛−𝐖𝐖^′𝐖⁻¹𝐖^′()=𝑍𝐘/𝑛−6_𝑅𝐽2RSS_𝐽/(),𝐹𝑛−6_𝐴𝐽=𝐘^′𝐗_𝐴𝐗^′_𝐴𝐗_𝐴⁻¹𝐗^′_𝐴−𝟏_𝑛𝟏^′_𝑛𝟏_𝑛⁻¹𝟏^′_𝑛𝐘𝐘^′𝐈_𝑛−𝐖𝐖^′𝐖⁻¹𝐖^′()=𝑍𝐘/𝑛−6_𝐴𝐽2RSS_𝐽/(),𝐹𝑛−6_𝐷𝐽=𝐘^′𝐗_𝐷𝐗^′_𝐷𝐗_𝐷⁻¹𝐗^′_𝐷−𝟏_𝑛𝟏^′_𝑛𝟏_𝑛⁻¹𝟏^′_𝑛𝐘𝐘^′𝐈_𝑛−𝐖𝐖^′𝐖⁻¹𝐖^′()=𝑍𝐘/𝑛−6_𝐷𝐽2RSS_𝐽/(),𝑛−6 (6) where 𝐗_𝑅=(𝟏_𝑛,𝐆_𝑅), 𝐗_𝐴=(𝟏_𝑛,𝐆_𝐴), 𝐗_𝐷=(𝟏_𝑛,𝐆_𝐷), and   𝐖=^{𝐗1𝐎𝑛1×3}_{𝐎𝑛2×3𝐗2}, 𝑍^𝑅_𝐽=𝑁₀+𝑁₁𝑁₂𝑛𝑦₀₁−𝑦₂,𝑍^𝐷_𝐽=𝑁₀𝑁₁+𝑁₂𝑛𝑦₀−𝑦₁₂,𝑍^𝐴_𝐽=𝑁₀𝑁₁+2𝑁₂𝑦₀−𝑁₁𝑁₀−𝑁₂𝑦₁−𝑁₂2𝑁₀+𝑁₁𝑦₂×𝑛𝑁₀𝑁₁+4𝑁₂+𝑁₁𝑁₂⁻¹,𝑦₀=1𝑁₀^𝑛10_𝑗=1𝑦_𝑗+^𝑛1+𝑛20_𝑗=𝑛1+1𝑦_𝑗,𝑦₁=1𝑁₁^{𝑛10+𝑛11}_{𝑗=𝑛10+1}𝑦_𝑗+^{𝑛1+𝑛20+𝑛21}_{𝑗=𝑛1+𝑛20+1}𝑦_𝑗,𝑦₂=1𝑁₂^𝑛1_{𝑗=𝑛10+𝑛11+1}𝑦_𝑗+^𝑛_{𝑗=𝑛1+𝑛20+𝑛21+1}𝑦_𝑗,𝑦₀₁=𝑁₀𝑦₀+𝑁₁𝑦₁𝑁₀+𝑁₁,𝑦₁₂=𝑁₁𝑦₁+𝑁₂𝑦₂𝑁₁+𝑁₂. (7) Furthermore, we propose the joint testing statistic as 𝐹^MAX_𝐽=max𝐹^𝑅_𝐽,𝐹^𝐴_𝐽,𝐹^𝐷_𝐽. (8)

In order to calculate the power of the proposed joint analysis, we have to get the thresholds, which is determined by the significance level. Denote the threshold for choosing the promising SNPs in Stage 1 by 𝑢₁, which is the solution of Pr_𝐻0𝐹^MAX₁>𝑢₁=𝛾. (9) Since the genome-wide significance level is 𝛼/𝑚, in order to control the false positive rate, we have Pr_𝐻0𝐹^MAX₁>𝑢₁,𝐹^MAX_𝐽>𝑢_𝐽=𝛼/𝑚, (10) where 𝑢_𝐽 is the cut-off point for the joint statistic. Once we have 𝑢₁ and 𝑢_𝐽, the power is calculated by Pr_𝐻1𝐹^MAX₁>𝑢₁,𝐹^MAX_𝐽>𝑢_𝐽. (11)

We now give the detail to calculate the cut-off point and power above. The left side of (10) can be further expressed as Pr_𝐻0𝐹^MAX₁>𝑢₁,𝐹^MAX_𝐽>𝑢_𝐽=1−Pr_𝐻0𝐹^MAX₁≤𝑢₁−Pr_𝐻0𝐹^MAX_𝐽≤𝑢_𝐽+Pr_𝐻0𝐹^MAX₁≤𝑢₁,𝐹^MAX_𝐽≤𝑢_𝐽. (12) For controlling the type I error rate and calculating the power, we need to know the distribution or the asymptotic distribution of (𝑍^𝑅₁,𝑍^𝐴₁,𝑍^𝐷₁,𝑍^𝑅_𝐽,𝑍^𝐴_𝐽,𝑍^𝐷_𝐽,RSS₁,RRS_𝐽)^′ under both 𝐻₀ and 𝐻₁.

Note that whether 𝐻₀ or 𝐻₁ holds, RSS₁ and RSS_𝐽 and (𝑍^𝑅₁,𝑍^𝐴₁,𝑍^𝐷₁,𝑍^𝑅_𝐽,𝑍^𝐴_𝐽,𝑍^𝐷_𝐽)^′ are mutually independent (the proof is given in Appendix A). Denote the correlation matrix of (𝑍^𝑅₁,𝑍^𝐴₁,𝑍^𝐷₁)^′ by 𝐕₁=(𝑣_𝑘𝑙)_3×3, whose entries are 𝑣₁₁=𝑣₂₂=𝑣₃₃=1, 𝑣₁₂=𝑣₂₁=√𝑛₁₂(2𝑛₁₀+𝑛₁₁)/(𝑛₁₀+𝑛₁₁)[𝑛₁₀(𝑛₁₁+4𝑛₁₂)+𝑛₁₁𝑛₁₂], 𝑣₁₃=𝑣₃₁=√𝑛₁₀𝑛₁₂/(𝑛₁₀+𝑛₁₁)(𝑛₁₁+𝑛₁₂), and 𝑣₂₃=𝑣₃₂=√𝑛₁₀(2𝑛₁₂+𝑛₁₁)/(𝑛₁₁+𝑛₁₂)[𝑛₁₀(𝑛₁₁+4𝑛₁₂)+𝑛₁₁𝑛₁₂], respectively. Similarly, let 𝐕_𝐽=(𝑣^∗_𝑘𝑙)_3×3 be the correlation matrix of (𝑍^𝑅_𝐽,𝑍^𝐴_𝐽,𝑍^𝐷_𝐽)^′ with 𝑣^∗₁₁=𝑣^∗₂₂=𝑣^∗₃₃=1, 𝑣^∗₁₂=𝑣^∗₂₁=𝑁₂(2𝑁₀+𝑁₁)/(𝑁₀+𝑁₁)[𝑁₀(𝑁₁+4𝑁₂)+𝑁₁𝑁₂], 𝑣^∗₁₃=𝑣^∗₃₁=𝑁₀𝑁₂/(𝑁₀+𝑁₁)(𝑁₁+𝑁₂), and 𝑣^∗₂₃=𝑣^∗₃₂=𝑁₀(2𝑁₂+𝑁₁)/(𝑁₁+𝑁₂)[𝑁₀(𝑁₁+4𝑁₂)+𝑁₁𝑁₂]. Then, we can derive that RRS₁/𝜎²∼𝜒²_𝑛1−3, RSS_𝐽/𝜎²∼𝜒²_𝑛−6, and 𝑍^𝑅₁,𝑍^𝐴₁,𝑍^𝐷₁,𝑍^𝑅_𝐽,𝑍^𝐴_𝐽,𝑍^𝐷_𝐽^′_𝐻0∼𝑁₆𝟎₆,𝜎²𝐕₁𝝆𝝆^′𝐕_𝐽, (13) where 𝝆=(𝜌_𝑘𝑙)_3×3 is the correlation matrix between (𝑍^𝑅₁,𝑍^𝐴₁,𝑍^𝐷₁)^′ and (𝑍^𝑅_𝐽,𝑍^𝐴_𝐽,𝑍^𝐷_𝐽)^′, with 𝜌₁₁=Corr𝑍^𝑅₁,𝑍^𝑅_𝐽=𝑛𝑛₁₀+𝑛₁₁𝑛₁₂𝑛₁𝑁₀+𝑁₁𝑁₂,𝜌₁₂=Corr𝑍^𝑅₁,𝑍^𝐴_𝐽=√𝑛𝑛₁₂2𝑛₁₀+𝑛₁₁𝑛₁𝑛₁₀+𝑛₁₁𝑁₀𝑁₁+4𝑁₂+𝑁₁𝑁₂,𝜌₁₃=Corr𝑍^𝑅₁,𝑍^𝐷_𝐽=√𝑛𝑛₁₂𝑛₁₀𝑛₁𝑛₁₀+𝑛₁₁𝑁₀𝑁₁+𝑁₂,𝜌₂₁=Corr𝑍^𝐴₁,𝑍^𝑅_𝐽=√𝑛𝑛₁₂2𝑛₁₀+𝑛₁₁𝑛₁𝑛₁₀𝑛₁₁+4𝑛₁₂+𝑛₁₁𝑛₁₂𝑁₀+𝑁₁𝑁₂,𝜌₂₂=Corr𝑍^𝐴₁,𝑍^𝐴_𝐽=𝑛𝑛₁₀𝑛₁₁+4𝑛₁₂+𝑛₁₁𝑛₁₂𝑛₁𝑁₀𝑁₁+4𝑁₂+𝑁₁𝑁₂,𝜌₂₃=Corr𝑍^𝐴₁,𝑍^𝐷_𝐽=√𝑛𝑛₁₀𝑛₁₁+2𝑛₁₂𝑛₁𝑛₁₀𝑛₁₁+4𝑛₁₂+𝑛₁₁𝑛₁₂𝑁₀𝑁₁+𝑁₂,𝜌₃₁=Corr𝑍^𝐷₁,𝑍^𝑅_𝐽=√𝑛𝑛₁₀𝑛₁₂𝑛₁𝑛₁₁+𝑛₁₂𝑁₀+𝑁₁𝑁₂,𝜌₃₂=Corr𝑍^𝐷₁,𝑍^𝐴_𝐽=√𝑛𝑛₁₀𝑛₁₁+2𝑛₁₂𝑛₁𝑛₁₁+𝑛₁₂𝑁₀𝑁₁+4𝑁₂+𝑁₁𝑁₂,𝜌₃₃=Corr𝑍^𝐷₁,𝑍^𝐷_𝐽=𝑛𝑛₁₀𝑛₁₁+𝑛₁₂𝑛₁𝑁₀𝑁₁+𝑁₂. (14)

Under 𝐻₁, for a given odds ratio OR=exp(𝛽₁) for subjects with two copies of risk allele corresponding to recessive model or one copy of risk allele corresponding to additive or dominant models, we have the following:(i) when the true genetic model is recessive, 𝑍^𝑅₁,𝑍^𝐴₁,𝑍^𝐷₁,𝑍^𝑅_𝐽,𝑍^𝐴_𝐽,𝑍^𝐷_𝐽^′_𝐻1∼𝑁₆𝝁^𝑅,𝜎²𝐕₁𝝆𝝆^′𝐕_𝐽, (15) where 𝝁^𝑅=(𝜇₁^𝑅𝑅,𝜇₁^𝑅𝐴,𝜇₁^𝑅𝐷,𝜇_𝐽^𝑅𝑅,𝜇_𝐽^𝑅𝐴,𝜇_𝐽^𝑅𝐷)^′ with 𝜇₁^𝑅𝑅=−𝑛₁₀+𝑛₁₁𝑛₁₂𝑛₁𝛽₁,𝜇₁^𝑅𝐴=−𝑛₁₂2𝑛₁₀+𝑛₁₁𝛽₁𝑛₁𝑛₁₀𝑛₁₁+4𝑛₁₂+𝑛₁₁𝑛₁₂,𝜇₁^𝑅𝐷=−√𝑛₁₀𝑛₁₂𝛽₁𝑛₁𝑛₁₁+𝑛₁₂,𝜇_𝐽^𝑅𝑅=−𝑁₀+𝑁₁𝑁₂𝑛𝛽₁,𝜇_𝐽^𝑅𝐴=−𝑁₂2𝑁₀+𝑁₁𝛽₁𝑛𝑁₀𝑁₁+4𝑁₂+𝑁₁𝑁₂,𝜇_𝐽^𝑅𝐷=−𝑁₀𝑁₂𝛽₁𝑛𝑁₁+𝑁₂, (16)(ii) when the true genetic model is additive, 𝑍^𝑅₁,𝑍^𝐴₁,𝑍^𝐷₁,𝑍^𝑅_𝐽,𝑍^𝐴_𝐽,𝑍^𝐷_𝐽^′_𝐻1∼𝑁₆𝝁^𝐴,𝜎²𝐕₁𝝆𝝆^′𝐕_𝐽, (17) where 𝝁^𝐴=(𝜇₁^𝐴𝑅,𝜇₁^𝐴𝐴,𝜇₁^𝐴𝐷,𝜇_𝐽^𝐴𝑅,𝜇_𝐽^𝐴𝐴,𝜇_𝐽^𝐴𝐷)^′ with 𝜇₁^𝐴𝑅=−√𝑛₁₂𝑛₁₁+2𝑛₁₀𝛽₁𝑛₁𝑛₁₀+𝑛₁₁,𝜇₁^𝐴𝐴=−𝑛₁₀𝑛₁₁+4𝑛₁₂+𝑛₁₁𝑛₁₂𝑛₁𝛽₁,𝜇₁^𝐴𝐷=−√𝑛₁₀𝑛₁₁+2𝑛₁₂𝛽₁𝑛₁𝑛₁₁+𝑛₁₂,𝜇_𝐽^𝐴𝑅=−𝑁₂𝑁₁+2𝑁₀𝛽₁𝑛𝑁₀+𝑁₁,𝜇_𝐽^𝐴𝐴=−𝑁₀𝑁₁+4𝑁₂+𝑁₁𝑁₂𝑛𝛽₁,𝜇_𝐽^𝐴𝐷=−𝑁₀𝑁₁+2𝑁₂𝛽₁𝑛𝑁₁+𝑁₂, (18)(iii) when the true genetic model is dominant, 𝑍^𝑅₁,𝑍^𝐴₁,𝑍^𝐷₁,𝑍^𝑅_𝐽,𝑍^𝐴_𝐽,𝑍^𝐷_𝐽^′_𝐻1∼𝑁₆𝝁^𝐷,𝜎²𝐕₁𝝆𝝆^′𝐕_𝐽, (19) where 𝝁^𝐷=(𝜇₁^𝐷𝑅,𝜇₁^𝐷𝐴,𝜇₁^𝐷𝐷,𝜇_𝐽^𝐷𝑅,𝜇_𝐽^𝐷𝐴,𝜇_𝐽^𝐷𝐷)^′ with 𝜇₁^𝐷𝑅=−√𝑛₁₂𝑛₁₀𝛽₁𝑛₁𝑛₁₀+𝑛₁₁,𝜇₁^𝐷𝐴=−𝑛₁₀𝑛₁₁+2𝑛₁₂𝛽₁𝑛₁𝑛₁₀𝑛₁₁+4𝑛₁₂+𝑛₁₁𝑛₁₂,𝜇₁^𝐷𝐷=−𝑛₁₀𝑛₁₁+𝑛₁₂𝑛₁𝛽₁,𝜇_𝐽^𝐷𝑅=−𝑁₂𝑁₀𝛽₁𝑛𝑁₀+𝑁₁,𝜇_𝐽^𝐷𝐴=−𝑁₀𝑁₁+2𝑁₂𝛽₁𝑛𝑁₀𝑁₁+4𝑁₂+𝑁₁𝑁₂,𝜇_𝐽^𝐷𝐷=−𝑁₀𝑁₁+𝑁₂𝑛𝛽₁. (20)

We develop a method for simplifying the calculations of Pr_𝐻0(𝐹^MAX₁≤𝑢₁) and Pr_𝐻0(𝐹^MAX_𝐽≤𝑢_𝐽) and Pr_𝐻0(𝐹^MAX₁≤𝑢₁,𝐹^MAX_𝐽≤𝑢_𝐽). The details are included in Appendix B, and the calculations of Pr_𝐻1(𝐹^MAX₁≤𝑢₁) and Pr_𝐻1(𝐹^MAX_𝐽≤𝑢_𝐽) and Pr_𝐻1(𝐹^MAX₁≤𝑢₁,𝐹^MAX_𝐽≤𝑢_𝐽) are essentially similar.

3. Results

3.1. Power Comparison

We conduct simulation studies to evaluate the performance of the proposed method under three commonly used genetic models (recessive, additive, and dominant models). We mainly compare the power of two approaches; one is the proposed method in this paper, and the other is the joint analysis based on the 𝐹-test statistics 𝐹^𝐴₁ and 𝐹^𝐴_𝐽. For convenience, we refer to the proposed method as MAXFJ and AFJ for the other one. We choose the sample size 𝑛=2000, and 𝑚=5×10⁵. The proportion of subjects genotyped in Stage 1 has three levels 𝜋=0.3,0.4,0.5. We set the genome-wide significance level as 𝛼=0.05 and that the significance level per SNP as 𝛼/𝑚=1×10⁻⁷. In Stage 1, the 𝑃-value threshold for SNPs selected for followup is set to be 1×10⁻⁴ and 2×10⁻⁴. We assume that the Hardy-Weinberg equilibrium holds in the general sample population, and then there are on average 𝑛×(1−MAF)², 2𝑛×MAF×(1−MAF), and 𝑛×MAF² individuals with genotype gg, Gg, and GG, respectively, where the minor allele frequency is set to be 0.15, 0.30 and 0.45. To make the power comparison more distinctly, we specify different genetic effect parameters 𝛽₁ under three genetic models as follows: 𝛽₁=0.5 for the recessive model, 𝛽₁=0.3 for the additive model, and 𝛽₁=0.4 for the dominant model.

The power results are displayed in Tables 1 and 2 for 𝛾=1×10⁻⁴ and 𝛾=2×10⁻⁴, respectively. They indicate that MAXFJ is more efficiency robust than AFJ across various inheritance models. As expected, AFJ is more powerful than MAXFJ under the additive model. However, MAFJ performs much more powerful than AFJ when the true genetic model is recessive. For instance, in Table 2, with 𝜋=0.4 and MAF = 0.3, the powers of AFJ and MAXFJ are 0.101 and 0.529, respectively. In summary, MAXFJ is substantially more powerful than AFJ in two-staged GWAS of quantitative traits, when the model for AFJ is misspecified.

Table 1 

Power comparison (𝑛=2000,𝛾=1×10⁻⁴,𝛼=0.05,  and 𝑚=5×10⁵).

Table 2 

Power comparison (𝑛=2000,𝛾=2×10⁻⁴,𝛼=0.05,  and 𝑚=5×10⁵).

3.2. An Illustration Example: Rheumatoid Arthritis

Rheumatoid arthritis (RA) is an autoimmune disease (resulting in a chronically systemic inflammatory disorder) which mainly attacks synovial joints. About 1% of the common adult population worldwide is affected by RA [16]. It has been pointed out that the genetic variants might play a major role in RA susceptibility [17]. Genetic Analysis Workshop 16 (GAW16) based on the North American Rheumatoid Arthritis Consortium (NARAC) is a GWAS testing association with RA using about 5×10⁵ SNPs [18–20]. It included 868 individuals who were RA positive (cases) and also had continuous trait anticyclic citrullinated peptide (anti-CCP) measures and 1194 controls sampled from the New York Cancer Project (NYCP) without RA which had no anti-CCP measures. Huizinga et al. [21] pointed out that a greater anti-CCP would be linked to better prediction of increased risk developing RA. Chen et al. [22] showed that SNP rs2476601 located in PTPN22 had the most significant association with RA. Here, we only focus on SNP rs2476601 and apply two joint analysis methods (AFJ and MAXFJ) to evaluate its statistical significance. The minimum of anti-CCP among 868 cases was affected to each control, and a log transformation of anti-CCP was applied in the analysis. Then, we considered 𝜋=0.3,0.4,0.5 three simulation circumstances. For 𝜋=0.3, thirty percent of individuals were randomly sampled from all cases and controls and were used as the data from Stage 1, and the rest of individuals were treated as the data of Stage 2. The 𝑃-values of AFJ and MAXFJ were calculated, respectively. We repeated the above procedure 1,000 times and saved the corresponding 𝑃-values. A base-10 logarithm transformation and an opposite transformation were successively applied to these 𝑃-values, and the histogram and density of these transformed data were obtained (Figure 1). Similarly, we conducted the simulation and calculation for 𝜋=0.4 and 0.5, and the corresponding histogram and density were presented in Figures 2 and 3. Examination of Figures 1–3 showed that the 𝑃-values of MAXFJ are more stable than those of AFJ and the estimated density curves of MAXFJ are more closer to the symmetrical normal distribution while the estimated density curves of AFJ are rather skewed, which indicated that MAXFJ possesses more robust performance when the real genetic models are unknown.

(a)

(b)

Figure 1 

The histogram and density of −log10𝑃 when 𝜋=0.3 (the left subgraph corresponds to MAXFJ while the right one for AFJ).

(a)

(b)

Figure 2 

The histogram and density of −log10𝑃 when 𝜋=0.4 (the left subgraph corresponds to MAXFJ while the right one for AFJ).

(a)

(b)

Figure 3 

The histogram and density of −log10𝑃 when 𝜋=0.5 (the left subgraph corresponds to MAXFJ while the right one for AFJ).

4. Discussion

We have developed a feasible two-stage design and the corresponding robust joint analysis approach for quantitative trait GWASs. The method is based on the 𝐹-statistics over three different genetic models. The denominator of the used 𝐹-statistic, which is constructed without assuming any genetic model, is different from the commonly used one. This adoption can reduce the computation intensity. Taking advantage of an ingenious design matrix, we successfully construct the common denominator of three 𝐹-test statistics for the joint analysis with combined raw data from both stages. The statistical significance (𝑃-value) for the proposed joint analysis method can be calculated with the derived analytic expressions on the basis of the asymptotic distributions, which greatly reduce the complexity and computational intensity compared with the resampling-type permutation and bootstrap procedures. Our numerical results demonstrate that this novel approach has the greater efficiency robustness for genetic model uncertainty than the 𝐹-statistic-based joint analysis which assumes the additive genetic model.

In this work, we did not investigate the power of joint analysis based on other existing robust association methods for quantitative traits such as So and Sham’s method. We find that it is very difficult to extend So and Sham’s method (score test-based MAX3) to two-staged GWASs with quantitative outcomes, since it is almost impossible to derive the joint distribution of score tests from two stages.

For simplicity, here we do not take into account the effects of covariates in the considered two-stage design. However, in real application, the proposed method can be easily applied to the situation including one or more covariates as shown by the original MAXF by Li et al. [10]. It is important to stress that we combine the raw data from two stages to construct the joint statistic, unlike the joint analysis for binary traits using the weighted sum of two statistics in Stages 1 and 2 [12]. Furthermore, one basic assumption in this paper is that the effect sizes of genetic variants between two stages are identical (i.e., no heterogeneity exists), which is the natural and reasonable precondition for the data fusion strategy. In addition, the population-based genetic association studies may be affected by the population stratification, and this needs future research to examine it.

Appendix

A. The Derivation of the Asymptotic Properties of 𝐹^𝑅_𝐽,𝐹^𝐴_𝐽,𝐹^𝐷_𝐽

Consider the linear model for the combined raw data from Stage 1 and Stage 2 as follows: 𝐘=𝐖𝝑+𝜺, (A.1) where 𝝑=(𝛽₀,𝜁₁,𝜁₂,𝛽^∗₀,𝜁^∗₁,𝜁^∗₂)^′ and 𝜺∼𝑁_𝑛(𝟎_𝑛,𝜎²𝐈_𝑛).

Denote 𝐂_𝑅=,𝐂100−1000100−1000100−1010000_𝐷=,𝐂100−1000100−1000100−101−1000_𝐴=,𝐂100−1000100−1000100−102−1000₀=.100−1000100−1000100−1010000001000 (A.2)

Based on the design matrix 𝐗𝐖=¹𝐎_𝑛1×3𝐎_𝑛2×3𝐗₂ (A.3) for the expanded full model above, we can get the ordinary least square estimator of 𝝑 by 𝝑=(𝐖^′𝐖)⁻¹𝐖^′𝐘, and the residual sum of square is given by RSS_𝐽=𝐘^′[𝐈_𝑛−𝐖(𝐖^′𝐖)⁻¹𝐖^′]𝐘. Furthermore, we denote the residual sum of squares under the following constraints: 𝐂_𝑅𝝑=𝟎₄, 𝐂_𝐷𝝑=𝟎₄, 𝐂_𝐴𝝑=𝟎₄, and 𝐂₀𝝑=𝟎₅, by RSS_𝑅,RSS_𝐷,RSS_𝐴, and RSS₀, respectively. After some algebras, we can obtain 𝐹^𝑅_𝐽=RSS₀−RSS_𝑅RSS_𝐽/()=𝑛−6RSS₀−RSS_𝐽−RSS_𝑅−RSS_𝐽RSS_𝐽/(),𝐹𝑛−6_𝐴𝐽=RSS₀−RSS_𝐴RSS_𝐽/()=𝑛−6RSS₀−RSS_𝐽−RSS_𝐴−RSS_𝐽RSS_𝐽/(),𝐹𝑛−6_𝐷𝐽=RSS₀−RSS_𝐷RSS_𝐽/()=𝑛−6RSS₀−RSS_𝐽−RSS_𝐷−RSS_𝐽RSS_𝐽/().𝑛−6 (A.4)

On the one hand, according to 𝐖^′𝐗𝐖=^′1𝐗₁𝐎_3×3𝐎_3×3𝐗^′₂𝐗₂, (A.5) it follows that 𝐖^′𝐖⁻¹=𝐗^′₁𝐗₁⁻¹𝐎_3×3𝐎_3×3𝐗^′₂𝐗₂⁻¹,𝐖𝐖^′𝐖⁻¹𝐖^′=𝐗₁𝐗^′₁𝐗₁⁻¹𝐗^′₁𝐎_𝑛1×𝑛2𝐎_𝑛2×𝑛1𝐗₂𝐗^′₂𝐗₂⁻¹𝐗^′₂. (A.6) So, we can get that 𝐘^′𝐈_𝑛−𝐖𝐖^′𝐖⁻¹𝐖^′𝐘=𝐘^′₁𝐈_𝑛1−𝐗₁𝐗^′₁𝐗₁⁻¹𝐗^′₁𝐘₁+𝐘^′₂𝐈_𝑛2−𝐗₂𝐗^′₂𝐗₂⁻¹𝐗^′₂𝐘₂. (A.7) That is, RSS_𝐽=RSS₁+RSS₂. Furthermore, RSS_𝐽/(𝑛−6) is also the unbiased estimator of the variance of the residual 𝜎² based on the independence and unbiasedness of RSS₁ and RSS₂.

On the other hand, RSS₀−RSS_𝐽 and RSS_𝑅−RSS_𝐽 are both independent of RSS_𝐽, so (𝑍^𝑅_𝐽)²=RSS₀−RSS_𝑅=(RSS₀−RSS_𝐽)−(RSS_𝑅−RSS_𝐽) is also independent of RSS_𝐽, and (RSS₀−RSS_𝑅)/𝜎²∼𝜒²₁, RSS_𝐽/𝜎²∼𝜒²_𝑛−6. Consequently, we have 𝐹^𝑅_𝐽∼𝐹_1,𝑛−6. Similarly, we can get 𝐹^𝐴_𝐽∼𝐹_1,𝑛−6 and 𝐹^𝐷_𝐽∼𝐹_1,𝑛−6.

B. The Detailed Calculation of Pr_𝐻0(𝐹^MAX₁≤𝑢₁) and Pr_𝐻0(𝐹^MAX_𝐽≤𝑢_𝐽) and Pr_𝐻0(𝐹^MAX₁≤𝑢₁,𝐹^MAX_𝐽≤𝑢_𝐽)

Denote 𝑤₀=(𝑛₁₀+𝑛₁₁)𝑛₁₂/(𝑛₁₀(𝑛₁₁+4𝑛₁₂)+𝑛₁₁𝑛₁₂) and 𝑤₁=𝑛₁₀(𝑛₁₁+𝑛₁₂)/(𝑛₁₀(𝑛₁₁+4𝑛₁₂)+𝑛₁₁𝑛₁₂). For a given 𝑐>0, Pr_𝐻0𝑍^𝑅₁𝜎𝑍≤𝑐,^𝐴1𝜎𝑍≤𝑐,^𝐷1𝜎=≤𝑐_Ω0𝑓𝑧₀,𝑧₁;𝚺₀𝑑𝑧₀𝑑𝑧₁, (B.1) where Ω₀={(𝑧₀,𝑧₁)∶|𝑧₀|≤𝑐,|𝑤₀𝑧₀+𝑤₁𝑧₁|≤𝑐,|𝑧₁|≤𝑐} and 𝑓(𝑧₀,𝑧₁;𝚺₀) is the bivariate normal density function for (𝑍^𝑅₁/𝜎,𝑍^𝐷₁/𝜎)^′ with mean 𝟎₂ and variance-covariance matrix 𝚺₀=^1𝑣13_𝑣131. Taking advantage of the symmetry of bivariate normal distribution, the above twofold integration can be only calculated at the right half space of Ω₀ and then multiplied by 2, which is _Ω0𝑓𝑧₀,𝑧₁;𝚺₀𝑑𝑧₀𝑑𝑧₁=2^{𝑐(1−𝑤1)/𝑤0}₀𝑑𝑧₀^𝑐_−𝑐𝑓𝑧₀,𝑧₁;𝚺₀𝑑𝑧₁+^𝑐_{𝑐1−𝑤1/𝑤0}𝑑𝑧₀^{(𝑐−𝑤0𝑧0)/𝑤1}_−𝑐𝑓𝑧₀,𝑧₁;𝚺₀𝑑𝑧₁. (B.2)

Based on the property of conditional distributions of the multivariate normal distribution, we have 𝑓𝑧₀,𝑧₁;𝚺₀=𝜙𝑧₀𝑓𝑧₁∣𝑧₀;𝑣₁₃, (B.3) where 𝜙(𝑧₀) is the probability density function of 𝑁(0,1) and 𝑓(𝑧₁∣𝑧₀;𝑣₁₃) is the density function of the conditional normal distribution 𝑁(𝑣₁₃𝑧₀,1−𝑣²₁₃). That is, 𝑓𝑧₁∣𝑧₀;𝑣₁₃=11−𝑣²₁₃𝜙𝑧₁−𝑣₁₃𝑧₀1−𝑣²₁₃. (B.4) Then, it follows that ^{𝑐(1−𝑤1)/𝑤0}₀𝑑𝑧₀^𝑐_−𝑐𝑓𝑧₀,𝑧₁;𝚺₀𝑑𝑧₁=^{𝑐(1−𝑤1)/𝑤0}₀Φ𝑐−𝑣₁₃𝑧₀1−𝑣²₁₃−Φ−𝑐−𝑣₁₃𝑧₀1−𝑣²₁₃×𝜙𝑧₀𝑑𝑧₀,_{𝑐𝑐1−𝑤1/𝑤0}𝑑𝑧₀^{(𝑐−𝑤0𝑧0)/𝑤1}_−𝑐𝑓𝑧₀,𝑧₁;𝚺₀𝑑𝑧₁=^𝑐_{𝑐(1−𝑤1)/𝑤0}Φ𝑐−𝑤₀𝑧₀/𝑤₁−𝑣₁₃𝑧₀1−𝑣²₁₃−Φ−𝑐−𝑣₁₃𝑧₀1−𝑣²₁₃𝜙𝑧₀𝑑𝑧₀. (B.5) Thus, _Ω0𝑓𝑧₀,𝑧₁;𝚺₀𝑑𝑧₀𝑑𝑧₁=2^{𝑐(1−𝑤1)/𝑤0}₀Φ𝑐−𝑣₁₃𝑧₀1−𝑣²₁₃𝜙𝑧₀𝑑𝑧₀+^𝑐_{𝑐1−𝑤1/𝑤0}Φ𝑐−𝑤₀𝑧₀/𝑤₂−𝑣₁₃𝑧₀1−𝑣²₁₃×𝜙𝑧₀𝑑𝑧₀−^𝑐₀Φ−𝑐−𝑣₁₃𝑧₀1−𝑣²₁₃𝜙𝑧₀𝑑𝑧₀. (B.6)

Denote 𝑤^∗₀=(𝑁₀+𝑁₁)𝑁₂/(𝑁₀(𝑁₁+4𝑁₂)+𝑁₁𝑁₂) and 𝑤^∗₁=𝑁₀(𝑁₁+𝑁₂)/(𝑁₀(𝑁₁+4𝑁₂)+𝑁₁𝑁₂). For any given 𝑐₁,𝑐₂>0, Pr_𝐻0𝑍^𝑅₁𝜎≤𝑐₁,𝑍^𝐴₁𝜎≤𝑐₁,𝑍^𝐷₁𝜎≤𝑐₁,𝑍^𝑅_𝐽𝜎≤𝑐₂,𝑍^𝐴_𝐽𝜎≤𝑐₂,𝑍^𝐷_𝐽𝜎≤𝑐₂=_Ω1𝑓𝑧₀,𝑧₁,𝑧₂,𝑧₃;𝚺₁𝑑𝑧₀𝑑𝑧₁𝑑𝑧₂𝑑𝑧₃, (B.7) where Ω₁=𝑧₀,𝑧₁,𝑧₂,𝑧₃∶𝑧₀≤𝑐₁,𝑤₀𝑧₀+𝑤₁𝑧₁≤𝑐₁,𝑧₁≤𝑐₁,𝑧₂≤𝑐₂,𝑤^∗₀𝑧₂+𝑤^∗₁𝑧₃≤𝑐₂,𝑧₃≤𝑐₂ (B.8) and 𝑓(𝑧₀,𝑧₁,𝑧₂,𝑧₃;𝚺₁) is the multivariate normal density function for (𝑍^𝑅₁/𝜎,𝑍^𝐷₁/𝜎,𝑍^𝑅_𝐽/𝜎,𝑍^𝐷_𝐽/𝜎)^′ with mean 𝟎₄ and variance-covariance matrix 𝚺₁=1𝑣₁₃𝜌₁₁𝜌₁₃𝑣₁₃1𝜌₃₁𝜌₃₃𝜌₁₁𝜌₃₁1𝑣^∗₁₃𝜌₁₃𝜌₃₃𝑣^∗₁₃1. (B.9)

Note that ∮_Ω1𝑓(𝑧₀,𝑧₁,𝑧₂,𝑧₃;𝚺₁)𝑑𝑧₀𝑑𝑧₁𝑑𝑧₂𝑑𝑧₃ is the sum of nine integrations as follows: 𝐿₁=^{𝑐1(1−𝑤1)/𝑤0}_{−𝑐1(1−𝑤1)/𝑤0}𝑑𝑧₀^𝑐1_−𝑐1𝑑𝑧₁^{𝑐2(1−𝑤∗1)/𝑤∗0}_{−𝑐21−𝑤∗1/𝑤∗0}𝑑𝑧₂^𝑐2_−𝑐2𝑓𝑧₀,𝑧₁,𝑧₂,𝑧₃;𝚺₁𝑑𝑧₃,𝐿₂=^{𝑐1(1−𝑤1)/𝑤0}_{−𝑐1(1−𝑤1)/𝑤0}𝑑𝑧₀^𝑐1_−𝑐1𝑑𝑧₁^{−𝑐2(1−𝑤∗1)/𝑤∗0}_−𝑐2𝑑𝑧₂^𝑐2_{−(𝑐2+𝑤∗0𝑧2)/𝑤∗1}𝑓𝑧₀,𝑧₁,𝑧₂,𝑧₃;𝚺₁𝑑𝑧₃,𝐿₃=^{𝑐1(1−𝑤1)/𝑤0}_{−𝑐1(1−𝑤1)/𝑤0}𝑑𝑧₀^𝑐1_−𝑐1𝑑𝑧₁^𝑐2_{𝑐21−𝑤∗1/𝑤∗0}𝑑𝑧₂^{(𝑐2−𝑤∗0𝑧2)/𝑤∗1}_−𝑐2𝑓𝑧₀,𝑧₁,𝑧₂,𝑧₃;𝚺₁𝑑𝑧₃,𝐿₄=^𝑐1_{𝑐1(1−𝑤1)/𝑤0}𝑑𝑧₀^{(𝑐1−𝑤0𝑧0)/𝑤1}_−𝑐1𝑑𝑧₁^{𝑐2(1−𝑤∗1)/𝑤∗0}_{−𝑐21−𝑤∗1/𝑤∗0}𝑑𝑧₂^𝑐2_−𝑐2𝑓𝑧₀,𝑧₁,𝑧₂,𝑧₃;𝚺₁𝑑𝑧₃,𝐿₅=^𝑐1_{𝑐1(1−𝑤1)/𝑤0}𝑑𝑧₀^{(𝑐1−𝑤0𝑧0)/𝑤1}_−𝑐1𝑑𝑧₁^{−𝑐2(1−𝑤∗1)/𝑤∗0}_−𝑐2𝑑𝑧₂^𝑐2_{−𝑐2+𝑤∗0𝑧2/𝑤∗1}𝑓𝑧₀,𝑧₁,𝑧₂,𝑧₃;𝚺₁𝑑𝑧₃,𝐿₆=^𝑐1_{𝑐1(1−𝑤1)/𝑤0}𝑑𝑧₀^{(𝑐1−𝑤0𝑧0)/𝑤1}_−𝑐1𝑑𝑧₁^𝑐2_{𝑐21−𝑤∗1/𝑤∗0}𝑑𝑧₂^{(𝑐2−𝑤∗0𝑧2)/𝑤∗1}_−𝑐2𝑓𝑧₀,𝑧₁,𝑧₂,𝑧₃;𝚺₁𝑑𝑧₃,𝐿₇=^{−𝑐1(1−𝑤1)/𝑤0}_−𝑐1𝑑𝑧₀^𝑐1_{−(𝑐1+𝑤0𝑧0)/𝑤1}𝑑𝑧₁^{𝑐2(1−𝑤∗1)/𝑤∗0}_{−𝑐21−𝑤∗1/𝑤∗0}𝑑𝑧₂^𝑐2_−𝑐2𝑓𝑧₀,𝑧₁,𝑧₂,𝑧₃;𝚺₁𝑑𝑧₃,𝐿₈=^{−𝑐1(1−𝑤1)/𝑤0}_−𝑐1𝑑𝑧₀^𝑐1_{−(𝑐1+𝑤0𝑧0)/𝑤1}𝑑𝑧₁^{−𝑐2(1−𝑤∗1)/𝑤∗0}_−𝑐2𝑑𝑧₂^𝑐2_{−𝑐2+𝑤∗0𝑧2/𝑤∗1}𝑓𝑧₀,𝑧₁,𝑧₂,𝑧₃;𝚺₁𝑑𝑧₃,𝐿₉=^{−𝑐1(1−𝑤1)/𝑤0}_−𝑐1𝑑𝑧₀^𝑐1_{−(𝑐1+𝑤0𝑧0)/𝑤1}𝑑𝑧₁^𝑐2_{𝑐21−𝑤∗1/𝑤∗0}𝑑𝑧₂^{(𝑐2−𝑤∗0𝑧2)/𝑤∗1}_−𝑐2𝑓𝑧₀,𝑧₁,𝑧₂,𝑧₃;𝚺₁𝑑𝑧₃. (B.10) Moreover, we can obtain that 𝐿₂=𝐿₃, 𝐿₄=𝐿₇, 𝐿₅=𝐿₉, and 𝐿₆=𝐿₈ based on the symmetry of the integration domain for (𝑧₀,𝑧₁) and (𝑧₂,𝑧₃), respectively.

We have 𝑓𝑧₀,𝑧₁,𝑧₂,𝑧₃;𝚺₁=𝜙𝑧₀𝑓𝑧₁∣𝑧₀;𝑣₁₃×𝑓𝑧₂∣𝑧₀,𝑧₁;𝑣₁₃,𝜌₁₁,𝜌₃₁×𝑓𝑧₃∣𝑧₀,𝑧₁,𝑧₂;𝑣₁₃,𝑣^∗₁₃,𝜌₁₁,𝜌₁₃,𝜌₃₁,𝜌₃₃, (B.11) where 𝑓(𝑧₂∣𝑧₀,𝑧₁;𝑣₁₃,𝜌₁₁,𝜌₃₁) is the conditional normal density function as 𝑓𝑧₂∣𝑧₀,𝑧₁;𝑣₁₃,𝜌₁₁,𝜌₃₁=11−𝜌²₁₁−2𝜌₁₁𝜌₃₁𝑣₁₃+𝜌²₃₁/1−𝑣²₁₃×𝜙𝑧₂−𝑧₀𝜌₁₁−𝑣₁₃𝜌₃₁+𝑧₁𝜌₃₁−𝜌₁₁𝑣₁₃/1−𝑣²₁₃×1−𝜌²₁₁+2𝜌₁₁𝜌₃₁𝑣₁₃+𝜌²₃₁/1−𝑣²₁₃⁻¹, (B.12) and 𝑓(𝑧₃∣𝑧₀,𝑧₁,𝑧₂;𝑣₁₃,𝑣^∗₁₃,𝜌₁₁,𝜌₁₃,𝜌₃₁,𝜌₃₃) is the conditional normal density function given by 𝑓𝑧₃∣𝑧₀,𝑧₁,𝑧₂;𝑣₁₃,𝑣^∗₁₃,𝜌₁₁,𝜌₁₃,𝜌₃₁,𝜌₃₃=11−𝜌₁₃,𝜌₃₃,𝑣^∗₁₃Γ⁻¹𝜌₁₃,𝜌₃₃,𝑣^∗₁₃^′𝑧×𝜙³−𝜌₁₃,𝜌₃₃,𝑣^∗₁₃Γ⁻¹𝑧₀,𝑧₁,𝑧₂^′1−𝜌₁₃,𝜌₃₃,𝑣^∗₁₃Γ⁻¹𝜌₁₃,𝜌₃₃,𝑣^∗₁₃^′ (B.13) with Γ the submatrix of 𝚺₁ formed by first three rows and three columns.

Denote (𝜎^∗)²=(𝜌₁₃,𝜌₃₃,𝑣^∗₁₃)Γ⁻¹(𝜌₁₃,𝜌₃₃,𝑣^∗₁₃)^′. Then, we have𝐿₁=∫^{𝑐1(1−𝑤1)/𝑤0}_{−𝑐11−𝑤1/𝑤0}𝜙𝑧₀𝑑𝑧₀∫^𝑐1_−𝑐111−𝑣²₁₃𝜙𝑧₁−𝑣₁₃𝑧₀1−𝑣²₁₃𝑑𝑧₁∫^{𝑐2(1−𝑤∗1)/𝑤∗0}_{−𝑐2(1−𝑤∗1)/𝑤∗0}Φ𝑐₂−𝜌₁₃,𝜌₃₃,𝑣^∗₁₃Γ⁻¹𝑧₀,𝑧₁,𝑧₂^′1−(𝜎^∗)²⁻¹−Φ−𝑐₂−𝜌₁₃,𝜌₃₃,𝑣^∗₁₃Γ⁻¹𝑧₀,𝑧₁,𝑧₂^′1−(𝜎^∗)²⁻¹𝑑𝑧₂,𝐿₂=∫^{𝑐1(1−𝑤1)/𝑤0}_{−𝑐11−𝑤1/𝑤0}𝜙𝑧₀𝑑𝑧₀∫^𝑐1_−𝑐111−𝑣²₁₃𝜙𝑧₁−𝑣₁₃𝑧₀1−𝑣²₁₃𝑑𝑧₁∫^{−𝑐2(1−𝑤∗1)/𝑤∗0}_−𝑐2Φ𝑐₂−𝜌₁₃,𝜌₃₃,𝑣^∗₁₃Γ⁻¹𝑧₀,𝑧₁,𝑧₂^′1−(𝜎^∗)²⁻¹−𝑐−Φ²+𝑤^∗₀𝑧₂𝑤^∗₁−𝜌₁₃,𝜌₃₃,𝑣^∗₁₃Γ⁻¹𝑧₀,𝑧₁,𝑧₂^′1−(𝜎^∗)²⁻¹𝑑𝑧₂,𝐿₄=∫^𝑐1_{𝑐11−𝑤1/𝑤0}𝜙𝑧₀𝑑𝑧₀∫^{(𝑐1−𝑤0𝑧0)/𝑤1}_−𝑐111−𝑣²₁₃𝜙𝑧₁−𝑣₁₃𝑧₀1−𝑣²₁₃𝑑𝑧₁∫^{𝑐2(1−𝑤∗1)/𝑤∗0}_{−𝑐21−𝑤∗1/𝑤∗0}Φ𝑐₂−𝜌₁₃,𝜌₃₃,𝑣^∗₁₃Γ⁻¹𝑧₀,𝑧₁,𝑧₂^′1−(𝜎^∗)²⁻¹−Φ−𝑐₂−𝜌₁₃,𝜌₃₃,𝑣^∗₁₃Γ⁻¹𝑧₀,𝑧₁,𝑧₂^′1−(𝜎^∗)²⁻¹𝑑𝑧₂,𝐿₅=∫^𝑐1_{𝑐11−𝑤1/𝑤0}𝜙𝑧₀𝑑𝑧₀∫^{(𝑐1−𝑤0𝑧0)/𝑤1}_−𝑐111−𝑣²₁₃𝜙𝑧₁−𝑣₁₃𝑧₀1−𝑣²₁₃𝑑𝑧₁∫^{−𝑐2(1−𝑤∗1)/𝑤∗0}_−𝑐2Φ𝑐₂−𝜌₁₃,𝜌₃₃,𝑣^∗₁₃Γ⁻¹𝑧₀,𝑧₁,𝑧₂^′1−(𝜎^∗)²⁻¹−𝑐−Φ²+𝑤^∗₀𝑧₂𝑤^∗₁−𝜌₁₃,𝜌₃₃,𝑣^∗₁₃Γ⁻¹𝑧₀,𝑧₁,𝑧₂^′1−(𝜎^∗)²⁻¹𝑑𝑧₂,𝐿₆=∫^𝑐1_{𝑐11−𝑤1/𝑤0}𝜙𝑧₀𝑑𝑧₀∫^{(𝑐1−𝑤0𝑧0)/𝑤1}_−𝑐111−𝑣²₁₃𝜙𝑧₁−𝑣₁₃𝑧₀1−𝑣²₁₃𝑑𝑧₁∫^𝑐2_{𝑐21−𝑤∗1/𝑤∗0}Φ𝑐₂−𝑤^∗₀𝑧₂𝑤^∗₁−𝜌₁₃,𝜌₃₃,𝑣^∗₁₃Γ⁻¹𝑧₀,𝑧₁,𝑧₂^′1−(𝜎^∗)²⁻¹−Φ−𝑐₂−𝜌₁₃,𝜌₃₃,𝑣^∗₁₃Γ⁻¹𝑧₀,𝑧₁,𝑧₂^′1−(𝜎^∗)²⁻¹𝑑𝑧₂. (B.14) Finally, _Ω1𝑓𝑧₀,𝑧₁,𝑧₂,𝑧₃;𝚺₁𝑑𝑧₀𝑑𝑧₁𝑑𝑧₂𝑑𝑧₃=𝐿₁+2𝐿₂+𝐿₄+𝐿₅+𝐿₆. (B.15)

Acknowledgments

The work was partially supported by the National Natural Science Foundation of China (no. 11225103; 11161054; 11171293; 81072386); the Key Fund of Yunnan Province (no. 2010CC003); The Youth Program of Applied Basic Research Programs of Yunnan Province (2013FD001); the Foundation of Yunnan University (no. 2012CG018).