Assortative mating biases marker-based heritability estimators

Abstract

Many traits are subject to assortative mating, with recent molecular genetic findings confirming longstanding theoretical predictions that assortative mating induces long range dependence across causal variants. However, all marker-based heritability estimators implicitly assume mating is random. We provide mathematical and simulation-based evidence demonstrating that both method-of-moments and likelihood-based estimators are biased in the presence of assortative mating and derive corrected heritability estimators for traits subject to assortment. Finally, we demonstrate that the empirical patterns of estimates across methods and sample sizes for real traits subject to assortative mating are congruent with expected assortative mating-induced biases. For example, marker-based heritability estimates for height are 14% - 23% higher than corrected estimates using UK Biobank data.

Fig. 1. Theoretical and empirical behavior of existing and corrected estimators.

HE regression (${\widehat{h}}_{\mathrm{HE}}^2$) and REML (${\widehat{h}}_{\mathrm{REML}}^2$) estimates for varying phenotypic mate correlations ($r$) in simulated datasets ($n$ = 64,000) assuming that $h_0^2=.5$, that all causal variants are measured, and that AM has reached equilibrium. Values of ${\widehat{h}}_{\mathrm{HE}}^2$ are consistent with our closed-form approximation ($\mathbb{E}\left[{\widehat{h}}_{\mathrm{HE}}^2\right]$) under the assumption of exchangeable loci. Our corrected estimators, ${\widehat{h}}_{\mathrm{\infty }}^2$and ${\widehat{h}}_0^2$, which are based on ${\widehat{h}}_{\mathrm{HE}}^2$ and $r$ and assume equilibrium, recover the true equilibrium ($h_{\infty }^2$) and panmictic ($h_0^2$) heritabilities. Further, for a given observation of ${\widehat{h}}_{\mathrm{HE}}^2$ for a trait under disequilibrium AM, the current generation heritability ($h_t^2$) and $h_0^2$ are probabilistically bounded in expectation (teal and red regions respectively; see main text).

Fig. 2. REML and HE estimates across varying sample sizes in simulated data.

a Comparison of HE regression and REML heritability estimates as functions of sample size for varying phenotypic mating correlation ($r_{\mathrm{pheno}}$) and fixed panmictic heritability ($h_0^2=0.5$) in simulated data. We computed multiple estimates per sample size for each estimator and parameter combination by applying estimators to independent sub-samples. Whereas HE regression estimates are upwardly biased independent of sample size, REML estimates slowly converge to the panmictic heritability as sample sizes increase. b Extended simulations demonstrating high-dimensional behavior of the REML estimator as a function of sample size. Forward time simulations required a larger population size ($N_{\mathrm{sim}}=3\times {10}^6$) to obtain samples of up to $n=648,000$ unrelated individuals. Obtaining REML estimates for samples larger than this was not computationally feasible, but the dashed red line shows predicted values for larger sample sizes extrapolated from a regression model including first and second order log-linear components. Results are consistent with theoretical predictions that the REML estimator converges to the panmictic heritability in very large samples.

Fig. 3. Naïve approaches to addressing AM induced bias and the impact of missing data.

a Simulations employing the same parameters described in Fig. 1 demonstrate that neither partitioned nor principal component adjusted approaches mitigate the impact of AM on HE (${\widehat{h}}_{\mathrm{HE}}^2$) or REML (${\widehat{h}}_{\mathrm{REML}}^2$) heritability estimates. Additionally, simulations confirm that LDSC is subject to equivalent biases. Single component: standard single genomic variance component models. Single comp. + 10 PCs: included the first ten within-sample PCs as covariates. Partitioned: included four annotation-based variance components generated by median splits of within-sample minor allele frequencies and LD scores. b Simulations demonstrate that conclusions regarding estimator bias do not change when some of the influence of causal variants is not captured by measured SNPs. Simulations employed the same parameters as above except that 0, 50, or 75% of randomly selected SNPs (both causal and non-causal) were dropped. As expected, estimates were attenuated when SNPs were dropped but overall patterns remained consistent.

Fig. 4. REML and HE estimates across varying sample sizes in UK Biobank data.

Comparison of HE (${\widehat{h}}_{\mathrm{HE}}^2$) and REML (${\widehat{h}}_{\mathrm{REML}}^2$) heritability estimates as a function of sample size for real traits in a sample of unrelated European ancestry UK Biobank participants. Points connected by thin lines represent estimates derived from pairs of non-overlapping subsamples of size 16,000 and $n$−16,000, whereas thick lines reflect average log-linear trends. Two traits with previous evidence for AM (height and years of education) and two negative control traits (body mass index and bone mineral density) were selected for analysis a priori. Consistent with theoretical predictions, height and years of education demonstrated significant estimator divergence with increasing sample size whereas body mass index and bone mineral density did not.