Interaction-LD score regression (i-LDSC) is a statistical framework for estimating heritability (i.e. the proportion of trait variance attributable to genetic variance). Here, we will give an overview of the i-LDSC method and its corresponding software, as well as detail how its underlying model differs from that of LDSC (Bulik-Sullivan et al., 2015b). We will assume that we are analyzing a GWAS dats set $𝒟=\{\mathbf{𝐗},\mathbf{𝐲}\}$ where $\mathbf{𝐗}$ is an $N\times J$ matrix of genotypes with $J$ denoting the number of SNPs (each of which is encoded as {0, 1, 2} copies of a reference allele at each locus $j$) and $\mathbf{𝐲}$ is an $N$-dimensional vector of measurements of a quantitative trait. The i-LDSC framework only requires summary statistics of individual-level data: namely, marginal effect size estimates for each SNP $\widehat{\mathit{𝜷}}$ and a sample LD matrix $\mathbf{𝐑}$ (which can be provided via reference panel data).

We begin by considering the following generative linear model for complex traits

where $b_0$ is an intercept term; $\mathit{𝜷}=({\beta }_1,\mathrm{\dots },{\beta }_J)$ is a $J$-dimensional vector containing the true additive effect sizes for an additional copy of the reference allele at each locus on ${\displaystyle \mathbf{y}}$; ${\displaystyle \mathbf{W}}$ is an $N\times M$ matrix of (pairwise) cis-acting SNP-by-SNP statistical interactions between some subset of causal SNPs, where columns of this matrix are assumed to be the Hadamard (element-wise) product between genotypic vectors of the form ${\mathbf{𝐱}}_j\circ {\mathbf{𝐱}}_k$ for the $j$-th and $k$-th variants; $\mathit{𝜽}=({\theta }_1,\mathrm{\dots },{\theta }_M)$ is an $M$-dimensional vector containing the interaction effect sizes; $\mathit{𝜺}$ is a normally distributed error term with mean zero and variance scaled according to the proportion of phenotypic variation not explained by genetic effects (Bulik-Sullivan et al., 2015b), which we will refer to as the broad-sense heritability of the trait denoted by $H^2$; and $\mathbf{𝐈}$ denotes an $N\times N$ identity matrix. For convenience, we will assume that the genotype matrix (column-wise) and the trait of interest have been mean-centered and standardized (Strandén and Christensen, 2011; de Los Campos et al., 2013; Zhou et al., 2013). Lastly, we will let the intercept term b₀ be a fixed parameter and we will assume that the effect sizes are each normally distributed with variances proportional to their individual contributions to trait heritability (Yang et al., 2010; Wu et al., 2011; Zhou et al., 2013; Crawford et al., 2017)

Effectively, we say that $𝕍[\mathbf{𝐗}\mathit{𝜷}]={\phi }_{\beta }^2$ is the proportion of phenotypic variation contributed by additive SNP effects under the generative model, while $𝕍[\mathbf{𝐖}\mathit{𝜽}]={\phi }_{\theta }^2$ makes up the proportion of phenotypic variation contributed by genetic interactions. While the appropriateness of treating genetic effects as random variables in analytical derivations has been questioned (de Los Campos et al., 2015), later, we will justify the theory presented here with simulation results showing that i-LDSC accurately recovers non-additive genetic variance in Equation 1 under a broad range of conditions.

There are two key takeaways from the generative model specified above. First, Equation 2 implies that the additive and non-additive components in Equation 1 are orthogonal to each other. In other words, $𝔼[{\mathit{𝜷}}^{⊺}{\mathbf{𝐗}}^{⊺}\mathbf{𝐖}\mathit{𝜽}]=𝔼[{\mathit{𝜷}}^{⊺}]{\mathbf{𝐗}}^{⊺}\mathbf{𝐖}𝔼[\mathit{𝜽}]=\mathrm{𝟎}$. This is important because it means that there is a unique partitioning of genetic variance when studying a trait of interest. The second key takeaway is that the genotype matrix $\mathbf{𝐗}$ and the matrix of genetic interactions $\mathbf{𝐖}$ themselves are correlated despite being linearly independent (see Materials and methods). This property stems from the fact that the pairwise interaction between two SNPs is encoded as the Hadamard product of two genotypic vectors in the form ${\mathbf{𝐰}}_m={\mathbf{𝐱}}_j\circ {\mathbf{𝐱}}_k$ (which is a nonlinear function of the genotypes).

A central objective in GWAS studies is to infer how much phenotypic variation can be explained by genetic effects. To achieve that objective, a key consideration involves incorporating the possibility of non-additive sources of genetic variation to be explained by additive effect size estimates obtained from GWAS analyses (Hill et al., 2008). If we assume that the genotype and interaction matrices are correlated, then ${\displaystyle \mathbf{X}}$ and $\mathbf{𝐖}$ are not completely orthogonal (i.e. such that ${\displaystyle {\mathbf{X}}^{⊺}\mathbf{W}\ne 0}$) and the following relationship between the moment matrix ${\displaystyle {\mathbf{X}}^{⊺}\mathbf{y}}$, the observed marginal GWAS summary statistics ${\displaystyle \hat{\beta }}$, and the true coefficient values ${\displaystyle \beta }$ from the generative model in Equation 1 holds in expectation (see Materials and methods)

where $\mathbf{𝐑}$ is a sample estimate of the LD matrix, and $\mathbf{𝐕}$ represents a sample estimate of the correlation between the individual-level genotypes $\mathbf{𝐗}$ and the span of genetic interactions between causal SNPs in $\mathbf{𝐖}$. Intuitively, the term ${\displaystyle \mathbf{V}\theta }$ can be interpreted as the subset of pairwise interaction effects that are tagged by the additive effect estimates from the GWAS study. Note that, when (i) non-additive genetic effects do not contribute to the overall architecture of a trait (i.e. such that ${\displaystyle \mathit{\theta }=\mathbf{0}}$) or (ii) the genotype and interaction matrices $\mathbf{𝐗}$ and $\mathbf{𝐖}$ are uncorrelated, the equation above simplifies to a relationship between LD and summary statistics that is assumed in many GWAS studies and methods (Hormozdiari et al., 2014; Nakka et al., 2016; Zhu and Stephens, 2017; Zhang et al., 2018; Zhu and Stephens, 2018; Cheng et al., 2020; Demetci et al., 2021).

The goal of i-LDSC is to increase estimates of genetic variance by accounting for sources of non-additive genetic effects that can be explained by additive GWAS summary statistics. To do this, we extend the LD score regression framework and the corresponding LDSC software (Bulik-Sullivan et al., 2015b). Here, according to Equation 3, we note that ${\displaystyle \hat{\mathit{\beta }}\sim \mathcal{N}(\mathbf{R}\mathit{\beta }+\mathbf{V}\mathit{\theta },\lambda \mathbf{R})}$ where $\lambda$ is a scale variance term due to uncontrolled confounding effects (Guan and Stephens, 2011; Song et al., 2022). Next, we condition on ${\displaystyle \mathbf{\Theta }=(\mathit{\beta },\mathit{\theta })}$ and take the expectation of chi-square statistics ${\displaystyle {\mathit{\chi }}^2=N\hat{\mathit{\beta }}{\hat{\mathit{\beta }}}^{⊺}}$ to yield

We define ${\displaystyle {\ell }_j=\sum _k{r}_{jk}^2}$ as the LD score for the additive effect of the $j$-th variant (Bulik-Sullivan et al., 2015b), and ${\displaystyle f_j=\sum _m{v}_{jm}^2}$ represents the ‘cis-interaction’ LD score which encodes the pairwise interaction between the $j$-th variant and all other variants within a genomic window that is a pre-specified number of SNPs wide (Crawford et al., 2017), respectively. By considering only the diagonal elements of LD matrix in the first term, similar to the original LDSC approach (Bulik-Sullivan et al., 2015b; Song et al., 2022), we get the following simplified regression model

where ${\displaystyle {\mathit{\chi }}^2=({\chi }_1^2,\dots ,{\chi }_J^2)}$ is a $J$-dimensional vector of chi-square summary statistics, and ${\displaystyle \mathit{\ell }=({\ell }_1,\dots ,{\ell }_J)}$ and ${\displaystyle \mathit{f}=(f_1,\dots ,f_J)}$ are $J$-dimensional vectors of additive and cis-interaction LD scores, respectively. Furthermore, we define the variance components ${\displaystyle \tau =N{\phi }_{\beta }^2/J}$ and ${\displaystyle \vartheta =N{\phi }_{\theta }^2/M}$ as the additive and non-additive regression coefficients of the model, and 1 is the intercept meant to model the bias factor due to uncontrolled confounding effects (e.g. cryptic relatedness structure). In practice, we efficiently compute the cis-interaction LD scores by considering only a subset of interactions between each $j$-th focal SNP and SNPs within a cis-proximal window around the $j$-th SNP. In our validation studies and applications, we base the width of this window on the observation that LD decays outside of a window of 1 centimorgan (cM); therefore, SNPs outside the 1 cM window centered on the $j$-th SNP will not significantly contribute to its LD scores. Note that the width of this window can be relaxed in the i-LDSC software when appropriate. We fit the i-LDSC model using weighted least squares to estimate regression parameters and derive p-values for identifying traits that have significant statistical evidence of tagged cis-interaction effects by testing the null hypothesis ${\displaystyle H_0:\vartheta =0}$. Importantly, under the null model of a trait being generated by only additive effects, the i-LDSC model in Equation 5 reduces to an infinitesimal model (Fisher, 1999) or, in the case some variants have no effect on the trait, a polygenic model.

Lastly, we want to note the empirical observation that the additive ($\mathbf{\ell }$) and interaction $(\mathit{𝒇}$) LD scores are lowly correlated. This is important because it indicates that the presence of cis-interaction LD scores in the model specified in Equation 5 has little-to-no influence over the estimate for the additive coefficient $\tau$. Instead, the inclusion of $\mathit{𝒇}$ creates a multivariate model that can identify the proportion of variance explained by both additive and non-additive effects in summary statistics. In other words, we can interpret $\widehat{\vartheta }$ as an estimate of the phenotypic variation explained by tagged cis-acting interaction effects. The concept of additive genetic effects partially explaining non-additive variation has also described in various studies from quantitative genetics (Hill et al., 2008; Hivert et al., 2021; Mäki-Tanila and Hill, 2014). Under Hardy-Weinberg equilibrium, it can be shown that the additive variance explained by $J$ SNPs takes on the following form (Materials and methods) (Falconer and Mackay, 1983)

The expression for the additive variance ${\sigma }_A^2$ in Equation 6 is important because it represents the theoretical upper bound on the proportion of total phenotypic variance that can be recovered from GWAS summary statistics using the i-LDSC framework. As a result, we use the sum of coefficient estimates $\widehat{\tau }+\widehat{\vartheta }\le {\sigma }_A^2$ to construct i-LDSC heritability estimates. A full derivation of the cis-interaction regression framework and details about its corresponding implementation in our software i-LDSC can be found in Materials and Methods.

We illustrate the power of i-LDSC across different genetic trait architectures via extensive simulation studies (Materials and methods). We generate synthetic phenotypes using real genome-wide genotype data from individuals of self-identified European ancestry in the UK Biobank. To do so, we first assume that traits have a polygenic architecture where all SNPs have a nonzero additive effect. Next, we randomly select a set of causal cis-interaction variants and divide them into two interacting groups (Materials and methods). One may interpret the SNPs in group #1 as being the ‘hubs’ in an interaction map (Crawford et al., 2017), whereas SNPs in group #2 are selected to be variants within some kilobase (kb) window around each SNP in group #1. We assume a wide range of simulation scenarios by varying the following parameters:

• heritability: H² = 0.3 and 0.6;

• proportion of phenotypic variation that is generated by additive effects: $\rho =$ 0.5, 0.8, and 1;

• percentage of SNPs selected to be in group #1: 1%, 5%, and 10%;

• genomic window used to assign SNPs to group #2: ± 10 and ± 100 kb.

We also varied the correlation between SNP effect size and minor allele frequency (MAF; as discussed in Schoech et al., 2019). All results presented in this section are based on 100 different simulated phenotypes for each parameter combination.

Figure 1 demonstrates that i-LDSC robustly detects significant tagged non-additive genetic variance, regardless of the total number of causal interactions genome-wide. Instead, the power of i-LDSC depends on the proportion of phenotypic variation that is generated by additive versus interaction effects (ρ), and its power tends to scale with the window size used to compute the cis-interaction LD scores (see Materials and methods). i-LDSC shows a similar performance for detecting tagged cis-interaction effects when the effect sizes of causal SNPs depend on their minor allele frequency and when we varied the number of SNPs assigned to be in group #2 within 10 kb and 100 kb windows, respectively (Figure 1—figure supplements 1–5).

Figure 1 with 5 supplements see all

Download asset Open asset

Importantly, i-LDSC does not falsely identify putative non-additive genetic effects in GWAS summary statistics when the synthetic phenotype was generated by only additive effects ($\rho =1$). Figure 2 illustrates the performance of i-LDSC under the null hypothesis $H_0:\vartheta =0$, with the type I error rates for different estimation window sizes of the cis-interaction LD scores highlighted in panel A. Here, we also show that, when no genetic interaction effects are present, i-LDSC unbiasedly estimates the cis-interaction coefficient in the regression model to be $\widehat{\vartheta }=0$ (Figure 2B), robustly estimates the heritability (Figure 2C), and provides well-calibrated p-values when assessed over many traits (Figure 2D). This behavior is consistent across different MAF-dependent effect size distributions, and p-value calibration is not sensitive to misspecification of the estimation windows used to generate the cis-interaction LD scores (Figure 2—figure supplements 1–2).

Figure 2 with 2 supplements see all

Download asset Open asset

One of the innovations that i-LDSC offers over the traditional LDSC framework is increased heritability estimates after the identification of non-additive genetic effects that are tagged by GWAS summary statistics. Here, we applied both methods to the same set of simulations in order to understand how LDSC behaves for traits generated with cis-interaction effects. Figure 3 depicts boxplots of the heritability estimates for each approach and shows that, across an array of different synthetic phenotype architectures, LDSC captures less of phenotypic variance explained by all genetic effects. It is important to note that i-LDSC can yield upwardly biased heritability estimates when the cis-interaction scores are computed over genomic window sizes that are too small; however, these estimates become more accurate for larger window size choices (Figure 3—figure supplement 1). In contrast to LDSC, which aims to capture phenotypic variance attributable to the additive effects of genotyped SNPs, i-LDSC accurately partitions genetic effects into additive versus cis-interacting components, which in turn generally leads the ability of i-LDSC to capture more genetic variance. The mean absolute error between the true generative heritability and heritability estimates produced by i-LDSC and LDSC are shown in Supplementary files 1 and 2, respectively. Generally, the error in heritability estimates is higher for LDSC than it is for i-LDSC across each of the scenarios that we consider.

Figure 3 with 13 supplements see all

Download asset Open asset

Next, we perform an additional set of simulations where we explore other common generative models for complex trait architecture that involve non-additive genetic effects. Specifically, we compare heritability estimates from LDSC and i-LDSC in the presence of additive effects, cis-acting interactions, and a third source of genetic variance stemming from either gene-by-environment (G×E) or or gene-by-ancestry (G×Ancestry) effects. Details on how these components were generated can be found in Materials and Methods. In general, i-LDSC underestimates overall heritability when additive effects and cis-acting interactions are present alongside G×E (Figure 3—figure supplement 2) and/or G×Ancestry effects when PCs are included as covariates (Figure 3—figure supplement 3). Notably, when PCs are not included to correct for residual stratification, both LDSC and i-LDSC can yield unbounded heritability estimates greater than 1 (Figure 3—figure supplement 4). Also interestingly, when we omit cis-interactions from the generative model (i.e. the genetic architecture of simulated traits is only made up of additive and G×E or G×Ancestry effects), i-LDSC will still estimate a nonzero genetic variance component with the cis-interaction LD scores (Figure 3—figure supplements 5–7). Collectively, these results empirically show the important point that cis-interaction scores are not enough to recover missing genetic variation for all types of trait architectures; however, they are helpful in recovering phenotypic variation explained by statistical interaction effects. Recall that the linear relationship between (expected) ${\chi }^2$ test statistics and LD scores proposed by the LDSC framework holds when complex traits are generated under the polygenic model where all causal variants have the same expected contribution to phenotypic variation. When cis-interactions affect genetic architecture (e.g. in our earlier simulations in Figure 3), these assumptions are violated in LDSC, but the inclusion of the additional nonlinear scores in i-LDSC help recover the relationship between the expectation of ${\chi }^2$ test statistics and LD.

As a further demonstration of how i-LDSC performs when assumptions of the original LD score model are violated, we also generated synthetic phenotypes with sparse architectures using the spike-and-slab model (Zhou et al., 2013). Here, traits were simulated with solely additive effects, but this time only variants with the top or bottom $\{1,5,10,25,50,100\}$ percentile of LD scores were given nonzero effects (see Materials and methods). Breaking the relationship assumed under the LDSC framework between LD scores and chi-squared statistics (i.e. that they are generally positively correlated) led to unbounded estimates of heritability in all but the (polygenic) scenario when 100% of SNPs contributed to the phenotypic variation (Figure 3—figure supplement 8).

Finally, we performed a set of polygenic simulations to assess if i-LDSC estimates of non-additive genetic variance could be spuriously inflated due to either (i) unobserved additive effects (see, for example, Hemani et al., 2014), (ii) unobserved SNPs that are involved in genetic interactions, or by (iii) nonzero correlation between the additive and interaction effect sizes in the generative model (i.e. breaking the independence assumption in Equation 2). In the first setting, we observed that, across a range of both minor allele frequencies and effect sizes, the omission of causal haplotypes had a negligible effect on the estimated value of the coefficients in i-LDSC (Figure 3—figure supplement 9). We hypothesize this is due to the fact that the simulations were done for polygenic architectures where all SNPs have at least an additive effect. As a result, not observing a small subset of SNPs does not hinder the ability of i-LDSC to estimate genetic variance because the effect size of each SNP is small. If these simulations were conducted for sparse architectures, we would have likely seen a greater impact on i-LDSC; although, we have already shown the LD score regression framework to be uncalibrated for traits with sparse genetic architectures (again see Figure 3—figure supplement 8). In the second setting, we observed that the i-LDSC framework protects against the false discovery of non-additive genetic effects and underestimates the variance component $\vartheta$ when causal variants involved in pairwise interactions were unobserved (Figure 3—figure supplements 10 and 11). As a direct comparison, estimates of the additive variance component $\tau$ in i-LDSC were not affected by the unobserved interacting variants. Lastly, in the third setting, we observed that the mean estimate of the genetic variance in both LDSC and i-LDSC had a slight upward bias as the correlation between additive and interaction effect sizes in the generative model increased; however, the median of these bias estimates was still near zero across all simulated scenarios and their corresponding replicates (Figure 3—figure supplements 12 and 13).

To assess whether pairwise interaction genetic effects are significantly affecting estimates of heritability in empirical biobank data, we applied i-LDSC to 25 continuous quantitative traits from the UK Biobank and BioBank Japan (Supplementary file 3). Protocols for computing GWAS summary statistics for the UK Biobank are described in the Materials and methods; while pre-computed summary statistics for BioBank Japan were downloaded directly from the consortium website (https://pheweb.jp/downloads). We release the cis-acting SNP-by-SNP interaction LD scores used in our analyses on the i-LDSC GitHub repository from two reference groups in the 1000 Genomes: 489 individuals from the European superpopulation (EUR) and 504 individuals from the East Asian (EAS) superpopulation (see also Supplementary files 4 and 5).

In each of the 25 traits, we analyzed in the UK Biobank, we detected significant proportions of estimated genetic variation stemming from tagged pairwise cis-interactions (Table 1). This includes many canonical traits of interest in heritability analyses: height, cholesterol levels, urate levels, and both systolic and diastolic blood pressure. Our findings in Table 1 are supported by multiple published studies identifying evidence of non-additive effects playing a role in the architectures of different traits of interest. For example, Li et al., 2020 found evidence for genetic interactions that contributed to the pathogenesis of coronary artery disease. It was also recently shown that non-additive genetic effects plays a significant role in body mass index (Song et al., 2022). Generally, we find that the traditional LDSC produces lower estimates of trait heritability because it does not consider the additional sources of genetic signal that i-LDSC does (Table 1). In BioBank Japan, 23 of the 25 traits analyzed had a significant nonlinear component detected by i-LDSC — with HDL and triglyceride levels being the only exceptions.

For each of the 25 traits that we analyzed, we found that the i-LDSC heritability estimates are significantly correlated with corresponding estimates from LDSC in both the UK Biobank ($r^2=0.988$, $P=5.936\times {10}^{-24}$) and BioBank Japan ($r^2=0.849$, $P=6.061\times {10}^{-11}$) as shown in Figure 4A. Additionally, we found that the heritability estimates for the same traits between the two biobanks are highly correlated according to both LDSC ($r^2=0.848$, $P=7.166\times {10}^{-11}$) and i-LDSC ($r^2=0.666$, $P=6.551\times {10}^{-7}$) analyses as shown in Figure 4B. After comparing the i-LDSC heritability estimates to LDSC, we then assessed whether there was significant difference in the amount of phenotypic variation explained by the non-additive genetic effect component in the GWAS summary statistics derived from the the UK Biobank and BioBank Japan (i.e. comparing the estimates of $\vartheta$; see Figure 4—figure supplement 1A). We show that, while heterogeneous between traits, the phenotypic variation explained by genetic interactions is relatively of the same magnitude for both biobanks ($r^2=0.372$, $P=0.0119$). Notably, the trait with the most significant evidence of tagged cis-interaction effects in GWAS summary statistics is height which is known to have a highly polygenic architecture.

Figure 4 with 1 supplement see all

Download asset Open asset

The intercepts estimated by LDSC and i-LDSC are also highly correlated in both the UK Biobank and the BioBank Japan (Figure 4—figure supplement 1B). Recall that these intercept estimates represent the confounding factor due to uncontrolled effects. For LDSC, this does include phenotypic variation that is due to unaccounted for pairwise statistical genetic interactions. The i-LDSC intercept estimates tend to be correlated with, but are generally different than, those computed with LDSC — empirically indicating that non-additive genetic variation is partitioned away and is missed when using the standard LD score alone. This result shows similar patterns in both the UK Biobank ($r^2=0.888$, $P=1.962\times {10}^{-12}$) and BioBank Japan ($r^2=0.813$, $P=7.814\times {10}^{-10}$).

Lastly, we performed an additional analysis in the UK Biobank where the cis-interaction scores are included as an annotation alongside 97 other functional categories in the stratified-LD score regression framework and its software s-LDSC (Gazal et al., 2017; Materials and methods). Here, s-LDSC heritability estimates still showed an increase with the interaction scores versus when the publicly available functional categories were analyzed alone, but albeit at a much smaller magnitude (Table 2). The contributions from the pairwise interaction component to the overall estimate of genetic variance ranged from 0.005 for MCHC ($P=0.373$) to 0.055 for HDL ($P=0.575$; Figure 4C and D). Furthermore, in this analysis, the estimates of the non-additive components were no longer statistically significant for any of the traits in the UK Biobank (Table 2). Despite this, these results highlight the ability of the i-LDSC framework to identify sources of ‘missing’ phenotypic variance explained in heritability estimation. Importantly, moving forward, we suggest using the cis-interaction scores with additional annotations whenever they are available as it provides more conservative estimates of the role of non-additive effects on trait architecture.

Our goal in this study is to reanalyze summary statistics from genome-wide association studies (GWAS) and estimate heritability while accounting for both additive genetic associations and tagged interaction effects. We begin by assuming the following generative linear model for complex traits which can be seen as an extended view of Equation 1 in the main text

where $\mathbf{𝐲}$ denotes an $N$-dimensional vector of phenotypic states for a quantitative trait of interest measured in $N$ individuals; $b_0$ is an intercept term; $\mathbf{𝐗}$ is an $N\times J$ matrix of genotypes, with $J$ denoting the number of single nucleotide polymorphism (SNPs) encoded as $\{0,1,2\}$ copies of a reference allele at each locus; $\mathit{𝜷}=({\beta }_1,\mathrm{\dots },{\beta }_J)$ is a $J$-dimensional vector containing the true additive effect sizes for an additional copy of the reference allele at each locus on ${\displaystyle \mathbf{y}}$; ${\displaystyle X_D}$ is an $N\times J$ matrix that represents the dominance for each genotype encoded as $\{0,1,1\}$ with corresponding effect sizes ${\displaystyle \omega }$; ${\displaystyle \mathbf{W}}$ is an $N\times M$ matrix of genetic interactions; $\mathit{𝜽}=({\theta }_1,\mathrm{\dots },{\theta }_M)$ is an $M$-dimensional vector containing the interaction effect sizes; $\mathit{𝜺}$ is a normally distributed error term with mean zero and variance scaled according to the proportion of phenotypic variation not explained by the broad-sense heritability of the trait, denoted by $H^2$; and $\mathbf{𝐈}$ denotes an $N\times N$ identity matrix. Note that the encoding for dominance in ${\mathbf{𝐗}}_D$ was chosen because it imposes orthogonality with the genotype encoding in $\mathbf{𝐗}$ (Purcell et al., 2007; Vitezica et al., 2017; Palmer et al., 2023).

For convenience, we will assume that the genotype matrix (column-wise), the dominance matrix (also column-wise), and trait of interest have all been standardized (Strandén and Christensen, 2011; de Los Campos et al., 2013; Zhou et al., 2013). Furthermore, while the matrix $\mathbf{𝐖}$ could encode any source of non-additive genetic interactions (e.g. gene-by-environmental effects) in theory, we limit our focus in this study to trait architectures that have been generated with contributions stemming from cis-acting statistical SNP-by-SNP (or pairwise) interactions. To that end, we assume that the columns of $\mathbf{𝐖}$ are the Hadamard (element-wise) product between genotypic vectors of the form ${\mathbf{𝐱}}_j\circ {\mathbf{𝐱}}_k$ for the $j$-th and $k$-th variants. We also want to point out that the generative formulation of Equation 7 can also be easily extended to accommodate other fixed effects (e.g. age, sex, or genotype principal components), as well as other random effects terms that can be used to account for sample non-independence due to other environmental factors.

As a final set of assumptions, we will let the intercept term $b_0$ be a fixed parameter while allowing the other coefficients to follow independent Gaussian distributions with variances proportional to their individual contributions to the trait heritability (Yang et al., 2010; Wu et al., 2011; Zhou et al., 2013; Jiang and Reif, 2015; Crawford et al., 2017)

for $j=1,\mathrm{\dots },J$ and $m=1,\mathrm{\dots },M$. The broad-sense heritability of the trait is defined as $H^2={\phi }_{\beta }^2+{\phi }_{\omega }^2+{\phi }_{\theta }^2$. Under the generative model in Equation 7, we then say that $𝕍[\mathbf{𝐗}\mathit{𝜷}]={\phi }_{\beta }^2$ is the proportion of phenotypic variation contributed by additive SNP effects, $𝕍[{\mathbf{𝐗}}_D\mathit{𝝎}]={\phi }_{\omega }^2$ is the proportion of phenotypic variation contributed by dominance effects, and the set of interactions involving some subset of causal SNPs contribute the remaining proportion to the heritability $𝕍[\mathbf{𝐖}\mathit{𝜽}]={\phi }_{\theta }^2$. As we mentioned in the main text, we recognize that the appropriateness of treating genetic effects as random variables in analytical derivations has been questioned (de Los Campos et al., 2015), but our simulation studies show that i-LDSC accurately recovers non-additive genetic variance in Equation 7 under a broad range of conditions.

Assuming that the effect sizes $\{\mathit{𝜷},\mathit{𝝎},\mathit{𝜽}\}$ in Equation 8 follow independent and zero mean Gaussian distributions leads to orthogonality between the additive and non-additive components in Equation 7. Since the genotypes $\mathbf{𝐗}$ and the dominance values ${\mathbf{𝐗}}_D$ are fixed orthogonal matrices, it is straightforward to show that $\text{Cov}[\mathbf{𝐗}\mathit{𝜷},{\mathbf{𝐗}}_D\mathit{𝝎}]=0$ (Vitezica et al., 2017; Palmer et al., 2023). The same relationship can be shown for the additive and the pairwise interaction genetic effects where

with ${\mathbf{𝐱}}_j$ and ${\mathbf{𝐰}}_m$ denoting the $j$-th and $m$-th column of the individual-level genotype matrix $\mathbf{𝐗}$ and the interaction matrix $\mathbf{𝐖}$, respectively. Note that a similar derivation to Equation 9 can also be done for the dominance and pairwise genetic interaction effects. This concept of orthogonality is important because we want to preserve a unique partitioning of genetic variance when modeling a trait of interest.

The design matrices $\mathbf{𝐗}$ and $\mathbf{𝐖}$ in Equation 7 are not linearly dependent because the pairwise interactions between two SNPs are encoded as the Hadamard product of two genotypic vectors in the form ${\mathbf{𝐱}}_j\circ {\mathbf{𝐱}}_k$ (which is a nonlinear function). Linear dependence would have implied that one could find a transformation between a SNP and an interaction term in the form ${\mathbf{𝐰}}_m=c\times {\mathbf{𝐱}}_j$ for some constant $c$. However, despite their linear independence, $\mathbf{𝐗}$ and $\mathbf{𝐖}$ are themselves not orthogonal and still have a nonzero correlation. This implies that the inner product between genotypes and their interactions is nonzero ${\mathbf{𝐗}}^{⊺}\mathbf{𝐖}\ne \mathrm{𝟎}$. To see this, we focus on a focal SNP ${\mathbf{𝐱}}_j$ and consider three different types of interactions:

• Scenario I: Interaction between a focal SNP with itself (${\mathbf{𝐱}}_j\circ {\mathbf{𝐱}}_j$).

• Scenario II: Interaction between a focal SNP with a different SNP (${\mathbf{𝐱}}_j\circ {\mathbf{𝐱}}_k$).

• Scenario III: Interaction between a focal SNP with a pair of different SNPs (${\mathbf{𝐱}}_k\circ {\mathbf{𝐱}}_l$).

The following derivations rely on the fact that: (1) we assume that genotypes have been mean-centered and scaled to have unit variance, and (2) under Hardy-Weinberg equilibrium, SNPs marginally follow a binomial distribution ${\mathbf{𝐱}}_j\sim \text{Bin}(2,p)$ where $p$ represents the minor allele frequency (MAF) (Wray et al., 2007; Lippert et al., 2013).

As previously mentioned, the key to this work is that SNP-level GWAS summary statistics can also tag non-additive genetic effects when there is a nonzero correlation between individual-level genotypes and their interactions (as defined in Equation 7). Throughout the rest of this section, we will use ${\mathbf{𝐗}}^{⊺}\mathbf{𝐗}/N$ to denote the LD or pairwise correlation matrix between SNPs. We will then let $\mathbf{𝐑}$ represent an LD matrix empirically estimated from external data (e.g. directly from GWAS study data, or using a pairwise LD map from a population that is representative of the samples analyzed in the GWAS study). The important property here is the following

where the term $r_{jk}$ is again defined as the Pearson correlation coefficient between the $j$-th and $k$-th SNPs, respectively.

In traditional GWAS studies, summary statistics of the true additive effects $\mathit{𝜷}={({\mathbf{𝐗}}^{⊺}\mathbf{𝐗})}^{-1}{\mathbf{𝐗}}^{⊺}\mathbf{𝐲}$ in Equation 7 are typically derived by computing a marginal least squares estimate with the observed data

There are two key identities that may be taken from Equation 11. The first uses Equation 10 and is the approximate relationship (in expectation) between the moment matrix ${\mathbf{𝐗}}^{⊺}\mathbf{𝐲}$ and the linear effect size estimates $\widehat{\mathit{𝜷}}$:

The second key point combines Equations 10 and 12 to describe the asymptotic relationship between the observed marginal GWAS summary statistics $\widehat{\mathit{𝜷}}$ and the joint coefficient values $\mathit{𝜷}$ where (in expectation)

After some algebra, the above mirrors a high-dimensional regression model (in expectation) where $\widehat{\mathit{𝜷}}=\mathbf{𝐑}\mathit{𝜷}$ with the estimated summary statistics as the response variables and the empirically estimated LD matrix acting as the design matrix (Hormozdiari et al., 2014; Hormozdiari et al., 2016; Zhang et al., 2018; Cheng et al., 2020; Demetci et al., 2021). Theoretically, the resulting coefficients output from this high-dimensional model are the desired true effect size estimates used to generate the phenotype of interest.

When interactions contribute to the architecture of complex traits (i.e. $\mathit{𝜽}\ne \mathrm{𝟎}$), the marginal GWAS summary statistics derived using least squares in Equation 11 will also explain non-additive variation when there is a nonzero correlation between genotypes and their interactions. To see this, we use the concept of ‘omitted variable bias’ (Barreto and Howland, 2005) where the fitted model aims to estimate the true additive coefficients $\mathit{𝜷}$ but does not account for contributions from the non-additive components which also contribute to trait architecture. In this case, we get the following

Since we assume that the genotypes are orthogonal to both the dominance effects in Equation 7, we know that ${\mathbf{𝐗}}^{⊺}{\mathbf{𝐗}}_D=\mathrm{𝟎}$. This simplifies the above to be the following

where the matrix ${\mathbf{𝐗}}^{⊺}\mathbf{𝐖}$(which we showed to be nonzero) can be interpreted as the sample correlation between individual-level genotypes and the cis-interactions between causal SNPs. By taking the expectation using Equations 10 and 12, we get the following alternative (approximate) relationship between the observed marginal GWAS summary statistics $\widehat{\mathit{𝜷}}$ and the true coefficient values $\mathit{𝜷}$

which results from our initial assumption that the residuals are normally distributed with mean zero $𝔼[\mathit{𝜺}]=\mathrm{𝟎}$ in Equation 7. Here, we define $\mathbf{𝐕}$ to represent a sample estimate of the correlation between the individual-level genotypes and the non-additive genetic interaction matrix such that $𝔼[{\mathbf{𝐗}}^{⊺}\mathbf{𝐖}]\approx N\mathbf{𝐕}$. Similar to the LD matrix $\mathbf{𝐑}$, the correlation matrix $\mathbf{𝐕}$ is also assumed to be computed from reference panel data. Intuitively, when $\mathit{𝜽}\ne \mathrm{𝟎}$ there is additional phenotypic variation contributed by pairwise interactions that can be explained by GWAS effect size estimates. Moreover, when $\mathbf{𝐕}\mathit{𝜽}=\mathrm{𝟎}$, then the relationship in Equation 16 converges onto the conventional asymptotic assumption (in expectation) between GWAS summary statistics and the true additive coefficients in Equation 13; Hormozdiari et al., 2014; Hormozdiari et al., 2016; Zhang et al., 2018; Cheng et al., 2020; Demetci et al., 2021.

The concept of additive genetic effects partially explaining non-additive variation has also described in classical quantitative genetics (Hill et al., 2008; Hivert et al., 2021; Mäki-Tanila and Hill, 2014). Consider an individual genotyped at $J$ loci each with major and minor alleles A and B, respectively. Let $p_j$ be the allele frequency of A at the $j$-th locus, $a_j$ denote the additive effect, and ${[aa]}_{jk}$ be the additive-by-additive (pairwise) interaction effect between loci $j$ and $k$, and ${[aaa]}_{jkl}$ represent a third order interaction between loci $j$, $k$, and $l$. For simplicity in presentation, assume that dominance only makes a small contribution to the genetic variance (Palmer et al., 2023; Pazokitoroudi et al., 2021; Zhu et al., 2015). The population mean is given as the following

We follow the assumption that the genetic variation in human complex traits can predominately be explained by additive effects, with the remainder variation being mostly explained by additive-by-additive effects (Weinreich et al., 2018; Jiang and Reif, 2015; Fisher, 1919; Lynch and Walsh, 1998). As a result, we will ignore the higher order interaction terms in Equation 17. Under Hardy-Weinberg equilibrium, we can find the average effect by taking the first derivative of the population mean with respect to the frequency of the increasing allele (Mäki-Tanila and Hill, 2014; Hivert et al., 2021). For the $j$-th SNP, the average effect (including terms up to second-order interaction) is given by the following

which notably contains both the additive effect and a summation of additive-by-additive interactions between pairs of loci. The additive genetic variance for the $j$-th SNP takes on the following form

which is the product of the square of the average effect in Equation 18 and the heterozygosity at $j$-th locus $𝕍[{\mathbf{𝐱}}_j]=2p_j(1-p_j)$ (again assuming that SNPs marginally follow a binomial distribution). The total additive variance is then obtained by summing over the $J$ loci such that ${\sigma }_A^2={\sum }_j{\sigma }_A^2(j)$ (Falconer and Mackay, 1983).

We can derive a parallel construction for additive genetic variance using the generative random effect model presented in Equation 7; Hivert et al., 2021. Here, we will leverage that with genotype data taken for $N$ individuals, ${\sum }_i{x}_{ij}/N=2p_j$. Ignoring the assumed small contributions from dominance effects, the population mean for a quantitative trait $\mathbf{𝐲}$ can be written as the following

To find the average effect for the $j$-th locus, we this time take the first derivative of the population mean in Equation 20 with respect to the allele frequency such that

which, similar to the theoretical form in quantitative genetics, also contains both the additive effect of the $j$-th SNP and additional terms encoding the interaction effect between the $j$-th SNP and all other variants in the data. Once again, under Hardy-Weinberg equilibrium, the additive variance for the $j$-th SNP is found as taking on the following form

where we can explicitly draw connections between the two frameworks by setting ${\beta }_j=a_j$ and ${\theta }_{jk}={[aa]}_{jk}$. Note that when there no non-additive effects (such that $\mathit{𝜽}=\mathrm{𝟎}$), the above reduces to ${\sigma }_A^2={\sum }_{j}2p_j(1-p_j){\beta }_j^2$ which resembles the classical form for the additive genetic variance (Lynch and Walsh, 1998).

In order to derive the interaction LD score (i-LDSC) regression framework, recall that our goal is to recover missing heritability from GWAS summary statistics by incorporating an additional score that measures the non-additive genetic variation that is tagged by genotyped SNPs. To do this, we build upon the LD score regression framework and the LDSC software (Bulik-Sullivan et al., 2015b). Here, we assume nonzero contributions from cis-acting pairwise interaction effects in the generative model of complex traits as in Equation 16, and we use the observed least squares estimates from Equation 11 to compute chi-square statistics ${\chi }_j^2=N{\widehat{\beta }}_j^2$ for every $j=1,\mathrm{\dots },J$ variant in the data. Taking the expectation of these statistics yields

We can simplify Equation 23 in two steps. First, by combining the prior assumption in Equation 8 and the asymptotic approximation in Equation 16, we can show that marginal expectation (i.e. when not conditioning on the true coefficients) $𝔼[{\widehat{\beta }}_j]=0$ for all variants. Second, by conditioning on the generative model from Equation 7, we can use the law of total variance to simplify $𝕍[{\widehat{\beta }}_j]$ where

since ${\mathbf{𝐱}}_j^{⊺}{\mathbf{𝐗}}_D=\mathrm{𝟎}$. Using the same logic from the original LDSC regression framework (Bulik-Sullivan et al., 2015b), we can use Isserlis’ theorem Isserlis, 1918 to write the above in terms of more familiar quantities based on sample correlations

where ${\tilde{r}}_{jk}$ is used to denote the sample correlation between additively-coded genotypes at the $j$-th and $k$-th variants, and ${\tilde{v}}_{jm}$ is used to denote the sample correlation between the genotype of the $j$-th variant and the $m$-th genetic interaction on the phenotype of interest (again see Equation 16). Furthermore, we can use the delta method (only displaying terms up to $𝒪(1/N^2)$) to show that (in expectation)

Next, we can then approximate the quantities in Equation 24 via the following

where ${\mathrm{\ell }}_j$ is the corresponding LD score for the additive effect of the $j$-th variant and f_j represents the “interaction” LD score between the $j$-th SNP and all other variants in the data set (Crawford et al., 2017), respectively. Altogether, this leads to the specification of the univariate framework with the $j$-th SNP

where we define $\tau =N{\phi }_{\beta }^2/J$ as estimates of the additive genetic signal, the coefficient $\vartheta =N{\phi }_{\theta }^2/M$ as an estimate of the proportion of phenotypic variation explained by tagged pairwise interaction effects, and 1 is the intercept meant to model the misestimation due to uncontrolled confounding effects (e.g. cryptic relatedness and population stratification). Similar to the original LDSC formulation, an intercept greater than one means significant bias. Note that the simplification for many of the terms above such as $(1-H^2)/N\approx 1/N$ results from our assumption that the number of individuals in our study is large. For example, the sample sizes for each biobank-scale study considered in the analyses of this manuscript are at least on the order of $N\ge {10}^4$ observations (see Supplementary file 5). Altogether, we can jointly express Equation 27 in multivariate form as

where ${\mathit{𝝌}}^2=({\chi }_1^2,\mathrm{\dots },{\chi }_J^2)$ is a $J$-dimensional vector of chi-square summary statistics, and $\mathbf{\ell }=({\mathrm{\ell }}_1,\mathrm{\dots },{\mathrm{\ell }}_J)$ and $\mathit{𝒇}=(f_1,\mathrm{\dots },f_J)$ are $J$-dimensional vectors of additive and cis-interaction LD scores, respectively. It is important to note that, while ${\mathit{𝝌}}^2$ must be recomputed for each trait of interest, both vectors $\mathbf{\ell }$ and $\mathit{𝒇}$ only need to be constructed once per reference panel or individual-level genotypes (see next section for efficient computational strategies).

To identify summary statistics that have significant tagged interaction effects, we test the null hypothesis $H_0:\vartheta =0$. The i-LDSC software package implements the same model fitting strategy as LDSC. Here, we use weighted least squares to fit the joint regression in Equation 28 such that

where $\mathbf{𝚿}$ is a $J\times J$ diagonal weight matrix with nonzero elements set to values inversely proportional to the conditional variance $𝕍[{\chi }_j^2|{\mathrm{\ell }}_j,f_j]={\psi }_{jj}^{-1}$ to adjust for both heteroscedasticity and over-estimation of the summary statistics for each SNP (Bulik-Sullivan et al., 2015b). Standard errors for each coefficient estimate are derived via a jackknife over blocks of SNPs in the data (Finucane et al., 2015), and we then use those standard errors to derive p-values with a two-sided test (i.e. testing the alternative hypothesis $H_A:\vartheta \ne 0$). It is worth noting that the block-jackknife approach tends to be conservative and yield larger standard errors for hypothesis testing (Efron, 1982). As an alternative, we could first run i-LDSC using the block-jackknife procedure over all traits in a study and then use the average of the standard errors to calculate the statistical significance of coefficient estimates; but we do not explore this strategy here and leave that for future work. The quantitative genetics expression for the additive variance ${\sigma }_A^2$ in Equation 22 is important because it represents the theoretical upper bound on the proportion of phenotypic variance that can be explained from GWAS summary statistics via i-LDSC. Using this relationship, we can write the following (approximate) inequality

For all analyses in this paper, we estimate proportion of phenotypic variance explained by genetic effects using a sum of the coefficients $\widehat{\tau }+\widehat{\vartheta }$ (i.e. the estimated additive component plus the additional genetic variance explained by the tagged pairwise interaction effects).

In practice, cis-interaction LD scores in i-LDSC can be computed efficiently through realizing two key opportunities for optimization. First, given $J$ SNPs, the full matrix of genome-wide interaction effects $\mathbf{𝐖}$ contains on the order of $J(J-1)/2$ total pairwise interactions. However, to compute the cis-interaction score for each SNP, we simply can replace the full $\mathbf{𝐖}$ matrix with a subsetted matrix ${\mathbf{𝐖}}_j$ which includes only interactions involving the $j$-th SNP. Analogous to the original LDSC formulation (Bulik-Sullivan et al., 2015b), we consider only interactive SNPs within a cis-window proximal to the focal $j$-th SNP for which we are computing the i-LDSC score. In the original LDSC model, this is based on the observation that LD decays outside of a window of 1 centimorgan (cM) (Bulik-Sullivan et al., 2015b); therefore, SNPs outside the 1 cM window centered on the $j$-th SNP $j$ will not significantly contribute to its LD score. The second opportunity for optimization comes from the fact that the matrix of interaction effects for any focal SNP, ${\mathbf{𝐖}}_j$, does not need to be explicitly generated. Referencing Equation 24, the i-LDSC scores are defined as ${\mathbf{𝐱}}_j^{⊺}{\mathbf{𝐖}}_j{\mathbf{𝐖}}_j^{⊺}{\mathbf{𝐱}}_j/N^2$. This can be re-written as ${\mathbf{𝐱}}_j^{⊺}({\mathbf{𝐃}}_j{\mathbf{𝐗}}^{(j)}){({\mathbf{𝐃}}_j{\mathbf{𝐗}}^{(j)})}^{⊺}{\mathbf{𝐱}}_j$, where ${\mathbf{𝐃}}_j=\text{diag}({\mathbf{𝐱}}_j)$ is a diagonal matrix with the $j$-th genotype as its nonzero elements (Crawford et al., 2017) and ${\mathbf{𝐗}}^{(j)}$ denotes the subset SNPs within a cis-window proximal to the focal $j$-th SNP. This means that the i-LDSC score for the $j$-th SNP can be simply computed as the following

With these simplifications, the computational complexity of generating i-LDSC scores reduces to that of computing LD scores — modulo a vector-by-vector Hadamard product which, for each SNP, is constant factor of $N$ (i.e. the number of genotyped individuals).

When computing cis-interaction LD scores, the most important decision is choosing the number of interacting SNPs to include in ${\mathbf{𝐗}}^{(j)}$ (or equivalently ${\mathbf{𝐖}}_j$ for each $j$-th focal SNP in the calculation of f_j in Equation 31). The i-LDSC framework considers different estimating windows to account for our lack of a priori knowledge about the ‘correct’ non-additive genetic architecture of traits. Theoretically, one could follow previous work Guan and Stephens, 2011; Carbonetto and Stephens, 2012; Zhou et al., 2013; Zhu and Stephens, 2017; Zhu and Stephens, 2018; Demetci et al., 2021 by considering an $L$-valued grid of possible SNP interaction window sizes. After fitting a series of i-LDSC regressions with cis-interaction LD scores ${\mathit{𝒇}}^{(l)}$ generated under the $L$-different window sizes, we could compute normalized importance weights using their maximized likelihoods via the following

As a final step in the model fitting procedure, we could then compute averaged estimates of the coefficients $\tau$ and $\vartheta$ by marginalizing (or averaging) over the $L$-different grid combinations of estimating windows

This final step can be viewed as an analogy to model averaging where marginal estimates are computed via a weighted average using the importance weights (Hoeting et al., 1999). In the current study, we explore the utility of cis-interaction LD scores generated with different window sizes ± 5, ± 10, ± 25, and ± 50 SNPs around each $j$-th focal SNP. In practice, we find that cis-interaction LD scores that are calculated using larger windows lead to the most robust estimates of heritability while also not over representing the total phenotypic variation explained by tagged non-additive genetic effects (see Figure 3—figure supplement 1). Therefore, unless otherwise stated, we use cis-interaction LD scores calculated with a ± 50 SNP interaction window for all simulations and real data analyses conducted in this work. For a direct comparison between choosing a single window size versus the model averaging strategy described above, see Supplementary files 1 and 2.

The LDSC software computes LD scores using annotations over equally spaced minor allele frequency (MAF) bins. These annotations enable the per trait relationship between the MAF and the effect size of each variant in the genome to vary based on the discrete category (or MAF bin) it is placed into. This additional flexibility is intended to help LDSC be more robust when estimating heritability. The relationship between MAF and effect size is already implicitly encoded in the LDSC formulation since we assume genotypes are normalized. When normalizing by the variance of each SNP (or equivalently its MAF), we make the assumption that rare variants inherently have larger effect sizes. There exists a true functional relationship between MAF and effect size which is likely to be somewhere between the two extremes of (i) normalizing each SNP by its MAF and (ii) allowing the variance per SNP to be dictated by its MAF.

Recent approaches have proposed using a single parameter $\alpha$ to better represent the nonlinear relationship between MAF and variant effect size. The main idea is that this $\alpha$ not only provides the same additional flexibility to LDSC as the MAF-based discrete annotations, but it also empirically yields even more precise heritability estimates (Zabad et al., 2021). Namely, we use

where $a_c(k)$ is the annotation value for the $c$-th categorical bin. The α parameter is unknown in practice and needs to be estimated for any given trait. While standard ranges for α can be used for heritability estimates, we use a restricted maximum likelihood (REML) based method which was recently developed (Schoech et al., 2019).

In the i-LDSC software, we use this α construction to handle the relationship between MAF and variant effect size for two specific reasons. First, by constructing the LD scores using α, we more accurately capture the variation in chi-square test statistics due to additive effects (Zabad et al., 2021). Second, we note that there is correlation between MAF and (i) LD scores, (ii) cis-interaction LD scores, and (iii) trait architecture. To that end, if we do not properly condition on MAF, there becomes additional bias, and we may falsely attribute some amount of variation in the chi-square test statistics to LD or the tagged interaction effects. Therefore, in our formulation, we include an α term on the LD scores to condition on this effect. We demonstrate in simulations that this removes the bias introduced by the relationship between MAF and trait architecture, and it mitigates potential inflation of type I error rates in the i-LDSC test.

In the main text, we analyzed 25 complex traits in both the UK Biobank and BioBank Japan data sets. In order to account for minor allele frequency (MAF) dependent trait architecture, we calculated $\alpha$ values for each trait that had not been analyzed by previous studies (Schoech et al., 2019). The α estimates for each of the 25 traits analyzed in this study are shown in Supplementary file 4. Intuitively, $\alpha$ parameterizes the weighting of the effects of each individual variant given its frequency in the study cohort and can take on values in the range of [–1,0]. More negative values of $\alpha$ indicate that lower frequency variants contribute more to the observed variation in a trait of interest, whereas values of α closer to zero indicate that common variants contribute a greater amount of variation to observed trait values.

We took α values for 11 traits (again see Supplementary file 4) that had previously been calculated from Schoech et al. For the remaining 14 traits analyzed in this study, we followed the estimation protocol described in the same manuscript. Specifically, using the variants passing the quality control step in our pipeline for 25,000 randomly selected individuals in the UK Biobank cohort, we constructed MAF-dependent genetic relatedness matrices for values of $\alpha =\{-1,-0.95,-0.9,\mathrm{\dots },0\}$ using the GRM-MAF-LD software (Schoech, 2018). We then used the GCTA software (Yang et al., 2011) to obtain heritability and likelihood estimates using REML for each $\alpha$-trait pairing. We then fit a trait-specific profile likelihood across the range of α values and estimate the maximum likelihood value of $\alpha$ using a natural cubic spline.

We used a simulation scheme to generate synthetic quantitative traits and SNP-level summary statistics under multiple genetic architectures using real genome-wide data from individuals of self-identified European ancestry in the UK Biobank. Here, we consider phenotypes that have some combination of additive effects, cis-acting interactions, and a third source of genetic variance stemming from either gene-by-environment (G×E) or gene-by-ancestry (G×Ancestry) effects. For each scenario, we select some set of SNPs to be causal and assume that complex traits are generated via the following general linear model

where $\mathbf{𝐲}$ is an $N$-dimensional vector containing all the phenotypes; $\mathbf{𝐗}$ is an $N\times J$ matrix of genotypes encoded as 0, 1, or 2 copies of a reference allele; β is a $J$-dimensional vector of additive effect sizes for each SNP; $\mathbf{𝐖}$ is an $N\times M$ matrix which holds all pairwise interactions between the randomly selected subset of the interacting SNPs with corresponding effects θ is an $N\times K$ matrix of either G×E or G×Ancestry interactions with coefficients $\mathit{𝜸}$; and $\mathit{𝜺}$ is an $N$-dimensional vector of environmental noise. The phenotypic variation is assumed to be $𝕍[\mathbf{𝐲}]=1$. All additive and interaction effect sizes for SNPs are randomly drawn from independent standard Gaussian distributions and then rescaled so that they explain a fixed proportion of the phenotypic variance $𝕍[\mathbf{𝐗}\mathit{𝜷}]+𝕍[\mathbf{𝐖}\mathit{𝜽}]+𝕍[\mathbf{𝐙}\mathit{𝜸}]=H^2$. Note that we do not assume any specific correlation structure between the effect sizes β, θ, and $\mathit{𝜸}$. We then rescale the random error term such that $𝕍[\mathit{𝜺}]=(1-H^2)$. In the main text, we compare the traditional LDSC to its direct extension in i-LDSC. For each method, GWAS summary statistics are computed by fitting a single-SNP univariate linear model via least squares where ${\widehat{\beta }}_j={({\mathbf{𝐱}}_j^{⊺}{\mathbf{𝐱}}_j)}^{-1}{\mathbf{𝐱}}_j^{⊺}\mathbf{𝐲}$ for every $j=1,\mathrm{\dots },J$ SNP in the data. These effect size estimates are used to derive the chi-square test statistics ${\chi }_j^2=N{\widehat{\beta }}_j^2$. We implement both LDSC and i-LDSC with the LD matrix $\mathbf{𝐑}={\mathbf{𝐗}}^{⊺}\mathbf{𝐗}/N$ and the cis-interaction correlation matrix $\mathbf{𝐕}={\mathbf{𝐗}}^{⊺}\mathbf{𝐖}/N$ being computed using a reference panel of 489 individuals from the European superpopulation (EUR) of the 1000 Genomes Project (https://mathgen.stats.ox.ac.uk/impute/data_download_1000G_phase1_integrated.html). The resulting matrices $\mathbf{𝐑}$ and $\mathbf{𝐕}$ are used to compute the additive and cis-interaction LD scores, respectively.