Linear Mixed-Effect Models Through the Lens of Hardy–Weinberg Disequilibrium

For genetic association studies with related individuals, the linear mixed-effect model is the most commonly used method. In this report, we show that contrary to the popular belief, this standard method can be sensitive to departure from Hardy–Weinberg equilibrium (i.e., Hardy–Weinberg disequilibrium) at the causal SNPs in two ways. First, when the trait heritability is treated as a nuisance parameter, although the association test has correct type I error control, the resulting heritability estimate can be biased, often upward, in the presence of Hardy–Weinberg disequilibrium. Second, if the true heritability is used in the linear mixed-effect model, then the corresponding association test can be biased in the presence of Hardy–Weinberg disequilibrium. We provide some analytical insights along with supporting empirical results from simulation and application studies.

1 Introduction

Genetic association tests are often derived from a regression model, regressing the phenotypic data of a complex trait (Y) on the genotypic data of a single-nucleotide polymorphism (SNP; G), as well as on the covariate data of important environmental factors (Z). When individuals in a sample are genetically related with each other, the linear mixed-effect model (LMM) is the most commonly used method for genome-wide association studies (GWAS) (Eu-Ahsunthornwattana et al., 2014). The variance–covariance matrix of the regression model is partitioned into a weighted sum of the genetic correlation matrix and the correlation matrix due to shared environmental effects. The genetic correlation matrix is typically represented by the kinship coefficient matrix, which is either inferred from the (correctly) known pedigree structure or estimated based on the available genome-wide genetic data (Yang et al., 2011; Dimitromanolakis et al., 2019). The weight for the genetic correlation matrix is referred to as the heritability of the trait (Visscher et al., 2006; Visscher et al., 2008); Falconer (1985) gave a theoretical modeling of the variance partition, which sets the foundation for heritability.

It is commonly assumed that these regression-based association tests are robust to departure from Hardy–Weinberg equilibrium (HWE) (Sasieni, 1997). HWE states that the two alleles in a genotype are independent draws from the same Bernoulli distribution, or, equivalently, genotype frequencies depend solely on the allele frequencies (Hardy et al., 1908; Weinberg, 1908). For a biallelic SNP with two possible alleles A and a, let p and 1 − p be the population allele frequencies, respectively. Under HWE, p_aa = (1 − p)², p_Aa = 2p(1 − p), and p_AA = p², where p_aa, p_Aa, and p_AA are the population genotype frequencies of genotypes aa, Aa, and AA, respectively. To quantify the departure from HWE or the amount of Hardy–Weinberg disequilibrium (HWD),

$\delta =p_{AA}-p^2(1)$

is a widely used measure (Weir, 1996), and δ = 0 indicates HWE holds. We note that a) HWE is also known as Hardy–Weinberg proportion and b) δ is also known as p(1 − p)F, where F is the inbreeding coefficient (Powell et al., 2010). Equivalently, instead of quantifying the genotype frequencies as p_aa = (1 − p)² + δ, p_Aa = 2p(1 − p) − 2δ, and p_AA = p² + δ based on δ (Weir, 1996), we can define them based on F as p_aa = (1 − p)² + p(1 − p)F, p_Aa = 2p(1 − p)(1 − F), and p_AA = p² + p(1 − p)F (Powell et al., 2010). As the classical Pearson χ² HWE testing is based on comparing the observed genotype counts with the expected under HWE (Zhang and Sun, 2021); we thus chose δ for this work to be consistent with the GWAS literature.

A truly associated or causal SNP can be out of HWE (Wittke-Thompson et al., 2005; Ryckman and Williams, 2008; Turner et al., 2011), which is often overlooked but an important consideration when studying a method’s robustness to HWD. Note that the HWD attributed to true association is typically not as extreme as the HWD caused by genotyping errors (Zhang and Sun, 2020). Thus, true HWD can remain in a “cleaned” dataset after applying the standard HWD-based quality control screening using a stringent p-value threshold [e.g., 10^–12 for an application of the UK Biobank data by Bycroft et al. (2018)]. With a sample of independent individuals, both theoretical and empirical results support that genotype-based association tests are robust to HWD (Sasieni, 1997; Schaid and Jacobsen, 1999; Zhang, 2021). However, in the presence of sample dependency, little has been discussed.

In this report, we first provide some analytical insights on why the standard LMM can be sensitive to HWD in pedigree data in contrast to when analyzing a sample of unrelated individuals. We then demonstrate with a simple sib-pair design that 1) when the heritability is estimated from the data as in practice, although the empirical type I error rate of the LMM is well controlled, the estimated heritability is biased, often upward biased; 2) when the true heritability is known and used, the empirical type I error rate of the LMM is then inflated when δ > 0, and deflated if δ < 0. The result of 2) is novel, but it is mostly of an academic interest as the true heritability of a trait is often unknown in practice. On the other hand, the result of 1) has important practical implications because if the estimate of a trait heritability is larger than the true value, then it helps explain some of the “missing heritability” (Manolio et al., 2009); the insightful work of Chen (2014) “discuss[es] the circumstances in which the HE [Haseman‐Elston] regression and the mixed linear model are equivalent.”

2 Methods

2.1 Traditional Y ∼ G Model With Independent Samples, T _Indep, Is Robust to HWD

Let Y be a (continuous) trait of interest, and G = 0, 1, and 2, respectively, for the genotypes aa, Aa, and AA of a SNP. Additionally, for notation simplicity but without loss of generality, we assume that there is only one additional covariate, denoted by Z. With a sample of n unrelated individuals, the traditional genotype-based association analysis assumes that

$y={\alpha }^{\ast }1+{\beta }^{\ast }\text{g}+{\gamma }^{\ast }\text{z}+{ϵ}^{\ast },{ϵ}^{\ast }\sim N\left(\text{0},{\sigma }^{\ast 2}I\right),(2)$

where y = (y₁, y₂,…, y_n) is a n × 1 vector for the phenotypic values, 1 is a n × 1 vector of 1’s, g = (g₁, g₂,…, g_n) is a n × 1 vector for the genotypes of the SNP, z = (z₁, z₂,…, z_n) is a n × 1 vector for the covariate values, ϵ* is the error term with variance σ*², and I is the identity matrix.

Score-based tests are often used for genetic association analyses (Derkach et al., 2015). In this case, the score statistic of testing H₀: β* = 0 can be easily derived as

$T_{\text{\hspace{0.17em}indep}}=n\cdot \frac{{\left\{{\left(g-\bar{g}1\right)}^T\left(\text{y}-\bar{y}1\right)-\frac{{\left(g-\bar{g}1\right)}^T\left(\text{z}-\bar{z}1\right){\left(\text{y}-\bar{y}1\right)}^T\left(\text{z}-\bar{z}1\right)}{{\left(\text{z}-\bar{z}1\right)}^T\left(\text{z}-\bar{z}1\right)}\right\}}^2}{\left[{\left(g-\bar{g}1\right)}^T\left(g-\bar{g}1\right)-\frac{{\left\{{\left(g-\bar{g}1\right)}^T\left(\text{z}-\bar{z}1\right)\right\}}^2}{{\left(\text{z}-\bar{z}1\right)}^T\left(\text{z}-\bar{z}1\right)}\right]\left[{\left(\text{y}-\bar{y}1\right)}^T\left(\text{y}-\bar{y}1\right)-\frac{{\left\{{\left(\text{y}-\bar{y}1\right)}^T\left(\text{z}-\bar{z}1\right)\right\}}^2}{{\left(\text{z}-\bar{z}1\right)}^T\left(\text{z}-\bar{z}1\right)}\right]}.(3)$

To observe T _indep’s connection with Hardy–Weinberg disequilibrium, it is instructive to employ some algebraic tricks and show that

$\frac{1}{n}{\left(g-\bar{g}1\right)}^T\left(g-\bar{g}1\right)=\widehat{\text{var}}\left(G\right)=2\left(\widehat{p}\left(1-\widehat{p}\right)+\left({\widehat{p}}_{AA}-{\widehat{p}}^2\right)\right)=2\left(\widehat{p}\left(1-\widehat{p}\right)+\widehat{\delta }\right).$

Because $\widehat{\delta }={\widehat{p}}_{AA}-{\widehat{p}}^2$ measures the amount of HWD present in the data (Weir, 1996), T _indep inherently adjusts for departure from HWE through $\widehat{\text{var}}(G)=2(\widehat{p}(1-\widehat{p})+\widehat{\delta })$. As a result, the traditional genotype-based association test is robust to HWD in independent samples.

When Y is binary, the classic logistic regression is commonly used. However, Chen (1983) showed that under some regularity conditions, the score test statistics have an identical form for the exponential family in independent samples, which was recently validated by Zhang and Sun (2021) for genetic association studies. Additionally, Derkach et al. (2015) showed that for Y-dependent sampling, “the score statistics are identical for conditional and full likelihood approaches, and are of the same form as those for ordinary random sampling.” Thus, in terms of association testing (not genetic effect estimation), we can conclude that genotype-based association studies of binary traits in independent samples are also robust to HWD.

2.2 Linear Mixed-Effect Model With Dependent Samples, T _LMM, Can Be Sensitive to HWD

Although a pedigree-based study design is rare for genome-wide association studies, individuals can be (cryptically) related with each other even in population-based GWAS (Sun et al., 2017). Omitting related individuals simplifies the association analysis but reduces the sample size and thus power. Instead, Σ_Φ, the kinship coefficient matrix, can be estimated using the available genome-wide data to capture the sample relatedness between the n individuals (Visscher et al., 2006; Yang et al., 2011). The association analysis using the full sample can be conducted using the linear mixed-effect model.

$\text{y}={\alpha }^{\ast }1+{\beta }^{\ast }\text{g}+{\gamma }^{\ast }\text{z}+{ϵ}^{\ast },\text{\hspace{0.17em}where\hspace{0.17em}}{ϵ}^{\ast }\sim N\left(\text{0},{\sigma }_y^2{\mathrm{\Sigma }}_y\right)\text{and}{\mathrm{\Sigma }}_y=h^2{\mathrm{\Sigma }}_{\mathrm{\Phi }}+\left(1-h^2\right)I.(4)$

Compared with the linear model used for independent samples, var(ϵ*) = σ*²I in Eq. 2 is replaced by ${\sigma }_y^2{\mathrm{\Sigma }}_y$ to reflect the sample dependence. The matrix Σ_y is a weighted average of two components, where Σ_Φ reflects the sample relatedness; naturally, the model is reduced to the linear model of Eq. 2 for independent samples when Σ_Φ = I. The weight h² is interpreted as the heritability of the trait (Visscher et al., 2008), $h^2{\sigma }_y^2$ as the phenotypic variation due to (additive) genetic variation, and $(1-h^2){\sigma }_y^2={\sigma }_e^2$ as the phenotypic variation due to environmental variation. The matrix Σ_Φ is the kinship matrix, where Σ_Φ(i, j) = 2ϕ_i,j and ϕ_i,j is the kinship coefficient between the ith and jth samples.

By convention, h² is defined as

$h^2=\frac{\sum \underset{k}{\overset{2}{\beta }}\text{var}\left(G_k\right)}{\text{var}\left(Y\right)}=\frac{\sum \underset{k}{\overset{2}{\beta }}2p_k\left(1-p_k\right)}{\sum \underset{k}{\overset{2}{\beta }}2p_k\left(1-p_k\right)+{\sigma }_e^2},$

where there could be multiple causal SNPs, k = 1,…, S. In reality, h² is estimated by the correlation between phenotypes of related individuals. Consider the simple case of sibling pairs, and let Y₁ and Y₂ be the phenotypes for sib 1 and sib 2, respectively. Allowing for HWD and adjusting for the kinship coefficient ϕ, the estimated h² is

${\widehat{h}}^2=\widehat{\text{corr}}\left(Y_1,Y_2\right)/2\varphi ,$

where corr(Y₁, Y₂) depends on the correlation between G_1k and G_2k between the siblings; see Zhang and Sun (2021) for the derivation of corr(G_1k, G_2k) accounting for kinship coefficient and HWD. Thus,

$\frac{E\left({\widehat{h}}^2\right)}{h^2}=\frac{\sum \underset{k}{\overset{2}{\beta }}2\left(p_k\left(1-p_k\right)+{\delta }_k\right)}{\sum \underset{k}{\overset{2}{\beta }}2p_k\left(1-p_k\right)},$

and the bias of the h² estimate is

$E\left({\widehat{h}}^2\right)-h^2=h^2\cdot \frac{\sum \underset{k}{\overset{2}{\beta }}{\delta }_k}{\sum \underset{k}{\overset{2}{\beta }}{p}_k\left(1-p_k\right)}.(5)$

Under the simple case of one causal SNP, the bias is simplified to h² ⋅ δ/(p(1 − p)).

Given the analytical insights provided so far, we then briefly examine the empirical properties of T _LMM through both application and simulation studies.

3 Results

3.1 Cystic Fibrosis Sib-Pair Data Application: T_LMM Has Correct Type I Error but h² Appears to Be Overestimated

We extracted 65 sibling pairs from a cystic fibrosis (CF) gene modifier study (Wright et al., 2011; Sun et al., 2012). The phenotype Y of interest is the lung function measurements of the 130 related individuals with CF. In total, there were 570,539 SNPs genotyped using the Illumina 610-Quad Beadchip after applying the standard quality control, including minor allele frequency (MAF) greater than 2%. To stabilize the variance estimation, we additionally required SNPs to have MAF greater than 5%. We then applied T _LMM to the remaining 505,172 SNPs. In the application, we treated h² as unknown and estimated it based on the linear mixed-effect model of Eq. 4 as in convention.

When h² was estimated from the data, our association testing based on T _LMM had good type I error control (results not shown), consistent with the empirical observations in the GWAS literature. However, the estimated h², obtained using the 65-pair sibling data, is ${\widehat{h}}^2=0.82$. This value is substantially greater than 0.5, the commonly believed “true” heritability of lung function in CF obtained from the classic monozygous (MZ) vs. dizygous (DZ) twin-based estimation method (Vanscoy et al., 2007).

To verify if the large heritability estimate from the LMM method in our application was due to chance, we conducted a proof-of-principle simulation study. We assumed that only one causal SNP, G_causal with MAF of 0.2, affects Y with h² = 0.5. Genotype and phenotype values for 65 sibling pairs were then simulated under the assumption of HWE (i.e., without HWD). Among the 100,000 independently simulated replicates, only 4.24% of the heritability estimates were greater than ${\widehat{h}}^2=0.82$. This suggests that ${\widehat{h}}^2=0.82$, the value that was observed in the CF data application, was unlikely if the true heritability was 0.5 and without HWD at the causal SNP.

To verify if HWD at the causal SNP can lead to a biased heritability estimate, we then conducted additional simulation studies, following the same sib-pair design as mentioned previously. Our goal is to demonstrate that 1) when h² is treated as a nuisance parameter, its estimate based on model (Eq. 4) cross-reference can be biased in the presence of HWD; and 2) assuming the true h² is known, the empirical type I error rate of LMM (Eq. 4) cross-reference inflates when δ > 0 and deflates when δ < 0.

3.2 Simulated Sib-Pair Data in the Presence of HWD: h² Estimate Is Biased

Consider a continuous trait Y with h² = 0.5 and influenced by one causal SNP, G_causal, with minor allele frequency of 0.2 and with HWD factor, δ_causal, ranging from −0.04 to 0.16. A non-associated SNP, G_tested, also has an MAF of 0.2 but with its own δ_tested, which may not be the same as δ _causal in a specific simulation study. The sample size was 65 sibling pairs, chosen to match with the sample size of the cystic fibrosis application study in Section 3.1.

Most practical implementations of the linear mixed-effect model (Eq. 4) cross-reference treat h² as a nuisance parameter, and no type I error issue has been reported. Indeed, when h² was estimated in our simulation study conducted in Section 3.3, the test size of T _LMM was correct at the nominal level (black squares in Figure 3 shown in Section 3.3) even if δ_tested ≠ 0 (i.e., out of HWE) and across the range of δ_causal values (from −0.04 to 0.16).

However, in this situation, when h² is treated as unknown, the impact of HWD is on the estimation of h². Specifically, Figure 1 shows that ${\widehat{h}}^2$ is downward biased when δ_causal < 0, and upward biased if δ_causal > 0. The bias can be substantial. For example, when δ _causal = 0.10, the estimated heritability ${\widehat{h}}^2$ is centered at 0.80 as compared to the true value of 0.5, with a bias of 0.30. Indeed, based on our theoretical insight in Section 2.2, the expected bias is h² ⋅ δ/(p(1 − p)) = 0.5 ⋅ 0.1/(0.2(1–0.2)) = 0.31.

FIGURE 1

FIGURE 1. Box plots of ${\widehat{h}}^2$, estimated from the linear mixed-effect model (Eq. 4) against δ_causal. The true heritability of the phenotype is h² = 0.5. The minor allele frequencies p_causal = p_tested = 0.2 and 10,000 independent replicates of phenotypes and genotypes for 65 sibling pairs were simulated for each δ_causal value. The empirical type 1 error rates are shown in Figure 3 as black squares.

In Figure 1, it is notable that ${\widehat{h}}^2$ can be greater than one. Since h² is the proportion of variance in Y explained by additive genetic variation, 0 ≤ h² ≤ 1 by definition. However, if δ_causal ≠ 0, ${\widehat{h}}^2$ based on the LMM, without additional truncation, is a biased estimate of h² with a bias of h² ⋅ δ/(p(1 − p)) for this sib-pair design as shown in Section 2.2; the bias is 0 (i.e., no bias) under HWE when δ = 0.

Additionally, although a larger sample that consists of 5,000 sibling pairs shrinks the variance of the h² estimate as expected, it does not shrink the bias, as shown in Figure 2. However, we also note that, in practice, it is unlikely to have so many sibling pairs.

FIGURE 2

FIGURE 2. Box plots of ${\widehat{h}}^2$, estimated from the linear mixed-effect model (Eq. 4) against δ_causal using 5,000 sibling pairs. The other parameter values are the same as in Figure 1, where the true heritability of the phenotype is h² = 0.5, the minor allele frequencies p_causal = p_tested = 0.2, and 10,000 independent replicates were simulated for each δ_causal value.

3.3 Simulated Sib-Pair Data in the Presence of HWD: When Using the True h² Value T _LMM Has Incorrect Test Size

Here, we conducted the association analysis between Y and the non-associated SNP, G_tested, using the LMM model of Eq. 4 but assuming h² = 0.5 is known.

Figure 3A plots the empirical type I error rates (blue circles) of T _LMM using the true h² = 0.5, for a nominal level of 0.05, estimated from independently simulated 10,000 replicates for each δ_causal value. (An empirical type I error greater than 0.05 + 3 ⋅ 0.002 = 0.056 can be considered inflated as the standard error of the empirical type I error rate can be estimated as $\sqrt{0.05\cdot 0.95/10000}=0.002$.) In Figure 3A, the trend of type I error inflation is clear as δ_causal increases.

FIGURE 3

FIGURE 3. Empirical type I error rate of T _LMM based on the linear mixed-effect model (Eq. 4) against δ_causal. (A) When G_tested of tested SNPs is in HWD with δ_tested = 0.06. (B) When G_tested of tested SNPs is in HWE with δ_tested = 0. The true heritability of the phenotype is h² = 0.5, the minor allele frequencies p_causal = p_tested = 0.2, and 10,000 independent replicates of phenotypes and genotypes for 65 sibling pairs were simulated for each δ_causal value. The blue circles are for T _LMM using the true heritability h² = 0.5, and the black squares are for T _LMM while estimating h² (results of ${\widehat{h}}^2$ are shown in Figure 1).

In Figure 3A, we set δ_tested = 0.06, but we note that the main cause of the type I error issue is δ_causal ≠ 0 when using the LMM of Eq. 4 with h² = 0.5 plugged in. Indeed, Figure 3B shows that even if G_tested is in HWE (i.e., δ_tested = 0), the problem remains, albeit less severe, as long as δ_causal ≠ 0.

4 Discussion

We used a sib-pair design to demonstrate that the linear mixed-effect model can be problematic in the presence of Hardy–Weinberg disequilibrium at the causal SNP(s). To demonstrate that the LMM-based heritability estimate can be biased, as a proof-of-principle, our simulation study assumed that the phenotype Y has only one causal SNP, which is unrealistic for complex traits. However, the analytical insight shown in Eq. 5 (i.e., bias expected to be $h^2\cdot \sum {\beta }_k^2{\delta }_k/\sum {\beta }_k^2{p}_k(1-p_k)$) suggests that the issue discussed here remains relevant in the case of multiple causal SNPs as $\sum {\beta }_k^2{\delta }_k$ is unlikely to be zero, even while allowing the signs of δ_k to differ.

Assuming the true heritability h² is known, we also demonstrated the potential type I error issue of the LMM in the presence of HWD using data that consist of related individuals only. In practice, this issue diminishes if the sample includes a large number of independent individuals or the magnitude of HWD at the causal SNP is small. Additionally, in practice, h² is treated as unknown, in which case, the type I error rate of the LMM is well controlled; indeed, no increased false positives of the LMM due to HWD have been reported in the literature to the best of our knowledge. However, the estimate of h² can be upward biased and upwardly so if δ_causal > 0, as demonstrated in the simulation study in Section 3.2 and seen in the cystic fibrosis application study in Section 3.1. This new observation offers a possible complementary explanation of the “missing heritability” discussed extensively in Maher (2008).

In practice, SNPs out of HWE are typically not analyzed due to concerns for low genotyping quality (Wellcome Trust Case Control Consortium, 2007; Bycroft et al., 2018; Marees et al., 2018). However, the observation made here remains relevant as the heritability estimates in LMM-based models are biased when the causal SNPs are in HWD (which is unknown in practice) but not the tested SNPs. This is also supported by Figure 3B. When there was HWD at the causal SNP (e.g., δ_causal = 0.10 on the X-axis), there was a type I error issue even if there was no HWD at the tested SNP (i.e., δ_tested = 0). Conversely, Figure 3A shows that if there was no HWD at the causal SNP (i.e., δ_causal = 0 on the X-axis), then the test is accurate even if there was HWD at the tested SNP (δ_tested = 0.06).

Additionally, the HWE-based screening practice itself can be called into question because a truly associated SNP is often in HWD (Wittke-Thompson et al., 2005; Ryckman and Williams, 2008; Turner et al., 2011). The potential of leveraging the HWD expected at a causal SNP to increase the power of association testing has been explored by several groups (Song and Elston, 2006; Wang and Shete, 2008; Zhang and Sun, 2020).

We have not examined the implication of HWD combined with linkage disequilibrium (LD) (Weir, 2008) on the LMM, which is an important future research question. Additionally, recent work has shown that dominant genetic effect could complicate the LD measure and interpretation (Palmer et al., 2021), which in turn could affect our examination of the effect of HWD on the LMM.

Although the linear mixed-effect model is a popular and powerful method for GWAS, conceptually, the use of kinship coefficient matrix (i.e., Σ_Φ), derived from G, as part of the variance–covariance matrix (i.e., Σ_y) of the LMM can be problematic because the response variable Y is the phenotype of interest. An alternative approach is to reverse the roles of Y and G in the regression model. Indeed, O’Reilly et al. (2012) proposed MultiPhen, a method that treats the genotype G of an SNP as the response variable and phenotype values Y of multiple traits as predictors, and uses an ordinal logistic regression applicable to independent samples. More recently, Zhang (2021) (Chapter 2) proposed a generalized reverse (or retrospective) regression model that can be applied to dependent samples, which takes the form of

$g=\alpha 1+\beta \text{y}+\gamma \text{z}+ϵ,ϵ\sim N\left(\text{0},{\sigma }^2{\mathrm{\Sigma }}_g\right),{\sigma }^2{\mathrm{\Sigma }}_g={\sigma }^2{\mathrm{\Sigma }}_{\mathrm{\Phi }}+{\mathrm{\Sigma }}_{\delta },(6)$

where Σ_Φ is the kinship coefficient matrix as defined earlier and Σ_δ is a function of δ that explicitly models the amount of HWD; the use of a linear model for the discrete genotype data G is motivated by the work of Chen (1983).

Interestingly, if the variance and covariance matrices in Eqs 4, 6 of the LMM were the same, the resulting score test statistics are also the same. However, conceptually, the model Eq. 6 correctly uses the kinship coefficient matrix to model the response variable G, in contrast to the LMM model of Eq. 4. Specifically, at a tested SNP, as the reverse regression is conditional on Y, the variance–covariance matrix only concerns G_tested, that is, Σ_g. The modeling and estimation of Σ_g can account for potential HWD through Σ_δ, in addition to the genetic correlation captured by the kinship coefficient matrix of Σ_Φ, resulting in a more robust association test for related individuals. Indeed, when the method was applied to the same simulated sib-pair data in Section 3.3, it had correct type I error control [results shown in Figure 2.2 of Chapter 2 of Zhang (2021)]. However, how to model gene–environment interaction through the reverse regression framework remains an open question.

Data Availability Statement

The data analyzed in this study are subject to the following licenses/restrictions: The CF application data are available by application to the Cystic Fibrosis Canada National data registry for researchers who meet the criteria for access to confidential clinical data for the purpose of CF research. Requests to access these datasets should be directed to cfregistry@cysticfibrosis.ca.

Ethics Statement

The studies involving human participants were reviewed and approved by the Canadian Gene Modifier Study (CGMS), the Research Ethics Board of the Hospital for Sick Children (# 0020020214 from 2012–2019 and #1000065760 from 2019-present), and all participating subsites. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.

Author Contributions

LZ and LS proposed the method and wrote the manuscript. LZ performed the analysis. LS obtained the application data and funding.

Funding

This research was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC; RGPIN-04934 and RGPAS-522594).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Dr. Lisa J. Strug and acknowledge her laboratory for providing the cystic fibrosis genotype data. LZ was a trainee and funding recipient of the CANSSI Ontario STAGE (Strategic Training for Advanced Genetic Epidemiology) program at the University of Toronto.

References

Bycroft, C., Freeman, C., Petkova, D., Band, G., Elliott, L. T., Sharp, K., et al. (2018). The uk Biobank Resource with Deep Phenotyping and Genomic Data. Nature 562, 203–209. doi:10.1038/s41586-018-0579-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, C.-F. (1983). Score Tests for Regression Models. J. Am. Stat. Assoc. 78, 158–161. doi:10.1080/01621459.1983.10477945

CrossRef Full Text | Google Scholar

Chen, G.-B. (2014). Estimating Heritability of Complex Traits from Genome-wide Association Studies Using IBS-Based Hasemanâ€"Elston Regression. Front. Genet. 5, 107. doi:10.3389/fgene.2014.00107

PubMed Abstract | CrossRef Full Text | Google Scholar

Derkach, A., Lawless, J. F., and Sun, L. (2015). Score Tests for Association under Response-dependent Sampling Designs for Expensive Covariates. Biometrika 102, 988–994. doi:10.1093/biomet/asv038

CrossRef Full Text | Google Scholar

Dimitromanolakis, A., Paterson, A. D., and Sun, L. (2019). Fast and Accurate Shared Segment Detection and Relatedness Estimation in Un-phased Genetic Data via Truffle. Am. J. Hum. Genet. 105, 78–88. doi:10.1016/j.ajhg.2019.05.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Eu-Ahsunthornwattana, J., Miller, E. N., Fakiola, M., Jeronimo, S. M. B., Blackwell, J. M., Cordell, H. J., et al. (2014). Comparison of Methods to Account for Relatedness in Genome-wide Association Studies with Family-Based Data. Plos Genet. 10, e1004445. doi:10.1371/journal.pgen.1004445

PubMed Abstract | CrossRef Full Text | Google Scholar

Falconer, D. S. (1985). A Note on Fisher's 'average Effect' and 'average Excess'. Genet. Res. 46, 337–347. doi:10.1017/s0016672300022825

PubMed Abstract | CrossRef Full Text | Google Scholar

Hardy, G. H. (1908). Mendelian Proportions in a Mixed Population. Science 28, 49–50. doi:10.1126/science.28.706.49

PubMed Abstract | CrossRef Full Text | Google Scholar

Maher, B. (2008). Personal Genomes: The Case of the Missing Heritability. Nature 456, 18–21. doi:10.1038/456018a

PubMed Abstract | CrossRef Full Text | Google Scholar

Manolio, T. A., Collins, F. S., Cox, N. J., Goldstein, D. B., Hindorff, L. A., Hunter, D. J., et al. (2009). Finding the Missing Heritability of Complex Diseases. Nature 461, 747–753. doi:10.1038/nature08494

PubMed Abstract | CrossRef Full Text | Google Scholar

Marees, A. T., de Kluiver, H., Stringer, S., Vorspan, F., Curis, E., Marie-Claire, C., et al. (2018). A Tutorial on Conducting Genome-wide Association Studies: Quality Control and Statistical Analysis. Int. J. Methods Psychiatr. Res. 27, e1608. doi:10.1002/mpr.1608

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Reilly, P. F., Hoggart, C. J., Pomyen, Y., Calboli, F. C., Elliott, P., Jarvelin, M. R., et al. (2012). Multiphen: Joint Model of Multiple Phenotypes Can Increase Discovery in Gwas. PLoS One 7, e34861. doi:10.1371/journal.pone.0034861

PubMed Abstract | CrossRef Full Text | Google Scholar

Palmer, D. S., Zhou, W., Abbott, L., Baya, N., Churchhouse, C., Seed, C., et al. (2021). Analysis of Genetic Dominance in the uk Biobank. bioRxiv

Google Scholar

Powell, J. E., Visscher, P. M., and Goddard, M. E. (2010). Reconciling the Analysis of Ibd and Ibs in Complex Trait Studies. Nat. Rev. Genet. 11, 800–805. doi:10.1038/nrg2865

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryckman, K., and Williams, S. M. (2008). Calculation and Use of the hardy-weinberg Model in Association Studies. Curr. Protoc. Hum. Genet. Chapter 1, Unit–18. doi:10.1002/0471142905.hg0118s57

CrossRef Full Text | Google Scholar

Sasieni, P. D. (1997). From Genotypes to Genes: Doubling the Sample Size. Biometrics 53, 1253–1261. doi:10.2307/2533494

PubMed Abstract | CrossRef Full Text | Google Scholar

Schaid, D. J., and Jacobsen, S. J. (1999). Blased Tests of Association: Comparisons of Allele Frequencies when Departing from Hardy-Weinberg Proportions. Am. J. Epidemiol. 149, 706–711. doi:10.1093/oxfordjournals.aje.a009878

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, K., and Elston, R. C. (2006). A Powerful Method of Combining Measures of Association and Hardy-Weinberg Disequilibrium for fine-mapping in Case-Control Studies. Statist. Med. 25, 105–126. doi:10.1002/sim.2350

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, L., Dimitromanolakis, A., and Chen, W.-M. (2017). “Identifying Cryptic Relationships,” in Statistical Human Genetics (Springer), 45–60. doi:10.1007/978-1-4939-7274-6_4

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, L., Rommens, J. M., Corvol, H., Li, W., Li, X., Chiang, T. A., et al. (2012). Multiple Apical Plasma Membrane Constituents Are Associated with Susceptibility to Meconium Ileus in Individuals with Cystic Fibrosis. Nat. Genet. 44, 562–569. doi:10.1038/ng.2221

PubMed Abstract | CrossRef Full Text | Google Scholar

Turner, S., Armstrong, L. L., Bradford, Y., Carlson, C. S., Crawford, D. C., Crenshaw, A. T., et al. (2011). Quality Control Procedures for Genome-wide Association Studies. Curr. Protoc. Hum. Genet. Chapter 1, Unit1–19. doi:10.1002/0471142905.hg0119s68

CrossRef Full Text | Google Scholar

Vanscoy, L. L., Blackman, S. M., Collaco, J. M., Bowers, A., Lai, T., Naughton, K., et al. (2007). Heritability of Lung Disease Severity in Cystic Fibrosis. Am. J. Respir. Crit. Care Med. 175, 1036–1043. doi:10.1164/rccm.200608-1164oc

PubMed Abstract | CrossRef Full Text | Google Scholar

Visscher, P. M., Hill, W. G., and Wray, N. R. (2008). Heritability in the Genomics Era - Concepts and Misconceptions. Nat. Rev. Genet. 9, 255–266. doi:10.1038/nrg2322

PubMed Abstract | CrossRef Full Text | Google Scholar

Visscher, P. M., Medland, S. E., Ferreira, M. A. R., Morley, K. I., Zhu, G., Cornes, B. K., et al. (2006). Assumption-free Estimation of Heritability from Genome-wide Identity-By-Descent Sharing between Full Siblings. Plos Genet. 2, e41. doi:10.1371/journal.pgen.0020041

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., and Shete, S. (2008). A Test for Genetic Association that Incorporates Information about Deviation from hardy-weinberg Proportions in Cases. Am. J. Hum. Genet. 83, 53–63. doi:10.1016/j.ajhg.2008.06.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Weinberg, W. (1908). On the Demonstration of Heredity in Man. in (1963) Papers on Human Genetics.

Google Scholar

Weir, B. (1996). Genetic Data Analysis II: Methods for Discrete Population Genetic Data. Sunderland, Massachusetts: Sinauer Series (Sinauer).

Google Scholar

Weir, B. S. (2008). Linkage Disequilibrium and Association Mapping. Annu. Rev. Genom. Hum. Genet. 9, 129–142. doi:10.1146/annurev.genom.9.081307.164347

PubMed Abstract | CrossRef Full Text | Google Scholar

Wellcome Trust Case Control Consortium (2007). Genome-wide Association Study of 14,000 Cases of Seven Common Diseases and 3,000 Shared Controls. Nature 447, 661–678. doi:10.1038/nature05911

PubMed Abstract | CrossRef Full Text | Google Scholar

Wittke-Thompson, J. K., Pluzhnikov, A., and Cox, N. J. (2005). Rational Inferences about Departures from hardy-weinberg Equilibrium. Am. J. Hum. Genet. 76, 967–986. doi:10.1086/430507

PubMed Abstract | CrossRef Full Text | Google Scholar

Wright, F. A., Strug, L. J., Doshi, V. K., Commander, C. W., Blackman, S. M., Sun, L., et al. (2011). Genome-wide Association and Linkage Identify Modifier Loci of Lung Disease Severity in Cystic Fibrosis at 11p13 and 20q13.2. Nat. Genet. 43, 539–546. doi:10.1038/ng.838

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J., Lee, S. H., Goddard, M. E., and Visscher, P. M. (2011). Gcta: a Tool for Genome-wide Complex Trait Analysis. Am. J. Hum. Genet. 88, 76–82. doi:10.1016/j.ajhg.2010.11.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L. (2021). A General Study of Genetic Association Tests and the Test of Hardy–Weinberg Equilibrium. Ph.D. thesis. University of Toronto.

Google Scholar

Zhang, L., and Sun, L. (2021). A Generalized Robust Allele-Based Genetic Association Test. Oxford, United Kingdom: Biometrics.

Google Scholar

Zhang, L., and Sun, L. (2020). Leveraging hardy-weinberg Disequilibrium for Association Testing in Case-Control Studies. bioRxiv.

Google Scholar