Testing for homogeneity of gametic disequilibrium across strata

Background Assessing the non-random associations of alleles at different loci, or gametic disequilibrium, can provide clues about aspects of population histories and mating behavior and can be useful in locating disease genes. For gametic data which are available from several strata with different allele probabilities, it is necessary to verify that the strata are homogeneous in terms of gametic disequilibrium. Results Using the likelihood score theory generalized to nuisance parameters we derive a score test for homogeneity of gametic disequilibrium across several independent populations. Simulation results demonstrate that the empirical type I error rates of our score homogeneity test perform satisfactorily in the sense that they are close to the pre-chosen 0.05 nominal level. The associated power and sample size formulae are derived. We illustrate our test with a data set from a study of the cystic fibrosis transmembrane conductance regulator gene. Conclusion We propose a large-sample homogeneity test on gametic disequilibrium across several independent populations based on the likelihood score theory generalized to nuisance parameters. Our simulation results show that our test is more reliable than the traditional test based on the Fisher's test of homogeneity among correlation coefficients.


Background
Measuring gametic disequilibrium can provide important information about aspects of population histories and mating behavior [1] and can be useful in locating disease genes [2]. The term gametic disequilibrium is used in this article instead of the traditional term linkage disequilibrium to measure the extent of non-random association because such non-random association may be present between unlinked loci [3]. Various measures of gametic disequilibrium have been proposed [4][5][6], ranging from pairs of diallelic loci model to multiple multiallelic loci model. In this article, we consider the gametic disequilib-rium which is defined as the difference between the gametic probability and its expected probability under the assumption of no statistical association of alleles, and the gametic disequilibrium calculations are based on twoallele, two-locus model [7].
Consider two loci, A and B, each having two possible alleles (A 0 , A 1 ) and (B 0 , B 1 ), respectively. With two loci and two alleles, there are four possible gametes, namely, A 0 B 0 , A 0 B 1 , A 1 B 0 and A 1 B 1 . The gametic disequilibrium between the two loci is defined by where and denote the allele probabilities of A i and B j , denotes the gamete probability of A i B j , i, j = 0, 1. Suppose that the gametic data are available from K strata and let p ijk denote the gametic probability of array of A i B j for the k-th stratum, i, j = 0, 1; k = 1,...,K, ∑ i,j p ijk = 1 for each k. According to the relationship between allelic probability and gametic probability, the allele probabilities of A 0 , A 1 , B 0 and B 1 are derived as p 0+k , p 1+k , p +0k and p +1k , respectively. Here "+" denote the summation over 0 and 1, for example, p 0+k = p 00k + p 01k . For stratum k (k = 1,...,K), the gametic disequilibrium is calculated as It is easy to show that D k is bounded by where D k,min = -min{p 1+k p +1k , p 0+k p +0k }, D k,max = min{p 1+k p +0k , p 0+k p +1k }. Testing for the homogeneity of gametic disequilibrium among strata can be informative in discriminating among the evolutionary agents generating them in natural population [8]. Detecting gametic disequilibrium can be informative in mapping gene and providing meaningful clues of population evolution.
Combining the evidence of gametic disequilibrium across several strata may be more sufficient to support the clues, in contrast to analysis with each strata. In this case, it is crucial to test the homogeneity of gametic disequilibrium across strata before combining the data. For this purpose, it is interesting to consider the following hypothesis H 0 : D 1 = ʜ = D K versus H 1 : D i ≠ D j for at least a pair i ≠ j.
( 1 ) Weir [9] recommended a homogeneity test on gametic disequilibrium, based on Fisher's test of homogeneity among correlation coefficients [10]. In his method, the gametic disequilibrium D k is first transformed to a correlation coefficient r k by r k = D k / , r k is then transformed to a normal variable z k by Fisher's z transformation, and a weighted sum of squares of the z values which has χ 2 distribution with K -1 degrees of freedom is finally proposed for testing homogeneity of gametic disequilibrium. As pointed out by Zapata and Alvarez [8], this test is actually for homogeneity of r values instead of D values. They may not be equivalent when the allele probabilities are different across strata. Instead, Zapata and Alvarez [8] suggested the use of the normalized difference D' [11]. Specifically, is the ratio of D k to D k,max when D k > 0, or the ratio of D k to -D k,min when D k < 0. Zapata and Alvarez obtained the bias-corrected confidence interval for each D' value across strata via the bootstrap method. Hence, acceptance or rejection of homogeneity of D' values can be determined by evaluating the obtained confidence intervals. For the example considered in Zapata and Alvarez [8], there is no intersection for the confidence intervals obtained from all strata. Hence, one has evidence to reject the null hypothesis of homogeneity. Unfortunately, Zapata and Alvarez [8] did not discuss the decision rules for cases such as intersections exist but the extent are different. Hence, no rigorous rule based on this confidence interval approach was proposed and this makes their method less practicable. However, no rigorous rule based on this confidence interval approach was proposed and this makes their method less practicable. It should be noted that the homogeneity test of either r values or D' values is not equivalent to the homogeneity test of D values. In particular, transformation D' only guarantees that the range of D' is [-1, 1]. However, there remains difficulties in interpreting the value of D'. Lewontin [11] noted that values of D' at different loci and in different populations tend to vary with the values of the allele probabilities, so that the problem of cross-locus and cross-population comparisons is not fully overcome by the use of D'. In this article, without doing any transformation, we develop an asymptotic homogeneity test directly based on D values via score method.

Homogeneity test
Let x ijk (i, j = 0, 1 and k = 1,ʜ,K) be the number of the gamete A i B j in the k-th stratum with the total gametes being n k = x 00k + x 01k + x 10k + s 11k . Let M(n k , {p ijk }) denote the quadrinomial distribution with parameter vector (p 00k , p 01k , p 10k , p 11k )'. Thus, we have {x ijk : i, j = 0, 1} ~ M(n k , {p ijk }) for k = 1,...,K. The homogeneity hypothesis in (1) is of interest in this article. Here, we assume that K is fixed and n k is sufficiently large for k = 1, 2,...,K. Noticing that p 00k = p 0+k p +0k + D k , p 01k = p 0+k p +1k -D k , p 10k = p 1+k p +0k -D k , p 11k = p 1+k p +1k + D k , the log-likelihood for the k-th stratum can be expressed in terms of D k , p 1+k and p +1k (k = 1,....,K). That is, , x k x k Here, D* is analogous to the well-known Mantel-Haenszel estimator [13]. It is a consistent estimator to D. In general, it is not an efficient estimator to D. The proof of consistency and the conditions for achieving asymptotic efficiency for D* is presented in Appendix. We notice that the calculation of in (2) is quite tedious. Nonetheless, it is easy to show that is simply given by (see Appendix for the proof). It can be shown that X 2* has an asymptotic chi-square distribution with K -1 degrees of freedom under H 0 . Therefore, the homogeneity hypothe-

Asymptotic power and sample size
We will present the asymptotic power and sample size formulae based on X 2* [14]. For this purpose, we assume n k = na k for some n and a k > 0. Let , and be the true parameter values for D k , p 1+k and p +1k under H 1 , where k = 1, 2,ʜ,K and for at least a pair k ≠ j. Thus, the asymptotic power for the homogeneity score test X 2* at α level is given by },

Availability and requirements
We have implemented the test procedures for computing our score statistic X 2* in a Matlab project. Project name: gametic disequilibrium homogeneity score test (GDHST); Project home page: http://math.nenu.edu.cn/jhguo/pro gram.htm; Operating system: Windows XP; Programming language: Matlab 6.1; Licence: GNU GPL.

Simulation results
To evaluate the performance of our proposed homogeneity score test, we include the homogeneity test recommended by Weir [9] in our comparison study. The corresponding test statistic for homogeneity is given by where K is the total number of strata, n k is the total gamete number in stratum k, is the Fisher's z transformation with and (x 00k , x 01k , x 10k , x 11k )' being the number of the gamete array in the k-th stratum, and is the average of the z k values.
We investigate the performance of X 2* and T 2 in terms of type I error rate and power. For type I error rates, we consider both equal and unequal allele probabilities varying from 0.1 to 0.  Table 1, 2, 3, 4. Table 1 shows the performance of empirical type I error rates for X 2* and T 2 with equal allele probabilities across K = 3 strata. We observe the following.
1. When D is large (i.e., Dmax), both tests generally appear to be quite liberal (e.g., empirical size being 10 times of the nominal level), especially for small sample size (e.g., nk = 50) and small allele probability (e.g., p1+ = p+1 = (0.1, 0.1, 0.1)'). Such liberty in empirical size is more severe in T2 than in our asymptotic homogeneity test X2* and is significantly alleviated in X2* when sample size increases. However, sample size increase does not alleviate the liberty of T2 much. In fact, even for nk = 3200 for k = 1, 2, 3, T2 is still very liberal for D = 0.045 with empirical type I errors rate being 0.456 (data are not shown).
2. For other settings, both tests perform quite satisfactorily in the sense that their empirical sizes are well controlled around the pre-chosen nominal level. In general, the larger the sample size, the closer the empirical type I error rate to the pre-chosen nominal level. Table 2 reports the empirical size performance of X 2* and T 2 for unequal allele probabilities across K = 3 strata. We observe similar phenomena above. However, our asymptotic homogeneity test X 2* performs quite well in all settings under consideration for moderate to large sample sizes (i.e., n k = 100 and 200) while it is not the case for T 2 .
For T 2 , the resultant empirical type I error rate can be extremely inflated even for large sample design (e.g., more than 17 times of the nominal level when n k = 200 (for k = 1, 2, 3), p 1+ = p +1 = (0.5, 0.3, 0.1)', and D = 0.045). Since many type I error rates for X 2* and T 2 are liberal in Tables 1 to 4. The two-sided t-test is conducted to determined if an empirical type I error rate is significantly different from the nominal lever of 0.05. The t-test statistics is where m = 5000 and W represents the empirical type I error rate of X 2* or T 2 . Here, the t-test is almost identical to the z-test for the sample size is very large. Those empirical type I error rates which are significantly different from the nominal level of 0.05 are underlined in Tables 1 to 4. In Table 1, the total number of significant difference from the nominal level of 0.05 for X 2* and T 2 is 28 and 38, respectively. The pair (28, 38) can be further decomposed to (14,14), (8,13) and (6,11) according to n = 50, 100 and 200. The decreasing rate of the number of empirical type I error rates which is significant different from the nominal level of 0.05 for X 2* is 14/18-6/18 = 44.4% as n increases from 50 to 200. While the corresponding decreasing rate for T 2 is 14/18-11/18 = 16.7%. It is easy to see that our X 2* is less liberal than T 2 as sample size increases.   The empirical type I error rate which is significant for the two-sided t-test at the nominal level of 0.05 is marked with underline.  The empirical type I error rate which is significant for the two-sided t-test at the nominal level of 0.05 is marked with underline.  The empirical type I error rate which is significant for the two-sided t-test at the nominal level of 0.05 is marked with underline.  The empirical type I error rate which is significant for the two-sided t-test at the nominal level of 0.05 is marked with underline.
For Table 2, the total number of significant difference from the nominal level of 0.05 for X 2* and T 2 is 17 and 33, respectively. The pair (17, 33) can again be decomposed to (14,14), (8,13) and (6,11) according to n = 50, 100 and 200. The decreasing rate of the number of empirical type I error rates which is significant different from the nominal level of 0.05 for X 2* is 10/18-1/18 = 50.0% as n increases from 50 to 200. While the decreasing rate for T 2 is 12/18-10/18 = 11.1%. The decreasing rate of our X 2* is again more significant than that of T 2 .
In Table 3 to 4, the strata increases from 3 to 5. However, the decreasing rates of the number of empirical type I error rates which is significant different from the nominal level of 0.05 for Tables 3 and 4 is very close to that of Tables 1 and 2, respectively. Therefore, we have reason to believe that this decreasing rate is not greatly affected by the number of strata. Table 5 summarizes the empirical powers for X 2* and T 2 .
In view of the above results, we prefer the proposed homogeneity test X 2* to the traditional T 2 which is based on the Fisher's test of homogeneity among correlation coefficient.

Real and hypothetical examples
It is reported that mutations at the cystic fibrosis transmembrane conductance regulator gene (CFTR) cause cystic fibrosis, the most prevalent severe genetic disorder  in individuals of European descent. Mateu [15] conducted a worldwide genetic analysis of the CFTR region and analyzed normal allele and haplotype variation at two singlenucleotide polymorphisms (SNPs), namely the T854/ AvaII (2694 T/G) and TUB20/PVUII (4006-200 G/A). The T854 and TUB20 markers can be used to define the core haplotypes since they are diallelic, have presumably much lower mutation rates than the other polymorphisms and the ancestral state can be inferred for them.
Mateu [15] reported the T854-TUB20 haplotype frequencies by 18 populations. After communicating with one of their coauthors (Prof. Kenneth, pers. comm. 1996), it was found that their reported gametic frequencies were actually the maximum likelihood estimates of the gametic probabilities obtained from HAPLO, a software which can be applicable to missing data. In other words, all individuals with results for at least one of the two markers were included to estimate the gametic frequencies and no actual gametic counts were available. To create the gametic counts for each population, we first estimate the total number of participants in each population by the number of individuals who yielded results for at least one of the two markers. The reported gametic frequencies of each population given in Mateu [15] are multiplied to the estimated number of participants of this population and the closest integers are then taken to be the estimated gametic counts. The estimated gametic counts across the 18 populations are reported in Table 6, which is adopted as the real data in all subsequent analysis.
It is noticed that the gametic counts for the populations of Japanese (14th) and Surui (18th) are (0, 32, 0, 12)' and (0, 7, 0, 35)', respectively and their estimated gametic disequilibrium D k , D k,min and D k,max are all equal to zero. Therefore, we will exclude these two populations for subsequent homogeneity testings. We consider the following scenarios.
(i) Homogeneity of gametic disequilibrium among the 16 populations (i.e., excluding Japanese and Surui). The statistic value of our proposed X 2* is 121.35 with p-value being less than 0.0001 while that of T 2 yields 99.64 with p-value being less than 0.0001. In this case, both tests reject the homogeneity hypothesis at the 0.05 nominal level.
(ii) Homogeneity of gametic disequilibrium among those populations with the same numbers of participants for both markers T854 and TUB20 (i.e., Mbuti, Yemenites, Druze, Adygei, Catalans, Basques, Chinese, and Nasioi). To end this section, we analyze the hypothetical example of gametic disequilibrium between tow loci (A, B) in ten populations described in Zapata and Alvarez [8]. Here, the gametic counts are simply set by multiplying the haplotype frequencies given in Zapata and Alvarez [8] by 1000. The data are reproduced in Table 7. Obviously, the r values are homogeneous across the ten populations. For D' values, Zapata and Alvarez [8] utilized the bias-corrected nonparametic bootstrap method to obtain the 95% confidence interval for each D' values. Observing that the resultant confidence intervals have no intersection, they concluded that D' are heterogeneous. They suggested tests for homogeneity of gametic disequilibrium should be based on D', whose range is allele probability independent, rather than r. Although, the D values in Table 7 seem to be homogeneous, our homogeneity score test yields X 2* = 33.44 with p-value being less than 0.0001. Therefore, our test procedure also suggests the rejection of the homogeneity of gametic disequilibrium across the ten populations. In this case, our test reaches the same conclusion drawn by Zapata and Alvarez [8].

Discussion
Verification of the homogeneity assumption of gametic disequilibrium across several populations is crucial in gametic disequilibrium analysis. We note that traditional homogeneity test on gametic disequilibrium is based on the Fisher's test of homogeneity among correlation coefficients. However, our simulations demonstrate that this traditional test may not perform satisfactorily. Specifically, it can be very conservative or liberal, for almost all the cases in which the common true gametic disequilibrium D is bounded away from zero. Most importantly, these kinds of conservativeness and liberty can not effectively alleviated with increased sample sizes.
Our proposed large-sample homogeneity score test on gametic disequilibrium across several independent populations requires the count of haplotypes as input. In practice, only genotype data can be obtained in most situations. To employ our method, one can use some haplotyping software, such as PHASE, HAPLOTYPER, to resolve the genotype data as haplotype data. In this way, it separates haplotype phasing and gametic disequilibrium homogeneity test. Naturally, it is more promising to extend our method which can directly handle the genotype data. In this sense, model assumptions are based on genotype data. However, the haplotype phase uncertainty for the double heterozygotes makes the definition of D p 1+ p +1 Table 7: Hypothetical example of gametic disequilibrium between two loci (A, B) with twoalleles (A 0 , A 1 and B 0 , B 1 , respectively) across ten populations

Gamete counts Allele frequencies Disequilibrium Estimation
Population gametic disequilibrium can not be directly expressed by the genotype data even assuming Hardy-Weinberg equilibrium holds. It may severely affect the further derivation of the corresponding score test. Thus, extending our method to handle genotype data is an avenue we intend to explore future.

Conclusion
In this article, we propose a large-sample homogeneity test on gametic disequilibrium across several independent populations based on the likelihood score theory generalized to nuisance parameters. Our simulation results show that our test is more reliable than the traditional test based on the Fisher's test of homogeneity among correlation coefficients. Although our test may also demonstrate conservativeness and liberty in some cases, unlike the traditional test these issues can be effectively resolved by increasing sample sizes. For design purpose, sample size formula that controls power is derived.
where and are the MLEs of D k , p 1+k and p +1k , respectively. Hence, the asymptotic variance of is .
On the contrary, by the Central Limit Theorem, follows an asymptotic normal distribution N(0, Σ k ). By δ method, we immediately get that follows an asymptotic normal distribution N(0, w k (D k , p 1+k , p +1k )). Therefore, we can obtain the exact expression (D k , p 1+k , p +1k ) = n k /w k (D k , p 1+k , p +1k ).
Naturally, the expression of (D, p 1+k , p +1k ) is just (D k , p 1+k , p +1k ) by substituting D for D k .