HEDE: Heritability estimation in high dimensions
by Ensembling Debiased Estimators

Accurate heritability estimation presents a significant challenge in statistical genetics. Distinguishing the contributions of genetic versus environmental factors is crucial for understanding how our genetic makeup influences the development of complex diseases. Importantly, genetic research over the past decade has highlighted a significant discrepancy between heritability estimates from twin studies and those from Genome-Wide Association Studies (GWAS) (Manolio et al.,, 2009; Yang et al.,, 2015) This discrepancy, known as the “missing heritability problem”, underscores the importance of developing reliable heritability estimation methods to understand the extent of missing heritability. In this paper, we propose an ensemble method for heritability estimation, using high-dimensional Single Nucleotide Polymorphisms (SNPs) commonly collected in GWAS, that exhibits superior performance across a wide variety of settings and significantly outperforms existing methods in practical scenarios of interest. We focus on the narrow-sense heritability yang2011gcta, which refers to the proportion of phenotypic variance explained by additive genotypic effects. For conciseness, we henceforth refer to this simply as the heritability.

GWAS data typically contains more SNPs than individuals. This presents a challenge when estimating heritability. Addressing this challenge necessitates modeling assumptions. Previous work in this area can be broadly categorized into two classes based on these assumptions: random effects models and fixed effects models.

Random effects models treat the underlying regression coefficients as random variables and have been widely used for heritability estimation. Arguably, the most popular approach in this vein is GCTA, introduced by yang2011gcta. GCTA employs the classical restricted maximum likelihood method (REML) (also called genomic REML/GREML in their work) within a linear mixed effects model. It assumes that the regression coefficients are i.i.d. draws from a Gaussian distribution with zero mean and a constant variance. Since the genetic effect sizes/regression coefficients might depend on the design matrix (the genotypes), subsequent work provides various improvements that relax this assumption, including GREML-MS (lee2013estimation), GREML-LDMS (yang2015genetic), BOLT-LMM (loh2015efficient), among others. Specifically, GREML-MS allows the signal variance to depend on minor allele frequencies (MAFs), while GREML-LDMS allows it to depende on linkage disequilibrium (LD) levels/quantiles. Prior statistical literature has studied the theoretical properties of several random effects methods (jiang2016high; bonnet2015heritability; dicker2016maximum; hu2022misspecification).

Despite this rich literature, a major limitation of random effects model based heritability estimation method is that its consistency often relies on the correct parametric assumptions of regression coefficients. Specifically, these methods assume that regression coefficients are either i.i.d, or depend on the design matrix in specific ways. When these assumptions are violated, random effects methods can produce significantly biased heritablity estimates. The issue is magnified when SNPs are in linkage disequilibrium (LD) or when there is additional heterogeneity across the LD blocks not captured by the random effects assumptions. We elaborate this issue further in Figure 1, where we demonstrate that GREML, GREML-MS and GREML-LDMS all incur significant biases when their respective random effects assumptions are violated.

These challenges associated with random effects methods have spurred a distinct strand of methods by positing regression coefficients to be fixed while treating the design matrix to be random. Referred to as fixed effects methods, these approaches often exhibit greater robustness to diverse realizations of the underlying signal without the need to make parametric assumptions on regression coefficients (Figure 1). Consequently, our focus lies on estimating the heritability within this fixed effects model framework.

Several fixed effects methods have emerged for heritability estimation. However, a substantial body of the proposed work requires the underlying signal to be sparse (fan2012variance; sun2012scaled; tony2020semisupervised; guo2019optimal), and they are subject to loss of efficiency when these sparsity assumptions are violated (dicker2014variance; janson2017eigenprism; bayati2013estimating), see Sections 5.2, 3.3, 2, respectively. The methods proposed by several authors (chen2022statistical; dicker2014variance; janson2017eigenprism) operate under less stringent assumptions regarding the signal, but they depend on ridge-regression-analogous concepts. Therefore they are expected to have good efficiency when the signal is dense, and are subject to loss of efficiency when signals are sparse. We illustrate these issues further in Section 5.2. In practice, the underlying genetic architecture for a given trait is unknown, i.e., the signal could be sparse or dense. Hence, it is desirable to develop a robust method across a broad spectrum of signal regimes from sparse to dense. In addition, fixed effects methods often rely on assuming the covariate distribution to be Gaussian (bayati2013estimating; dicker2014variance; janson2017eigenprism; verzelen2018adaptive). This assumption fails to capture genetic data where genotypes are discrete, taking values in ${0,1,2}$.

To address these issues, we propose a robust ensemble heritability estimation method HEDE (Heritability Estimation via Debiased Ensemble), which ensembles the debiased Lasso and ridge estimators, with degrees-of-freedom corrections (bellec2019biasing; javanmard2014hypothesis). Debiased estimators, originally proposed to correct the bias for Lasso/ridge estimators (zhang2014confidence; van2014asymptotically; javanmard2014confidence), have been employed for tasks such as confidence interval construction and hypothesis testing, and have nice properties such as optimal testing properties javanmard2014confidence. Unlike these methods, HEDE utilizes debiasing for heritability estimation. To develop HEDE, we first demonstrate that a reliable estimator for heritability can be formulated using any linear combination of these debiased estimators, in which we devise a consistent bias correction strategy. Subsequently, we refine the linear combination based on this bias correction strategy. This methodology enables the consistent estimation of heritability using any linear combination of the debiased Lasso and the debiased Ridge. HEDE then uses a data-driven approach to selecting an optimal ensemble from the array of possible linear combinations. Finally, we introduce an adaptive method for determining the underlying regularization parameters. We develop both adaptive strategies—ensemble selection and tuning parameter selection for the Lasso and the Ridge—by minimizing a certain mean square error. This results in enhanced statistical performance for heritability estimation.

Underpinning our method is a rigorous statistical theory that ensures consistency under minimal assumptions. Specifically, we operate under sub-Gaussian assumptions, which are consistent with a broad class of designs observed in GWAS that motivates this paper. We show that the consistency of heritability estimation using HEDE guarantees remain intact despite employing adaptive strategies for selecting the tuning parameters and the debiased Lasso and ridge ensemble, all while accommodating sub-Gaussian designs.

We perform extensive simulation studies to evaluate the finite sample performance of HEDE compared with the existing methods in a range of settings mimicking GWAS data. We demonstrate that HEDE overcomes the limitations in random effects methods, remaining unbiased across a wide variety of signal structures, specifically signals with different amplitudes across various combinations of MAF and LD levels, as well as LD blocks. our simulation also illustrates that HEDE consistently achieves the lowest mean square error when compared to popular fixed effects methods across a wide spectrum of sparse and dense signals, as well as low and high heritability values. We compare HEDE’s performance with both random and fixed effects methods on the problems of estimating height and BMI heritability using the UK Biobank data.

The rest of this paper is organized as follows. Section 2 introduces the problem setup. Section 3 presents our HEDE while Section 4 discusses its theoretical properties. Section 5 presents extensive simulation studies while Section 6 complements these numerical investigations with a real-data application on the UK Biobank data (sudlow2015uk). Finally, Section 7 concludes with discussions.

In this section, we formally introduce our problem setup. We consider a high-dimensional regime where the covariate dimension grows with the sample size. To elaborate, we consider a sequence of problem instances $\{\bm{y}(n),\bm{X}(n),\bm{\beta}(n),\bm{\epsilon}(n)\}_{n\geq 1}$ that satisfies a linear model

where $\bm{y}(n)=(y_{1},\ldots,y_{n})^{\top}\in\mathbb{R}^{n}$ represents observations of a real-valued phenotype of interest, such as height, BMI, blood pressure, for $n$ unrelated individuals from a population. The unknown regression coefficient vector, $\bm{\beta}(n)=(\beta_{1},\ldots,\beta_{p(n)})^{\top}\in\mathbb{R}^{p(n)}$, captures the effects of $p(n)$ SNPs and is assumed to be fixed parameters. The environmental noise, $\bm{\epsilon}(n)=(\epsilon_{1},\ldots,\epsilon_{n})\in\mathbb{R}^{n}$, is independent of $\bm{X}(n)$, with each entry also independent. Lastly, $\bm{X}(n)\in\mathbb{R}^{n\times p(n)}$ signifies the matrix of observed genotypes with i.i.d. rows $\bm{x}_{i\bullet}$, which are normalized as follows

where $G_{ij}=0,1,2$ is the genotype containing the minor allele count (the number of copies of the less frequent allele) at SNP $j$ for individual $i$ and $\bar{G}_{j}=\sum_{i=1}^{n}G_{ij}/n$ is the minor allele frequency of SNP $j$ in the sample.

We assume the covariates have sub-Gaussian distributions. Since $\bm{X}$ is random and $\bm{\beta}$ is deterministic, we have placed ourselves in the fixed effects framework. Moving forward, we suppress the dependence on $n$ when the meaning is clear from the context. Our parameter of interest is the heritability, which is the proportion of variance in $y_{i}$ explained by features in $\bm{x}_{i\bullet}$ defined as

We emphasize the following remark regarding the aforementioned definition. We assume that $n$ and $p$ diverge with $p/n\rightarrow\delta>0$. This setting–also known as the proportional asymptotics regime–has gained significant recent attention in high-dimensional statistics. A major advantage of this regime is that theoretical results derived under such asymptotics demonstrate attractive finite sample performance (sur2019likelihood; sur2019modern; candes2020phase; zhao2022asymptotic; liang2022precise; jiang2022new). Furthermore, prior research in heritability estimation (jiang2016high; dicker2016maximum; janson2017eigenprism) utilized this framewrok to characterize theoretical properties of GREML and its variants. In practical situations, we observe a single value for $n$ and $p$, so we may substitute the corresponding ratio $p/n$ into our theory in place of $\delta$ to execute data analyses. Formal assumptions on the covariate and noise distributions, and the signal are deferred to Section 4.

We begin by describing the overarching theme underlying our method. For simplicity, we discuss the case where $\bm{X}$ has independent columns, deferring the extension to correlated columns to Section 3.4. The denominator of $h^{2}$ in (3) can be consistently estimated using the sample variance of $\bm{y}$. The numerator of $h^{2}$ in (3) simplifies to $\|\bm{\beta}\|_{2}^{2}$ in the setting of independent columns. Consequently, we seek to develop an accurate estimator for this norm.

To this end, we turn to the debiased Lasso estimator and the debiased ridge estimator. Denote the Lasso and the ridge by $\hat{\bm{\beta}}_{\textrm{L}}$ and $\hat{\bm{\beta}}_{\textrm{R}}$ respectively. These penalized estimators are known to shrink estimated coefficients to $0$, therefore inducing bias. The existing literature provides debiased versions of these estimators, denoted as $\hat{\bm{\beta}}_{\textrm{L}}^{d}$ and $\hat{\bm{\beta}}_{\textrm{R}}^{d}$. These debiased estimators roughly satisfy the following property: $\hat{\bm{\beta}}_{\textrm{L}}^{d}\approx\bm{\beta}+\sigma_{\textrm{L}}Z_{1}$, $\hat{\bm{\beta}}_{\textrm{R}}^{d}\approx\bm{\beta}+\sigma_{\textrm{R}}Z_{2}$ meaning that each entry of $\hat{\bm{\beta}}_{\textrm{L}}^{d},\hat{\bm{\beta}}_{\textrm{R}}^{d}$ is approximately centered around the corresponding true signal coefficient with Gaussian fluctuations and variances given by $\sigma_{L},\sigma_{R}$ respectively. By taking norms on both sides, we obtain that $\|\hat{\bm{\beta}}_{k}^{d}\|_{2}^{2}-p\cdot\sigma_{k}^{2}\approx\|\bm{\beta}\|_{2}^{2}$, where $k=\textrm{L}\,\,\text{or}\,\,\textrm{R}$ depending on whether we consider the Lasso or the ridge. Thus, precise estimates for the variances $\sigma_{k}^{2}$’s produce accurate estimates for $h^{2}$ starting from each debiased estimator.

This intuition extends to any linear combination of the debiased Lasso and ridge in the form of $\tilde{\bm{\beta}}_{C}^{d}:=\alpha_{\textrm{L}}\hat{\bm{\beta}}_{\textrm{L}}^{d}+(1-\alpha_{\textrm{L}})\hat{\bm{\beta}}_{\textrm{R}}^{d}$. A similar relationship holds for the norm: $\|\tilde{\bm{\beta}}_{C}^{d}\|_{2}^{2}-p\cdot\tilde{\tau}_{C}^{2}\approx\|\bm{\beta}\|_{2}^{2}$. We establish that for any $\alpha_{L}\in[0,1]$, we can derive a consistent estimator for the ensemble variance $\tilde{\tau}_{C}^{2}$. Among all possible ensemble choices with different ensembling parameter $\alpha_{\textrm{L}}$ and different tuning parameters $\lambda_{\textrm{L}},\lambda_{\textrm{R}}$ for calculating the Lasso and the ridge, we adaptively choose the combination that minimizes the mean square error of the ensemble estimator. This method, dubbed HEDE, harnesses the strengths of $L_{1}$ and $L_{2}$ regression and demonstrates strong practical performance, since we design our adaptive tuning process to particularly enhance statistical accuracy. We provide further elaboration on HEDE in the subsequent subsections: Section 3.1, discusses the specific forms of the debiased estimators utilized in our framework; Section 3.2 describes our ensembling approach, and Section 3.3 discusses the final HEDE method with our adaptive tuning strategy.

We turn to discuss the key theoretical properties of HEDE. We first focus on the case of independent covariates. Recall that throughout, we assume an asymptotic framework where $p,n\rightarrow\infty$ with $\lim_{n\rightarrow\infty}p/n\rightarrow\delta\in(0,\infty)$. We work with the following assumptions on the covariate, signal, and noise.

It is worth noting that the uniform sub-Gaussianity in Assumption 1(i) implies that the $\bar{G}_{j}$’s from (2) are bounded. This is satisfied in GWAS, as $\bar{G}_{j}$ is minor allele frequency (MAF) (psychiatric2009genomewide). The remaining Assumptions 1(ii) and (iii) are all weak and natural to impose. With Assumption 1 in mind, our primary goal is to establish that HEDE consistently estimates the heritability. The following theorem establishes this result under mild additional conditions on the minimum singular values of submatrices of $\bm{X}$ (Assumption 2, Supplementary Materials A.1). The proof is provided in Supplementary Materials F.3.

Theorem 4.1 states that HEDE is consistent when $\lambda_{L},\lambda_{R}$ and $\alpha_{\textrm{L}}$ are arbitrarily chosen in $[\lambda_{\min},\lambda_{\max}]$ and $[0,1]$, respectively. It follows that HEDE is consistent using our proposed choices $\hat{\alpha}_{\textrm{L}},\hat{\lambda}_{\textrm{L}},\hat{\lambda}_{\textrm{R}}$. Here $[\lambda_{\min},\lambda_{\max}]$ denotes the range of $\lambda_{L},\lambda_{R}$ that we use. See details in Supplementary Materials A.2.

Proving the consistency for given fixed choices of these hyperparameters is relatively straightforward, However, establishing the uniform gurantee is non-trivial. It requires analyzing the debiased Lasso and ridge estimators jointly, and furthermore, pinning down this joint behavior uniformly over all hyperparameter choices. We achieve both these goals under relatively mild assumptions on the covariate and noise distributions. Once we establish the uniform control, Part (ii) follows trivially since the uniform result provides consistency for the particular choices $\hat{\alpha}_{L},\hat{\lambda}_{L},\hat{\lambda}_{R}$ produced by our algorithm. Below, we present the key theorem necessary for proving Theorem 4.1.

We now discuss how Proposition 4.2 leads to our previous Theorem 4.1. Proposition 4.2 equivalently proposes a consistent estimator for $\|\bm{\beta}\|_{2}^{2}$, which is directly related to $h^{2}$ by a factor of the population variance of $y_{i}$. Since the population variance of $y_{i}$ can be consistently estimated by the sample variance of the observed outcomes, we can safely replace the former by the latter in the denominator. Theorem 4.1(i) then follows on using the final fact that clipping does not affect consistency.

Proposition is 4.2 proved in Supplementary materials F.2. Here we briefly outline our main idea. Recall from (8) that $\tilde{\bm{\beta}}_{C}^{d}$ is a linear combination of $\hat{\bm{\beta}}_{\textrm{L}}^{d}$ and $\hat{\bm{\beta}}_{\textrm{R}}^{d}$. To develop uniform guarantees for $\tilde{\bm{\beta}}_{C}^{d}$, it therefore suffices to establish the corresponding uniform guarantees for the debiased Lasso and ridge jointly, which we establish in the following theorem.

Theorem 4.3 is proved in Supplementary materials D. It states that the joint “Gaussianity” described in (6) occurs at the population level in terms of the value of Lipschitz functions of the estimators, in addition to a per-coordinate basis. While results of this nature have appeared in the recent literature on exact high-dimensional asymptotics, majority of these works focus on a single regularized estimator (or its debiased version) with a fixed tuning parameter value. Theorem 4.3 provides the first joint characterization of multiple debiased estimators that is uniform over a range of tuning parameter values.

Among related works, celentano2021cad established joint characterization of two estimators and their debiased versions for given fixed tuning parameter choices. Meanwhile, miolane2021distribution characterized certain properties of the Lasso and the debiased Lasso uniformly across tuning parameter values. However, developing a heritability estimator that performs well in practice across sparse to dense signals requires the stronger result of the form in Theorem 4.3. Dealing with the debiased lasso and ridge jointly, while requiring uniform characterization over tuning parameters as well as allowing sub-Gaussian covariates, poses unique challenges. We overcome these challenges in our work, building upon miolane2021distribution, celentano2020lasso, celentano2021cad and universality techniques proposed in han2022universality, which in turn rely on the Convex Gaussian Minmax Theorem (GCMT) (thrampoulidis2015regularized). We detail the connections and differences in the Supplementary Materials.

Note that Theorem 4.3 involves the variance-covariance parameters $\tau_{L}^{2},\tau_{LR},\tau_{R}^{2}$ (formal definitions deferred to later (23) (24) in the interest of space. We establish in Supplementary Materials, Theorem F.1 that $\hat{\tau}_{\textrm{L}}^{2},\hat{\tau}_{\textrm{R}}^{2},\hat{\tau}_{\textrm{L}\textrm{R}}$, defined in (7), consistently estimate these parameters. Furthermore, we show that this consistency holds uniformly over the required range of tuning parameter values. This result, coupled with Theorem 4.3, yields our key Proposition 4.2.

Finally we turn to the case where the observed features are correlated. Recall from Section 3.4 that here we consider the design matrix to take the form $\bm{X}=\bm{Z}\bm{\Sigma}^{1/2}$, where $\bm{Z}$ now has independent entries and satisfies the conditions in Assumption 1(i). Thus, our previous theorems apply to $\bm{Z}$ (which is unobserved), but not to the observed $\bm{X}$. Applying whitening with estimated $\hat{\bm{\Sigma}}$ results in final features of the form $\bm{X}\hat{\bm{\Sigma}}^{-1/2}=\bm{Z}\bm{\Sigma}^{1/2}\hat{\bm{\Sigma}}^{-1/2}$. When $\hat{\bm{\Sigma}}^{1/2}$ approximates $\bm{\Sigma}$ accurately, $\bm{\Sigma}^{1/2}\hat{\bm{\Sigma}}^{-1/2}$ behaves as the identity matrix $\bm{I}$. In this case, $\bm{X}\hat{\bm{\Sigma}}^{-1/2}$ is close to $\bm{Z}$ which satisfies our previous Assumption 1. HEDE therefore continues to enjoy uniform consistency guarantees, as formalized below.

Although we present this theorem for a general sequence of covariance matrices $\bm{\Sigma}$ that allow accurate estimators in operator norm, we use it only for block diagonal covariances relevant for our application, where the block sizes are relatively small compared to the sample size. In Supplementary Materials, Proposition A.1, we establish that for such covariance matrices, the block-wise sample covariance matrix provides such an acccurate estimator. Theorem 4.4 is proved in Supplementary Materials G. The key idea lies in establishing that the estimation error incurred while replacing $\bm{\Sigma}$ by $\hat{\bm{\Sigma}}$ does not affect the uniform consistency properties of HEDE.

We evaluate the performance of HEDE through simulation studies. In Section 5.1, we compare HEDE to some representative random effects methods for estimating the heritability under different signal structures. In Section 5.2, we compare HEDE against several fixed effect based methods for estimating the heritability.

We next describe our real data experiments. We apply HEDE on the unrelated white British individuals from the UK biobank GWAS data which contain common variants (sudlow2015uk) to estimate the heritability for two commonly studied phenotypes: height and BMI. For the design matrix $\bm{X}$, we follow the preprocessing steps outlined in Section 5.1, but include all $22$ chromosomes, unlike Section 5.1 which considers only chromosome $22$. For each phenotype, we exclude individuals with missing values. We center the outcomes and denote the resultant vector by $\tilde{\bm{y}}$. We account for non-genetic factors by including covariates for age, age², sex as well as $40$ ancestral principal components. We derive the final response vector $\bm{y}$ by regressing out the influence of these non-genetic covariates from $\tilde{\bm{y}}$ using Ordinary Least Squares (OLS). The large sample size of approximately $330,000$ individuals compared to only $43$ non-genetic covariate, ensures that OLS can be confidently applied.

We apply HEDE to preprocessed $\bm{X}$ and $\bm{y}$. For comparison, we also apply GREML-SC and GREML-LDMS (see Section 5.1 for details) from the random effects methods family, as well as MM and AMP from the fixed effects methods family. We did not include Eigenprism and EstEqn since their runtimes are both $O(n^{3})$ (assuming $p/n\rightarrow\delta$), so they could not be completed within a reasonable time frame. We calculate standard errors using the standard deviation of point estimates from $15$ disjoint random subsets, each containing $20,000$ individuals.²²2In our high-dimensional setting, it is known that traditional resampling techniques for estimating standard errors such as the bootstrap, subsampling, etc. fail (el2018can; clarte2024analysis; bellec2024asymptotics). Corrections have been proposed to address issues with using these for estimating standard errors while inferring individual regression vector coefficients and closely related problems. However, such strategies for heritability estimation are yet to be developed. To circumvent this issue, we use disjoint subsets of observed data to mimic independent training data draws from the population. Additionally, due to memory constraints, many methods (including HEDE) could not process all $22$ chromosomes simultaneously. Therefore, for all methods, we estimate the total heritability by summing the heritabilities calculated for each chromosome separately. This approach is expected to introduce minimal error, given the widely accepted assumption of independence across chromosomes.

We summarize the results in Table 2. First, we observe that MM and AMP estimates differ more from HEDE compared to other methods. This is consistent with Figure 2, where we observe that AMP often has larger MSE than other fixed effects methods. Thus, we expect AMP estimates to differ significantly from HEDE, and MM. For MM, the situation is more nuanced. Prior work suggests the true height heritability should be on the higher end of the range shown in Figure 2. In this region, MM has higher MSE than HEDE across all sparsity levels. Therefore, MM estimates of height heritability may be less accurate than HEDE. For BMI, prior work provides heritability estimates within the range $0.285$ to $0.436$ using various methods. Referring to Figure 2, we see that MM’s relative MSE compared to HEDE depends on the sparsity level in this BMI heritability range. Without knowing the exact sparsity, it is difficult to draw conclusions. However, we still observe the general trend of MM producing larger estimates than other methods for BMI, similar to height.

Comparing HEDE to the two random effect methods, we observe that HEDE point estimates fall between GREML-SC and GREML-LDMS ones. For both height and BMI, HEDE estimates lie between GREML-SC and GREML-LDMS,, with the latter always undershooting. This is consistent with the trend observed in the last column of Figure 1. This difference may indicate a heterogeneous underlying genetic signal, especially in a way that differentiates the SC and LDMS approaches within the GREML framework.

In this paper, we introduce HEDE, a new method for heritability estimation that harnesses the strengths of $\ell_{1}$ and $\ell_{2}$ regression in high-dimensional settings through a sophisticated ensemble approach. We develop data-driven techniques for selecting an optimal ensemble and fine-tuning key hyperparameters, aiming to enhance statistical performance in heritability estimation. As a result, we present a competitive estimator that surpasses existing fixed effects heritability estimators across diverse signal structures and heritability values. Notably, our approach circumvents bias issues inherent in random effects methods. In summary, our contribution offers a dependable heritability estimator tailored for high-dimensional genetic problems, effectively accommodating underlying signal variations. Importantly, our method maintains consistency guarantees even with adaptive tuning and under minimal assumptions on the covariate distribution. We validate the efficacy of our approach through comprehensive simulations and real-data analysis.

HEDE’s computational complexity is primarily driven by the Lasso/ridge solutions, making it more efficient than many existing heritability estimators (e.g. Eigenprism (janson2017eigenprism), where an SVD is required, and EstEqn (chen2022statistical), where chained matrix inversions are involved). HEDE is potentially scalable to much larger dimensions via snpnet (rad2020scalable). Furthermore, with proper Lasso/ridge solving techniques, HEDE only needs knowledge of summary statistics $\bm{X}^{\top}\bm{X}$ and $\bm{X}^{\top}\bm{y}$. This makes it an attractive option in situations where individual data sharing is difficult due to privacy concerns. However, it is important to note that these discussions assume the availiability of an accurate estimate of the LD blocks and strength. This is a common limitation that other fixed effects methods also face. Any potential errors in covariance estimation could naturally lead to additional errors in heritability estimation. Therefore, improving and scaling covariance estimation in high dimensions could independently benefit heritability estimation. We manage to circumvent this issue by leveraging the hypothesis of independence across chromosomes, which significantly aids us computationally.

Our theoretical framework relied on two pivotal assumptions: the independence of observed samples and a sub-Gaussian tailed distribution for covariates. Natural settings occur where these assumptions are violated. Within the genetics context, familial relationships and repeated measures in longitudinal studies among observed samples can introduce dependence. For the broader problem of signal-to-noise ratio estimation in various scientific fields, both assumptions might be violated. For instance, applications such as wireless communications (eldar2003asymptotic) and magnetic resonance imaging (benjamini2013shuffle) may introduce challenges such as temporal dependence or heavier-tailed covariate distributions. Exploring analogues of HEDE under such conditions presents an intriguing avenue for future study. Recent advancements in high-dimensional regression, which extend debiasing methodology and proportional asymptotics theory for regularized estimators to these complex scenarios (li2023spectrum; lahiry2023universality), offer valuable insights for further exploration. We defer these investigations to future research.

We focus on continuous traits in this paper. Discrete disease traits occur commonly in genetics. Such situations are typically modeled using logistic, probit, or binomial regression models, depending on the number of possible categorical values of the response. Adapting our current methodology to these scenarios is a crucial avenue for future research. Debiasing methodologies for generalized linear models exist in the literature. A compelling question would be to understand how various debiasing techniques can be ensembled to optimize their benefits for heritability estimation, similar to our approach in this work for the linear model.