New Equivalence Tests for Hardy–Weinberg Equilibrium and Multiple Alleles

: We consider testing equivalence to Hardy–Weinberg Equilibrium in case of multiple alleles. Two different test statistics are proposed for this test problem. The asymptotic distribution of the test statistics is derived. The corresponding tests can be carried out using asymptotic approximation. Alternatively, the variance of the test statistics can be estimated by the bootstrap method. The proposed tests are applied to three real data sets. The ﬁnite sample performance of the tests is studied by simulations, which are inspired by the real data sets.


Introduction
Hardy-Weinberg Equilibrium (HWE) plays an important role in the field of the population genetics and related scientific domains. HWE is a common assumption in many areas of research so that assessing the compatibility of observed genotype frequencies with HWE is a basic step of a complete statistical analysis. There are two main approaches to this undertaking: goodness of fit tests and equivalence tests.
A vast amount of literature exists on the goodness of fit tests for HWE, which includes application of the asymptotic χ 2 and likelihood ratio tests. The specific exact goodness of fit tests for HWE are developed in [1][2][3][4][5] among others. The null hypothesis of all these tests is that the underlying population is exactly in HWE. Hence, the goodness of fit tests are tailored to establish lack of compatibility with HWE.
The equivalence tests are appropriate to establish sufficiently good agreement of the observed genotype frequencies with HWE. The exact and approximate equivalence tests for the biallelic case are developed recently in [6][7][8]. To our best knowledge, there are not any equivalence tests for HWE and multiple alleles. Two different equivalence tests are developed in this paper for the case of multiple alleles. The tests can be carried out using the asymptotic approximation or bootstrap method.
A distribution of diploid genotypes at a k-allele locus can be represented as a lower triangular matrix p, where p (i, j) is the probability of the genotypes with alleles i and j. Let a (p) denote the allele distribution under p. The probability of the allele i under p can be calculated as a (p, i) = 1 2 ∑ k j=1 (p (i, j) + p (j, i)). If the population is in HWE, then the genotype distribution fulfills the conditions p(i, j) = 2a(p, i)a(p, j) for i < j and p(i, i) = a(p, i) 2 . Let e (p) denote the genotype distribution under the assumption of HWE, which is implied by the allele distribution a (p).
Euclidean distance l 2 (p, e(p)) can be considered a conditional distance between the genotype distributions p and e (p) under the joint allele distribution a (p). The equivalence test problem is then defined by H 0 = {l 2 (p, e (p)) ≥ ε} and H 1 = {l 2 (p, e (p)) < ε} (1) where ε is a tolerance parameter.

35
Let M denote the family of all possible genotype distributions at HWE. The minimum distance between p and M is defined by d (p, M) = min q∈M l 2 (p, q). The corresponding equivalence test problem is given by We observe the genotype frequencies p n of the sample size n. The natural test statistic for (1) is which can be easily computed. The appropriate test statistic for (2) is which requires optimization for the calculation of d (p n , M). The test statistic T c (p n ) can be considered a numerically efficient approximation to T m (p n ) because of l 2 (p n , e (p n )) ≥ d (p n , M). The subscript * will be used instead of c and m in the reminder of the paper, if statements are appropriate for both cases.
If Hypothesis (1) or (2) of the non-equivalence can be rejected for some appropriate value of ε then the true underlying genotype distribution is close to HWE with the probability greater than 1 − α, where α is the nominal level of the test. The appropriate value of ε depends on the application and the available sample size. The value of the parameter ε can be found by simulation as shown in Subsection 3.2. Alternatively, the minimum tolerance parameter ε, for which H 0 can be rejected, can be computed and reported, see Section 2 for details.

Equivalence tests
In this section, we derive the asymptotic distributions of the test statistics T c (p n ) and T m (p n ). We provide also an algorithm for the asymptotic and bootstrap-based tests.
Let v be the usual bijective mapping of the matrix p to the vector (p(1, 1), p (1, 2) , . . . , p (k, k)). Letd c denote the derivative of the function q → l 2 2 v −1 (q) , e v −1 (q) , where q is a vector of length k 2 . The derivatived c can be derived using the chain rule. Let p 0 ∈ H 0 fulfill the boundary condition l 2 (p 0 , e (p 0 )) = ε and let q 0 = v (p 0 ). Then the asymptotic distribution of T c (p n ) under p 0 is Gaussian with mean zero and variance σ 2 where Σ (q 0 ) = D q − qq t is a covariance matrix and D q is a square diagonal matrix, whose diagonal entries are elements of q. The proof of the statement can be found in [9].
The test statistic T m (p n ) converges weakly under the assumption that there exists a continuous function h on an open neighborhood of p 0 such that h (p) ∈ M and d (p, M) = l 2 (p, h (p)). The existence of a continuous minimizer h is also an important requirement for the numerical computation of d (p, M). We assume the existence of a continuous minimizer h on an open neighborhood of p 0 for the reminder of the paper. Letd m denote the derivative of the function q → l 2 2 v −1 (q) , h (p 0 ) . Then the asymptotic distribution of T m (p n ) under p 0 is Gaussian with mean zero and variance [10] for details.
The asymptotic variance σ 2 * (p 0 ) is unknown and can be estimated by σ 2 * (p n ). The asymptotic test can be carried out as follows: 1. Given are the genotype frequencies p n , the tolerance parameter ε and the significance level α. 2. Compute the tests statistic T * (p n ). 3. Estimate the asymptotic variance by σ 2 * (p n ).
The minimum tolerance parameter ε, for which the asymptotic test can reject H 0 , can be computed as To improve the finite sample performance of the proposed tests, the bootstrap method is applied to estimate the variance of T * (p n ), see [11], Section 6 for details. The estimator σ * (p n ) is then replaced by the bootstrap estimator of the variance. Otherwise, everything stays the same.

Simulation study
The proposed tests are implemented in R and are freely available on GitHub under https://github. com/TestingEquivalence/HardyWeinbergEquilibriumR. All simulations are performed in R-Studio on a usual scientific workstation.

Real data sets
The equivalence tests are applied to the following data sets, which are already analyzed in the literature on the goodness of fit tests: 1. from rheumatoid arthritis study, [12]; 2. from the documentation included with the GENEPOP software package, [13]; 3. genotype frequency data at Rhesus locus, [14]. The genotype distributions of the data sets are given in Table 1.
The minimum tolerance parameters ε, for which the tests can reject the corresponding H 0 , are displayed in Table 2. The distances d (p n , M) and l 2 (p n , e (p n )) are close to each other in all cases so that l 2 (p n , e (p n )) provides a good approximation to d (p n , M). The test results are also similar for T c and T m . The bootstrap tests are slightly more conservative than the asymptotic tests in all cases.
It could not be shown that data sets 1 and 2 are close to HWE. All goodness of fit tests in [15] reject also the null hypothesis of HWE for data set 2 at the nominal level 0.05. Data set 3 is very close to HWE. This observation corresponds to the results of the goodness of fit tests in [5,15]. Table 1. The data sets: 1) from rheumatoid arthritis study, [12]; 2) from the documentation included with the GENEPOP software package, [13]; 3) genotype frequency data at Rhesus locus, [14].

Test power
In this subsection we study the test power at HWE. We restrict yourself to the genotype distributions at HWE, which are implied by the real data sets from Subsection 3.1, because the family M is very large. To shed some light on the appropriate values of the tolerance parameter ε, the test power is computed for different values of ε, see Table 3. The value of ε may be considered appropriate if the test power is approximately 0.9. Hence, the appropriate value of ε is 0.1 for data set 1, 0.1 for data set 2 and 0.018 for data set 3.
The observed genotype frequencies p n are subjected to the sampling error. It is important for the test efficiency that the sampling error has a small influence on the test power at HWE. The test power is computed at 100 random genotype frequencies, where the corresponding random samples of size n are drawn from the implied genotype distribution e (p n ). The simulation results are summarized in Table 4. The power of all considered tests varies little from point to point. Hence, the impact of the sampling error on the test power at HWE is very small. Table 3. Simulated rejection probability of the equivalence tests at the nominal level 0.05. The rejection probability is simulated for different values of the tolerance parameter ε at the HWE distributions e (p n ), which are implied by data sets 1, 2 and 3. The sample size equals the size of the corresponding data set. The number of replications is 1000 for each experiment. A stands for the asymptotic test and B stands for the bootstrap test.

Type I error
We study the type I error rates of the proposed tests in this subsection. The boundary of H 0 is so complex that it is very difficult to find boundary points, which have the largest rejection probability. We consider therefore randomly selected boundary points of H 0 , which are based on the three real data sets from Subsection 3.1. The boundary points are generated using the following algorithm: 1. Given are p n and ε. 2. Draw a sample of size n from p n and compute the sample genotype frequency p n . 3. If T * ( p n ) < 0 then reject p n and repeat step 2. Otherwise accept p n . 4. Consider the linear combination a p n + (1 − a) e (p n ) for a ∈ [0, 1]. Find a n ∈ [0, 1] such that T * (a n p n + (1 − a n ) e (p n )) = 0. The value of a n can be found using any line search method. 5. Return a n p n + (1 − a n ) e (p n ), which is a random boundary point of H 0 .
The tolerance parameter ε is close to l 2 (p n , e (p n )) for each data set under consideration so that a n is usually not far from 1. The corresponding random boundary point is then close to p n . Hence, we explore the boundary of H 0 in the neighborhood of the given data set. The test power at 100 random boundary points is summarized in Table 5. The test power varies considerable from point to point. The asymptotic test based on T c is not conservative for all three data sets. The asymptotic test based on T m shows some anti-conservative tendencies for data sets 2 and 3. The bootstrap test based on T c is conservative for all three data sets. The bootstrap test based on T m shows slight non conservative tendencies.
The power at the boundary points is larger than the nominal level due to the following reasons. If the number of observations np n (i, j) is too small for some i and j then the distribution of T * (p n ) may be far away from the normal approximation and also may have considerable jumps. The critical values of the asymptotic and bootstrap tests are then incorrect. If the vectord * (v (p n )) contains zero elements then the power of the asymptotic tests tends to be above the nominal level α. The power of the bootstrap tests is closer to α in this case because the vectord * (v (p n )) is not used for the variance estimation by the bootstrap method.

Summary
Two different test statistics are proposed to establish equivalence of the genotype distributions to HWE. The critical values of the tests are calculated using the asymptotic approximation by the normal distribution. The variance of the test statistic is estimated asymptotically or by the bootstrap method. The minimum tolerance parameter ε, for which H 0 can be rejected, is derived. The tests are successfully applied to three real data sets, which are frequently considered in the literature. The test power at HWE and the type I error rates are studied at a large number of points, which are inspired by the real data sets. The asymptotic tests have anti-conservative tendencies and should be used with caution. The bootstrap-based tests are sufficiently conservative for the most practical situations. If more conservative tests are required then the nominal level may be halved or the tolerance parameter ε may be reduced. We recommend to perform all proposed tests in any case and compare the results. The appropriate value of ε depends on the application and the available sample size. The reasonable values of the parameter ε can be found by simulation as shown in Subsection 3.2. Additionally, the rejection probabilities at the close random boundary points may be studied as shown in Subsection 3.3.
Funding: This research received no external funding. The author declare no conflict of interest.