Population Stratification Sensitivity in SNP Heritability Estimates

Why Ancestry Structure Matters

SNP heritability describes a fraction of phenotypic variation attributable to the genetic effects tagged by a specified SNP set in a specified population. It is not a property of a phenotype that is automatically constant across samples. Ancestry composition can matter because allele frequencies and linkage disequilibrium patterns differ across ancestry groups, and because ancestry can also be associated with environmental, social, or measurement differences in the phenotype.

The central statistical concern is that genetic similarity may reflect ancestry as well as the genetic relationships relevant to the intended heritability interpretation. If individuals with similar ancestry also share environmental conditions that shift their phenotype, a genetic similarity matrix can align with that non-genetic phenotype pattern. GREML may then attribute some of the pattern to its genetic random effect. In GWAS summary statistics, ancestry-related differences can similarly contribute to test-statistic inflation; how much appears in the LD score regression slope or intercept depends on how the confounding relates to LD scores.

Definition: Population stratification is structure in a sample associated with ancestry or population membership. It becomes a concern for SNP heritability when that structure induces phenotype differences or associations that are not adequately represented by the genetic model and its adjustment. An ancestry-sensitive estimate is one that changes materially under reasonable adjustments or analyses targeting more homogeneous ancestry groups.

Sensitivity is not itself proof that an estimate is biased: a change may also reflect a changed target population, lower precision, or differences in allele frequencies and LD. Conversely, a stable estimate does not prove that all relevant stratification has been removed. The goal is to understand which sources of structure could affect the estimate, assess whether diagnostics reveal them, and report the scope and limits of the analysis.

How Structure Can Enter GREML

A basic GREML model for a quantitative phenotype is

\mathbf{y}=\mathbf{X}\boldsymbol{\beta}+\mathbf{g}+\boldsymbol{\varepsilon}, \qquad \mathbf{g}\sim N(\mathbf{0},\sigma_g^2\mathbf{K}), \qquad \boldsymbol{\varepsilon}\sim N(\mathbf{0},\sigma_e^2\mathbf{I}).

Here \(\mathbf{X}\) contains fixed-effect covariates, \(\mathbf{K}\) is the genetic relationship matrix (GRM) constructed from the chosen SNP set, and the genetic and residual effects are assumed independent. The model attributes covariance between individuals to the genetic component in proportion to their GRM similarity, after accounting for the fixed effects. The estimated SNP heritability is typically a ratio of the genetic variance to the modeled phenotypic variance, with the precise scale depending on the phenotype and model.

Suppose two ancestry groups differ in their average phenotype because of an environmental exposure. If members of the same group also tend to have more similar genotypes, entries of \(\mathbf{K}\) can reflect group membership. A phenotype pattern aligned with that similarity can resemble genetic covariance, even if the group mean difference is not caused by the SNP effects being estimated. This is a form of model misspecification: the assumed genetic covariance does not adequately separate the genetic signal of interest from ancestry-associated structure.

Including ancestry principal components (PCs), or other justified ancestry covariates, in \(\mathbf{X}\) can adjust for phenotype variation associated with those measured axes. It is not a guarantee against stratification. PCs may not capture all relevant structure; a linear adjustment may not capture nonlinear or localized patterns; and within-ancestry environmental or technical confounding may remain. Covariate choice also has to be considered alongside the study design and phenotype measurement.

Worked Example: A Group Difference and Its Potential Confounding

Consider hypothetical data from 600 individuals, split evenly into two ancestry groups. Suppose the phenotype has within-group variance 1.00 in each group, with group means 0.00 and 0.40 standardized units. Let \(Z=0\) for the first group and \(Z=1\) for the second.

The overall mean is \((0.00+0.40)/2=0.20\). The variance of \(Z\) is \(0.5(0-0.5)^2+0.5(1-0.5)^2=0.25\). The covariance of \(Z\) and the phenotype is \(0.5(0-0.5)(0.00-0.20)+0.5(1-0.5)(0.40-0.20)=0.10\).

The total phenotype variance is the within-group variance plus the between-group contribution: \(1.00+0.5(0.00-0.20)^2+0.5(0.40-0.20)^2=1.04\). Thus the group indicator–phenotype correlation is \(0.10/\sqrt{0.25\times1.04}=0.1961\). If genetic similarity also distinguishes the groups, this pattern could affect a GREML estimate unless the relevant structure is adequately modeled. The calculation demonstrates a possible alignment; it does not establish that the group difference is genetic.

What Ancestry Covariates and Within-Ancestry Analyses Do

Ancestry PCs summarize major axes of genetic variation in the analyzed genotypes. Including selected PCs as fixed effects can reduce phenotype mean differences associated with those axes before estimating the residual genetic covariance. The adjustment is useful when the PCs capture structure relevant to the phenotype, but it is not a universal correction: results can depend on which PCs and other covariates are included, and residual confounding can persist.

A within-ancestry analysis restricts the sample to a group or a more genetically homogeneous subset and estimates heritability there. This can reduce between-group structure, but does not ensure that all structure within the subset is absent. It also changes the analyzed population and usually reduces sample size, increasing uncertainty. Differences in allele frequencies, LD, phenotype distributions, or SNP coverage can mean that the within-group estimate targets a different tagged genetic variance from the pooled estimate.

Consequently, a pooled estimate with PCs and a within-ancestry estimate are informative sensitivity analyses, not automatically two estimates of precisely the same quantity. Their interpretation should state the population, SNP set, phenotype scale, and adjustment used. Comparing them is most useful when the analysis design is otherwise as consistent as possible.

Worked Example: Interpreting a GREML Sensitivity Pattern

Suppose a hypothetical quantitative-trait analysis reports an estimate of \(0.32\) with standard error \(0.06\) in the pooled sample without ancestry PCs, \(0.24\) with standard error \(0.06\) after adding PCs, and \(0.20\) with standard error \(0.09\) in a within-ancestry subset. All estimates are on the same stated phenotype scale and use the same SNP set.

The first change is \(0.24-0.32=-0.08\): the PC-adjusted estimate is lower by 0.08. The subset estimate is \(0.20-0.32=-0.12\) relative to the unadjusted pooled estimate. These directions are consistent with the possibility that ancestry-related structure affected the unadjusted estimate, but the numbers alone do not establish the cause.

The within-ancestry estimate has a larger standard error, \(0.09\) rather than \(0.06\), as can occur with a smaller sample. It also targets a subset population, so its difference from the pooled estimate may reflect both reduced structure and changed genetic or phenotypic composition. The defensible conclusion is that the estimated heritability is sensitive to these analysis choices and that the sensitivity merits investigation—not that one number is proven to be the unbiased value.

Population Structure in LD Score Regression

For a GWAS SNP \(j\), a standard LD score regression expectation can be written approximately as

E[\chi_j^2]\approx 1+a+\frac{N h_{\mathrm{SNP}}^2}{M}l_j, \qquad l_j=\sum_k r_{jk}^2.

Here \(N\) is the GWAS sample size, \(M\) is the number of SNPs in the modeled set, \(l_j\) is the LD score for SNP \(j\), and \(h_{\mathrm{SNP}}^2\) is the SNP heritability represented in the regression under its assumptions. The slope is \(N h_{\mathrm{SNP}}^2/M\); the intercept is \(1+a\). In the standard interpretation, an excess component of inflation approximately unrelated to LD score contributes to the intercept, whereas polygenic signal related to LD score contributes to the slope.

Population stratification does not necessarily appear only in the intercept. If its association with test statistics is correlated with LD score, some of its contribution can affect the slope. In that case the slope-based heritability estimate can be biased even when the intercept is not strikingly elevated. The intercept is therefore a calibration diagnostic, not a complete test that rules out ancestry-related bias.

LD scores should also be appropriate to the GWAS ancestry. Differences in LD across populations can make a reference panel poorly matched to the sample, undermining the expected relationship between association statistics and LD score. A reasonable interpretation therefore considers ancestry matching, the GWAS covariates and sample, and other diagnostics together rather than relying on one regression output.

Worked Example: When the Intercept Does Not Capture All Inflation

For illustration, suppose a hypothetical LDSC analysis has \(N=100{,}000\) individuals and \(M=1{,}000{,}000\) SNPs. If the target SNP heritability is \(0.20\), the expected polygenic slope under the standard model is \(N h_{\mathrm{SNP}}^2/M=(100{,}000\times0.20)/1{,}000{,}000=0.020\).

Now suppose stratification contributes a constant 0.10 to expected chi-square statistics, independently of LD score. At LD scores 5, 10, and 15, the expected values are \(1+0.10+0.020(5)=1.20\), \(1+0.10+0.020(10)=1.30\), and \(1+0.10+0.020(15)=1.40\). A line through these values has intercept 1.10 and slope 0.020. In this simplified case, the intercept reflects the constant excess, and the slope remains consistent with \(h_{\mathrm{SNP}}^2=0.020M/N=0.20\).

Instead, suppose the expected values are 1.20, 1.35, and 1.50. The slope is \((1.50-1.20)/(15-5)=0.030\), and the intercept is \(1.20-0.030(5)=1.05\). Converting the slope gives \(0.030M/N=0.30\). If the extra inflation increased with LD score because of confounding, this slope could overstate the target heritability; the intercept of 1.05 would not reveal the full problem. These constructed values illustrate the logic, not a universal numerical correction.

Diagnostics and a Practical Sensitivity Workflow

No single diagnostic settles whether population stratification has materially affected a heritability estimate. A useful assessment asks whether measured ancestry structure is related to the phenotype, whether adjustment changes the estimate, and whether the analysis is appropriate to the ancestry composition and genetic data.

Describe the analyzed sample.
Report relevant ancestry composition, exclusions, phenotype scale, SNP set, and the covariates used. Clarify whether the target is a pooled or within-ancestry population.

Inspect measured structure.
Examine ancestry PCs and other available population indicators. Check whether the phenotype varies with those axes and consider whether the pattern could reflect environmental or technical differences.

Fit reasonable adjusted analyses.
Compare estimates with justified ancestry covariates and, where scientifically appropriate, within-ancestry estimates. Keep other modeling choices consistent when possible and report uncertainty for each estimate.

Check LDSC assumptions and calibration.
For summary-statistic analyses, assess the ancestry match of the LD reference and interpret slope and intercept jointly with other diagnostics. Do not treat an intercept near one as proof that stratification is absent.

State the limits of the comparison.
Explain changes in sample size or target population and distinguish a sensitivity pattern from proof of a particular bias mechanism.

Worked Example: A Phenotype–PC Diagnostic

In a hypothetical pooled sample, a regression of a standardized phenotype on PC1 gives a coefficient of \(0.18\) phenotype standard deviations per one standard deviation of PC1, with standard error \(0.06\). The test statistic is \(0.18/0.06=3.00\). Under a two-sided standard normal approximation, the p-value is about \(0.0027\). This is evidence that the phenotype mean varies along the measured ancestry axis in this sample.

After adding justified ancestry covariates, suppose a corresponding residual diagnostic gives a coefficient of \(0.03\) with standard error \(0.06\). Its test statistic is \(0.03/0.06=0.50\), with a two-sided normal-approximation p-value of about \(0.617\). The diagnostic no longer shows a clear association along that measured axis. This does not prove that all stratification is removed: it only assesses the measured relationship under this model, and the PC association itself does not establish whether the cause is genetic, environmental, or technical.

Common Mistake and Key Interpretation

Wrong reasoning: “The LDSC intercept is close to one, so population stratification cannot bias the heritability estimate.”
Correction: An intercept is most informative about inflation that is approximately independent of LD score. Confounding correlated with LD score can affect the slope, so ancestry matching and other diagnostics remain important.
Wrong reasoning: “The estimate fell after adding PCs, so the original estimate was entirely stratification bias.”
Correction: A change indicates sensitivity to the model adjustment. It can be consistent with reduced confounding, but changes in the modeled target or other sample characteristics may also matter.

Key takeaway: Assess ancestry sensitivity by combining the GREML model and its covariates, phenotype–ancestry diagnostics, pooled and within-ancestry analyses, and LDSC calibration and reference matching. Interpret changes as evidence about robustness and scope, not as automatic proof of bias or its absence.

Check Your Understanding

Use the model assumptions and examples above to answer the following questions.

How can a phenotype difference between ancestry groups contribute to a GREML estimate if genetic similarity also reflects group membership?
Why does adding ancestry PCs not guarantee that all population stratification has been removed?
In LD score regression, what type of inflation is the intercept intended to capture, and why can stratification still bias the slope?
Why might a within-ancestry estimate differ from a pooled, PC-adjusted estimate even if both are correctly calculated?
What does a phenotype association with a PC establish, and what does it not establish?