To appreciate the limitations of past work on methods of testing HWE across strata, we briefly review some of the key aspects, because some points provide a deeper understanding of the issues and some developments are useful for our exact methods. For notation, we use A and B to represent the rare and common alleles, respectively, of a locus, with respective allele frequencies p and q=1-p (pq). As explained by Weir,⁴ the frequencies of the three genotypes can be expressed in terms of the allele frequencies and a measure of departure from HWE (coefficient of disequilibrium D):

Haldane⁵ was the first to develop a stratified test for HWE. He did this by recognizing that D is expected to be zero when HWE is true. When HWE holds, the allele counts are sufficient statistics, and so the probability of the genotype counts, conditional on the allele counts, allows one to compute the mean and variance of the parameter of interest. Let N_AA, N_AB, and N_BB denote the counts of the genotypes. To estimate D from a sample, it may be tempting to plug in the sample estimates, , , and . However, because genotype counts are negatively correlated, this would lead to a biased estimate for small samples; the bias diminishes as the total sample size N increases. As emphasized by Smith,⁶ this may not be an important bias for large samples, but, when adding contributions across strata, each of small sample size, the bias can be amplified. Hence, Haldane used an unbiased estimate of D,

To combine the h values over strata, Haldane first standardized each stratum’s h by its SE, , and then summed these standardized terms over the K strata to compute the combined statistic , which has an approximate standard normal distribution when HWE is true. A problem with this approach is first standardizing and then summing. A more powerful approach would be to first sum and then standardize, much like the way the Mantel-Haenszel test is constructed for testing a common odds ratio over strata⁷ or the way NPL statistics can be optimally combined across pedigrees.⁸ Hence, we propose the statistic , which also has an approximate standard normal distribution. Positive values of T imply an excess of homozygotes, and negative values an excess of heterozygotes, making it simple to interpret significant departures from HWE.

In contrast to Haldane’s method, Smith⁶ computed a weighted sum of the h_k values, using weights proportional to the inverse of the variance. However, his derivations were a bit odd, because he assumed that the allele frequencies are constant over strata, which is counter to what we wish to assume.

The methods by Haldane and Smith are appropriate if the h_k values are all in the same direction (positive or negative) over strata, but they can cancel each other if this is not the case, which will weaken power. For this reason, others have assumed that the ratio θ=P²_AB/4P_AAP_BB is constant over strata, much like the assumption of a constant odds ratio across stratified 2×2 tables in epidemiological studies.^9–11 Nonetheless, the resulting test for HWE across strata, derived by Olson,^9,10 is also based on a weighted sum of h_k values. Nam¹¹ derived score statistics based on likelihoods that depend on the θ parameter and showed that his combined tests for HWE had properties similar to the test proposed by Olson. Hence, the variety of proposed tests for HWE that combine information across strata are all based on the stratum-specific h_k values, with merely slightly different ways of weighting the contribution from each stratum. Simulations (not shown) suggest that the type I error rate and power of the different methods are similar, and our simple T statistic would provide a powerful test for HWE when the sample sizes of the strata are not too small and departures from HWE are all in the same direction.

Because summing h_k values over strata can cancel each other when they differ in sign, Troendle and Yu¹² proposed a statistic that is analogous to summing h²_k/Var(h_k) over strata. The resulting statistic has a χ² distribution with K df. Although this method can have greater power when the h_k values differ in sign, it is likely to have weak power in general, because of the many df. An alternative approach is to test whether the θ_k values significantly differ over the strata, because significant heterogeneity implies departure from HWE. To compute this type of statistic under the null hypothesis of homogeneity (yet allowing departure from HWE), one needs to estimate a common θ parameter. Using estimating equations, Olson and Foley¹⁰ derived a consistent estimator for θ, whereas Nam¹¹ used an iterative maximum-likelihood method. Both of these approaches, however, can run into undefined parameter estimates; Olson’s θ is undefined when there are no AA homozygotes across all strata (or no BB homozygotes); similar problems occur for the maximum-likelihood estimator. In these cases, the test for homogeneity breaks down.

Because of the above complications with large-sample statistical tests for HWE across strata or for homogeneity of departures from HWE across strata, exact methods are appealing. Instead of summing a measure of departure from HWE, an exact test evaluates the combined evidence over strata by considering the probability of genotype configurations when the null hypothesis is true; extreme departures in different directions are rare under HWE, giving a small P value, yet a sum of measures of departure from HWE could, in fact, completely miss this situation. Exact tests also avoid numerical problems (e.g., division by zero), and they provide appropriate control of the type I error rates. To date, only Olson and Foley¹⁰ considered exact methods. However, because their methods allowed for an arbitrary number of alleles, exact computations were not feasible. Rather, they needed to rely on Markov Chain–Monte Carlo methods. Because of the broad use of SNPs, we present efficient computational methods to compute exact tests both for HWE over strata and for homogeneity of departures of HWE over strata. We demonstrate our methods by applying them to the SNP genotype data from ^Perlegen,2 and ^HapMap.³ The results illustrate the advantages of our exact stratified test for HWE over the minimum exact-test P value or other approximate methods. Furthermore, simulations confirm our empirical findings.

Now, to extend these ideas to strata, let be the vector of observed counts of AB heterozygotes for the different strata, and let be a vector containing a configuration of possible values of heterozygotes for the different strata. The probability of under HWE is the product of expression (2) over the K strata, . This allows us to compute an exact stratified P value by

A naive approach to compute the exact P value is to evaluate all possible configurations, which is inefficient. Rather, we first compute P(x) for all possible values of x within each stratum. This avoids having to recompute P(x) many times. Using recursion makes this fast, and using log-probabilities avoids numerical imprecision. Because it is of critical importance, we order the log-probabilities such that we can stop the summation for the P values as soon as possible. To do this, we sort the log-probabilities into increasing order within each stratum, using quick sort. Then, we begin to evaluate different possible configurations by summing, across strata, the log-probabilities for values of x in the vector. If this sum is less than or equal to the log-probability of the observed data, then we exponentiate it and add it to the running sum for the P value. The prior sorting of the log-probabilities within each stratum allows us to skip over computations that would generate and, hence, would not contribute to the sum for the P value. This is explained in greater detail by use of an example in appendix A. An advantage of our approach is that, when small P values are used as a quality-control (QC) filter (e.g., P<.001 as used elsewhere^2,3), we can stop computations early when the computed P value exceeds a specified threshold.

Under the assumption of constant θ, the probability of an configuration is

To efficiently compute the exact P value, we first enumerate possible x_k values for each stratum (for now, ignoring the constraint that the x_k values must sum to the total number of observed heterozygotes) and use recursion to determine the contribution of the x_k values to . To see this, let q(x_k) denote the ratio of factorials for stratum k (i.e., ). It is easy to verify that

We applied our exact methods to SNPs in the ^Perlegen,2 and ^HapMap,3 data sets. For the Perlegen data, we used 1,585,674 SNPs from all chromosomes. For the X chromosome, we used males and females for the pseudoautosomal regions and only females for other regions on X. The Perlegen data have a total of 71 subjects from three ethnic groups: 23 African Americans, 24 European Americans, and 24 Han Chinese. Furthermore, the Perlegen data were “cleaned” by a number of criteria, including exact tests for HWE within each of the ethnic groups. SNPs were given a poor quality score if the smallest P value across the three ethnic groups was <.001. Hence, the Perlegen data are useful to evaluate whether combining information across strata detects significant departure from HWE that was missed by using the minimum P value.

To provide a more complete comparison of using the minimum P value versus the exact stratified P value, we also applied our exact methods to SNPs in the HapMap data, using only the autosomes. These data have a total of 210 independent subjects from four ethnic groups: 60 Yoruba from Ibadan, Nigeria; 60 U.S. residents with northern and western European ancestry (CEPH samples); 45 Han Chinese; and 45 Japanese. Note that the offspring of “trios” were not used in our tests of HWE.

For the HapMap data, instead of using the cleaned data, we used the “redundant-unfiltered” genotype data. This allows us to evaluate our methods on data that did not have genotypes removed because of prior tests of HWE within strata. For this data, various QC flags were used to indicate reasons why SNPs failed the QC criteria, including an exact test for HWE within each of the strata; P<.0001 in any of the strata was flagged as a failure. For our analyses, we did not eliminate SNPs that failed for this reason. Rather, we eliminated SNPs that failed the QC criteria for any reason not indicated by a HWE failure. For the duplicate samples, we coded a genotype as “missing” if the duplicates did not agree and then removed the duplicates for analyses. This resulted in the examination of 3,798,286 SNPs.

Figure 1. Perlegen data. Exact stratified test versus the minimum within-strata exact test; plot of .

Because the test for homogeneity of departures from HWE could also be used as a test for HWE, we plot in figure 2 the exact P values for the exact stratified test for HWE and the test for homogeneity. This illustrates that the lgP values for the exact stratified test tend to be larger than the lgP values for the homogeneity test, which suggests that using the homogeneity test as a way to test for HWE is not likely to be as powerful as using the stratified test. Olson and Foley¹⁰ also noted the weak power of the homogeneity test.

Figure 2. Perlegen data. Exact stratified test versus exact test of homogeneity; plot of .

We also applied our adaptation of Haldane’s stratified test for HWE, as well as Olson’s⁹ stratified test, and contrasted these with the exact stratified test. The results in figure 3 illustrate that Haldane’s and Olson’s stratified tests give nearly identical P values, suggesting that the different ways to weight the h values over strata has little impact in real applications. In figure 4, we compare the exact stratified test with our version of Haldane’s stratified test. This figure illustrates that the exact test and Haldane’s test can have large discrepancies. At one extreme, where the exact test gives large lgP values and Haldane’s test gives lgP values near 0, the cause was typically that the h values differed in sign over strata, causing them to cancel each other, to result in Haldane’s summary statistic to be near zero. In contrast, the exact test was able to detect these types of departure from HWE. At the other extreme, where the exact stratified test gave lgP values near 0 and Haldane’s test gave lgP values near 1, the summary Haldane statistics were typically negative, implying too few homozygotes. In these cases, only one type of homozygote was observed over all strata, yet the observed and HWE expected genotype counts were quite close. Furthermore, Haldane’s test tended to have more extreme values of lgP than did the exact test; 109 instances with lgP>12. In all of these cases, the strata had no heterozygotes and either one or two rare homozygotes. These results empirically emphasize the inadequacy of the normal distribution for Haldane’s stratified statistic when there are sparse data, leading to P values that are likely much too small.

Figure 3. Perlegen data. Haldane’s versus Olson’s stratified tests; plot of .

Figure 4. Perlegen data. Exact versus Haldane’s stratified tests: plot of .

Figure 5. HapMap data. Exact stratified test versus the minimum within-strata exact test; plot of .

Like the Perlegen results displayed in figure 2, we found for the HapMap data that the lgP values for the exact stratified test tended to be larger than the lgP values for the homogeneity test, emphasizing that using the homogeneity test as a way to test for HWE is not likely to be useful (results not shown). Also, for the HapMap data, Haldane’s and Olson’s stratified tests gave similar results (results not shown).

Table1. Simulation Type I Error Rates

Simulations for power were conducted for only P=.20 because all tests have the correct type I error rate for this situation, again restricted to 25 subjects per stratum. Results for power are presented in table 2, for when the departure from HWE is in the same direction and magnitude across all strata. In this case, the Haldane and Olson statistics had the greatest power, as expected, because they were derived under the assumption of constant departure from HWE across strata. However, the decreases in power for the other tests were generally small. Furthermore, the power was approximately the same for the exact test and the χ² statistic of Troendle and Yu.¹² Results for power when the departure from HWE differed over strata are presented in table 3. In this case, the exact test and Troendle and Yu’s statistic had similar power that was greater than that for Haldane’s and Olson’s tests.

Table2. Simulation Power When Departure from HWE is the Same for All Strata

Table3. Simulation Power When Departure from HWE Differs over Strata

Table4. Timing of hweStrata

Application of our exact stratified test for HWE to the ^Perlegen and ^HapMap data sets provides an empirical comparison of our new methods with the common approach that uses the minimum exact-test P value. Both data sets emphasize the greater power of the exact stratified test, which makes intuitive sense because it simultaneously evaluates HWE over all strata rather than independently testing each stratum. The exact stratified test also accounts for testing multiple strata; a Bonferroni correction would be needed to control the type I error rate when using the minimum P value over the strata.

Although a number of approximate stratified tests of HWE have been proposed, our applications illustrate that our version of Haldane’s test gives results nearly identical to those of Olson’s stratified test for HWE, suggesting that the different ways of weighting the contribution from each stratum do not have major influences on the tests. By comparing results from Haldane’s stratified test with those from the exact stratified test, we found that using the standard normal distribution to approximate the distribution of Haldane’s test can give exceptionally small P values when there are sparse genotype counts, which suggests that the normal distribution is not adequate and that the exact stratified test is more reliable. Finally, both the Perlegen and HapMap data illustrated that the exact stratified test for HWE is much more powerful to detect departures from HWE than the exact test for homogeneity, echoing the simulation results of Olson and Foley.¹⁰

Our simulation results confirmed that the exact stratified test provides the correct type I error rate, whereas the tests proposed by Haldane, Olson, and Troendle and Yu can have inflated type I error rates in the presence of sparse data (i.e., small strata and rare alleles). Furthermore, our simulations confirmed that the exact test provides the greatest power when the departure from HWE is in different directions across strata. Finally, the simulations suggest that when the strata sizes are not small and the frequency of the rare allele is at least 5%, the omnibus test of Troendle and Yu would be a good substitute for the exact stratified test.

Although our work was motivated by the relatively small ethnic groups within the Perlegen and HapMap data sets and by the potential for using these data sets for planning large-scale genome association studies in large “homogeneous” ethnic groups, our exact tests should prove useful for many genetic studies. Some examples are follow-up studies in multiple ethnic groups, each of which may not be large (note that association analyses would need to account for the different ethnic groups, such as a Mantel-Haenszel stratified analysis), population genetic studies in multiple ethnic groups from a geographic region, or studies in which apparently homogeneous ethnic groups can be clustered into smaller ethnic subsets on the basis of many measured markers.¹⁴

In conclusion, our simulations and the application of exact and approximate stratified tests for HWE to more than five million SNPs, with strata sizes ranging from 23 to 60 subjects, provide convincing results that the exact stratified test provides the most-robust and most-powerful results. Furthermore, efficient computational algorithms for SNP genotype data, which we developed in the C programming language, allow the exact stratified test to be computed within reasonable computing time for sample sizes on the order of the HapMap data (e.g., strata sizes ranging from 45 to 60 subjects over four strata). The C source code for our software, called “hweStrata,” is available from our Web site.

This work was supported by U.S. Public Health Service, National Institutes of Health, contract grants GM065450 and GM61388 (The Pharmacogenetics Research Network). The constructive criticisms from two anonymous reviewers helped to improve the evaluation and presentation of the proposed exact stratified test.

To illustrate how sorting the log-probabilities within each stratum leads to an efficient way to compute an exact P value, we present in table A1 a hypothetical example of genotypes for three strata. For this example, the total log-probability of the observed genotypes is −24.3649. Although there are 288 possible configurations for the possible number of heterozygotes over the strata, configurations that have log-probabilities larger than the observed value of −24.3649 do not contribute to the P value, so they can be skipped. The log-probabilities, sorted within each stratum, are illustrated in table A2. Note that summing the first log-probability across strata gives the log-probability of the most rare configuration. The indices for these log-probabilities are 1, 1, 1, with a resulting sum of −30.0905 (see table A3). Because this sum is smaller than the summed log-probabilities for the observed data, this configuration contributes to the P value. Increasing the index for the third stratum, up to its maximum value of 4, gives summed log-probabilities that are less than those for the observed data, so all these configurations contribute to the P value. When the index for the second stratum is increased to 2, we find that the indices 1, 2, 2, lead to a summed log-probability of −22.3117 (configuration 6 in table A3), which is larger than the observed value. Increasing the index for stratum 3 to 3 or 4 would only lead to larger summed log-probabilities, so these values can be “jumped” over. The next array of indices is 1, 3, 1, which contributes to the P value, but the following array, 1, 3, 2, is another “jump” array. The arrays that indicate the jumps are given in table A3. For this example, we need to evaluate only 13 configurations of the total 288 configurations.

To illustrate how we enumerate configurations that meet the constraints (1) that each x_k is restricted to a range of values determined by N_A,k (0–N_A,k if N_A,k is even and 1–N_A,k otherwise) and (2) that the elements in sum to the total number of observed heterozygotes, we present a simple example in table A4. This table is a slightly modified version of the data in table A1. For this example, the observed total number of heterozygotes is seven, so the maximum number of heterozygotes N_A,k cannot be achieved for any of the strata.

To enumerate possible configurations, we determine the minimum and maximum possible values of x for strata 1 and 2. The value for x₃ is x₃=7-x₁-x₂. We then decrease x₃ by increments of 2 and increase x₂ by increments of 2, keeping x₁ fixed. Once x₂ achieves its maximum value, we increment x₁ by 2, set x₂ to its minimum value, determine x₃=7-x₁-x₂, and then again decrement x₃ and increment x₂. This process continues until x₁ achieves its maximum value. This pattern is easily viewed below with the example data of possible configurations that sum to 7, the observed total number of heterozygotes.

TableA1.

Demonstration Genotype Data for Exact Stratified Test of HWE

Stratum	NAA	NAB	NBB	Log-Probability	NA	No. of Possible Heterozygotes
1	10	2	13	−11.0767	22	12
2	5	0	8	−8.3254	10	6
3	3	0	5	−4.9628	6	4
Total	18	2	26	−24.3649	…	…

TableA2.

Sorted Values of Log-Probabilities for Each Stratum

TableA3.

Sum of Log-Probabilities over Strata Configurations

TableA4.

Demonstration Genotype Data for Exact Test of Homogeneity of HWD

Stratum	NAA	NAB	NBB	NA	Range of Possible Heterozygotes
1	10	2	13	22	0–22
2	5	2	8	12	0–12
3	3	3	5	9	1–9
Total	18	7	26	…	…