Estimating Broad Sense Heritability via Kernel Ridge Regression

Genetic heritability, which measures the proportion of phenotypic variance explained by genetic factors, is a crucial quantity that can enhance our understanding of the genetic mechanisms underlying complex phenotypes [9]. As noted by Zhu and Zhou [21], two types of heritability can be estimated. The first is broad-sense heritability ($H^{2}$), which evaluates the proportion of phenotypic variance explained by all genetic factors, including additive effects, dominance effects, and epistatic interactions. The second type, narrow-sense heritability ($h^{2}$), represents the proportion of phenotypic variance explained solely by additive genetic effects.

Since the success of Genome-Wide Association Studies (GWAS), hundreds of single nucleotide polymorphisms (SNPs) associated with complex traits have been discovered [4, 2]. However, SNPs identified through GWAS explain only a small fraction of the genetic variation for any given trait, a well-known issue referred to as ”missing heritability” [10, 11]. One possible explanation is that each individual SNP contributes only a tiny effect to the trait; thus, many SNPs remain undetected as they fail to meet the stringent significance threshold in GWAS [20]. As a result, existing literature on narrow-sense heritability estimation often relies on modeling the joint effects of multiple SNPs, with linear mixed-effects models (LMMs) playing a key role in constructing narrow-sense heritability estimators. For example, [20] applied Restricted Maximum Likelihood (REML) multi-variance component analysis on large sample sizes to estimate heritability and uncover genetic architecture. This method demonstrated significant reductions in computation time and memory use compared to previous approaches. The BOLT-REML algorithm revealed that SNP heritability increases with higher-frequency SNPs, while rarer SNPs have a larger effect per allele but explain less variance per SNP. The study further characterized schizophrenia as having a highly polygenic architecture.

In line with this trend, numerous advancements have been made to improve LMM-based narrow-sense heritability estimators. For instance, [9] addressed three main challenges—scale, ascertainment, and SNP quality control—and demonstrated that genetic variation in disease liability, linked to common SNPs in linkage disequilibrium (LD), could be estimated from genome-wide association studies (GWAS) data. First, variance explained on the observed risk scale was transformed to an underlying liability scale. Second, the method accounted for ascertainment, which led to a much higher proportion of cases in the analyzed sample than in the general population. Third, the new approach minimized inflated heritability estimates by preventing artificial differences between case and control genotypes through rigorous GWAS quality control. This method was tested using simulated data and real GWAS data from three Wellcome Trust Case Control Consortium (WTCCC) datasets: Crohn’s disease, bipolar disorder, and type I diabetes. Additionally, Speed et al. [15] proposed a modified kinship matrix that weights SNPs based on local LD, significantly reducing bias in heritability estimates and enhancing estimation precision.

Several additional approaches have been developed to improve narrow-sense heritability estimation. For example, with the reduced cost of whole-genome sequencing and advances in data storage, biobank data has become accessible. [3] explored an estimator derived under a model that minimizes assumptions about genetic architecture, acknowledging that models with numerous assumptions risk bias due to the unknown nature of genetic architecture. They proposed a generalized random effects (GRE) model that makes minimal assumptions about genetic architecture. They found that the GRE estimator is nearly unbiased across all architectures compared to existing methods, with a maximum bias of around 2%. As anticipated, the GRE model demonstrated robust performance across architectures, with a runtime of O($N{p_{k}}^{2}$) for chromosome $k$ with $p_{k}$ SNPs, taking approximately 3 hours and 60 GB of memory for 50,000 SNPs. Another direction involves estimating narrow-sense heritability using only summary data, such as the LD matrix. Although individual-level data can provide more precise estimates, its accessibility is often limited by data-sharing restrictions. Recently, Li et al. [8] proposed an iterative procedure called HEELS, which uses summary statistics to solve REML-score equations by adapting Henderson’s algorithm for variance component estimation in LMMs. This approach improves the statistical efficiency of heritability estimators based on summary statistics, achieving performance comparable to that of REML.

Although accurately estimating narrow-sense heritability has received considerable attention in the literature, the true relationship between genetic components and phenotypic traits is unknown, and focusing solely on additive (linear) effects may result in a loss of information. Therefore, it is also important to estimate broad-sense heritability accurately. Mathematically, suppose the relationship between SNPs ($\boldsymbol{Z}$) and a phenotypic trait ($Y$) is given by $Y=g(\boldsymbol{Z})+\epsilon$, where $\epsilon$ is a random error. The broad-sense heritability can then be expressed as

If the function $g$ were known and i.i.d. data pairs $(\boldsymbol{Z}_{1},Y_{1}),(\boldsymbol{Z}_{2},Y_{2}),\ldots,(\boldsymbol{Z}_{n},Y_{n})$ were observed, a natural estimator for $\textrm{Var}[g(\boldsymbol{Z})]$ would be the sample variance $(n-1)^{-1}\sum_{i=1}^{n}(g(\boldsymbol{Z}_{i})-\bar{g})^{2}$, with $\bar{g}=n^{-1}\sum_{i=1}^{n}g(\boldsymbol{Z}_{i})$. Similarly, a natural estimator for $\textrm{Var}[\epsilon]$ is its sample counterpart, $n^{-1}\sum_{i=1}^{n}(Y_{i}-g(\boldsymbol{Z}_{i}))^{2}$. An estimator of $H^{2}$ can be obtained by replacing the variances with their respective estimators. However, since $g$ is unknown, it must also be estimated. In this paper, this estimation will be achieved through kernel ridge regression.

The remainder of the paper is organized as follows: Section 2 provides the details on our proposed estimator for the broad-sense heritability as well as some theoretical properties of the estimator, followed by some simulation studies whose results are presented in Section 3. A real data analysis of the proposed method for the data from the Alzheimer’s Neuroimaging Initiative (ADNI) is given in Section 4 and Section 5 provides a discussion of the proposed method and future work.

Notations: Throughout the remainder of the paper, we use bold lower-cased alphabetic letters and bold lower-cased Greek letters to denote non-random vectors. Capitalized alphabetic letters and Greek letters will be used to denote random variables, while bold capitalized alphabetic letters and Greek letters will be used to denote matrices or random vectors. For a matrix $\boldsymbol{A}$, $\left\|{\boldsymbol{A}}\right\|_{op}$ is used to denote the operator norm, which is the largest singular value of $\boldsymbol{A}$. For a vector $\boldsymbol{a}$, $\left\|{\boldsymbol{a}}\right\|$ is used to denote the regular Euclidean norm. For a square matrix $\boldsymbol{A}$, $l_{1}(\boldsymbol{A})\geq l_{2}(\boldsymbol{A})\geq\cdots\geq l_{n}(\boldsymbol{A})$ are used to denote the eigenvalues of $\boldsymbol{A}$ and $\left\|{\boldsymbol{A}}\right\|_{F}$ and $\left\|{\boldsymbol{A}}\right\|_{op}$ denote the Frobenius norm and the operator norm of $\boldsymbol{A}$ respectively.

Consider the following semi-parametric model for the population:

where $\boldsymbol{Z}\in\mathbb{R}^{p}$ is a vector containing $p$ SNPs. Suppose that $g\in\mathcal{H}_{K}$ with $\mathcal{H}_{K}$ being a reproducing kernel Hilbert space associated with the kernel function $K(\cdot,\cdot)$. Now, given the observed data pairs $\{(Y_{i},\boldsymbol{Z}_{i})\}_{i=1}^{n}$ from a random sample, the sample version of model (2) is given by

Here we assume that $\epsilon_{i}$ are i.i.d with mean zero and variance $\sigma_{\epsilon}^{2}$. Estimation of $g\in\mathcal{H}_{K}$ can be accomplished via minimizing the penalized least squares loss function

Thanks to the Representer Theorem [6, 7], the general solution for the nonparametric function $g(\cdot)$ in (4) can be expressed as

so that model (3), in a matrix-vector form, can be expressed as

where $\boldsymbol{K}=[K(\boldsymbol{Z}_{i},\boldsymbol{Z}_{j})]_{i,j=1}^{n}$ is a kernel matrix constructed based on the observed data $\boldsymbol{Z}_{1},\ldots,\boldsymbol{Z}_{n}$ and the optimization problem (4) becomes

Set $Q(\boldsymbol{\alpha})=\frac{1}{n}(\boldsymbol{Y}-\boldsymbol{K\alpha})^{T}(\boldsymbol{Y}-\boldsymbol{K\alpha})$. Taking derivatives with respect to $\boldsymbol{\alpha}$ gives

By letting it equal to $\boldsymbol{0}$ and the solution to the equation $\frac{\partial Q}{\partial\boldsymbol{\alpha}}=\boldsymbol{0}$ is given by

We estimate the broad sense heritability by using $\hat{\boldsymbol{\alpha}}$ as given above and using $\hat{\sigma}_{\epsilon}^{2}=\frac{1}{n}(\boldsymbol{Y}-\boldsymbol{K}\hat{\boldsymbol{\alpha}})^{T}(\boldsymbol{Y}-\boldsymbol{K}\hat{\boldsymbol{\alpha}})$. Then the broad sense heritability can be estimated via

Now, we will look at the theoretical properties of the proposed heritability estimator $\hat{H}^{2}$ by providing asymptotic upper bound and lower bound. For simplicity, we denote $\boldsymbol{g}=[g(\boldsymbol{Z}_{1}),\ldots,g(\boldsymbol{Z}_{n})]^{T}$, $\bar{\boldsymbol{J}}_{n}=\frac{1}{n}\boldsymbol{J}_{n}$ and $\hat{\sigma}_{g}^{2}=\frac{1}{n-1}\hat{\boldsymbol{\alpha}}\boldsymbol{K}(\boldsymbol{I}_{n}-\bar{\boldsymbol{J}}_{n})\boldsymbol{K}\hat{\boldsymbol{\alpha}}$. To start with, let us decompose $\hat{\sigma}_{g}^{2}$ and $\hat{\sigma}_{\epsilon}^{2}$ as follows:

where

Similarly,

where

Our first observation is that the cross-product term $\frac{1}{n-1}I_{3}^{g}$ and $\frac{1}{n}I_{e}^{\epsilon}$ are asymptotically negligible under conditions (C1) and (C2). (C1) assumes that the random errors are sub-Gaussian random variables and (C2) requires that the underlying function $g$ is bounded. These assumptions are common in literature.

• (C1)

  (Sub-Gaussian Errors) $\epsilon_{1},\ldots,\epsilon_{n}$ are i.i.d. sub-gaussian random variables with $\max_{i}\left\|{\epsilon_{i}}\right\|_{\psi_{2}}\leq C$ for some constant $C>0$.

• (C2)

  (Bounded Signals) $\left\|{g}\right\|_{\infty}\leq M$ for some $M>0$.

Since $I_{3}^{g}$ and $I_{3}^{\epsilon}$ are asymptotically negligible, the asymptotic performance of $\hat{H}^{2}$ is solely determined by the asymptoci behavior of $I_{i}^{g}$ and $I_{i}^{\epsilon}$, $i\in\{1,2\}$. We first focus on the asymptotic performance of $I_{2}^{g}$ and $I_{2}^{\epsilon}$.

For $I_{1}^{g}$ and $I_{1}^{\epsilon}$, since both terms are quadratic form, so they naturally are naturally lower bounded from 0. To provide nontrivial bounds for these quadratic forms, we further impose the following assumptions. Condition (C3) requires that $\boldsymbol{g}$ and $\boldsymbol{v}_{1}$ are not orthogonal so that $\boldsymbol{g}$ does not lie in the null space of $\boldsymbol{K}$. Similarly, (C3) also requires that neither $\boldsymbol{v}_{1}$ nor $\boldsymbol{g}$ lies in the null space of $\boldsymbol{J}_{n}$. Condition (C4) requires that there is sufficient cap between the largest eigenvalue and the second largest eigenvalue of the kernel matrix $\boldsymbol{K}$.

• (C3)

  There exist constants $c\in(0,1]$ such that

  	$\displaystyle\left|{\boldsymbol{v}_{1}^{T}\boldsymbol{1}_{n}}\right|$	$\displaystyle\geq c\sqrt{n},$
  	$\displaystyle\left|{\boldsymbol{v}_{1}^{T}\boldsymbol{g}}\right|$	$\displaystyle\geq c\left\|{\boldsymbol{g}}\right\|\textrm{ a.s.},$
  	$\displaystyle\left|{\boldsymbol{1}_{n}^{T}\boldsymbol{g}}\right|$	$\displaystyle\geq c\sqrt{n}\left\|{\boldsymbol{g}}\right\|\textrm{a.s.}$

  where $\boldsymbol{v}_{1}$ is the eigenvector corresponding to the largest eigenvalue of $\boldsymbol{K}$.

• (C4)

  (Spectral Gap) There exists constant $\alpha\in(\max(2c^{2}-1,0),1]$ such that

  $$\frac{l_{1}(\boldsymbol{K})}{l_{2}(\boldsymbol{K})}>\frac{\alpha+1-c^{2}}{c^{2}},\textrm{ a.s.}$$

  where $l_{2}(\boldsymbol{K})$ is the second largest eigenvalue of the kernel matrix $\boldsymbol{K}$.

As one can see from Proposition 3, to obtain a nontrivial lower bound, besides the conditions (C3) and (C4), one also requires that the regularization parameter $\lambda_{n}$ should be chosen properly as well. In addition, as a consequence of Proposition 3, the quantity $\frac{1}{n}(\boldsymbol{1}_{n}^{T}\boldsymbol{K}(\boldsymbol{K}+n\lambda_{n})^{-1}\boldsymbol{g})^{2}$ behaves similarly as $\frac{1}{n}(\boldsymbol{1}_{n}^{T}\boldsymbol{g})^{2}$. This observation is given in the following proposition.

As a result, the bounds for $\hat{\sigma}_{\epsilon}^{2}$ and $\hat{\sigma}_{g}^{2}$ further provide an upper bound and a lower bound for $\hat{\sigma}_{\epsilon}^{2}/\hat{\sigma}_{g}^{2}$, which is an important quantity for $\hat{H}^{2}$ since $\hat{H}^{2}=1/(1+\hat{\sigma}_{\epsilon}^{2}/\hat{\sigma}_{g}^{2})$. In addition, by careful choice of the regularization parameter $\lambda_{n}$, it can be shown that the bounds provided covers the true value $\sigma_{\epsilon}^{2}/\sigma_{g}^{2}$.

To evaluate the performance of the proposed estimators, we started by creating populations, from which response variables were generated based on the following model:

where $N$ is the population size specified in the following subsections; $\epsilon_{i}$ are i.i.d. normal random variables with mean 0 and standard deviation 0.5. The genetic variants $\boldsymbol{Z}_{i}\in\mathbb{R}^{p}$ were generated under two different scenarios: (1). The genetic variants simulated based on Hardy-Weinberg equilibrium (HWE) and (2). The genetic variants sampled from the 1000-Genome Project [1]. Three underlying functions were used in the simulation studies:

1. 1.

   Linear function: $g(\boldsymbol{Z}_{i})=2\boldsymbol{Z}_{i}^{T}\boldsymbol{\beta}+5$.

2. 2.

   Quadratic function: $g(\boldsymbol{Z}_{i})=(\boldsymbol{Z}_{i}^{T}\boldsymbol{\beta})^{2}$.

3. 3.

   Trigonometric function: $g(x\boldsymbol{Z}_{i})=\sin(\boldsymbol{Z}_{i}^{T}\boldsymbol{\beta})+2\boldsymbol{Z}_{i}^{T}\boldsymbol{\beta}$.

The genetic effect $\boldsymbol{\beta}\sim\mathcal{N}_{p}(\boldsymbol{0},\sigma_{g}^{2}\boldsymbol{I}_{n})$ and the variance component $\sigma_{g}$ is chosen so that the population broad-sense heritability lies between 0.2 and 0.8. In each scenario, we considered two cases: the low-dimensional case ($N>p$) and the high-dimensional case ($N<p$). We now explain in detail how the $\boldsymbol{Z}_{i}$’s were simulated in these two scenarios.

Alzheimer’s disease (AD) is the most common neurodegenerative disease and has substantial genetic components [5, 16]. AD primarily targets the hippocampus region of the brain, which plays an important role in learning and memory. In the early stages of AD, the hippocampus suffers from volume loss [17, 12]. Additionally, studies have shown that all hippocampal subregion volumes are highly heritable [19]. Accurately estimating genetic heritability can help us understand more about the genetic inheritance of AD and support the development of effective genetic treatments for the disease.

We applied our proposed genetic heritability estimator to whole genome sequencing data from the Alzheimer’s Disease Neuroimaging Initiative (ADNI). A total of 808 samples from the screening and baseline of the ADNI1 and ADNI2 studies included whole genome sequencing data. After retaining SNPs that were called or imputed with $r^{2}>0.9$, we performed quality control by removing all SNPs with MAF ¡0.05 and any SNPs that did not pass the Hardy-Weinberg equilibrium test ($p$-value $<1e-6$). Combining individuals having phenotype information, 576 individuals remain for the analysis.

We began by regressing the natural logarithm of hippocampal volume on several common clinical covariates, including age, gender, education status, and the number of alleles in APOE$\varepsilon$4. The residuals were then used as the response variable in the kernel ridge regression. The three kernel matrices used in the simulation studies—linear kernel, polynomial kernel of degree 2, and Gaussian kernel-were constructed based on the SNP matrix.

After applying our method to estimate the heritability based on the residuals using $n\lambda_{n}=0.8$, we obtained $\hat{H}^{2}=0.028$ for the linear kernel; $\hat{H}^{2}=0.115$ for the polynomial kernel and $\hat{H}^{2}=0.024$ for the Gaussian kernel. As in the simulation studies, the polynomial kernel gives the best performance among the three kernel matrices. We hypothesize that the reason of seeing these small numbers is that the age predictor has explained much of the variations in the response.

In this paper, we proposed an estimator for the broad-sense genetic heritability of continuous phenotypes using kernel ridge regression. Theoretically, we provided an upper bound and a lower bound for the heritability estimator under certain regularity conditions. Consequently, we demonstrated that with an appropriate choice of the regularization parameter $\lambda_{n}$ in the kernel ridge regression, the upper and lower bounds for the estimator encompass the underlying heritability. Simulation studies across various scenarios suggest that the proposed estimator can indeed capture the underlying heritability.

In the methodology section, we provided a general result on the bounds for our heritability estimator. It is worth noting that in some special cases, the bounds could be simpler. For instance, when $r:=\textrm{rank}(\boldsymbol{K})=o(n)$, it then follows from the von Neumann trace inequality [13] that as $n\to\infty$,

Such a case occurs when the kernel matrix is rank deficient, which is generally true in the low-dimensional case. Consequently, $\hat{\sigma}_{g}^{2}=\frac{1}{n-1}I_{1}^{g}+o_{p}(1)$. In other words, the integral term for the empirical spectral distribution of $\boldsymbol{K}$ will vanish.

Despite the popularity of the Gaussian kernel in practice, the results of the simulation studies and real data analysis suggest that the polynomial kernel generally provides the best performance. In reality, for a given dataset, there is little prior knowledge about which kernel matrix to use. As a future project, we will explore using machine learning techniques, such as neural networks, to learn from the data and self-determine the best kernel matrix. Additionally, estimators obtained from kernel ridge regression are usually biased, and this bias may deteriorate the performance of the proposed estimator. Whether debiasing techniques can improve the performance of heritability estimation will also be part of our future work.

Acknowledgements

The authors would like to thank MathWorks at Texas State University for supporting this research.

• •

  Funding: No funding was received for conducting this study.

• •

  Competing interests: The authors declare they have no financial interests. In terms of non-financial interests, Xiaoxi Shen has served as an editorial review board member for Frontiers in Systems Biology.

• •

  Ethics approval and consent to participate: Not applicable.

• •

  Consent for publication: Not applicable.

• •

  Data availability: The data used in the simulation studies can be found at https://github.com/SxxMichael/KRR-Heritability

• •

  Materials availability: Not applicable

• •

  Code availability: The code used to implement the methods can be found at https://github.com/SxxMichael/KRR-Heritability

• •

  Author contribution: Conceptualization: Olivia Bley, Elizabeth Lei, Andy Zhou; Methodology: Olivia Bley, Elizabeth Lei, Andy Zhou, Xiaoxi Shen; Formal analysis and investigation: Xiaoxi Shen; Writing - original draft preparation: Olivia Bley, Elizabeth Lei, Andy Zhou; Writing - review and editing: Olivia Bley, Elizabeth Lei, Andy Zhou, Xiaoxi Shen; Supervision: Xiaoxi Shen.