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 = α ∗ 1 + β ∗ g + γ ∗ z + ϵ ∗ , ϵ ∗ ∼ N 0 , σ ∗ 2 I , (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  indep = n ⋅ g − g ̄ 1 T y − y ̄ 1 − g − g ̄ 1 T z − z ̄ 1 y − y ̄ 1 T z − z ̄ 1 z − z ̄ 1 T z − z ̄ 1 2 g − g ̄ 1 T g − g ̄ 1 − g − g ̄ 1 T z − z ̄ 1 2 z − z ̄ 1 T z − z ̄ 1 y − y ̄ 1 T y − y ̄ 1 − y − y ̄ 1 T z − z ̄ 1 2 z − z ̄ 1 T z − z ̄ 1 . (3) To observe T _indep’s connection with Hardy–Weinberg disequilibrium, it is instructive to employ some algebraic tricks and show that 1 n g − g ̄ 1 T g − g ̄ 1 = var ̂ G = 2 p ̂ 1 − p ̂ + p ̂ A A − p ̂ 2 = 2 p ̂ 1 − p ̂ + δ ̂ .

Because δ ̂ = p ̂ A A − p ̂ 2 measures the amount of HWD present in the data (Weir, 1996), T _indep inherently adjusts for departure from HWE through var ̂ ( G ) = 2 ( p ̂ ( 1 − p ̂ ) + δ ̂ ) . 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.

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. y = α ∗ 1 + β ∗ g + γ ∗ z + ϵ ∗ ,  where  ϵ ∗ ∼ N 0 , σ y 2 Σ y and Σ y = h 2 Σ Φ + 1 − h 2 I . (4)

Compared with the linear model used for independent samples, var(ϵ*) = σ*² I in Eq. 2 is replaced by σ y 2 Σ 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 σ y 2 as the phenotypic variation due to (additive) genetic variation, and ( 1 − h 2 ) σ y 2 = σ 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 = ∑ β k 2 var G k var Y = ∑ β k 2 2 p k 1 − p k ∑ β k 2 2 p k 1 − p k + σ 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 h ̂ 2 = corr ̂ Y 1 , Y 2 / 2 ϕ , 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, E h ̂ 2 h 2 = ∑ β k 2 2 p k 1 − p k + δ k ∑ β k 2 2 p k 1 − p k , and the bias of the h ² estimate is E h ̂ 2 − h 2 = h 2 ⋅ ∑ β k 2 δ k ∑ β k 2 p k 1 − p k . (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.

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 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 h ̂ 2 = 0.82 . This suggests that 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.

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 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 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.

In Figure 1, it is notable that 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, 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.

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 0.05 ⋅ 0.95 / 10000 = 0.002 .) In Figure 3A, the trend of type I error inflation is clear as δ _causal increases.

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.