Imputation-Based
Analysis of Association Studies:
Candidate Regions and Quantitative Traits
Bertrand Servin¤a*, Matthew Stephens¤b
Department of Statistics, University of Washington, Seattle, Washington, United States of America

We introduce a new framework for the analysis of association studies, designed to allow untyped variants to be more
effectively and directly tested for association with a phenotype. The idea is to combine knowledge on patterns ofcorrelation among SNPs (e.g., from the International HapMap project or resequencing data in a candidate region ofinterest) with genotype data at tag SNPs collected on a phenotyped study sample, to estimate (‘‘impute’’) unmeasured
genotypes, and then assess association between the phenotype and these estimated genotypes. Compared withstandard single-SNP tests, this approach results in increased power to detect association, even in cases in which thecausal variant is typed, with the greatest gain occurring when multiple causal variants are present. It also provides
more interpretable explanations for observed associations, including assessing, for each SNP, the strength of theevidence that it (rather than another correlated SNP) is causal. Although we focus on association studies withquantitative phenotype and a relatively restricted region (e.g., a candidate gene), the framework is applicable andcomputationally practical for whole genome association studies. Methods described here are implemented in asoftware package, Bim-Bam, available from the Stephens Lab website http://stephenslab.uchicago.edu/software.html.
Citation: Servin B, Stephens M (2007) Imputation-based analysis of association studies: Candidate regions and quantitative traits. PLoS Genet 3(7): e114. doi:10.1371/journal.
pgen.0030114

and a Bayesian regression approach to perform the test for
association, allowing for potential errors in the imputed
genotypes. Although we focus speciﬁcally on methods foranalyzing quantitative phenotypes in candidate gene studies,
the same general framework can also be applied to discrete
traits, and/or genome-wide scans.
These imputation-based methods can be viewed as anatural analysis complement to the ‘‘tag SNP’’ design strategy
for association studies, which attempts to choose SNPs thatare highly correlated with, and hence good predictors of,untyped SNPs. We are simply directly exploiting thisproperty, together with recently developed statistical models
for multi-locus LD ([4,5]) to infer the untyped SNP genotypes.
Our approach is also somewhat analogous to multipoint
approaches to linkage mapping (e.g., [6]), in which observed

Introduction
Although the development of cheap high-throughput
genotyping assays have made large-scale association studies
a reality, most ongoing association studies genotype only asmall proportion of SNPs in the region of study (be that thewhole genome, or a set of candidate regions). Because ofcorrelation (linkage disequilibrium, LD) among nearby
markers, many untyped SNPs in a region will be highly
correlated with one or more nearby typed SNPs. Thus,
intuitively, testing typed SNPs for association with aphenotype will also have some power to pick up associations
between the phenotype and untyped SNPs. In practice,typical analyses involve testing each typed SNP individually,
and in some cases combinations of typed SNPs jointly (e.g.,
haplotypes), for association with phenotype, and hoping thatthese tests will indirectly pick up associations due to untyped
SNPs. Here, we present a framework for more directly andeffectively interrogating untyped variation.
In outline, our approach improves on standard analyses by
exploiting available information on LD among untyped andtyped SNPs. Partial information on this is generally available
from the International HapMap project [1]; in some cases
more detailed information (e.g., resequencing data) may also
be available, either through public databases (e.g., SeattleSNPs [2]), or through data collected as a part of theassociation study design (e.g., [3]). Our approach combines
this background knowledge of LD with genotypes collected at
typed SNPs in the association study, to explicitly predict(‘‘impute’’) genotypes in the study sample at untyped SNPs,
and then tests for association between imputed genotypes
and phenotype. We use statistical models for multi-marker
LD to perform the genotype imputation, with uncertainty,
PLoS Genetics | www.plosgenetics.org

Editor: David B. Allison, University of Alabama at Birmingham, United States ofAmerica
Received December 14, 2006; Accepted May 30, 2007; Published July 27, 2007
A previous version of this article appeared as an Early Online Release on May 30,
2007 (doi:10.1371/journal.pgen.0030114.eor).
Copyright: Ó 2007 Servin and Stephens. This is an open-access article distributed
under the terms of the Creative Commons Attribution License, which permits
unrestricted use, distribution, and reproduction in any medium, provided theoriginal author and source are credited.
Abbreviations: BF, Bayes factor; df, degree of freedom; IBD, identity by descent;LD, linkage disequilibrium; MCMC, Markov Chain Monte Carlo; QTN, quantitative
trait nucleotide* To whom correspondence should be addressed. E-mail: bservin@toulouse.inra.fr¤a Current address: INRA Laboratoire de Génétique Cellulaire, Castanet-Tolosan,
France¤b Current address: Departments of Statistics and Human Genetics, University of
Chicago. Chicago, Illinois, United States of America

1296

July 2007 | Volume 3 | Issue 7 | e114

Imputation-Based Analysis of Association Studies

Author Summary

random sample from the population, and consider application to other designs in the discussion.
Our strategy is to use patterns of LD in the panel, together
with the tag SNP genotypes in the cohort, to explicitly predictthe genotypes at all markers for members of the cohort, andthen analyze the data as if the cohort had been genotyped at all
markers. There are thus two components to our approach: (i)predicting (‘‘imputing’’) cohort genotypes, and (ii) analyzing
association between cohort genotypes and phenotypes. For(i), we use existing models for population genetic variation
across multiple markers [4,5], which perform well at estimating missing genotypes, and provide a quantitative assessmentof the uncertainty in these estimates [5]. For (ii), we introduce
a new approach based on Bayesian regression, and describe
how this approach can yield not only standard Bayesianinference, but also p-values for testing the null hypothesis ofno genotype–phenotype association. We chose to take a
Bayesian approach partly because it provides a natural way toconsider uncertainty in estimated genotypes. However, theBayesian approach has other advantages; in particular, it
provides a measure of the strength of the evidence for an
association (the Bayes factor, BF) that is, in some respects,
superior to conventional p-values. Furthermore, in our
simulations, p-values from our Bayesian approach provide
more powerful tests than standard tests, even if the cohort isactually genotyped at all markers (including all causal
variants).

Ongoing association studies are evaluating the influence of genetic
variation on phenotypes of interest (hereditary traits and susceptibility to disease) in large patient samples. However, althoughgenotyping is relatively cheap, most association studies genotype
only a small proportion of SNPs in the region of study, with many
SNPs remaining untyped. Here, we present methods for assessing
whether these untyped SNPs are associated with the phenotype ofinterest. The methods exploit information on patterns of multimarker correlation (‘‘linkage disequilibrium’’) from publically available databases, such as the International HapMap project or the
SeattleSNPs resequencing studies, to estimate (‘‘impute’’) patient
genotypes at untyped SNPs, and assess the estimated genotypes for
association with phenotype. We show that, particularly for common
causal variants, these methods are highly effective. Compared withstandard methods, they provide both greater power to detectassociations between genetic variation and phenotypes, and also
better explanations of detected associations, in many cases closely
approximating results that would have been obtained by genotyping all SNPs.

genotypes at multiple markers predict patterns of identity bydescent (IBD) at nearby positions without markers, and test
for correlation between these patterns of IBD and observed
phenotypes. In the association context, we are predicting
identity by state rather than IBD, and the methods ofpredicting identity by state versus IBD differ greatly, but
the approaches share the idea of using multipoint information to predict single-point information, and, at least in theirsimplest form, subsequently assessing correlation with phenotype at the single-point level. This strategy provides a clean
and rigorous way to avoid the ‘‘curse of dimensionality’’ thatcan plague haplotype-based analyses, without making ad hoc
decisions such as pooling rare haplotypes into a single class.
Although our methods are developed in a Bayesian framework, they can also be used to compute p-values assessing
signiﬁcance of observed genotype–phenotype associations.
Our approach should therefore be of interest to practitioners
whether or not they favor Bayesian procedures in general. It
has two main advantages over more standard approaches.
First, it provides greater power to detect associations. Part ofthis increased power comes from incorporating extra
information (knowledge on patterns of LD among typed
and untyped SNPs), but, unexpectedly, we also found anincreased power of our Bayesian approach even when all
SNPs were actually typed. Second, and perhaps more
importantly, it provides more interpretable explanations forpotential associations. Speciﬁcally, for each SNP (typed anduntyped), it provides a probability that it is causal. This
contrasts with standard single-SNP tests, which provide a pvalue for each SNP, but no clear way to decide which SNPs
with small p-values might be causal.

Bayesian Regression Approach
We now provide further details of our Bayesian regression
approach. The literature on Bayesian regression methods is
too large to review here, but papers particularly relevant to
our work include [7–9].
For simplicity, we focus on the situation where cohort
genotypes are known at all SNPs (tag and non-tag). Extension
to the situation, where the cohort is genotyped only at tag
SNPs and other genotypes are imputed using sampling-basedalgorithms such as PHASE [10,11] or fastPHASE [5], isrelatively straightforward (see below).
Let G denote the cohort genotypes for all n individuals in
the cohort, and y ¼ (y1, . . ., yn) denote the corresponding(univariate, quantitative) phenotypes. We model the phenotypes by a standard linear regression:Xxij bj þ ei ;ð1Þ
yi ¼ l þjwhere yi is the phenotype measurement for individual i, l is
the phenotype mean of individuals carrying the ‘‘reference’’
genotype, the xijs are the elements of a design matrix X(which depends on the genotype data; see below), the bjs are
the corresponding regression coefﬁcients, and ei is a residual.
We assume eis are independent and identically distributed;N(0,1/s), where s denotes the inverse of the variance,
usually referred to as the precision (we choose this parameterization to simplifyX
notation in later derivations). Thus
xij bj ; 1=sÞ and:yi jl; xi ; b; s ; Nðl þ

Methods
We focus on an association study design in which genotype
data are available for a dense set of SNPs on a panel ofindividuals, and genotypes are available for a subset of these
SNPs (which for convenience we refer to as ‘‘tag SNPs’’) on acohort of individuals who have been phenotyped for a
univariate quantitative trait. We assume the cohort to be a
PLoS Genetics | www.plosgenetics.org

j

X
pﬃﬃﬃ
xij bj ÞÞ2 :Pðyi jxi ; l; b; sÞ} sexp½0:5sðyi  ðl þð2Þj

We assume a genetic model where the genetic effect is
additive across SNPs (i.e., no interactions) and where the1297

July 2007 | Volume 3 | Issue 7 | e114

Imputation-Based Analysis of Association Studies

priors are sensible. For prior D2, as we show in Protocol S1
the posteriors are proper, and the BF tends to a sensible limit.
For prior D1 we believe this to be true, although we have notproven it.
Prior on SNP effects. For brevity, we refer to SNPs that
affect phenotype as QTNs, for quantitative trait nucleotides.
Our prior on the SNP effects has two components: a prior onwhich SNPs are QTNs and a prior on the QTN effect sizes.
Prior on which SNPs are QTNs. We assume that with someprobability, p0, none of the SNPs is a QTN; that is, the ‘‘nullmodel’’ of no genotype–phenotype association holds. Otherwise, with probability (1  p0), we assume there are l QTNs,
where l has some distribution p(l) on f1, 2, . . ., nsg where ns
denotes the number of SNPs in the region. Given l, we assume
all subsets of l SNPs are equally likely. Both p0 and p(l) can becontext-dependent, and choice of appropriate values isdiscussed below.
Prior on effect sizes. If SNP j is a QTN, then its effect ismodeled by two parameters, aj and dj ¼ ajkj. The parameter aj
measures a deviation from the mean l and will depend on theunit of measurement of the phenotype. To reﬂect this, we
scale the prior on aj by thephenotypic standard deviationpﬃﬃﬃﬃﬃﬃﬃ
within each genotype class, 1=s. Speciﬁcally, our prior on aj
is Nð0; r2a =sÞ, where ra reﬂects the typical size of a QTN effectcompared with the phenotype standard deviation within eachgenotypic class. Choice of ra may be context-dependent, andis discussed below.
The parameter dj ¼ ajkj measures the dominance effect of a
QTN. If kj ¼ 0, then the QTN is additive: the heterozygotemean is exactly between the means of the two homozygotes. If
kj ¼ 1 (respectively, 1), allele 1 (respectively, 0) is dominant.
The case jkj j.1 corresponds to overdominance of allele 1 orallele 0. We investigate two different priors for thedominance effect:
Prior D1: We assume that kj is a priori independent of aj,
with kj ; Nð0; r2k Þ. We chose rk ¼ 0.5, which givesPðjkj j.1Þ ’ 0:05, reﬂecting a belief that overdominance isrelatively rare.
Prior D2: We assume that dj is a priori independent of aj,
with dj ; Nð0; r2d =sÞ, where we took rd ¼ 0.5ra. This prior on
dj induces a prior on kj in which kj is not independent of aj.
Prior D1 has the attractive property that the priorprobability of overdominance is independent of the QTN
additive effect aj. However, the posterior distributions of aj
and kj must be estimated via a computationally intensive
Markov Chain Monte Carlo (MCMC) scheme (see Protocol
S2). (An alternative, which we have not yet pursued, would be
to approximate BFs under prior D1 by numerical methods,
such as Laplace Approximation; e.g., [14]). Prior D2 is moreconvenient, as, when combined with the priors on l and s inEquation 4, posterior probabilities of interest can becomputed analytically (Protocol S1).
For both priors D1 and D2 we assume effect parameters fordifferent SNPs are, a priori, independent (given the otherparameters).
Choice of p0, p(l), and ra. The above priors include‘‘hyperparameters,’’ p0 and ra, and a distribution p(l) thatmust be speciﬁed. The hyperparameter p0 gives the priorprobability that the region contains no QTNs. While choiceof appropriate value is both subjective and context-speciﬁc,
for candidate regions we suggest p0 will typically fall in therange 102 to 0.5. If data on multiple regions are available,

three possible genotypes at each SNP (major allele homozygote, heterozygote, and minor allele homozygote) haveeffects 0, a þ ak and 2a, respectively [12]. We achieve this byincluding two columns in the design matrix for each SNP, one
column being the genotypes (coded as 0, 1, or 2 copies of theminor allele), and the other being indicators (0 or 1) for
whether the genotype is heterozygous. The effect of SNP j isthen determined by a pair of regression coefﬁcients (bj1, bj2),
which are, respectively, the SNP additive effect aj anddominance effect dj ¼ ajkj. While there are other ways tocode the correspondence between genotypes and the design
matrix, we chose this coding to aid specifying sensible priors(see below).
Priors for (b, l, s). Prior speciﬁcation is intrinsicallysubjective, and specifying priors that satisfy everyone isprobably a hopeless goal. Our aim is to specify ‘‘useful’’
priors, which avoid some potential pitfalls (discussed below),
facilitate computation, and have some appealing properties,
while leaving some room for context-speciﬁc subjective input.
In particular, we describe two priors below, which we refer toas prior D1 and D2, that were developed based on thefollowing considerations: (i) inference should not depend on
the units in which the phenotype is measured; (ii) even if thephenotype is affected by SNPs in this region, the majority of
SNPs will likely not be causal; (iii) for each causal variant
there should be some allowance for deviations from additiveeffects (i.e., dominant/recessive effects) without entirely
discarding additivity as a helpful parsimonious assumption;and (iv) computations should be sufﬁciently rapid to makeapplication to genome-wide studies practical (this lastconsideration refers to prior D2).
Priors on the phenotype mean and variance. The parameters l and s relate to the mean and variance of thephenotype, which depend on units of measurement. It seems
desirable that estimates (and, more generally, posterior
distributions) of these parameters scale appropriately withthe units of measurement, so, for example, multiplying all
phenotypes by 1,000 should also multiply estimates of l by
1,000. Motivated by this, for prior D1 we used Jeffreys’ priorfor these parameters:
Pðl; sÞ}1=s:ð3ÞThis prior is well known to have the desired scaling
properties in the simpler context where observed data are
assumed to be N(l, 1/s) [13], and we conjecture that our prior
D1 also possesses these desired scaling properties in the more
complex context considered here, although we have notproven this.
For prior D2 we used a slightly different prior, based onassuming a prior for (l, s) of the forms ; Cðj=2; k=2Þljs ; Nð0; r2l =sÞ:ð4ÞSpeciﬁcally, our prior D2 assumes the limiting form of thisprior as j,k ! 0 and r2l ! ‘. In Protocol S1 we show that theposterior distributions obtained using this limiting priorscale appropriately.
Both prior distributions above are ‘‘improper’’ (meaningthat the densities do not integrate to a ﬁnite value). Great
care is necessary before using improper priors, particularly
where one intends to compute BFs to compare models, as we
do here. However, we believe results obtained using these
PLoS Genetics | www.plosgenetics.org

1298

July 2007 | Volume 3 | Issue 7 | e114

Imputation-Based Analysis of Association Studies

BF ¼ PðyjG; H1 Þ=PðyjG; H0 Þ;

then it might be possible to estimate p0 from the data,
although we do not pursue this here. Instead, we mostly
sidestep the issue of specifying p0 by focusing on the BF(described below), which allows readers of an analysis to use
their own value for p0 when interpreting results.
In specifying the prior, p(l), for the number of QTNs, wesuggest concentrating most of the mass on models withrelatively few QTNs. Indeed, here we focus mainly on theextreme case in which p(l) is entirely concentrated on l ¼ 1:
that is, the ‘‘alternative’’ model is that the region contains asingle QTN. Although rather restrictive, this seems a goodstarting point in practice, particularly since our results showthat it can perform well even if multiple QTNs are present.
Nonetheless, there are advantages to considering models withmultiple QTNs, and so we also consider a prior where p(l)puts equal mass on l ¼ 1, 2, 3, or 4. This prior sufﬁces toillustrate the potential of our approach, although in practice
it would probably be preferable to place decreasing probabilities on larger numbers of QTNs (e.g., p(l ¼ 2) , p(l ¼ 1)). An
alternative would be to sidestep specifying p(l) by computing
BFs comparing, say, 4-QTN, 3-QTN, 2-QTN, and 1-QTN
models versus the ‘‘null’’ model. However, interpreting andacting on these BFs will inevitably correspond to implicit
assumptions about the relative prior plausibility of these
multi-QTN models.
Finally, speciﬁcation of the standard deviation of the effectsize, ra, involves subtle issues. Although it may seem tempting
to use ‘‘large’’ ra to reﬂect relative ‘‘ignorance’’ about effectsizes [15], we believe this is inadvisable. Although large rayields a ﬂat prior on effect sizes, this prior is far from
uninformative, in that it places almost all its mass on large
effect sizes. The result would essentially allow only zero
effects (i.e., the ‘‘null’’ model), or large effects (the ‘‘alternative’’ model). If in truth the causal SNPs have relatively
small effect, which is probably generally realistic, then (forrealistic sample sizes) the null model would be strongly
favored over the alternative, because the data would be moreconsistent with zero effects than with large effects. Choice ofra can thus strongly affect inference, particularly the BF,
which we use to summarize evidence for the region
containing any QTNs. Partly because of this, in practice wesuggest averaging results over several values for ra (equivalent
to placing a prior on ra). It may also be helpful to examine
sensitivity of results to ra. For example, if the BF is small for
all values of ra, then there is no evidence for any QTN in theregion; if the BF is large for some values and small for others,
then the evidence depends on the extent to which you believein large versus small effects. However, for simplicity, allresults in this paper were obtained using a ﬁxed value of ra ¼
0.5.

where H0 denotes the null hypothesis that none of the SNPs isa QTN (aj ¼ dj ¼ 0 for all j), and H1 denotes the complementary
event (i.e., at least one SNP is a QTN). Computing the BF
involves integrating out unknown parameters, as described inProtocols S1 and S2. In interpreting a BF, it is helpful to bearin mind the formula ‘‘posterior odds ¼ prior odds 3 BF,’’ so,for example, if the prior odds are 1:1 (i.e., p0 ¼ 0.5, so
association with genetic variation in the region is consideredequally plausible, a priori, as no association) then a BF of 10
gives posterior odds of 10:1, or ;91% probability of an
association.
In the special case where we allow at most one QTN,
Equation 5 reduces tonSXPðyjG; Hj Þ=PðyjG; H0 Þ
ð6Þ
BF ¼ ð1=nS Þj¼1

where Hj denotes the event that SNP j is the QTN. The jthterm in this sum corresponds to the BF for Hj versus the nullmodel, and involves the genotype data at SNP j only. We referto these terms as the ‘‘single-SNP’’ BFs, so in this special casethe overall BF is the mean of the single-SNP BFs. This natural
way for combining information across (potentially correlated) SNPs is an attractive property of BFs compared withsingle-SNP p-values. Furthermore, in terms of detecting a
genotype–phenotype association it can work well even if
multiple QTNs are present (see Results).
The Bayes/non-Bayes compromise. From a Bayesian viewpoint, the BF provides the measure of the strength of evidence
for genotype–phenotype association. That is, if one accepts
our prior distributions and modeling assumptions, then the
BF is all that is necessary to decide whether a genuine
association is present. However, given the potential for
debate over prior distributions, and for deviations frommodeling assumptions, it is helpful to note that a p-value fortesting H0 can be obtained from a BF through permutation.
Speciﬁcally, one can compute the BF for the observed data,
and for artiﬁcial data sets created by permuting observed
phenotypes among cohort individuals, and obtain a p-value as
the proportion of permuted data sets for which the BF
exceeds the BF for the observed data. Being based on
permutation, the resulting p-value is valid irrespective of whether
the model or priors are appropriate. This p-value also provides ahelpful way to compare our approach with standard tests ofassociation, and, as we show below, tests based on BF appear
to perform well in a wide variety of situations. Using BFs as
test statistics to obtain p-value is referred to as the ‘‘Bayes/
non-Bayes compromise’’ by Good [17].
Explaining and interpreting associations. To ‘‘explain’’observed associations we compute posterior distributions for
SNP effects (aj þ dj and 2aj for the heterozygote and minorallele homozygote, respectively), with particular focus on theposterior probability that each SNP is a QTN P(aj 6¼ 0). Here,our Bayesian regression approach has an important qualitative advantage over standard multiple regression. Specifically, if a genetic region contains multiple highly correlated
SNPs, each highly correlated with the phenotype, then thecorrect conclusion would be any of these SNPs could be
causal, without identifying which one. This will be reﬂected in
the posterior distribution of the effects: the overall proba-

Inference
We focus on two key inferential problems: (i) detecting
association between genotypes and phenotype, and (ii)explaining observed associations. In the model of Equation 1,
these translate to answering (i) are any bjs non-zero? and (ii)which bjs are non-zero and how big are they? We view theability to address both questions within a single framework to
be an advantage of our approach.
Detecting association. To measure the evidence for any
association between genotypes and phenotypes, we use theBF, [16] given by
PLoS Genetics | www.plosgenetics.org

ð5Þ

1299

July 2007 | Volume 3 | Issue 7 | e114

Imputation-Based Analysis of Association Studies

bility that at least one SNP is a QTN will be high, but (at leastin the simplest case where we assume at most one QTN) thisprobability will be spread out over the multiple correlated
SNPs. In contrast, if multiple highly correlated SNPs areincluded in a standard multiple regression it is possible thatno one of them will produce a signiﬁcant p-value.
We also argue that the imputation-based approach brings
us closer to being able to interpret estimated effects for each
SNP as actual causal effects, rather than simply associations.
Indeed, the key to making the leap from association tocausality is controlling for all potential confounding factors,
and by imputing genotypes at nearby SNPs, the imputationbased approach controls for one important set of confounding factors (the nearby SNPs), which would otherwise beignored. Thus, while functional studies provide the ultimate
route to convincingly demonstrating causal effects, ourapproach may help target such studies on the most plausible
candidate SNPs.

Imputing genotypes
In the tagSNP design, observed genotypes Gobs consist ofpanel genotypes at all SNPs and cohort genotypes at tagSNPs
only. To apply our methods in this situation, we use samplingbased algorithms (PHASE [10,11], or fastPHASE [5]) togenerate multiple imputations for the complete genotype
data (all individuals at all SNPs) by sampling from PðGjGobs Þ.
We then incorporate these imputations into our inference:for prior D1, this involves adding a step in the MCMC scheme
to sample the imputed genotypes from their posterior
distribution given all the data; for prior D2 it involves simply
averaging relevant calculations over imputations. Details are
given in Protocols S1 and S2.

Availability of Software
Methods described here are implemented in a software
package, Bim-Bam (Bayesian IMputation-Based Association
Mapping), available from the Stephens Lab website http://
stephenslab.uchicago.edu/software.html.

Results‘‘Power’’ and Comparisons with Other Approaches
We compared the power of our approach to other common
approaches via simulation. We simulated genotype andphenotype data (with l ¼ 0 and s ¼ 1) for genetic regions oflength 20 kb containing a single QTN, and genetic regions oflength 80 kb containing four QTNs, as follows:(1) Using a coalescent-based simulation program, msHOT[18], simulate 600 haplotypes from a constant-sized randommating population, under an ‘‘inﬁnite sites’’ mutation model,
with (population-scaled) mutation rate h ¼ 0.4/kb and ‘‘background’’ recombination rate q ¼ 0.8/kb, and a recombination
hotspot (width 1 kb; recombination rate 50q per kb) in thecenter of the region.
(2) Form genotypes for a ‘‘panel’’ of 100 individuals byrandomly pairing 200 haplotypes, and a ‘‘cohort’’ of 200
individuals by randomly pairing the other 400 haplotypes.
(3) Select tag SNPs from the panel data using the approachof Carlson et al. [19] with an r2 cutoff of 0.8. As in Carlson etal. [19], SNPs with panel minor allele frequency (MAF) ,0.1
were not tagged.
(4) Select which SNPs are QTNs, and their effect sizes, and
PLoS Genetics | www.plosgenetics.org

simulate phenotype data for each cohort individual according to Equation 1. We considered four scenarios: (A) a‘‘common’’ (MAF.0.1) QTN, with a range of effect sizes a ¼
0.2, 0.3, 0.4, 0.5 and ‘‘mild’’ dominance for the minor allele (d¼ 0.4a); (B) a common QTN, with a ¼ 0.3 and ‘‘strong’’
dominance for the major allele (d ¼ a); (C) a ‘‘rare’’ (MAF
0.01  0.05) QTN, with a ¼ 1 and no dominance (d ¼ 0); (D)four common, relatively uncorrelated, QTNs, each with a ¼
0.3 and d ¼ 0.4a. In each situation, we randomly chose a QTN
satisfying the relevant MAF requirements (in the 600 sampled
haplotypes), except under scenario (D) we ﬁrst chose four tag
SNP ‘‘bins’’ at random and then randomly chose a QTN
satisfying the MAF requirement in each bin, thereby ensuring
the four QTNs were relatively uncorrelated. (While real data
may contain multiple highly correlated QTNs, we did not
explicitly consider this case, since their effect would be
similar to a single QTN.)
We compared power of tests based on the BF (under prior
D2, allowing at most one QTN, using Equation 6 with fourother signiﬁcance tests:(1) Two tests based on pmin, the minimum p-value obtained
from testing each SNP individually (via standard ANOVAbased methods) for association with the phenotype. These
two tests differed in whether the single SNP p-values wereobtained using the 1 degree-of-freedom (df) ‘‘allelic’’ test,
which assumes an additive model where the mean phenotype
of heterozygotes lies midway between the two homozygotes(equivalent to linear regression of phenotype on genotype),or the 2 df ‘‘genotype’’ test, which treats the mean of theheterozygotes as a free parameter.
(2) A test based on preg, the global p-value obtained fromlinear regression of phenotype on all SNP genotypes (using
the standard F statistic, coding the genotypes as 0,1 and 2 ateach SNP, and assuming additivity across SNPs). See Chapman et al. [20] for example.
(3) A test based on BFmax, the maximum single-SNP BF. We
included this test for comparison with the mean single-SNP BF(Equation 6), to examine whether averaging information
across SNPs in Equation 6 improved power.
For each test, we analyzed each dataset in two ways: as if
data had been collected using (i) a ‘‘resequencing design’’ (i.e.,
all individuals were completely resequenced, so genotype
data are available at all SNPs in all individuals); and (ii) a ‘‘tag
SNP design’’ (i.e., in panel individuals genotype data are
available at all SNPs, but in cohort individuals genotype data
are available at tag SNPs only). For the tag SNP design, we
assumed haplotypic phase is known in the panel (as it is,mostly, for the HapMap data for example), but not in thecohort; however our approach can also deal with unknown
phase in the panel. For preg and pmin, tests were performed on
all SNPs for the resequencing design, and on tag SNPs only
for the tag SNP design. For BF and BFmax, single-SNP BFs
were computed for all SNPs in both designs (averaging over
imputed genotypes for non-tag SNPs in the tag SNP design).
For preg, we computed a p-value assessing signiﬁcance using
the standard asymptotic distribution for the F statistic; forthe other tests we found p-values by permutation, using 200–
500 random permutations of phenotypes assigned to cohort
individuals. (The relatively small number of permutations
limits the size of the smallest possible p-value, causing
discontinuities near the origin in Figure 1).
Figure 1 shows power of each test versus type I error under
1300

July 2007 | Volume 3 | Issue 7 | e114

Imputation-Based Analysis of Association Studies

PLoS Genetics | www.plosgenetics.org

1301

July 2007 | Volume 3 | Issue 7 | e114

Imputation-Based Analysis of Association Studies
Figure 1. Power Comparisons(A) single common variant, modest dominance; (B) single common variant, strong dominance for minor allele; (C) single rare variant, no dominance; (D)multiple common variants.
Each colored line shows power of test varying with significance threshold (type I error). Black: BF from our method (prior D2); Green: pmin (allelic test);
Red: pmin (genotype test); Blue: preg, multiple regression; Grey: BFmax. Each column of figures shows results for data analyzed under the ‘‘resequencing
design’’ (left) and the ‘‘tag SNP design’’ (right). Each row shows results for the four different simulation scenarios.
doi:10.1371/journal.pgen.0030114.g001

both resequencing and tag SNP designs. For Scenario (A) (asingle common QTN), the relative performances of methods
were similar for all four effect sizes examined (unpublisheddata), and so we pooled these results in the ﬁgure.
Comparing pmin and preg, the single-SNP tests (pmin) weremore powerful when all variants (including the causal
variant) were typed, or when the QTN was a common SNP
and therefore ‘‘tagged’’ by a tag SNP, while the regressionbased approach (preg) was more powerful when the QTN wasa rare SNP not ‘‘tagged’’ by any tag SNP. Among the two
single-SNP tests, the 1 df allelic test performed as well as, or
better than, the 2 df genotypic test, except in Scenario (B),where the major allele exhibits strong dominance. In
particular, for Scenario (A), where the causal variant
exhibits dominance, the allelic test (which assumes nodominance), performed better than the genotypic test. Thisis presumably because, with the effect and sample sizes
considered, the extra parameter estimated in the genotypic
test does not sufﬁciently improve model ﬁt. Althoughrelative performance of pmin and preg in the tag SNP design
could depend on tag SNP selection scheme (and the one weused, based on pairwise LD, would seem to favor pmin), it
seems reasonable to expect single-SNP tests to be effectiveat detecting ‘‘direct’’ associations between the phenotype
and a causal variant, or ‘‘near-direct’’ association between a
SNP that tags a causal variant, and the regression-based
approach to be better at detecting indirect associations
between a phenotype and a variant not ‘‘tagged’’ by a single
SNP (the intuition, from Chapman et al. [20], is that suchvariants can be highly correlated with linear combinationsof tag SNPs, and thus be detected by linear regression). In
principle, preg could also effectively capture ‘‘direct’’associations, but our empirical results suggest that it is lesseffective at this than the single SNP tests. (However, poor
performance of preg under the resequencing design may be
due in part to inadequacy of the asymptotic theory whenlarge numbers of correlated covariates are used. This might
be alleviated by assessing signiﬁcance of preg by permutation.)
Turning now to our approach, except for Scenario (B) in
the tag SNP design, the test based on the BF is as powerful ormore powerful than the other tests. Thus, unlike preg and pmin,
the BF performs well in detecting both ‘‘direct’’ and‘‘indirect’’ associations: if the QTN is typed, the BF detects
it using observed genotype data at that SNP; otherwise, it
detects it using the imputed genotype data at the QTN. In
Scenario (B), where the major allele exhibits strong dominance, our approach suffered slightly in power comparedwith the genotypic test, presumably because our prior places
relatively low weight on strong dominance. However, thepower loss was small compared with that of the allelic test.
Thus our prior ‘‘allows’’ for dominance without suffering thefull penalty incurred by the extra parameter in the genotypic
test when dominance is less strong (Scenario [A]).
PLoS Genetics | www.plosgenetics.org

In Scenario (D), which involved multiple QTNs, tests basedon the BF clearly outperformed other tests considered, eventhough the BF was computed allowing at most one QTN.
Our explanation is that the BF, being the average of singleSNP BFs, has greater opportunity to capture the presence ofmultiple QTNs than does the minimum p-value. This
explanation is supported by the fact that the maximum BF,
BFmax, performs less well than BF. To examine whether
power might be further increased by explicitly allowing for
multiple QTNs, we compared power for BFs computed using
1-QTN and 2-QTN models (in the 2-QTN model p(l ¼ 1) ¼ p(l¼ 2) ¼ 0.5). We found little difference in power, although BFs
for the 2-QTN model tended to be larger than BFs for the 1QTN model, so allowing for multiple QTNs may help if the
BF itself, rather than a p-value based on the BF, is used tomeasure the strength of evidence for association. In
addition, considering multiple-QTN models should haveadvantages when attempting to explain an association (seebelow).
A second, and perhaps more surprising, situation where
the BF outperforms other methods is when all SNPs are
typed and tested (i.e., Scenario (A), resequencing design).
Here, in contrast to Scenario (D), BFmax performs similarly tothe standard BF, suggesting that the power gain is due not toaveraging, but to an intrinsic property of single-SNP BFs
that makes them better measures of evidence than singleSNP p-values. Our explanation is that the BF tends to be lessinﬂuenced by less informative SNPs (e.g., those with verysmall MAF, of which there are many in the resequencing
design), whereas p-values tend to give equal weight to all
SNPs, regardless of information content. Speciﬁcally, BFs for
relatively uninformative SNPs will always lie close to 1, andshould not greatly inﬂuence either the maximum or the
average of the single-SNP BFs (or, more precisely, will not
greatly inﬂuence differences in these test statistics among
permutations of phenotypes). In contrast, p-values for each
SNP are forced, by deﬁnition, to have a uniform distribution
under H 0, and so p-values from a large number of
uninformative SNPs unassociated with the phenotype could
swamp any signal generated by a single informative SNP
associated with the phenotype. Although the resequencing
design is currently uncommon, this observation suggests thatit may generally be preferable to rank SNPs according totheir BFs, rather than by p-values (e.g., in genome scans). It
also highlights a general (rarely considered, and perhaps
underappreciated) drawback of p-values as a measure ofevidence: the strength of evidence of a given p-value depends
on the informativeness of the test being performed, or, more
speciﬁcally, on the distribution on the p-values under thealternative hypothesis, which is generally not known. Thus,
for example, a p-value of 105 in a study involving few
individuals may be less impressive than the same p-value in alarger study. In contrast, the interpretation of a BF does notdepend on study size or similar factors.
1302

July 2007 | Volume 3 | Issue 7 | e114

Imputation-Based Analysis of Association Studies

Figure 2. Comparison of Results for Resequencing Design (x-axis) and Tag SNP Design (y-axis)
Panels show: (a) errors in the estimates (posterior means) of the heterozygote effect (a þ d); (b) errors in the estimates (posterior means) of the main
effect (a); and (c) posterior probability of being a QTN (P((a, d) 6¼ (0, 0))) assigned to the causal variant.
doi:10.1371/journal.pgen.0030114.g002

Resequencing versus Tag SNP Designs

noticeable drop in power for the tag SNP design versus theresequencing design, and the BFs, posterior probabilities, andeffect size estimates under the two designs often differ
substantially (unpublished data). This may seem slightly
disappointing: one might have hoped that, even with tag
SNPs chosen to capture common variants, they might also
capture some rare variants. Indeed, this can happen: in somesimulated data sets the rare causal variant was clearly
identiﬁed by our approach, presumably because it was highly
correlated with a particular haplotype background, and could
thus be accurately predicted by tag SNPs. However, thisoccurred relatively rarely (just a few simulations out of 100).
We wondered whether a different tagging strategy, aimed
at capturing rare variants, might improve performance when
the causal variant is rare. The development of such strategies
lies outside the scope of this paper, but, to assess potential
gains that might be achieved, we analyzed rare-variant
simulations assuming that all SNPs except the causal variant
were typed in the cohort. Power from this approach (Figure3) gives a conservative upper bound on what could beachieved using a more effective tagging design, without
actually typing the causal variant. Although power was higher
than with the r2-based tag SNP selection, it remained
substantially lower than in the resequencing design, where
the causal variant is typed.
We also wondered whether a different approach to impute
missing genotypes (in the cohort at non-tag SNPs) might
improve performance. For results above, we used the software
fastPHASE [5] to impute the genotypes, so we re-ran theanalysis using a different imputation algorithm [10,11].
Results for these two approaches (Figure 4) show little
difference in terms of power, consistent with previous results[5] suggesting the two approaches have similar accuracy in
imputing missing genotypes.
In summary, imputation-based methods appear to increase
power of the tag SNP design to detect rare variants, butnevertheless remain notably less powerful than BFs based on
the complete resequencing data.

An important feature of Figure 1 is that, for Scenarios (A),(B), and (D), where the causal SNPs are common, power issimilar for the resequencing and tag SNP designs. Indeed, inthese cases most other aspects of inference are also similar.
For example, Figure 2 shows that, under Scenarios (A) and(B), estimated effect sizes, BFs, and posterior probability thatthe actual causal variant is a QTN, are typically similar forboth designs. Thus under these scenarios, our imputationbased approach effectively recreates results that would have beenobtained by resequencing all individuals.
In contrast, when the causal variant is rare, there is a

Figure 3. Examination of Potential Effect of Different Tag SNP Strategies
on Power, When the Causal Variant is Rare (0.01 , MAF , 0.05)
Solid line: Resequencing design; dashed line: tag SNP design, with tags
selected using method from [19]; and dotted line: tag SNP design, withall SNPs except the causal SNP as tags.
doi:10.1371/journal.pgen.0030114.g003
PLoS Genetics | www.plosgenetics.org

Comparison of Prior D1 and D2
Priors D1 and D2 differ in their assumed correlation
between the dominance effect (d ¼ ak) and main effect a: in1303

July 2007 | Volume 3 | Issue 7 | e114

Imputation-Based Analysis of Association Studies

Figure 4. Power of the Multipoint Approach in the Rare Variant Scenario
for Two Different Imputation Algorithms
doi:10.1371/journal.pgen.0030114.g004

Figure 5. Scatter Plot of Samples from Prior Distribution of a (x-axis) anda þ d (y-axis), for Priors D1 (Black) and D2 (Blue)The solid yellow line corresponds to d ¼ 0 (additivity). The dashed redlines are the limits above and below which a SNP exhibits overdominance.
doi:10.1371/journal.pgen.0030114.g005

D1 the prior probability of overdominance is independent ofa, whereas under D2 overdominance is more likely for small athan for large a (Figure 5). In this respect, D1 is perhaps moresensible than D2; however, D2 is computationally much
simpler. To examine the effects of these priors on inference,
we compared (i) the BF and (ii) the posterior probability
assigned to the actual causal variant under each prior for thedatasets from Scenarios (A) and (B). Results agreed quite
closely (Figure 6), suggesting prior D2 provides a reasonableapproximation to prior D1 in the scenarios considered. Thisis important, since prior D2 is computationally practical forcomputing BFs for very large datasets (e.g., genome-wide
association studies with hundreds of thousands of SNPs), forwhich sampling posterior distributions of parameters using
an MCMC scheme would be computationally daunting.

this possibility may provide substantially higher BFs. Figure 7
shows the marginal posterior probabilities for each SNP beinga QTN, under the one-QTN and multi-QTN models, conditional on at least one SNP in the region being a QTN.
(Summarising the more complex information on posterior
probabilities for combinations of SNPs is an important futurechallenge.) Under the one-QTN model, only one of the fourcausal SNPs has a large marginal posterior probability,
whereas under the multi-QTN model all four are moderately
large. Of course, other SNPs correlated with the four QTNs
were also associated with the phenotype, and so have elevatedposterior probabilities. This example illustrates the potential
for the multi-QTN model to provide fuller explanations for
associations.

Allowing for Multiple Causal Variants
When analyzing a candidate region, one would ideally likenot only to detect any association, but also to identify thecausal variants (QTNs). Since a candidate region could
contain multiple QTNs, we implemented an MCMC scheme(using prior D1) to ﬁt multi-QTN models where the numberof QTNs is estimated from the data; here, we consider amulti-QTN model with equal prior probabilities on 1, 2, 3, or
4 QTNs. (A similar MCMC scheme could also be implemented
for prior D2, and could exploit the analytical advantages ofthis prior to reduce computation. Indeed, for regions
containing a modest number of SNPs it would be possible
to examine all subsets of SNPs, and entirely avoid MCMC.)
We compare this multi-QTN model with a one-QTN model
on a dataset simulated with four QTNs (scenario [D]). The
estimated BF for a one-QTN model was ;6,000, while for themulti-QTN model it was .105 (we did not perform sufﬁcient
iterations to estimate how much bigger than 105). Thus, if aregion contains multiple causal variants, then allowing for
PLoS Genetics | www.plosgenetics.org

SCN1A Polymorphism and Maximum Dose of
Carbamazepine
We applied our method to data from association studies
involving the SCN1A gene and the maximum dose of
carbamazepine in epileptic patients [21,22]. For this analysis,the ‘‘panel’’ consisted of parents from 32 trios of European
descent from the CEPH Utah collection [21] and the ‘‘cohort’’consisted of 425 patients of European descent for whom the
maximum dose of carbamazepine had been determined [22].
Genetic data on the trios were available for 15 polymorphisms,comprising 14 SNPs and one indel, which corresponded to snps 1–15 and indel12 in Table 2 of Weale et al.
[21]. For cohort individuals, genotype data are available atfour tag SNPs: snp1 (rs590478), snp5 (rs8191987), snp7
(rs3812718), and snp9 (rs2126152). These SNPS were chosen
to summarize haplotype diversity at the 15 panel polymorphisms (for details, see Tate et al. [22]).
1304

July 2007 | Volume 3 | Issue 7 | e114

Imputation-Based Analysis of Association Studies

Figure 6. Comparison of Inferences using Prior D1 and D2 for the BF (Left) and the Posterior Probability Assigned to the Causal Locus Being a QTN(Right)Results shown are for all datasets for the common variant Scenario (A) and (B) and for both the resequencing design and the tag SNP design. The
discrepancy between the larger estimated BFs is caused by the fact that we used insufficient MCMC iterations to accurately estimate very large BFs(.106) under prior D1.
doi:10.1371/journal.pgen.0030114.g006

We ﬁrst estimated haplotypes in 64 parents using the trio
option in PHASE [23]. Since trio information allows
haplotypes to be accurately determined [23] we assumed
these estimated panel haplotypes were correct in subsequent
analyses. We then applied our method to compute a BF for
overall association between genetic data and the phenotype,
and to compute, for each SNP, the posterior probability that
it was a QTN. In applying our method we used PHASE to
impute the genotypes in the cohort at non-tag SNPs, andperformed analyses under priors D1 and D2.
BFs for priors D1 and D2 were, respectively, 3.15 and 2.33,
and the corresponding p-values (estimated using 1,000
permutations) were 0.006 and 0.019, respectively. We also
computed p-values using single SNP tests at tag SNPs andobtained 0.007 for the allelic test and 0.019 for the genotype
test. (These are essentially the two tests performed by Tate etal. [22], who reported the smallest p-values uncorrected for
multiple comparisons.) These BFs represent only modest
evidence for an association. If one were initially evensomewhat skeptical about SCN1A as a candidate for inﬂuencing this phenotype, one might remain somewhat skepticalafter analyzing these data. For example, with a 20% priorprobability on variation in SCN1A inﬂuencing phenotype, theposterior probability of association under either prior is,50%. (Prior probability of 0.2 gives prior odds of 0.2:(1–0.2),
or 1:4; a BF of 3 then gives posterior odds of 3:4, whichtranslates to a posterior probability of 3/7.) On the other
hand, SCN1A might be considered a relatively good candidate
for inﬂuencing response to carbamazepine, since it is the
drug’s direct target. And, depending on follow-up costs andpotential beneﬁts of ﬁnding a functional variant, posterior
probabilities of very much ,50% might be deemed worth
following-up.
PLoS Genetics | www.plosgenetics.org

Among the 15 SNPs analyzed, snp7 was assigned the highest
posterior probability of being a QTN (Figure 8). This SNP,
which is a tag SNP, was also implicated by the analysis in Tateet al. [22]. However, the posterior probability of this SNP
represents only 34 % of the posterior mass. Six additional
SNPs are needed to encompass 90% of the posterior mass:
snp6 (rs3812719), snp8 (rs490317), snp9 (rs2126152), snp10(rs7601520), snp11 (rs2298771) and snp13 (rs7571204). The
posterior distributions of the main effect, a, for each of theseseven SNPs, conditional on it being a QTN, are very similar(Figure 8).
In summary, these data provide modest evidence of
association between SCN1A and maximum dose of carbamazepine, and, among the SNPs analyzed, snp7 (rs3812718)
appears to be the best candidate for being causal. A recent
follow-up study appears to conﬁrm this variant as being
functionally important [24].

Discussion
We described a new approach for analysis of association
studies, with two important components: (i) it uses imputation of unknown genotypes, based on statistical modeling of
patterns of LD, to allow untyped SNPs to be directly assessed
for association with phenotype; (ii) it uses BFs, rather than pvalues, to assess genotype–phenotype association.
The idea of trying to ﬁnd associations between phenotypes
and untyped variants is old, and underlies many existing
methods for assessing association. In some cases this aim is
implicit (e.g., testing for association between haplotypes andphenotypes can be thought of as an attempt to indirectly test
untyped variants that may lie on a particular haplotype
background), and in others it is explicit (to give just one
example, Zöllner and Pritchard [25] place mutations on an
1305

July 2007 | Volume 3 | Issue 7 | e114

Imputation-Based Analysis of Association Studies

Figure 7. Illustration of How a Multi-QTN Model Can Provide Fuller Explanations Than a One-QTN Model for Observed Associations
The figure shows, for each SNP in a dataset simulated under Scenario (D), the estimated posterior probability that it is a QTN, conditional on an
association being observed. Left: Results from one-QTN model. Right: Results from multi-QTN model allowing up to four QTNs. The four actual QTNs areindicated with a star. Colors of the vertical lines indicate tag SNP ‘‘bins’’ (i.e., groups of SNPs tagged by the same variant).
doi:10.1371/journal.pgen.0030114.g007

of freedom incurred by the general 2 df hypothesis test; (ii)the BF better reﬂects the informativeness of each SNP, inparticular, that SNPs with small MAF are typically lessinformative than SNPs with larger MAF (this advantage
presumably being greatest for SNP panels containing many
SNPs with small MAF); (iii) it provides a principled way totake into account prior information on each SNP, e.g.,whether it lies in or near a gene whose function is believed
likely to inﬂuence the trait; and (iv) averaging single-SNP BFs
provides a convenient, and in some ways effective, approachto combining information across multiple SNPs in a region.
Perhaps the most important disadvantage of BFs comparedwith p-values is that a BF is strictly ‘‘valid’’ only under theassumption that both the prior and the model are ‘‘correct.’’
Since this is never the case in practice, BFs are never strictly
valid, Our hope is to make the prior and model sufﬁciently
accurate that resulting BFs are ‘‘useful.’’ (Note that p-values
may be valid but useless: e.g., p-values simulated from a
uniform distribution independent of phenotype and genotype data are valid, in that they are uniformly distributed
under the null hypothesis, but useless.) Here, it is helpful todistinguish two different uses of BFs: as test statistics tocompute permutation-based p-values, as in the powercomparisons in this paper, and as direct measures of evidence(e.g., in ‘‘posterior odds ¼ BF 3 prior odds’’). Our limited
experience is that p-values obtained from BFs are relatively
robust to prior and modeling assumptions, but that theabsolute values of BFs are substantially more sensitive. In
particular, BFs tend to be sensitive to both (i) choice of ra, rd;and (ii) the normality assumption in the phenotype model.
We now discuss each of these issues in turn.
Choice of ra, rd corresponds to quantifying prior beliefs
about likely additive and dominance effect sizes. In thispaper, we used (in prior D2) ra ¼ 0.5 and rd ¼ ra/2. We now
believe these values are likely larger than appropriate for

estimated tree, and test resulting genotypes for association
with phenotype). A key difference between our approach andthese existing methods is that we focus on testing variants
about which something is known (i.e., SNPs that are known toexist, and have documented patterns of LD), and exploiting
this information. This idea, which seems in many ways morecompelling than testing hypothetical untyped variants about
which nothing is known, has been recently developed byseveral groups [5,26–31]. While there are, no doubt, multiple
effective ways to implement the general strategy, attractive
key features of our approach include the use of ﬂexible
statistical models for multi-locus LD to estimate missing
genotypes, with uncertainty; and the use of Bayesian methods
to account for uncertainty in estimated genotypes.
While several papers have suggested Bayesian approaches
to association studies (e.g., [15,32,33]), our work includes somedistinctive contributions. First, our prior distributions for
single-SNP effects have a number of desirable properties: (i)they scale appropriately with changes in measurement unitsof the phenotype, (ii) they center on an additive model whileallowing for dominance, and (iii) they facilitate rapid
calculations. This last feature means that our work can form
the foundation of simple Bayesian analyses in genome-wide
association studies, e.g., computing a single-SNP BF for each
SNP, as a Bayesian analogue of single-SNP hypothesis tests.
This option is available in our software, but to further
facilitate its use by others, and to emphasize the simplicity of
the analytical calculations, we give R code for computing the
BF for typed SNPs under prior D2 (see Protocol S1). A second
distinctive contribution is that we compare our Bayesian
approach directly with standard p-value based approaches,
providing both qualitative insight and quantitative support
for several advantages of single-SNP BFs over single-SNP pvalues. These advantages include: (i) the BF allows for both
additive and dominant effects without the additional degree
PLoS Genetics | www.plosgenetics.org

1306

July 2007 | Volume 3 | Issue 7 | e114

Imputation-Based Analysis of Association Studies

Figure 8. Results for the SCN1A Dataset
Left panel shows the posterior probability assigned to each SNP being a QTN, with filled triangles denoting tag SNPs and open circles denoting non-tag
SNPs. The right panel shows (in gray) estimated posterior densities of the additive effect for each of the seven SNPs assigned the highest posterior
probabilities of non-zero effect (representing 90% of the posterior mass). The average of these curves is shown in black.
doi:10.1371/journal.pgen.0030114.g008

G, and to perform inference for the genetic effects
parameters, b, based on the conditional likelihood L(b) ¼P(yj fyg, G, b). However, this conditional likelihood does notappear to be analytically tractable, and so analysis of thisdesign may require development of computationally tractable approximations. Similarly, adapting our approach to
standard case-control designs will require development ofappropriate priors and computational algorithms, andrepresents an important area for future work.

most studies of complex phenotypes, placing too little weighton small, but realistic, effect sizes. Our current suggested‘‘default’’ procedure is to average BFs computed with ra ¼
0.05, 0.1, 0.2, and 0.4, and rd ¼ ra/4, which places more weighton smaller effect sizes, and less weight on overdominance. We
would expect to modify these values in the light of further
information about typical effect sizes for particular traits. It
could also be argued that, in addition to allowing acontinuum of deviations from the additive model, it may
make sense to specify prior probabilities for ‘‘pure’’ recessiveor dominant models (i.e., d ¼ a, a). BFs under these models
can be computed easily by simply replacing all heterozygous
genotypes with homozygous genotypes for the major orminor allele.
Regarding the normality assumption, following a suggestion by Mathew Barber (personal communication), in
practical applications, we are currently applying a normal
quantile transform to phenotypes (replacing the rth biggest ofn observations with the (r  0.5)/nth quantile of the standard
normal distribution) before applying our methods andcomputing BFs. Imposing normality on our phenotype inthis way is different from the normality assumption in our
phenotype model, which states that the residuals are normally
distributed. However, in this context, where effect sizes are
expected to be generally rather small, normality of phenotype
and normality of residuals are somewhat similar assumptions,suggesting that this transform may be effective.
Throughout this paper, we have assumed a ‘‘population’’sampling design in which phenotype and genotype data are
available on a random sample from a population, andperform analyses conditional on the observed genotype data.
An alternative common design involves collecting genotypes
only on individuals whose phenotypes lie in the tails of the
distribution [34]. To apply our methods to such designs, wesuggest conditioning on unordered observed phenotypes,
denoted fyg, in addition to conditioning on the genotypes
PLoS Genetics | www.plosgenetics.org

Supporting Information
Protocol S1. Analytical Computations for Prior D2
Found at doi:10.1371/journal.pgen.0030114.sd001 (100 KB PDF).
Protocol S2. MCMC Sampling for Prior D1
Found at doi:10.1371/journal.pgen.0030114.sd002 (66 KB PDF).
Accession Numbers
The National Center for Biotechnology Information (NCBI) Entrez(http://www.ncbi.nlm.nih.gov/gquery/gquery.fcgi) Gene ID of theSCN1A gene is 6323.

Acknowledgments
We thank D. Goldstein for access to the SCN1A data, and M. Weale
and S. Tate for providing the data in a convenient electronic form.
We thank N. Patterson for pointing us to the I. Good reference, J.
Marchini and P. Donnelly for helpful conversations, and J. Pritchard
and two anonymous referees for useful comments on earlier versionsof the manuscript. Computing support was provided by the
University of Washington Center for Study of Demography andEcology, High Performance Computing Cluster Cooperative.
Author contributions. BS and MS conceived and designed theexperiments, performed the experiments, analyzed the data, contributed reagents/materials/analysis tools, and wrote the paper.
Funding. This work was supported by NIH grant RO1 HG02585–01
to MS.
Competing interests. The authors have declared that no competinginterests exist.
1307

July 2007 | Volume 3 | Issue 7 | e114

Imputation-Based Analysis of Association Studies
References
1. The International HapMap Consortium (2005) A haplotype map of thehuman genome. Nature 437: 1299–1320.
2. SeattleSNPs. Seattle (Washington): NHLBI Program for Genomic Applications. Available: http://pga.gs.washington.edu. Accessed 12 June 2007.
3. Kraft P, Pharoah P, Chanock SJ, Albanes D, Kolonel LN, et al. (2005)
Genetic variation in the HSD17B1 gene and risk of prostate cancer. PLoS
Genet 1: e68. doi:10.1371/journal.pgen.0010068
4. Li N, Stephens M (2003) Modeling linkage disequilibrium and identifying
recombination hotspots using single-nucleotide polymorphism data.
Genetics 165: 2213–2233.
5. Scheet P, Stephens M (2006) A fast and ﬂexible statistical model for largescale population genotype data: Applications to inferring missing
genotypes and haplotypic phase. Am J Hum Genet 78: 629–644.
6. Almasy L, Blangero J (1998) Multipoint quantitative-trait linkage analysis in
general pedigrees. Am J Hum Genet 62: 1198–1211.
7. George EI, McCulloch RE (1997) Approaches to Bayesian variable selection.
Statistica Sinica 7: 339–373.
8. Nott JD, Green PJ (2004) Bayesian variable selection and the SwendsenWang algorithm. J Comput Graph Statist 13: 141–157.
9. Yang X, Belin TR, Boscardin WJ (2005) Imputation and variable selection in
linear regression models with missing covariates. Biometrics 61: 498.
10. Stephens M, Smith N, Donnelly P (2001) A new statistical method for
haplotype reconstruction from population data. Am J Hum Genet 68: 978–989.
11. Stephens M, Scheet P (2005) Accounting for decay of linkage disequilibrium in haplotype inference and missing-data imputation. Am J Hum
Genet 76: 449–462.
12. Lynch M, Walsh B (1998) Genetics and analysis of quantitative traits.
Sunderland (Massachussetts): Sinauer Associates. 980 p.
13. Jeffreys H (1946) An invariant form for the prior probability in estimation
problems. Proc R Soc Lond Ser A Math Phys Eng Sci 186: 453–461.
14. Raftery AE (1996) Approximate Bayes factors and accounting for model
uncertainty in generalised linear models. Biometrika 83: 251–266.
15. Lunn DJ, Whittaker JC, Best N (2006) A Bayesian toolkit for genetic
association studies. Genet Epidemiol 30: 231–47.
16. Kass RE, Raftery AE (1995) Bayes factors. J Amer Statist Assoc 90: 773–795.
17. Good I (1992) The Bayes/Non-Bayes compromise: A brief review. J Amer
Statist Assoc 87: 597–606.
18. Hellenthal G, Stephens M (2007) msHOT: Simulating crossover and geneconversion hotspots with Hudson’s (2002) ms simulator. Bioinformatics 23:
520–521.
19. Carlson CS, Eberle MA, Rieder MJ, Yi Q, Kruglyak L, et al. (2004) Selecting amaximally informative set of single-nucleotide polymorphisms for association analyses using linkage disequilibrium. Am J Hum Genet 74: 106–120.
20. Chapman JM, Cooper JD, Todd JA, Clayton DG (2003) Detecting disease

PLoS Genetics | www.plosgenetics.org

associations due to linkage disequilibrium using haplotype tags: A class oftests and the determinants of statistical power. Hum Hered 56: 18–31.
21. Weale ME, Depondt C, Macdonald SJ, Smith A, Lai PS, et al. (2003)
Selection and evaluation of tagging SNPs in the neuronal-sodium-channel
gene SCN1A: Implications for linkage-disequilibrium gene mapping. Am J
Hum Genet 73: 551–565.
22. Tate SK, Depondt C, Sisodiya SM, Cavalleri GL, Schorge S, et al. (2005)
Genetic predictors of the maximum doses patients receive during clinicaluse of the anti-epileptic drugs carbamazepine and phenytoin. Proc Natl
Acad Sci U S A 102: 5507–5512.
23. Marchini J, Cutler D, Patterson N, Stephens M, Eskin E, et al. (2006) A
comparison of phasing algorithms for trios and unrelated individuals. Am J
Hum Genet 78: 437–450.
24. Heinzen EL, Yoon W, Tate SK, Sen A, Wood NW, et al. (2007) Nova2
interacts with a cis-acting polymorphism to inﬂuence the proportions ofdrug reponsive splice variants of SCN1A. Am J Hum Genet. In press.
25. Zöllner S, Pritchard JK (2005) Coalescent-based association mapping andﬁne mapping of complex trait loci. Genetics 169: 1071–1092.
26. Soranzo N, Cavalleri GL, Weale ME, Wood NW, Depondt C, et al. (2004)
Identifying candidate causal variants responsible for altered activity of the
ABCB1 multidrug resistance gene. Genome Res 14: 1333–1344.
27. Marchini J, Myers S, McVean G, Donnelly P (2005) A novel Bayesian
approach to localising disease genes. The American Society of Human
Genetics 55th annual meeting. Available: http://www.ashg.org/genetics/
ashg05s/. Accessed 25 June 2007.
28. Pe’er I, De Bakker PIW, Maller J, Jones K, Altshuler MD, et al. (2005) Using
haplotype maps to improve power of whole genome association studies
with ﬁxed sets of markers. In: The American Society of Human Genetics
55th Annual Meeting. Available: http://www.ashg.org/genetics/ashg05s/.
Accessed 25 June 2007.
29. Servin B, Stephens M (2005) A Bayesian approach to the analysis ofcandidate genes association studies. In: The American Society of Human
Genetics 55th Annual Meeting. Available: http://www.ashg.org/genetics/
ashg05s/. Accessed 25 June 2007.
30. Nicolae DL (2006) Testing Untyped Alleles (TUNA)-applications togenome-wide association studies. Genet Epidemiol 30: 718–727.
31. Dai JY, Ruczinski I, LeBlanc M, Kooperberg C (2006) Imputation methods
to improve inference in SNP association studies. Genet Epidemiol 30: 690–
702.
32. Sillanpaa MJ, Bhattacharjee M (2005) Bayesian association-based ﬁne
mapping in small chromosomal segments. Genetics 169: 427–439.
33. Morris AP (2006) A ﬂexible Bayesian framework for modeling haplotype
association with disease, allowing for dominance effects of the underlying
causative variants. Am J Hum Genet 79: 679–694.
34. Cohen J, Pertsemlidis A, Kotowski IK, Graham R, Garcia CK, et al. (2005)
Low LDL cholesterol in individuals of African descent resulting fromfrequent nonsense mutations in PCSK9. Nat Genet 37: 161–165.

1308

July 2007 | Volume 3 | Issue 7 | e114