Simple Algorithms to Calculate Asymptotic Null Distributions of Robust Tests in Case-Control Genetic Association Studies in R

The case-control study is an important design for testing association between genetic markers and a disease. The Cochran-Armitage trend test (CATT) is one of the most commonly used statistics for the analysis of case-control genetic association studies. The asymptotically optimal CATT can be used when the underlying genetic model (mode of inheritance) is known. However, for most complex diseases, the underlying genetic models are unknown. Thus, tests robust to genetic model misspeciﬁcation are preferable to the model-dependant CATT. Two robust tests, MAX3 and the genetic model selection (GMS), were recently proposed. Their asymptotic null distributions are often obtained by Monte-Carlo simulations, because they either have not been fully studied or involve multiple integrations. In this article, we study how components of each robust statistic are correlated, and ﬁnd a linear dependence among the components. Using this new ﬁnding, we propose simple algorithms to calculate asymptotic null distributions for MAX3 and GMS, which greatly reduce the computing intensity. Furthermore, we have developed the R package Rassoc implementing the proposed algorithms to calculate the empirical and asymptotic p values for MAX3 and GMS as well as other commonly used tests in case-control association studies. For illustration, Rassoc is applied to the analysis of case-control data of 17 most signiﬁcant SNPs reported in four genome-wide association studies.


Introduction
The case-control association study is a useful design for testing genetic association (Risch and Merikangas 1996). In such a design, case-control samples for a diallelic marker are summarized in a 2×3 contingency table, where the rows correspond to case-control status and the columns to three genotypes. The null hypothesis of no association is equivalent to no association in the contingency table. Under the alternative hypothesis, if one of the two alleles confers a high risk of the disease, an individual's risk having the disease increases with the number of risk alleles in the genotype. In other words, the alternative is ordered or penetrances are ordered, where a penetrance is the probability of the disease given one of the three genotypes. Four genetic models named recessive (REC), multiplicative (MUL), additive (ADD) and dominant (DOM) models are commonly used (Sasieni 1997;Freidlin et al. 2002).
Two tests, Pearson's χ 2 test (Pearson's test) and the Cochran-Armitage trend test (CATT), are commonly used for the analysis of case-control genetic association studies (Sasieni 1997;Balding 2006). Pearson's test ignores the ordered alternative but the CATT takes the order into account. To apply the CATT, increasing scores are specified a priori for the three genotypes depending on the genetic model. The three CATTs optimal for REC, ADD/MUL and DOM models are obtained (Freidlin et al. 2002;Zheng et al. 2003). When the underlying genetic model can be approximately pre-specified, the CATT with the optimal score is always more powerful than Pearson's test. However, misspecification of the genetic model may result in a substantial loss of power for the CATT. In this case, Pearson's test is more robust (score independent). Therefore, there is a trade-off of robustness and efficiency between Pearson's test and the CATT (Yamada and Okada 2009;Zheng et al. 2009).
In real data analysis, the true genetic model is often unknown. Hence, robust tests are preferable to the CATT or Pearson's test (Freidlin et al. 2002). Robust tests for case-control genetic association studies include the maximin efficiency robust test (MERT, Gastwirth 1966, 1985 with the recently developed R package lawstat (Hui et al. 2008), MAX3 (Freidlin et al. 2002;Gonzalez et al. 2008), the constrained likelihood ratio test (CLRT, Wang and Sheffield 2005), the GMS (Zheng and Ng 2008), and the optimal dose-effect trend test (Yamada and Okada 2009). A detailed review of robust tests and their applications to genetic linkage and association studies can be found in Joo et al. (2009a). In the following we focus on commonly used robust tests for case-control studies.
Suppose a family of scientifically plausible models is defined. Corresponding to each model, an asymptotically optimal, normally distributed test is obtained. Hence, a family of normally distributed tests is formed. When the model is uncertain, a pre-specified test from this family is not fully efficient. The minimum efficiency (e.g., Pitman asymptotic relative efficiency; Noether (1955)) of each test can be obtained over the family of models. A test with higher minimum efficiency is more robust. The MERT achieves the maximum minimum efficiency (Gastwirth 1966). Under some conditions, the MERT can be written as the weighted average of two normally distributed tests with the minimum correlation (called the extreme pair) (Gastwirth 1985). In case-control association studies, the extreme pair corresponds to the CATTs under the REC and DOM models. On the other hand, MAX3 of Freidlin et al. (2002) takes the maximum of the absolute values of three CATTs respectively optimal for the REC, ADD and DOM models. The MERT is often less powerful than MAX3 for case-control studies (Freidlin et al. 2002). However, the MERT is easier to use because it asymptotically follows a normal distribution under the null hypothesis. The MAX3 of Gonzalez et al. (2008) is similar to that of Freidlin et al. (2002). Gonzalez et al. (2008) considered the maximum of the likelihood ratio tests (LRTs) for the three genetic models. A test with performance similar to MAX3 is the CLRT (Wang and Sheffield 2005). It is a LRT but restricts the alternative space to REC, ADD and DOM models. Yamada and Okada (2009) found Pearson's test is a special trend test with a data-driven score, which was also noticed by Zheng et al. (2009). Based on this finding, Yamada and Okada (2009) proposed an optimal dose-effect mode trend test where the genetic effect of the heterozygous genotype is restricted between two homozygous genotypes. The performance of the test of Yamada and Okada (2009) is similar to that of the CLRT. Another recent proposed robust procedure is the GMS (Zheng and Ng 2008), which is a two-phase adaptive test. In the first phase, the underlying genetic model is selected using the Hardy-Weinberg disequilibrium (HWD) trend test between cases and controls (Song and Elston 2006). In the second phase, the CATT with the selected genetic model is applied to test for association. Zheng and Ng (2008) studied an asymptotic null distribution of the GMS.
Among all the tests we considered above, MAX3 and GMS have greater efficiency robustness than Pearson's test and single CATT (Freidlin et al. 2002;Zheng and Ng 2008). On the other hand, MAX3, the CLRT, and the optimal dose-effect mode trend test have comparable power. The GMS seems to be slightly more powerful than MAX3 (Zheng and Ng 2008). However, there is no direct comparison between the GMS and the optimal dose-effect trend test. The asymptotic distributions of the CLRT, the MAX3 of Gonzalez et al. (2008), and the optimal dose-effect mode trend test have been studied. If we directly derive the asymptotic distributions for MAX3 of Freidlin et al. (2002) and GMS, our work would be similar to those of Wang and Sheffield (2005); Gonzalez et al. (2008); Zheng and Ng (2008) and Yamada and Okada (2009), who derived different asymptotic distributions for their robust tests. In addition, there would have no computation benefits. However, we identify a linear dependence structure of the CATTs and use this finding to derive the asymptotic distribution of MAX3 and provide algorithms to obtain the p values of MAX3 and GMS. By doing this way, we greatly simplify the computation and make computation of MAX3 and GMS more efficient. Our finding cannot be applied to the robust tests of Wang and Sheffield (2005); Gonzalez et al. (2008) and Yamada and Okada (2009). Thus, we only focus on MAX3 and GMS.
One common feature of MAX3 and GMS is that their asymptotic null distributions do not follow a standard normal distribution. In practice, a parametric bootstrap procedure may be used to obtain their empirical null distributions or approximate p values. In the parametric bootstrap approach, case-control data are resampled based on the observed case-control genotype counts under the null hypothesis. For each of the m replicates, the two robust test statistics are calculated. The m calculated statistics for each robust test based on the bootstrapped data form its empirical null distribution. In this article, we study correlations among the three trend tests and identify a new asymptotic linear dependence among them. Using this finding, we consider a simple Monte-Carlo approach to approximate null distributions for MAX3 and GMS. This approach does not need to generate case-control data, calculate trend tests and form a robust test. Instead, it only needs to generate bivariate normal distribution with a given correlation and form the robust test. Furthermore, using the same finding, we propose to use the asymptotic null distributions to obtain the p values for MAX3. We find that the proposed method gives a very good approximation to the results of the Monte-Carlo approaches, but the proposed algorithm is much more computationally efficient.
Finally, we have developed a package called Rassoc in the R system (R Development Core Team 2009) for the Monte-Carlo algorithms and the asymptotic algorithms of MAX3 and GMS as well as other commonly used tests in case-control association studies. For illustration, we apply the developed R package to genetic markers reported in four genome-wide association studies (Klein et al. 2005;Hunter et al. 2007;Yeager et al. 2007;The Wellcome Trust Case Control Consortium 2007). These data sets are also incorporated as a data frame in Rassoc.
The rest of this article is organized as follows. In Section 2, the two robust procedures, MAX3 and GMS, with new algorithms are presented. The R package is presented in Section 3 with illustrations. Conclusions are given in Section 4.

Two robust procedures: MAX3 and GMS
Consider a diallelic marker, e.g., a single nucleotide polymorphism (SNP) with alleles D and d and assume D is the allele conferring high risk of the disease. Denote the three genotypes by G 0 = dd, G 1 = Dd and G 2 = DD with genotype frequencies g i = P(G i ) for i = 0, 1, 2. Denote the allele frequencies by P(D) = p and P(d) = 1 − p = q. When Hardy-Weinberg equilibrium (HWE) proportions hold in the population, g i can be written as g 0 = q 2 , g 1 = 2pq and g 2 = p 2 . Denote the penetrance by f i = P(case|G i ), the disease prevalence by k = P(case) = 2 i=0 f i g i and the genotype frequencies in cases and controls by p i = P(G i |case) = f i g i /k and q i = P(G i |control) = (1 − f i )g i /(1 − k) respectively for i = 0, 1, 2. Define genotype relative risks (GRRs) as λ i = f i /f 0 for i = 1, 2 (f 0 > 0). A genetic model is REC, MUL, ADD and DOM if λ 1 = 1, λ 1 = √ λ 2 , λ 1 = (1 + λ 2 )/2 and λ 1 = λ 2 , respectively. Under the null hypothesis of no association H 0 : λ 1 = λ 2 = 1, i.e., p i = q i for i = 0, 1, 2. Under the alternative hypothesis H 1 with D conferring the high risk, the alternative hypothesis can be expressed as H 1 : λ 2 ≥ λ 1 ≥ 1 and λ 2 > 1. Thus, the alternative hypothesis is ordered.
Denote the genotype counts of (G 0 , G 1 , G 2 ) in r cases and s controls by (r 0 , r 1 , r 2 ) and (s 0 , s 1 , s 2 ), respectively. Thus, r = 2 i=0 r i and s = 2 i=0 s i . Denote n i = r i + s i and n = r + s. The case-control samples for a SNP are summarized in Table 1.
The CATT has been employed to test association for the data in Table 1 (Sasieni 1997), which can be written as where (x 0 , x 1 , x 2 ) = (0, x, 1) are scores for (G 0 , G 1 , G 2 ) and x is a real number between 0 and 1. Under H 0 , Z x asymptotically follows a standard normal distribution N (0, 1). Zheng et al. (2003) showed that Z 0 , Z 1/2 and Z 1 are asymptotically optimal for REC, MUL/ADD and DOM models, respectively.

MAX3
When the genetic model is unknown, MAX3, denoted by Z MAX3 = MAX |Z 0 |, |Z 1/2 |, |Z 1 | , has been proposed and used in practice (Freidlin et al. 2002;Sladek et al. 2007; dd Dd DD total Case r 0 r 1 r 2 r Control s 0 s 1 s 2 s total n 0 n 1 n 2 n Table 1: Genotype data of a single SNP. 2008b,a). Because MAX3 covers a wider range of genetic models, it is more robust than any CATT.
Parametric bootstrap is commonly used to approximate the null distribution of MAX3. Given the observed data, a bootstrapped data set is simulated m times. For each bootstrapped data set, MAX3 is calculated, denoted by Z MAX3,j for j = 1 . . . m. Finally, Z MAX3,i , i = 1, . . . , m, form an empirical null distribution for MAX3, which can be used to approximate its critical value or the p value. This bootstrap approach is denoted as boot.
The boot method is computation intensive, in particular for a large m, because the bootstrap case-control samples are sampled each time. The choice of m is at least 10,000 for genome-wide association studies (GWAS) (Sladek et al. 2007;Li et al. 2008a,b). An asymptotic method which reduces the computation is given as follows.

Asymptotic null distribution and p value
Before we present the asymptotic null distribution of MAX3, we first report the following result whose proof is given in Appendix A.
Lemma 1. Under H 0 , Z 0 , Z 1/2 and Z 1 are asymptotically linearly dependent. In addition, 0,1 ) and ρ i,j is an asymptotic null correlation between Z i and Z j for i, j = 0, 1/2, 1, which are given in Appendix A.
Lemma 1 can be used to derive the asymptotic null distribution of MAX3. Denote the joint density function of (Z 0 , Z 1 ) by f (z 0 , z 1 , Σ) where Σ is the variance-covariance matrix of (Z 0 , Z 1 ). Then Applying the above densities, we have where ρ 0,1 can be estimated under H 0 . We use the function integrate in R to approximate the integration numerically. According to the help documentation in R, the numerical solution is calculated based on a global adaptive interval subdivision in connection with extrapolation by the Epsilon algorithm (Piessens et al. 1983). Hence, if the observed statistical value t * is obtained, the p value is given by P . This approach is denoted by asy.
The asy method is more computation efficient because it does not require generating casecontrol samples. Note that Gonzalez et al. (2008) studied the asymptotic distribution of MAX3 of three LRTs under the null hypothesis and proposed formula to calculate the associated p value. They studied the correlations among three LRTs using the Delta method and employed three-fold integrations to calculate the asymptotic distribution. We identified and used the linear dependence of the CATTs and only require double integrations.
Using Lemma 1, a simplified algorithm to approximate the null distribution of MAX3 can also be obtained as follows. First, generate (Z 0 , Z 1 ) from the bivariate normal distribution and calculate Z 1/2 using lemma 1. Second, calculate Z MAX based on the three CATTs obtained. Last, repeat the above procedure m times to simulate an empirical null distribution of MAX3. Note that all correlations are estimated under H 0 . This method is denoted as bvn. Notice that both boot and bvn methods require simulations of m replications to approximate the p values. However, unlike the boot method, the bvn method does not require to generate case-control data.

Numerical results
A simulation is conducted to compare three methods: the existing boot method, and the bvn and asy methods. The empirical cumulative distribution functions based on the three methods are reported based on 1,000,000 replicates. In each replicate, r = 500 cases and s = 500 controls are used under HWE. The results are summarized in Figure 1 with different minor allele frequencies (MAFs). The plots show that the asymptotic null distribution (asy) approximates the two simulation-based empirical null distributions (boot and bvn) very well.   Furthermore, we are interested in the most significant p values in GWAS. Hence, we re-plot in Figure 2 the region of Figure 1 with critical values greater than 3.0 (the region corresponding to small p values). The results show that when MAF = 0.1, the asy and bvn method is slightly different from the boot methods. All three procedures match very well for other MAFs.
The critical values of MAX3 using the above three methods can also be obtained. The results are presented in Table 2. Overall, the critical values calculated using the above three methods match very well except when MAF = 0.1. Besides, the critical values are not sensitive to the change of MAFs.
We also investigate the accuracy of the asymptotic distribution of MAX3 for small smaple sizes (n = 100 and 200) by simulations with 1 million replicates. Table 3 reports the empirical type I error rates of MAX3 using simulations. The results show that the test MAX3 based on the asymptotic null distribution is conservative when the MAF and sample sizes are small, i.e., MAF = 0.1 and n = 100 or 200. In fact, with MAF = 0.1, the expected total counts n 0 in the 2 × 3 tables (Table 1) are 1 and 2 for n = 100 and 200 respectively, and the tables are highly unbalanced. Under this situation, the asymptotic theory of MAX3 may not work well when the sample size is small.

GMS
The GMS is another robust method which has two phases (Zheng and Ng 2008). In phase 1, the Hardy-Weinberg disequilibrium trend test (HWDTT) proposed by Song and Elston (2006) is used to detect the underlying genetic model. In phase 2, an optimal CATT corresponding  to the selected genetic model is used for testing association. Since the GMS needs the HWE proportions to hold in the population, we assume HWE proportions to hold throughout this section.
The Hardy-Weinberg disequilibrium (HWD) coefficients in cases and controls are denoted by ∆ P = P 2 − (P 2 + P 1 /2) 2 and ∆ Q = Q 2 − (Q 2 + Q 1 /2) 2 where P i = P(G i |case) and Q i = P(G i |control), i = 1, 2. Zheng and Ng (2008) showed that, under HWE and when D confers a high risk of the disease, ∆ P − ∆ Q > 0 under the REC model and ∆ P − ∆ Q < 0 under the DOM model. The HWDTT of Song and Elston (2006) is the standardized test of ∆ P − ∆ Q = { p 2 − ( p 2 + p 1 /2) 2 } − { q 2 − ( q 2 + q 1 /2) 2 }, where p i = r i /r and Q i = s i /s for i = 0, 1, 2. The HWDTT can be written as (Song and Elston 2006) Under the null hypothesis of no association, Z HWDTT ∼ N (0, 1). With a pre-specified threshold c > 0 (commonly set at Φ −1 (0.95) = 1.645), Zheng and Ng (2008) classified the underlying genetic model as REC if Z HWDTT > c, DOM if Z HWDTT < −c and MUL/ADD otherwise. Then, when the underlying genetic model is selected, Z x optimal for the selected genetic model is used for testing association in phase 2.
Note that we assume D conferring a high risk of the disease in the above discussion. Although the sign of Z HWDTT is independent of which allele confers a high risk, the optimal CATT for the REC or DOM models does depend on such an assumption. If D confers a high risk as we assume, Z 0 and Z 1 are optimal for the REC and DOM models respectively. On the other hand, if d confers a high risk, then Z 0 and Z 1 are optimal for DOM and REC models, respectively. In this case, the expected values of Z 0 and Z 1 are negative. Which allele confers a high risk can be detected using Z 1/2 (Joo et al. 2009b). That is, if Z 1/2 > 0, Z 0 , Z 1/2 and Z 1 are optimal for REC, MUL/ADD and DOM models. On the other hand, if Z 1/2 < 0, −Z 1 , −Z 1/2 and −Z 0 are optimal for REC, MUL/ADD and DOM models. Finally, the test statistic for GMS can be written as (Joo et al. 2009b) Z GMS does not follow N (0, 1) under H 0 .
Like MAX3, the previous bootstrap (boot) method can be applied similarly. The case-control data are resampled and the GMS is applied to each replicate. In addition, a similar asymptotic method can be proposed as follows.

Asymptotic null distribution and p value
Denote the joint density function (Z 0 , Z 1/2 , Z HWDTT ) and (Z 1 , Z 1/2 , -Z HWDTT ) by f (z, Σ 1 ) and f (z, Σ 2 ) respectively, where f is the density function of a trivariate normal distribution and Σ 1 and Σ 2 are the variance-covariance matrixes. From ( 1), follow Joo et al. (2009b), for any t ≥ 0, We use the function pmvnorm from function mvtnorm (Genz et al. 2009) in R to numerically calculate the normal probabilities. According to the R documentation, we adopt the GenzBretz algorithm to approximate the numerical solution. The methodology is described in Genz (1992Genz ( , 1993. Substituting the correlations by the estimates under H 0 , we can use this formula to calculate the critical values and p values of the GMS.
Besides, we discuss a similar simplified bivariate normal distribution simulation. First, we present the following lemma whose proof is given in Appendix B.
Lemma 2. Under H 0 , Z 0 , Z HWDTT and Z 1 are asymptotically linearly dependent. In addition, , where ρ 0,1 is an asymptotic null correlation between Z 0 and Z 1 and ρ i is an asymptotic null correlation between Z i and Z HWDTT for i = 0, 1, which are given in Appendix B.
Using Lemma 2, the algorithm bvn for MAX3 can be obtained similarly to calculate Z HWDTT and obtain empirical null distribution of the GMS.

Numerical results
Simulation is used to compare the above three methods for the GMS. Analogy to Figures 1  and 2   8e-4 1e-3 1e-4 7e-5 9e-5 1e-5 4e-6 9e-6 Table 5: Empirical type I error rates for GMS with small sample sizes n and equal numbers of cases and controls. The nominal level α is based on the asymptotic distribution. and 4 for different regions of probabilities. The parameter settings are the same as those in Figure 1 for MAX3. The plots in Figures 3 and 4 show that these three methods match very well except for small MAF (MAF = 0.1). For example, in Figure 4, when MAF = 0.1, the cumulative distribution functions for the GMS at 3.0 are around 0.994 using the asy and bvn methods while around 0.996 using the boot method.
Next we report the critical values of the GMS using the above three methods based on 1,000,000 replicates (Table 4). The results show that the critical values calculated from different methods match well except for MAF = 0.1. For example, in Table 4, when MAF = 0.1 and the significance level was 0.05, the critical value calculated from the boot method was 2.153 while those calculated from the asy and bvn methods were both 2.207. Finally, analogy to MAX3, we also report the empirical type I error rates of the GMS when the asymptotic null distribution is used. The results for small sample sizes n = 100 and 200 are summarized in Table 5. It is observed that, with a small sample size, the test GMS based on the asymptotic null distribution is generally conservative especially when MAF is small (MAF = 0.1). For example, when n = 100 and MAF = 0.1 with nominal level equal to 0.05 and 0.01, the empirical type I error rates were 0.016 and 0.002, respectively.

R package: Rassoc
In this section we describe the R package Rassoc developed for calculating the p values of MAX3 and GMS based on the three methods introduced above. This package also contains some other commonly used tests in case-control association studies. The package is available from the Comprehensive R Archive Network at http://CRAN.R-project.org/package= Rassoc. After it has been installed, the package may be loaded with:

R> library("Rassoc")
Rassoc contains the following functions: ABT(data): Calculate the statistic and associated p value of the allelic based test (ABT) in case-control association studies (Sasieni 1997). The ABT detects association by comparing the allele frequencies between cases and controls. Under the null hypothesis of no association, the ABT follows the standard normal distribution N (0, 1).
CATT(data, x): Calculate the statistic and associated p value of the Cochran-Armitage trend test (CATT) in case-control association studies for a given genetic model (x = 0, 1/2, or 1). Under the null hypothesis of no association, the CATT follows the standard normal distribution N (0, 1).

MERT(data):
Calculate the statistic and associated p value of the maximin efficiency robust test (MERT) in case-control association studies (Freidlin et al. 2002). Under the null hypothesis of no association, the MERT follows the standard normal distribution N (0, 1).

MAX3(data, method, m):
Calculate the statistic and associated p value of the MAX3 test in case-control association studies. The p value is calculated based on two empirical Monte-Carlo methods or the proposed asymptotic method.
GMS(data, method, m): Calculate the statistic and associated p value of the GMS test in case-control association studies. The p value is calculated based on two empirical Monte-Carlo methods or the proposed asymptotic method.
The arguments in the functions are defined below: data: The 2 × 3 contingency table for analysis. The rows represent disease status and the columns represent genotypes. See the data formats in Table 1 x: The score for the CATT. It can be any real number between 0 and 1.

method:
The methods for calculating the p values of MAX3 and the GMS. "boot" represents the parametric bootstrap method; "bvn" represents the simulation from the bivariate normal distribution; "asy" represents the asymptotic null distribution method.
m: The number of replicates for "boot" and "bvn". It can be any positive integer for "asy".
Since the ABT, CATT and MERT follow the standard normal distribution N (0, 1) under the null hypothesis, we focus on MAX3 and GMS in this package. For each variable that has been imputed, MAX3 or GMS creates a list containing the name of the method selected, the values of test statistics and their corresponding p values. Therefore, using these two functions, we can calculate the p values of MAX3 and GMS for a dataset. For example, a simulated case-control dataset called a under the null hypothesis is given below: R> p <-c(0.25, 0.5, 0.25) R> ca <-rmultinom(1, 500, p) R> co <-rmultinom(1, 500, p) R> a <-matrix(rbind(ca, co), nrow = 2, byrow = TRUE) R> a 139 249 112 136 244 120 We can apply our programs to this dataset to calculate the p values of MAX3 and GMS respectively as follows: R> MAX3(a, "boot", 100000) The MAX3 test using the boot method data: a statistic=0.5993 p-value=0.7907 R> MAX3(a, "bvn", 100000) The MAX3 test using the bvn method data: a statistic=0.5993 p-value=0.7935 R> MAX3(a, "asy", 1) The MAX3 test using the asy method data: a statistic=0.5993 p-value=0.7933 R> GMS(a, "boot", 100000) The GMS test using the boot method data: a statistic=0.4894 p-value=0.6608

R> GMS(a, "bvn", 100000)
The GMS test using the bvn method data: a statistic=0.4894 p-value=0.6609 R> GMS(a, "asy", 1) The GMS test using the asy method data: a statistic=0.4894 p-value=0.6621 From the above output, the test statistics for MAX3 and GMS are 0.5993 and 0.4894 respectively. The corresponding p values of MAX3 for the simulated dataset a using "boot", "bvn" and "asy" methods are 0.7907, 0.7935 and 0.7933 respectively. On the other hand, for the GMS, they are 0.6608, 0.6609 and 0.6621 respectively. Note that the number of replicates 1 in the functions MAX3 and GMS can be any positive integer by using "asy" method. The empirical p values are calculated based on 100,000 replicates.

Real data analysis
For the purpose of illustration, we apply Rassoc to SNPs reported from four GWAS with 100,000 to 500,000 SNPs for age-related macular degeneration (AMD) (Klein et al. 2005), two cancer studies (Hunter et al. 2007;Yeager et al. 2007), and a hypertension study (The Wellcome Trust Case Control Consortium 2007). The datasets are summarized in Table 6 which were also given in Li et al. (2008a).

SNP ID
Note that all the empirical p values are calculated based on 1,000,000 replicates due to the computational burden. The results are summarized in Table 7.
Results show that when the analytical p values are greater than 10 −5 , the analytical p values and simulated p values match well. However, when the analytical p values are of order of magnitude 10 −6 , they do not match too well. One reason may be due to the fact that one million replicates may not be enough for p values of that magnitude. In conclusion, the parametric bootstrap and bivariate normal distribution methods based on simulation can be easily used in candidate-gene association studies in which several markers are tested, while the asymptotic distribution method is preferred in large-scale association studies such as GWAS.
The computing times of the three approaches are also different. The asymptotic distribution method does not need any replication, so it takes least time to compute. The parametric bootstrap method takes more time to compute. The bivariate normal distribution method replicates directly on statistics while the parametric bootstrap method needs to generate case-control data in each replicate. For MAX3, using our Pentium(R) D CPU 3.00 GHZ, 1.00 GB of RAM computer, it took less than 1 second to calculate the asymptotic p value of the first row of the data in Table 6, about 1 minute using the bivariate normal distribution method, and nearly 18 minutes using the parametric bootstrap method. For the GMS, the three methods took about 1 second, 4 minutes, and 22 minutes, respectively. All empirical methods were based on 1,000,000 replicates.

Conclusion
Unlike the model-based CATT which assumes the underlying genetic model is known, MAX3 and GMS do not rely on any prior information of the underlying genetic model while perform robustly across a family of scientifically plausible genetic models. For many complex diseases, the underlying genetic models are unknown. Thus, MAX3 and GMS are preferable. However, their null distributions (asymptotic and empirical) have not been fully investigated. In this article, we studied the dependence structures among the CATTs for the three most common genetic models and identified an asymptotic linear relationship among them. Using this finding, we proposed to obtain asymptotic p values of the test of statistics MAX3 and GMS based on numerical integrations. The proposed method was compared to the bivariate normal simulation and the parametric bootstrap methods. Our simulation results demonstrated that the proposed asymptotic distribution method performs well. The proposed method is computationally efficient because no simulation is needed.
Moreover, we developed the R package Rassoc which incorporates all the methods discussed above as well as other commonly used tests in case-control association studies. We illustrated the use of Rassoc based on simulated datasets and some associated SNPs from four real GWAS datasets. In practice, empirical p values of MAX3 and GMS for a candidate SNP analysis or the corresponding analytical p values for GWAS can be obtained by simply setting the arguments in our R programs incorporated in the package. We also plan to update Rassoc by adding other recently developed tests in case-control genetic association studies in the future, including those of Gonzalez et al. (2008) and Yamada and Okada (2009). In this article, we only consider the null distributions of some robust tests. The asymptotic distributions under the alternative hypothesis are not considered. Li et al. (2009) studied the asymptotic power for MAX3 based on multivariate normal among the CATTs. Applying the linear dependence structure, their computation can be simplified. This may be also worth investigating in future.