Bayesian variable selection regression for genome-wide association studies and other large-scale problems

Yongtao Guan; Matthew Stephens

Source-linked AI summary

Bayesian variable selection regression for genome-wide association studies and other large-scale problems

Yongtao Guan, Matthew Stephens

arXiv:1110.6019v1Oct 27, 2011stat.AP

TL;DR

GWAS requires variable selection among hundreds of thousands or millions of SNPs, often with many relevant effects too small for complete confident identification. The paper applies BVSR with interpretable posterior inclusion probabilities and PVE-focused priors, showing that it can be applied to large problems and estimate overall signal despite incomplete variant identification.

Problem
GWAS seeks relevant genetic variants in extremely high-dimensional data, but many small effects can make the complete relevant subset impossible to identify confidently.
Method
The paper applies BVSR with MCMC inference, priors on sparsity and effect size, and a prior specification based on the implied proportion of variance explained.
Results
BVSR can be successfully applied to large problems and can estimate overall signal even when confidently identifying relevant covariates is not possible.
Takeaways & Limitations
BVSR provides data-driven hyperparameter estimates, interpretable covariate inclusion probabilities, and PVE estimates relevant to assessing measured genetic signal and missing heritability.
Takeaways & Limitations
The simple computational scheme can perform poorly for multiple correlated covariates that are far apart along a chromosome and require joint updates.

Abstract

from arXiv · show

We consider applying Bayesian Variable Selection Regression, or BVSR, to genome-wide association studies and similar large-scale regression problems. Currently, typical genome-wide association studies measure hundreds of thousands, or millions, of genetic variants (SNPs), in thousands or tens of thousands of individuals, and attempt to identify regions harboring SNPs that affect some phenotype or outcome of interest. This goal can naturally be cast as a variable selection regression problem, with the SNPs as the covariates in the regression. Characteristic features of genome-wide association studies include the following: (i) a focus primarily on identifying relevant variables, rather than on prediction; and (ii) many relevant covariates may have tiny effects, making it effectively impossible to confidently identify the complete "correct" subset of variables. Taken together, these factors put a premium on having interpretable measures of confidence for individual covariates being included in the model, which we argue is a strength of BVSR compared with alternatives such as penalized regression methods. Here we focus primarily on analysis of quantitative phenotypes, and on appropriate prior specification for BVSR in this setting, emphasizing the idea of considering what the priors imply about the total proportion of variance in outcome explained by relevant covariates. We also emphasize the potential for BVSR to estimate this proportion of variance explained, and hence shed light on the issue of "missing heritability" in genome-wide association studies.

1. Introduction.

The paper frames GWAS as a large-scale variable-selection problem and evaluates BVSR as a practical approach focused on identifying relevant variants. It emphasizes interpretable posterior inclusion probabilities, prior specification through explained variance, and PVE estimation when complete variable identification is infeasible.

GWAS measures hundreds of thousands or millions of SNPs in hundreds to tens of thousands of individuals to identify regions affecting phenotypes.
BVSR treats phenotypes as regression responses and SNPs as covariates, identifying genomic regions through SNPs with likely nonzero coefficients.
BVSR produces posterior probabilities that individual covariates have nonzero coefficients, offering interpretable confidence measures for GWAS identification goals.
The paper demonstrates that BVSR can be practical for problems with hundreds of thousands of covariates and thousands of observations.
The proposed prior-specification approach focuses on the total proportion of response variance explained by relevant covariates, or PVE.
BVSR may estimate PVE even when small effects prevent reliable identification of all relevant covariates.

2. Models and priors.

The model represents sparse associations with binary inclusion indicators and Bayesian priors on key hyperparameters. Its novel prior specification targets interpretable beliefs about the proportion of variance explained rather than coefficient scale alone.

The model uses standard normal linear regression with observed response y and covariate matrix X, while inferring µ, β, and residual precision τ.
Because GWAS often has p ≫ n with few associated covariates, binary indicators γ identify which regression coefficients are nonzero.
In GWAS, γj = 1 can usually indicate a causal or correlated functional SNP, assuming population structure has been corrected.
The hyperparameters π and σa represent model sparsity and typical nonzero-effect size, and are assigned priors so their values can be informed by data.
The prior on π is chosen so its limits correspond to expectations of 1 and M covariates, with M = 400 used in applications partly for computational reasons.
The exchangeable priors treat all covariates as equally plausible and ignore functional clustering, although the framework could incorporate such information in future work.
2.1. Novel prior on σ2: The proposed prior makes the induced prior on PVE relatively flat over (0,1) by placing a uniform prior on h independently of γ.

3. Computation and inference.

The computational framework uses MCMC over model and hyperparameter space, with analytical integration and targeted estimators for inclusion probabilities, PVE, and prediction. The authors note that the simple scheme can perform poorly for certain correlated-covariate configurations.

MCMC samples the posterior of (h, π, γ), while β and τ are integrated out for marginal likelihood calculations and then sampled conditionally.
The Metropolis–Hastings algorithm explores models by adding, removing, or switching individual covariates and uses proposal strategies to improve computational performance.
The simple computational scheme may perform poorly when multiple correlated covariates are far apart along a chromosome and require joint updates.
Rao–Blackwellization estimates each covariate’s posterior inclusion probability by averaging conditional probabilities across MCMC samples.
Posterior PVE samples are obtained by computing PVE(β(i),τ(i)) for each sampled pair from the MCMC procedure.
Prediction for a future individual uses observed covariates together with posterior information about regression coefficients.
The relative prediction gain normalizes prediction error against mean prediction and optimal prediction, indicating the proportion of extractable signal obtained.

4. Goals and expectations.

The paper aims to extract more information from GWAS signals than conventional single-SNP analyses, while avoiding claims about a single best model or nonadditive SNP combinations. BVSR also accounts for local SNP correlations and controls for identified SNPs when evaluating subsequent ones.

BVSR seeks to extract more information from marginal associations than standard single-SNP analyses.
Rao–Blackwellized posterior inclusion probabilities effectively evaluate each SNP while controlling for SNPs currently included in the model.
The method does not seek a single best model because small effects and p ≫n make any particular model unlikely and the correct model effectively impossible to identify.
BVSR does not target nonadditive SNP interactions whose effects appear only through combinations of variants.
Local SNP correlations create a highly banded X′X structure, a feature that may make BVSR more useful in GWAS than in other settings.

5. Simulations and comparisons with other methods.

The simulations mimic GWAS settings with many more covariates than individuals, sparse causal architecture, and mostly small effects. They use simulated and real genotype data at scales up to approximately 550K SNPs and compare BVSR with single-SNP analysis and LASSO.

The simulations reproduce GWAS characteristics with p ≫n, approximately 10,000–300,000 covariates, about 1,000 individuals, and extreme sparsity.Most simulations use approximately 30 covariates affecting the response.
Most relevant covariates individually explain less than 1% of the variance of y.
The study uses 10K simulated independent SNPs, approximately 317K real SNPs, and approximately 550K real SNPs to assess realistic scale and scalability.The 317K data contain 980 individuals; the 550K data contain 988 individuals after quality control.
Phenotypes are simulated by selecting 30 causal SNPs, drawing their effects from a specified distribution, and choosing residual variance to achieve a target PVE.
BVSR uses one million iterations for 10K data and two million iterations for 317K and 550K data, with runtimes ranging from minutes to about one day.
Comparisons include single-SNP association analysis and the penalized regression method LASSO.Single-SNP rankings use Bayes factors, which perform similarly to standard single-SNP p-value rankings.

5.3. Inference of PVE, and its relationship to heritability.

BVSR estimates PVE by Bayesian model averaging over uncertainty in which covariates are relevant, connecting this quantity to questions about heritability. Simulations show generally good calibration in ordinary settings, but substantially less precision when many tiny effects produce the same total PVE.

PVE is closely related to heritability but reflects the optimal predictive accuracy from measured common genetic variants rather than the full genetic contribution.
Bayesian model averaging lets BVSR estimate PVE while accounting for uncertainty about which covariates are relevant.This differs from single-SNP and standard penalized regression approaches that typically identify one covariate set.
Results from 50 phenotype simulations per genotype setting are qualitatively similar under normal and double-exponential effect-size distributions.The displayed results use the double-exponential distribution DE(1).
85% of 90% credible intervals for PVE covered the true value across simulated and real SNP data.Posterior uncertainty increased with more SNPs and smaller true PVE, although intervals still excluded large values when true PVE was small.
For 1,000 causal SNPs among 317K, PVE estimates were considerably less precise, though credible intervals differed systematically between PVE = 0.3 and PVE = 0.For PVE = 0.3, 6 intervals covered 0.3; for PVE = 0, 1 interval covered 0.3.
For 300 causal SNPs among 10K with PVE = 0.3, 7 of 10 simulations had 90% credible intervals including PVE = 0.3.The interquartile range spanned the true PVE in 5 of 10 simulations.

5.5. Identifying the causal SNPs.

Across simulations, multi-SNP methods improved causal-variant identification over single-SNP analysis at a given false-positive level, while BVSR and LASSO were similar for independent SNPs and BVSR was better for correlated genotypes. BVSR also showed higher predictive gain than optimal LASSO in most data sets, whereas LASSO refitting was often counterproductive.

Causal SNP identification: Multi-SNP methods always achieved at least as many true positives as single-SNP analysis for a given number of false positives.This held across the simulated comparisons of causal SNP identification.
Causal SNP identification: BVSR and LASSO performed similarly with independent SNPs, but BVSR performed better with real correlated genotypes.The authors note that the reason for this difference was unclear.
Causal SNP identification: Multi-SNP analyses gained measurable power even for independent SNPs, likely because finite samples create nonzero correlations among nominally independent SNPs.Controlling for identified causal variants can improve detection of subsequent variants, although the gain is expected to be small for nearly orthogonal SNPs.
Causal SNP identification: Most multi-SNP gains occurred when false positives were small but nontrivial, by elevating moderately difficult causal SNPs in the rankings.The authors suggest combining multi-SNP analysis with other data or analyses to distinguish true and false positives near the top of the rankings.
Prediction performance: BVSR’s RPG exceeded directly optimized LASSO without refitting in 82 of 100 data sets, with mean RPG 0.315 versus 0.261.The comparison used the optimal LASSO solution, effectively granting LASSO oracle knowledge of its tuning parameter.
Prediction performance: LASSO refitting often worsened prediction: RPG2 was frequently negative, indicating overfitting when the initial model contained false positives.Refitting helped when LASSO selected relevant covariates precisely but was counterproductive when selection included several irrelevant covariates.
Prediction performance: Sparse BVSR predictors using the top 10, 30, or 100 PIPs had average RPGs of 0.313, 0.315, and 0.315, essentially matching the nonsparse posterior mean.The results also suggest that Bayesian model averaging accounts for a moderate portion of BVSR’s predictive gain over RPG1.
Prediction performance: The authors identify comparisons with more flexible penalized regressions, such as Elastic Net or SCAD, as future work.These alternatives may separate shrinkage and sparsity more flexibly or avoid overshrinking large effects.

5.7. Calibration of the posterior inclusion probabilities.

BVSR produces posterior inclusion probabilities whose calibration depends on estimating prior hyperparameters from the data. Fixing π or σa to substantially incorrect values can make the resulting PIPs poorly calibrated.

Calibration: BVSR’s PIPs were reasonably calibrated under the priors described in the paper, with binned proportions of causal SNPs lying near y = x.The calibration plot bins PIPs into 20 intervals and compares mean PIP with the observed causal-SNP proportion.
Calibration: Fixing π or σa at values approximately five times too large or too small produced poorly calibrated PIPs.The other hyperparameter was held at its true value in each sensitivity analysis.
Application: BVSR was applied to a GWAS of plasma C-reactive protein concentration in the PARC study.CRP is associated with inflammation and predictive of future cardiovascular disease.

5.8. Real data analysis: PARC GWAS for C-reactive protein.

BVSR applied to PARC CRP data identifies genomic regions associated with CRP and quantifies uncertainty at both region and SNP levels. Its PVE estimate suggests measured SNPs explain more signal than stringent single-SNP significance alone, but not all missing heritability.

Data and computation: 1,682 individuals provided phenotype and genotype data for four independent MCMC runs, with three showing strong agreement and one showing mild convergence or mixing concerns.The runs used 2 or 4 million iterations; the longer runs required approximately 60 and 90 CPU hours.
Interpretation: Correlated SNPs can receive distributed PIPs, making regional summaries more appropriate than simply ranking individual SNPs.This reflects uncertainty about which correlated SNP is included even when evidence supports including at least one variant in the region.
Regional analysis: Overlapping 1 Mb regions summarized BVSR evidence using posterior expected SNP counts and probabilities of containing one, two, or more included SNPs.Adjacent regions overlapped by 0.5 Mb, and expected counts were computed by summing SNP-level PIPs.
Regional analysis: Two regions showed very strong association evidence with E > 0.95, while a third showed moderately strong evidence with E > 0.75.Multiple additional regions had modest signals of E = 0.1 to 0.5, most of which were unlikely to be truly associated at that signal level.
Biological findings: The strongest association regions contained CRP, HNF1A, and APOE/APOC, concordant with robust associations reported in other single-SNP GWAS.The largest single-SNP Bayes factors in these regions were 10^6.2, 10^5.5, and 10^4.9, respectively.
PVE and missing heritability: BVSR estimated that measured SNPs explained a posterior mean PVE of 0.14, with a symmetric 90% CI of [0.05,0.25].The lower interval limit was similar to the approximately 6% PVE estimated from significant SNPs, while the upper limit remained below previous heritability estimates.

6. Extension to binary phenotypes.

The paper extends BVSR to binary phenotypes through a probit model with latent Gaussian variables. Simulations indicate generally reliable latent-variable PVE estimation and results for analyses limited to thousands of SNPs.

Probit extension: A probit extension links binary outcomes to latent variables z through yi = 1 when zi > 0.Posterior inference integrates out z with MCMC and requires one additional update relative to quantitative-phenotype analysis.
Probit extension: The binary-trait extension reuses quantitative-outcome priors by applying them to the unobserved latent Gaussian variables.Overall signal is summarized by estimating the PVE of the latent variables.
Simulation design: Simulations used 10,000 independent SNPs and samples of n = 1,000 or n = 6,000 individuals, with binary outcomes formed by thresholding simulated latent phenotypes.Thirty causal SNPs generated latent phenotypes with prespecified PVE before the largest n/2 values were assigned y = 1.
Simulation results: The method generally produced reliable binary-trait results when analyses were limited to thousands of SNPs, including agreement across independent runs.The authors suggest this supports analyses of small genomic regions such as individual genes.

7. Discussion.

The discussion presents BVSR as a scalable, interpretable alternative for GWAS that estimates hyperparameters and PVE, while identifying effect-size, data-quality, prior-independence, and computational limitations.

Contributions: BVSR was successfully applied to large problems and can estimate hyperparameters and overall PVE that relevant covariates may explain.The discussion frames these as potential benefits compared with standard single-SNP analyses.
Contributions: The paper introduces an interpretable prior-specification approach and presents BVSR as a competitive alternative to LASSO.The prior perspective focuses on what parameter choices imply about total variance explained.
Limitations: Genotyping error can cause false-positive associations when error correlates with phenotype, particularly in case-control studies with differing DNA quality.This makes data-quality control an important practical concern for GWAS analysis.
Limitations: Normally distributed effect sizes are an important assumption, and phenotypes with unusually large regional effects may violate it.The authors suggest checking robustness by analyzing data both including and excluding known large-effect loci.
Future extensions: Independent, common inclusion priors may be overly skeptical of multiple functional variants within one gene because they ignore local spatial dependence.The authors propose incorporating genomic annotations or dependence among nearby variables in future priors.
Limitations: BVSR remains considerably more computationally demanding than penalized regression methods, motivating more efficient inference algorithms.The authors view improved scalability as an important direction for large applications.

APPENDIX A: DETAILS OF MCMC SCHEME

The MCMC scheme samples model inclusion, prior inclusion probability, and variance-explained parameters using Metropolis–Hastings proposals. It combines local add/remove/exchange moves with rank-based additions and occasional longer-range proposals.

Posterior sampling: MCMC samples (h, π, γ) from their posterior distribution over continuous hyperparameters and binary covariate-inclusion configurations.The state space is (0,1)×(0,1)×{0,1}^p.
Marginal likelihood: The method integrates out β and τ analytically to compute the marginal likelihood p(y|h,γ).This analytic marginalization is used in posterior sampling of the model and hyperparameters.
Conditional sampling: After sampling h and γ, the algorithm samples β and τ from their conditional posterior distributions given y, γ, and h.The regression coefficients and residual precision are handled conditionally after the model state is sampled.
Metropolis–Hastings updates: The Metropolis–Hastings proposal updates γ locally, proposes π from Beta(|γ′|, p−|γ′|+1), and then proposes a new h.The proposal jointly modifies model composition and associated hyperparameters.
Long-range proposals: Longer-range proposals compound a random number of local moves, with probability 0.3 per iteration, to improve the theoretical convergence rate.The number of compounded moves is uniform on (2,...,20), and the method is called a small-world proposal.
Rank-based proposals: Covariate additions use ranks based on single-covariate Bayes factors so proposals focus more on covariates likely to enter the model.The ranking compares each covariate-only model with the null model at σ_a = 1.
Model moves: The γ update uses add, remove, and exchange moves to modify the currently included covariate set.Each move begins from the current inclusion vector and changes one covariate or swaps covariates in and out.
Move probabilities: The implementation selects add, remove, and exchange moves with probabilities 0.45, 0.45, and 0.1, respectively.These probabilities govern the move type chosen at each update.

APPENDIX B: CALCULATIONS FOR RAO–BLACKWELLIZED ESTIMATES

The appendix derives calculations for terms in equation (3.2), including residual-based univariate regressions and posterior quantities for covariates. It also explains parameter dependence and computational reuse across covariates.

Posterior calculations: The appendix separately specifies calculations for the numerator and denominator of the first term.These steps provide the components needed to compute the equation (3.2) terms.
Prior-dependent terms: The second term accounts for dependence between β−j and γj because the prior variance σa depends on h and γ.The associated calculation uses β−j|γ,τ,h as i.i.d. N(0,σ2(h,γ)/τ).
Posterior calculations: The first-term numerator uses a design matrix X = (1, Xj) with an intercept and covariate j, followed by the posterior distribution of βj.The matrix is n × 2, and Ω is defined from XtX + ν−1.
Residual-based calculations: The calculations reduce to univariate regressions of residuals R against covariate j.Residuals are defined as R = y − Xγ−jβγ−j.
Residual-based calculations: For covariates j ∉ γ+, the same residuals can be reused, whereas residuals must be recomputed for j ∈ γ+.This distinguishes the computational treatment of included and excluded covariates.

Acronyms used in the paper.

The paper defines acronyms for its Bayesian modeling, association-study, variable-selection, sampling, and variance-explanation terminology.

Acronyms: The acronyms include BMA, BVSR, GWAS, LASSO, MCMC, PIP, PVE, RPG, SIS, and SNP.They expand to Bayesian model averaging; Bayesian variable selection regression; genome-wide association studies; least absolute shrinkage and selection operator; Markov chain Monte Carlo; posterior inclusion probability; proportion of variance explained; relative prediction gain; sure independence screen; and single nucleotide polymorphism.