Synthetic surrogates improve power for genome-wide association studies of partially missing phenotypes in population biobanks

Abstract

Within population biobanks, incomplete measurement of certain traits limits the power for genetic discovery. Machine learning is increasingly used to impute the missing values from the available data. However, performing genome-wide association studies (GWAS) on imputed traits can introduce spurious associations, identifying genetic variants that are not associated with the original trait. Here we introduce a new method, synthetic surrogate (SynSurr) analysis, which makes GWAS on imputed phenotypes robust to imputation errors. Rather than replacing missing values, SynSurr jointly analyzes the original and imputed traits. We show that SynSurr estimates the same genetic effect as standard GWAS and improves power in proportion to the quality of the imputations. SynSurr requires a commonly made missing-at-random assumption but relaxes the requirements of existing imputation methods by not requiring correct model specification. We present extensive simulations and ablation analyses to validate SynSurr and apply it to empower the GWAS of dual-energy X-ray absorptiometry traits within the UK Biobank.

Extended Data Figure 1:. Robustness and precision of SynSurr with an uninformative and informative synthetic surrogate.

In all cases, the number of subjects with observed phenotypes was $n={10}^3$. The number of subjects with missing phenotypes was varied to achieve the indicated level of missingness. The standard estimator utilizes the observed values of $Y$ only. In panel A, the synthetic surrogate has correlation $\rho =0.00$ with the target phenotype, and is in fact independent of the target phenotype. Use of the SynSurr estimator with this uninformative surrogate results in no loss of efficiency relative to the standard analysis. In panel B, the synthetic surrogate has correlation $\rho =0.75$ with the target phenotype. SynSurr becomes more efficient as the number of subjects with missing target outcomes increases. The center of the box plot is the median, the upper and lower bounds of the box are the 75th and 25th percentiles, and the whiskers extend from the minimum to the maximum. The number of simulation replicates is 5 × 10³.

Extended Data Figure 2:. Signal recovery of SynSurr relative to the oracle GWAS for height and FEV1.

A slope of 1.0 indicates that the estimated effect sizes are consistent with the oracle effect sizes. Note that although the slope deviates from 1.0 at 90% missigness, the slope approaches 1.0 as missingness declines. The following figure, which assesses signal recovery for standard GWAS, provides a point of comparison for the $R^2$ values.

Extended Data Figure 3:. Signal recovery of imputation-based approaches and SynSurr relative to the oracle GWAS for height with 50% missingness.

A slope of 1.0 indicates that the estimated effect sizes are consistent with the oracle effect sizes, whereas a slope deviating from 1.0 suggests the presence of bias.

Extended Data Figure 4:. Predicted vs. observed values of body composition phenotypes within the model-building and GWAS data sets.

A random forest was trained to predict each of the 6 body composition phenotypes, obtained via DEXA scan, using 4,584 subjects allocated to the model-building data set. The GWAS dataset consists of 29,577 unrelated subjects with body compositions measured via DEXA. Model inputs included age, sex, height, body weight, body mass index, and 5 impedance measures (whole body, left arm, right arm, left leg and right leg).

Extended Data Figure 5:. Distribution of predicted body masses comparing subjects with and without DEXA measurements.

The violin plot shows the kernel density estimation of the distribution of the data, with the tips of the violin indicating the maximum and minimum observed values among subjects. Sample sizes: $n=29,577$ independent subjects with DEXA measurements; $n=317,921$ subjects without DEXA measurements.

Extended Data Figure 6:. SynSurr remains unbiased and properly controls the type I error when the same data are utilized for model training and for GWAS.

The number of subjects with observed phenotypes was $n={10}^3$, while the number with missing phenotypes was varied to achieve the indicated level of missingness. The model that generated the synthetic surrogate was either trained in the GWAS data set or in an independent data set of size $n={10}^3$. Upper shows the distribution of effect sizes across 20 × 10³. The true genetic effect size is ${\beta }_G=0.1$. The center of the box plot is the median, the upper and lower bounds of the box are the 75th and 25th percentiles, and the whiskers extend from the 5th to the 95th percentile. Lower shows the average ${\chi }^2$ statistic under $H_0:{\beta }_G=0$ across 50 × 10³ simulation replicates, for which the expected value is 1.0. Error bars are 95% confidence intervals for the mean. Panel A (left) considers a “misspecified” ($k=2$) model that can only capture quadratic dependence of $Y$ on $X$, while Panel B (right) considers a “correctly specified” model ($k=3$) that can capture the cubic dependence. As seen, the validity of SynSurr is not contingent on correct specification of the surrogate model.

Extended Data Figure 7:. Survey of assumptions surrounding missing phenotypic data in GWAS.

The methods sections of all studies contributing summary statistics to the GWAS catalog between May 1st and November 1st, 2023, were manually reviewed. Among 47 studies, 24 did not address missing phenotypic data. Of the 23 remaining, 21 made an assumption of missing at random (MAR) or missing completely at random (MCAR).

Figure 1:. Graphical overview of a SynSurr GWAS.

A. The data set is first split into a fully labeled model-building data set, including the target phenotype and surrogates, and a partially labeled inference data set, which also includes genetics. B. Within the model-building data set, an imputation model is trained to predicted the target phenotype on the basis of surrogates. C. The imputation model is transferred to the partially labeled inference data set and applied to predict the target outcome for all subjects. The predicted value of the target outcome is referred to as the “synthetic surrogate”. Importantly, the synthetic surrogate is maintained as a separate and distinct outcome from the partially missing target phenotype. D. Finally, within the inference data set, the partially missing target phenotype and the fully observed synthetic surrogate are jointly regressed on genotype and covariates to identify genetic variants associated with the target outcome.

Figure 2:. Unlike imputation-based estimators, SynSurr is robust to misspecification of the imputation model.

The true value for the parameter of interest is ${\beta }_G=0.1$, corresponding to a variant with $h^2=1\%$. For each estimator, the sample size is $n={10}^4$, the mean value across 10³ simulations is shown by the point, and two 95% confidence intervals (CIs) are presented: the dotted CI is based on the analytical standard error (SE) while the solid CI is based on the empirical SE. The oracle estimator has access to the complete version of $Y$, before 25% of values were set to missing. The standard estimator has access to the observed values of $Y$ only. The imputation-based estimators impute the missing values of $Y$ using an imputation-model fit on an independent data set. The set of covariates used to fit the imputation model are shown as a tuple: the imputation model based on $G$ and $X$ is correctly specified, whereas that based on $G$ alone or $X$ alone is misspecified. The SynSurr estimator jointly analyzes the partially missing $Y$ with the synthetic surrogate $\widehat{Y}$, where $\widehat{Y}$ is generated for all subjects from the imputation model. The key observation is that SynSurr does not require a correctly specified generative model to yield unbiased estimation and valid inference.

Figure 3:. SynSurr controls type I error across missingness rates and target-surrogate correlations.

Type I error is the probability of incorrectly rejecting the null hypothesis $H_0:{\beta }_G=0$. In all cases, the number of subjects with observed phenotypes was $n={10}^3$. The number of subjects with missing phenotypes was varied to achieve the indicated level of missingness. The synthetic surrogate has correlation $\rho \in \{0.00,0.25,0.50,0.75\}$ with the target phenotype. The number of simulation replicates is 10⁶. P-values are two-sided and were calculated by SynSurr. Error bands (dashed black lines) represent 95% confidence intervals around the expected −log₁₀ (p-values) under the null hypothesis. Adherence to the diagonal (red line) indicates that the p-values are uniformly distributed under the null.

Figure 4:. Power of SynSurr across various missing rates, target-surrogate correlations, and SNP heritabilities.

Power is the probability of correctly rejecting the null hypothesis $H_0:{\beta }_G=0$. In all cases, the number of subjects with observed phenotypes was $n={10}^3$. The number of subjects with missing phenotypes was varied to achieve the indicated level of missingness. In each panel, the synthetic surrogate has correlation $\rho \in \{0.00,0.25,0.50,0.75\}$ with the target phenotype and the SNP heritability was varied from 0.1% to 1%. When there is no missingness, SynSurr is equivalent to the standard analysis and shows no variation across values of $\rho$. The power of SynSurr increases with increasing missingness and target-surrogate correlation. The number of simulation replicates is 10⁴.

Figure 5:. Comparing SynSurr and standard GWAS with respect to the number and significance of genome-wide significant associations for body composition traits.

A. Number of genome-wide significant (GWS) associations ($p<5\times {10}^{-8}$) with DEXA body composition traits for standard and synthetic surrogate (SynSurr) GWAS. P-values are two-sided and are calculated by linear regression (Standard) or SynSurr. B. Average ${\chi }^2$ statistic at the union of variants that reached genome-wide significance under either method. A greater expected ${\chi }^2$ statistic directly corresponds to greater power to detect an association. Error bars are 95% confidence intervals for the mean. The number of independent GWS variants averaged across is shown in A.

Figure 6:. External validation via overlap of genome-wide significant variants for body composition with associations from the GWAS catalog.

Variants from the GWAS catalog associated with body fat distribution, body fat percentage, fat body mass, and lean body mass were compiled. A study variant was considered overlapped if it fell within 250 kb of a GWAS catalog variant. Panels A and B show the counts and proportions of overlapped variants, respectively. Note that, with 1 exception, all variants identified by standard GWAS were also identified by SynSurr (Supplementary Table 21). The perfect overlap of the standard GWAS variants with known body composition associations in panel B is a direct consequence of the standard GWAS detecting very few genome-wide significant variants (8.3 on average), and indicates that all of these variants were previously known.