Hardy-Weinberg Equilibrium
Hardy-Weinberg equilibrium is the condition that the allele and genotypes frequencies in a population remain constant between generations, absent of other evolutionary influences. Another way to express this statement is that in a randomly mating (large) population with no selection or mutation, if the frequency of an allele \(D\) is the proportion \(\rho\) then the expected frequency of heterozygous individuals (those with a \(D\) and \(R\) allele) is \(2\rho(1-\rho)\). To check whether a population is in Hardy-Weinberg equilibrium, we can compare this expected allele proportion to the observed proportion with heterozygous alleles. Therefore, we want to compute \(\gamma = \pi - 2\rho(1-\rho)\), where \(\pi\) is the observed proportion. While \(\delta\) is easy to compute, the variance is less straightforward. In general, one needs to apply the delta method to combine the variance estimates for \(\rho\) and \(\pi\), which involves computing the derivatives. An alternative is estimating equations, which automates the delta method and provides a variance estimate without needing to do any by-hand calculations.
To illustrate how estimating equations simplify estimation of whether a population is in Hardy-Weinberg equilibrium, we provide an application using data collected by AS Wiener on the M-N blood group locus among white individuals in New York City white between 1931 and 1969. Here, we consider whether the M-N blood group is in Hardy-Weinberg equilibrium.
Setup
[2]:
import numpy as np
import scipy as sp
import pandas as pd
import delicatessen as deli
from delicatessen import MEstimator
print("Versions")
print("NumPy: ", np.__version__)
print("SciPy: ", sp.__version__)
print("Pandas: ", pd.__version__)
print("Delicatessen:", deli.__version__)
Versions
NumPy: 2.3.5
SciPy: 1.16.3
Pandas: 2.3.3
Delicatessen: 4.2
The data used here can be found in Table 4.3 of The Evaluation of Forensic DNA Evidence. This tabled data is reproduced below and converted into individual-level data
[12]:
# Table 4.3 — Wiener data
table = pd.DataFrame()
table['S'] = [1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6]
table['A'] = [1., 0.5, 0.] * 6
table['H'] = [0, 1, 0] * 6
table['n'] = [71, 116, 49, 132, 232, 97, 166, 289, 127,
1037, 1623, 608, 287, 481, 186, 158, 249, 93]
assert 6001 == np.sum(table['n'])
# Expanding to individual level
d = pd.DataFrame(np.repeat(table.values, table['n'], axis=0), columns=table.columns)
Overall
Using this data, we can then determine whether Hardy-Weinberg equilibrium holds. To begin, we will ignore that these data come from different convenience samples. The estimating equations for testing whether a population is in Hardy-Weinberg equilibrium are
where \(A_i\) is the allele frequency for individual \(i\) and \(H_i\) is an indicator whether an individual is heterozygous. The final estimating function is for the Hardy-Weinberg equilibrium comparison, which follows from some algebra. The code provided below implements these estimating functions with delicatessen
[8]:
alleles = np.asarray(d['A'])
heterozygote = np.asarray(d['H'])
[9]:
def psi_delta(theta):
# Dividing up input parameters
rho, pi, gamma = theta
# Estimating function for mean of MM group
ee_p = alleles - rho
# Estimating function for proportion of MN group
ee_h = heterozygote - pi
# Estimating function for Hardy-Weinberg equilibrium
ef_d = pi - 2*rho*(1-rho) - gamma * np.ones(alleles.shape[0])
# Stacking the estimating functions together
return np.vstack([ee_p, ee_h, ef_d])
[18]:
estr = MEstimator(psi_delta, init=[0.5, 0.5, 0.])
estr.estimate()
[19]:
estr.print_results(decimals=3)
==============================================================
Estimation Method: M-estimator
--------------------------------------------------------------
No. Observations: 6001 | No. Parameters: 3
Solving algorithm: lm | Max Iterations: 5000
Solving tolerance: 1e-09 | Allow P-Inverse: 1
Derivative Method: approx | Deriv Approx: 1e-09
Small N Correction: None | Distribution: Z-stat
==============================================================
Theta StdErr Z-score LCL UCL P-value S-value
--------------------------------------------------------------
0.558 0.005 123.599 0.549 0.566 0.000 inf
0.498 0.006 77.196 0.486 0.511 0.000 inf
0.005 0.006 0.766 -0.008 0.017 0.443 1.173
==============================================================
Here, the last coefficient is near zero which indicates that the population is in Hardy-Weinberg equilibrium (provided the other assumptions for this conceptual model all hold). From the output, we can see we get confidence intervals and a corresponding P-value for this directly. This example highlights how estimating equations simplify statistical inference for transformation of different parameters.
Strata-Specific
Given these results, we might also be interested in whether Hardy-Weinberg equilibrium holds for all the convenience samples. The following code applies a stratified version of the previous estimator. For stratification, we only need to multiply by an indicator of the strata. Effectively, this will zero out the contribution of any individuals outside of the strata. The following code applies this process
[22]:
def psi_delta_s(theta):
# Setup items for a loop through the estimating functions
start, end = 0, 3
est_funcs = []
# Loop over the specific subgroups
for s in [1, 2, 3, 4, 5, 6]:
# Indicator if in corresponding strata
strata = np.where(d['S'] == s, 1, 0)
# Estimating function for corresponding parameters and strata
theta_s = theta[start: end]
ef_s = psi_delta(theta=theta_s) * strata
# Storage and updating for next steps
est_funcs.append(ef_s)
start = start + 3
end = end + 3
# Stacking the estimating functions together
return np.vstack(est_funcs)
[23]:
estr = MEstimator(psi_delta_s, init=[0.5, 0.5, 0.]*6)
estr.estimate()
[25]:
estr.print_results(subset=[2, 5, 8, 11, 14, -1], decimals=3)
==============================================================
Estimation Method: M-estimator
--------------------------------------------------------------
No. Observations: 6001 | No. Parameters: 18
Solving algorithm: lm | Max Iterations: 5000
Solving tolerance: 1e-09 | Allow P-Inverse: 1
Derivative Method: approx | Deriv Approx: 1e-09
Small N Correction: None | Distribution: Z-stat
==============================================================
Theta StdErr Z-score LCL UCL P-value S-value
--------------------------------------------------------------
-0.004 0.032 -0.128 -0.067 0.059 0.898 0.155
0.006 0.023 0.265 -0.039 0.052 0.791 0.338
-0.001 0.021 -0.058 -0.042 0.039 0.954 0.068
0.005 0.009 0.611 -0.012 0.022 0.541 0.886
0.010 0.016 0.612 -0.022 0.041 0.540 0.888
0.006 0.022 0.294 -0.037 0.050 0.769 0.379
==============================================================
Again these strata-specific results indicate that Hardy-Weinberg equilibrium holds. There is a challenge here. Namely, that reporting the confidence intervals for this full collection of parameters is not quite appropriate. If we consider statistical inference for all six parameters, then the reported confidence intervals will be too narrow (i.e., coverage will be below the implied 95%). To correct for this, we can instead compute the confidence bands. The following code computes the corresponding confidence bands for these six parameters
[26]:
estr.confidence_bands(subset=[2, 5, 8, 11, 14, -1], seed=198041)
[26]:
array([[-0.08903632, 0.08077719],
[-0.0547735 , 0.06704523],
[-0.05548282, 0.05310035],
[-0.0173577 , 0.02785834],
[-0.03229401, 0.05188822],
[-0.05136222, 0.06426222]])
As expected, the confidence bands are wider. Overall our conclusions of Hardy-Weinberg equilibrium remain the same.
References
Hardy, G. H. (1908). Mendelian proportions in a mixed population. Science, 28(706), 49-50.
National Research Council (US) Committee on DNA Forensic Science: An Update. The Evaluation of Forensic DNA Evidence. Washington (DC): National Academies Press (US); 1996. 4, Population Genetics. Available from: https://www.ncbi.nlm.nih.gov/books/NBK232608/
Weir, B. S. (1996). Genetic data analysis II. Sunderland. MA: Sinauer Associates, 161-173.