A Fast and Flexible Statistical Model for Large-Scale Population Genotype
Data: Applications to Inferring Missing Genotypes and Haplotypic Phase
Paul Scheet and Matthew Stephens
Department of Statistics, University of Washington, Seattle

We present a statistical model for patterns of genetic variation in samples of unrelated individuals from natural
populations. This model is based on the idea that, over short regions, haplotypes in a population tend to cluster
into groups of similar haplotypes. To capture the fact that, because of recombination, this clustering tends to belocal in nature, our model allows cluster memberships to change continuously along the chromosome according
to a hidden Markov model. This approach is flexible, allowing for both “block-like” patterns of linkage disequilibrium (LD) and gradual decline in LD with distance. The resulting model is also fast and, as a result, is practicablefor large data sets (e.g., thousands of individuals typed at hundreds of thousands of markers). We illustrate theutility of the model by applying it to dense single-nucleotide–polymorphism genotype data for the tasks of imputing
missing genotypes and estimating haplotypic phase. For imputing missing genotypes, methods based on this model
are as accurate or more accurate than existing methods. For haplotype estimation, the point estimates are slightly
less accurate than those from the best existing methods (e.g., for unrelated Centre d’Etude du Polymorphisme
Humain individuals from the HapMap project, switch error was 0.055 for our method vs. 0.051 for PHASE) but
require a small fraction of the computational cost. In addition, we demonstrate that the model accurately reflects
uncertainty in its estimates, in that probabilities computed using the model are approximately well calibrated. The
methods described in this article are implemented in a software package, fastPHASE, which is available from theStephens Lab Web site.
With the advent of cheap, quick, and accurate genotyping technologies, there is a need for statistical models
that both are capable of capturing the complex patternsof correlation (i.e., linkage disequilibrium [LD]) that exist among dense markers in samples from natural populations and are computationally tractable for large data
sets. Here, we present such a model and assess its ability
to accurately capture patterns of variation by applying
it to estimate missing genotypes and to infer haplotypicphase from unphased genotype data.
The model is motivated by the observation that, overshort regions (say, a few kilobases in human genomes),
haplotypes tend to cluster into groups of similar haplotypes. This clustering tends to be local in nature because, as a result of recombination, those haplotypes thatare closely related to one another and therefore similar
will vary as one moves along a chromosome. To capture
this, we allow the cluster membership of observed haplotypes to change continuously along the genome according to a hidden Markov model (HMM). (This idea
has been proposed, independently of our work, by others, including Sun et al. [2004], Rastas et al. [2005], andKimmel and Shamir [2005a]; see the “Discussion” section.) Each cluster can be thought of as (locally) representing a common haplotype, or combination of al-leles, and the HMM assumption for cluster membershipresults in each observed haplotype being modeled as amosaic of a limited number of common haplotypes (fig.
1). This approach seems more flexible than “blockbased” cluster models, which divide the genome into
blocks (segments of high LD) and allow cluster membership to change only across block boundaries (e.g.,
Greenspan and Geiger 2004; Kimmel and Shamir
2005b). The hope is that this more flexible assumption
will allow the model to capture complex patterns of LD
that are not well captured by block-based alternatives
while continuing to capture any “block-like” patterns
that are present.
Another model that also aims to flexibly capture patterns of LD is the PAC model of Li and Stephens (2003),
which partially underlies the PHASE software for haplotype inference and estimation of recombination rates(Stephens et al. 2001; Stephens and Donnelly 2003; Stephens and Scheet 2005). One way to view the model we
present here is as an attempt to combine the computational convenience of cluster-based models with the flexibility of the PAC model. Indeed, in terms of computational convenience, our model is substantially more
attractive than the PAC model, both in that computation
increases only linearly with the number of individuals

Received September 30, 2005; accepted for publication January 25, 2006; electronically published February 17, 2006.
Address for correspondence and reprints: Mr. Paul Scheet, University of Washington, Department of Statistics, Box 354322, Seattle, WA
98195-4322. E-mail: pscheet@alum.wustl.edu
Am. J. Hum. Genet. 2006;78:629–644.  2006 by The American Society of Human Genetics. All rights reserved. 0002-9297/2006/7804-0010$15.00

www.ajhg.org

The American Journal of Human Genetics Volume 78 April 2006

629

Figure 1

Illustration of how our model allows cluster membership to change continuously along a chromosome. Each column represents
a SNP, with the two alleles indicated by open and crossed squares. Successive pairs of rows represent the estimated pair of haplotypes for
successive individuals. Colors represent estimated cluster membership of each allele, which changes as one moves along each haplotype. Locally,each cluster can be thought of as representing a (common) combination of alleles at tightly linked SNPs, and the figure illustrates how eachhaplotype is modeled as a mosaic of these common combinations. The figure was produced by fitting our model to the HapMap data from 60
unrelated CEPH individuals (see the “Results” section) and then taking a single sample of cluster memberships and haplotypes from theirconditional distribution, given the genotype data and parameter estimates (appendix B). For brevity, haplotypes from only 10 individuals areshown.

(vs. quadratically for the PAC) and in that the model
can be applied directly to unphased genotype data, withunknown haplotypic phases integrated out analyticallyrather than via a time-consuming and tedious-to-implement Markov chain–Monte Carlo scheme, such as thatused by PHASE.
The price we pay for this computational convenienceis that our model is purely predictive; in common withthe block-based models mentioned above but in contrast
to the PAC model, our model does not attempt to directly relate observed genetic variation to underlying demographic or evolutionary processes, such as population
size or recombination. As such, it is not directly suited
to drawing inferences about these processes. However,
it is suited to two other applications that we consider
here: inferring unknown (“missing”) genotypes and inferring haplotypes from unphased genotype data. These
two applications are important for at least two reasons.
First, many methods for analyzing population data (e.g.,
methods that aim to draw inferences regarding demo630

graphic and evolutionary processes) struggle to deal withmissing genotypes or unphased data. This may be because dealing with these factors can create an impracticalcomputational burden or simply because necessary additional computer programming has not been done.
Therefore, applying those methods in practice often involves initially estimating missing genotypes and haplotypes by some other method. Second (as we expand
on in the discussion), the ability to accurately impute
missing genotypes and infer haplotypes has implications
for the development of methods for association-based
mapping of variants involved in common complex diseases. Comparisons with existing approaches for theseapplications suggest that our model has something tooffer, in terms of both speed and accuracy.
Material and Methods
ModelsTo introduce notation and the basic concepts underlying ourmodel, we begin by describing a simple cluster-based model

The American Journal of Human Genetics Volume 78 April 2006

www.ajhg.org

for haplotypes sampled from a population, in which each haplotype is assumed to have arisen from a single cluster. We thendescribe a modification to this model that allows cluster membership to change along each haplotype, to capture the fact
that, although sampled haplotypes exhibit clusterlike patterns,
these patterns tend to be local in nature. Finally, we describe
the extension of this model for haplotype data to unphased
genotype data through the assumption of Hardy-Weinberg
equilibrium (HWE) and a further extension that allows forcertain types of deviation from HWE.
Cluster model for haplotypes.—Suppose we observe n haplotypes, h p (h1 , … ,hn ), each comprised of data at M markers.
Let him denote the allele in the ith haplotype at marker m, sothat hi p (hi1 , … ,hiM). Throughout, we assume the markers arebiallelic SNPs, with alleles labeled 0 and 1 (arbitrarily) at eachsite, although the model is easily extended to multiallelicmarkers.
A simple cluster model for haplotypes can be developed asfollows. Assume that each sampled haplotype originates fromone of K clusters (labeled 1, … ,K). For simplicity, we initiallyassume K is known, but we will relax this assumption later.
Let z i denote the (unknown) cluster of origin for hi, and let
ak denote the relative frequency of cluster k, so that p(z i pkFa) p ak, where a p (a1 , … ,aK).
We assume that, given the cluster of origin of each haplotype, alleles observed at each marker are independent draws
from cluster-specific (and marker-specific) allele frequencies.
Thus, if vkm denotes the frequency of allele 1 in cluster k atmarker m, and v denotes the matrix of these values, then写 (clusters may seem slightly counterintuitive in this setting,where one expects to observe strong dependence among markers. The following observations may aid intuition. Each cluster
corresponds to a single row of the v matrix, which is a vectorof numbers in the range [0,1], with one number per marker.
(In fact, we imposed the constraint 0.01 ⭐ vkm ⭐ 0.99 on theelements of v, motivated by the idea that this might make themodel more robust to factors such as genotyping error.) For
estimates of v obtained from real data sets that we have examined, a moderate proportion (approximately two-thirds for
the HapMap data considered below) are very close to either
0 or 1. As a result, each row of v tends to look like a haplotype(a string of zeros and ones), but with occasional “fuzziness”
indicating uncertainty about the alleles at some positions. Thatis, each cluster can be thought of as representing a particular
combination of alleles at a subset of the markers, thus capturing strong dependence among these positions, and the assumption of independence can be thought of as relating todeviations from this base combination.
Local clustering of haplotypes.—Although sampled haplotypes certainly exhibit cluster-like patterns, these patterns tend
to be local in nature (fig. 1). To capture this, we replace theassumption that each haplotype originates from one of K
clusters with the assumption that each allele originates fromone of the clusters, and we use an HMM to model the factthat alleles at nearby markers are likely to arise from thesame cluster. Specifically, if z im denotes the cluster of origin for
him, we assume z i p (z i1 , … ,z iM) forms a Markov chain on{1, … ,K}, with initial-state probabilitiesp(z i1 p k) p ak1

M

p(hFz
i
i p k,v) p

mp1

him
vkm
1 ⫺ vkm)1⫺him .

and transition probabilities pm(k r k  ) given by

Since the clusters of origin are actually unknown, the probability of hi is obtained by summing equation (1) over all possible values of z i and weighting by their probabilities:

冘
冘 写 (
K

p(hFa,v)pip(z i p kFa)p(hFzii p k,v)

kp1

M

K

p

kp1

ak

mp1

him
vkm
1 ⫺ vkm)1⫺him .

(2)

Finally, specification of the model for h p (h1 , … ,hn ) is completed by assuming that h1 , … ,hn are independent and identically distributed from equation (2).
This simple model is essentially a haploid version of a modelthat has been widely used to capture the clustering that canoccur among individuals typed at unlinked (or loosely linked)
markers because of population structure in natural populations(Smouse et al. 1990; Rannala and Mountain 1997; Pritchard
et al. 2000). However, in those applications, the clusters represent “populations,” whereas here the clusters represent
groups of closely related haplotypes. The idea of using thismodel to capture clustering of haplotypes at tightly linked
markers seems to originate with Koivisto et al. (2003), whoused it to model data at tightly linked SNPs within a haplotype
block (see also Kimmel and Shamir 2005b).
The assumption of independence across markers withinwww.ajhg.org

(3)

(1)

pm(k r k  ): p p(z im p k Fz i(m⫺1) p k,a,r):pe⫺rm d m ⫹ (1 ⫺ e⫺rm d m)ak m ,1 ⫺ e⫺rm d m)ak m ,{(k p kk ( k(4)

for m p 2, … ,M, where dm is the physical distance between
markers m ⫺ 1 and m (assumed to be known) and wherer p (r2 , … ,rM) and a p (akm ) are parameters to be estimated.
This Markov chain is a discretized version of a continuous
Markov jump process, with jump rate rm per bp between markers m ⫺ 1 and m and with transition probabilitiesp(z im p k Fz i(m⫺1) p k, jump occurs) p ak m .

(5)

Informally, we think of rm as being related to the recombination
rate between m ⫺ 1 and m, although simulation results (notshown) suggest that generally there may be little correspondence between actual recombination rate and estimates of r.
If the physical distances between markers were not known,
then the compound parameter rm dm in equation (4) could be
replaced by a single parameter without loss of information.
Indeed, this is true even if the physical distances are known,
unless some constraint is placed on r (e.g., constraining all
rm to be equal). All results presented here were based on theunconstrained model and thus do not actually use the physicaldistances between markers. However, the algorithmic deriva-

The American Journal of Human Genetics Volume 78 April 2006

631

tions in the appendixes can be used for both the constrained
and the unconstrained models.
Given the cluster of origin of each allele, we assume, as
before, that the alleles are independent draws from the relevant
cluster allele frequencies, so

写

mp1

p(himFz im ,v) ,(1 ⫺ vk1 m )(1 ⫺ vk 2 m ),p vk1 m(1 ⫺ vk 2 m ) ⫹ vk 2 m(1 ⫺ vk1 m ),vk1 m vk 2 m ,{Since z i is unknown, the probability of hi is obtained by summing equation (6) over all possible values of z i and weighting
by their probabilities:zip(z iFa,r)p(hFzii ,v) ,

(7)2{(a2a )a, , kk p( kk
k11

k11

k 21

1

2

1

2

(8)

and transition probabilitiespm({k1 ,k 2 } r {k1 ,k 2 })p{

pm(k1 r k1 )pm(k 2 r k 2 ) ⫹ pm(k1 r k 2 )pm(k 2 r k1 ) ,k1 ( k 2 and k1 ( k 2pm(k1 r k1)pm(k 2 r k 2 ), otherwise,(9)

where pm(k r k  ) is defined in equation (4). These expressions
come from pairing two independent Markov chains with transition probabilities given in equation (4).
Given the clusters of origin, z•i, we again assume that the632

gim p 0
gim p 1 .
gim p 2(10)Note that, if some gim are missing, this is easily dealt with byreplacing the corresponding p(gimFz•im p {k1 ,k 2 },v) with any
positive constant (e.g., 1.0); this corresponds to the assumptionthat the genotypes are missing at random.
Since z•i is unknown, the probability of gi is obtained bysumming over all possible values:p(giFa,v,r) p

where p(z iFa,r) is determined by equations (3) and (4). Naive
computation of this sum would require a sum over KM possible
values for z i, but the Markov assumption for z i allows the sum
to be computed much more efficiently (with computational
cost increasing linearly with KM) using standard methods for
HMMs (e.g., Rabiner 1989).
Extension to genotype data.—Now suppose that, instead ofobserving haplotypes, we observe unphased genotype datag p (g1 , … ,gn ) on n diploid individuals. Let gim denote thegenotype at marker m in individual i, which we will code as
the sum of its alleles, so gim has the value 0, 1, or 2. One
approach to extending the haplotype-based model above tounphased genotype data is to assume that the two haplotypes
that make up each multilocus genotype are independent andidentically distributed from equation (7)—that is, to assume
HWE. Under this assumption, if z•im denotes the (unordered)pair of clusters from which genotype gim originates, thenz•i p (z•i1 , … ,z•iM) form a Markov chain with initial-state
probabilitiesp(z•i1 p {k1 ,k 2 }) p

p gimFz•im ,v) ,p( gimFz•im p {k1 ,k 2 },v)him
p(himFz im p k,v) p vkm(1 ⫺ vkm )1⫺him .

冘

mp1

(6)

where

p(hFa,v,r)
p
i

写(
M

p( giFz•i ,v) p
where

M

p(hFz
i
i ,v) palleles are independent draws from the relevant cluster allele
frequencies, so冘z•ip(z•iFa,r)p(giFz•i ,v) ,(11)

where p(z•iFa,r) is determined by equations (8) and (9). As before, methods for HMMs allow this sum to be computed efficiently (with computational cost increasing linearly withK 2M) (see appendix A).
This model (11) is reminiscent of the “linkage” model of
Falush et al. (2003), who modeled genotype data at loosely
linked markers in structured populations. One difference between their model and ours is that they allowed a (q in theirnotation) to vary among individuals but fixed a across markers, whereas we allow a to vary across markers but assume it
to be fixed across individuals. The reason for this differenceis that the interpretation of these parameters is very different
in the two applications. In the model of Falush et al. (2003),
this parameter controls each individual’s proportion of ancestry in each subpopulation (which would be expected to differ
across individuals), whereas here it controls the relative frequency of the common haplotypes (which would be expected
to differ in different genomic regions). Falush et al. (2003) also
restricted r to be constant, whereas we allow it to vary in eachmarker interval.
Modeling deviations from HWE and incorporation of subpopulation labels.—Although the assumption of HWE will not
hold exactly for real populations, previous studies have consistently suggested that models based on HWE can perform
well at haplotype inference and missing-data imputation, evenwhen there are clear and substantial deviations from HWE(e.g., Fallin and Schork 2000; Stephens and Scheet 2005). Nevertheless, we examined the potential benefits of modifying theabove model to deal with a certain type of deviation from
HWE. Specifically, we consider the situation in which the sampled individuals might be considered to be sampled from S
distinct subpopulations. (For example, the SeattleSNPs dataconsidered below consist of samples from 24 African Americans and 23 individuals of European descent, and we treat
these as separate subpopulations.) Let si 苸 1, … ,S denote thesubpopulation of origin for individual i (which is assumed to

The American Journal of Human Genetics Volume 78 April 2006

www.ajhg.org

be known here, although extension to the case where si is
unknown may also be interesting). Because of shared ancestry,different subpopulations may share some of their haplotype
structure, although haplotype frequencies and levels of LD
would be expected to differ. To capture this, we let the v parameters be shared across subpopulations but allowed the aand r parameters to vary among the subpopulations. Thus,a p (a(1), … ,a(S)) and r p (r (1), … ,r (S)), where (a(j),r (j)) denotethe parameters relating to subpopulation j, and写(np( gFa,v,r,s) p

ip1

p giFa(si),v,r (si) ) ,(12)

where s p (s1 , … ,sn ) and p( giFa(si),v,r (si) ) is as defined in equation (11).
Computation and Parameter Estimation
In this section, we outline the methods used to fit the models(11) and (12) and to perform two applications: missinggenotype imputation and haplotype inference.
Parameter estimation.—We use an expectation-maximization (EM) algorithm (Dempster et al. 1977) to estimate theparameters n p (v,a,r) of our model (see appendix C for details). The computational complexity of the algorithm isO(nMK 2 ) and, in particular, is linear in the number of sampledindividuals and markers, which allows it to be fitted to large
data sets.
As with any EM algorithm, our algorithm will typically finda local maximum of the likelihood function L(n; g ) p p( gFn).
For realistic data sets, this likelihood surface will have manydifferent local maxima, and thus different starting points for
the EM algorithm will typically lead to different parameter
estimates. A standard approach to dealing with this problem
is to first apply the algorithm T times from T different starting
points, obtaining T estimates nˆ 1 , … ,nˆ T , and then select whichever of these estimates gives the highest value for the likelihood.
However, because our focus here is on using the model for
prediction and not on parameter estimation itself, it is notnecessary to settle on a single estimate for the parameters, and,in our tests, we found that we were able to obtain more accurate predictions by combining results across T estimates, as
described below. We also found that, using this strategy ofcombining across estimates, reliable performance can be obtained with relatively few iterations of the EM algorithm per
starting point—probably far fewer than would be required bymost methods of monitoring convergence. For the results presented here, we typically used T p 20 starts of the EM algorithm, with up to 25 iterations per start, although results ofexperiments suggest that these values could be reduced without
sacrificing accuracy (results not shown). For each initialization of the EM algorithm, we set r p 0.00001, chose vkm
to be independent and identically distributed uniform on[0.01,0.99], and let a7m ∼ Dirichlet(1, … ,1). Our methods appeared somewhat robust to deviations from these choices (e.g.,setting akm p 1/K, for all m and k, produced similar mean
accuracy), although we did not undertake a detailed study.
Missing-genotype imputation.—For any genotype gim that is
unobserved (“missing”), it is straightforward to compute thewww.ajhg.org

probability that gim p x (x p 0,1,2), given all observed genotypes g and parameter values n, by use of冘冘 (Kp( gim p xFg,n) pK

k1p1 k 2pk1

p gim p xFz•im p {k1 ,k 2 },n)# p(z•im p {k1 ,k 2 }Fgi ,n) .
The first term in this sum is given by equation (10), and thesecond term is the conditional distribution of the hidden variables in the HMM, which can be obtained using standard
methods for HMMs (appendix A).
A natural point estimate for gim is then obtained by choosing
the value of x that maximizes this expression. As noted above,
we have found it helpful to combine results over several setsof parameter estimates nˆ 1 , … ,nˆ T , obtained from T different applications of the EM algorithm using different starting points.
Specifically, we used the estimate

冘Tgˆ im p argmaxx苸{0,1,2}

1p( gim p xFg,nˆ t ) .
T tp1

This method imputes genotypes marginally and provides a“best guess” for each genotype. It is also straightforward tosample from the joint distribution of the missing genotypes
given observed data—for example, by sampling from the conditional distribution of the haplotypes for all individuals, as
described below.
Haplotype inference.—We consider two aspects of the haplotype inference problem: (1) sampling the pairs of haplotypesof all individuals from their joint distribution given the unphased genotype data—this provides a useful way to assess or
account for uncertainty in haplotype estimates—and (2) constructing point estimates of the haplotypes carried by eachindividual—this is how many haplotype-reconstruction methods are used in practice (as a prelude to subsequent analysisof the estimated haplotypes). It also provides a convenient basis
for comparison with other haplotype reconstruction methods.
Sampling haplotypes from their joint distribution.—We will
use the term “diplotype” to refer to a pair of haplotypes thatcomprise the genetic data for an individual. Let di denote thediplotype for individual i, and d p (d1 , … ,dn ). Conditional ona particular parameter value n, the diplotypes of different individuals are independent (i.e., p(dFg,n) p 写 i p(diFgi ,n)), andthus one can sample from p(dFg,n) by sampling independently
from p(diFgi ,n) for each i. A method for doing this is describedin appendix B.
To combine results from T initiations of the EM algorithm,
we obtain a sample d˜ of size T # B as the union of B independent samples from each of p(dFg,nˆ 1 ), … ,p(dFg,nˆ T ). When
constructing dˆ SW
below, we used T p 20 and B p 50; foridˆ INi , we used T p 20 and B p 200.
Point estimation.—We have implemented two different
methods for producing point estimates, dˆ i, of the diplotypesof individual i. Both are obtained by first creating a sampled˜ i of diplotypes for individual i, as described above. From thissample, we define the following estimates.
1. dˆ INi , which is the diplotype that appears most often in

The American Journal of Human Genetics Volume 78 April 2006633

d˜ i. This estimate is motivated by an attempt to maximize
the probability that the whole diplotype is correct or,equivalently, to minimize the “individual error rate”—
that is, the proportion of individuals whose haplotypes
are not determined by their genotype data (Stephens etal. 2001).
2. dˆ SWi , which is constructed as follows. Starting at one endof the genetic region, we move through the heterozygous
sites, phasing each site relative to the previous heterozygous site by selecting the two-site diplotype that occurs
most frequently (at that pair of sites) in d˜ i. This estimate
is motivated by an attempt to minimize the “switch error”—that is, the proportion of heterozygous sites thatare phased incorrectly relative to the previous heterozygous site (this is 1 minus the switch accuracy of Lin etal. [2002]).
Note that, for data sets containing a large number of markers, individuals may have a very large number of plausible
diplotype configurations, none of which are overwhelmingly
more probable than the others. In such cases, the most probable diplotype configuration is both difficult to reliably identify(requiring a very large sample d˜ i) and not especially interesting(in that it is very unlikely to be correct). Therefore, for largescale studies, we would tend to favor the use of dˆ SWover dˆ iIN.
i
Selecting K.—Selection of K is essentially a model-selection
problem, which is, in general, a tricky statistical problem. We
found that standard approaches to model selection, such as
Akaike information criterion (Akaike 1973) and Bayesian information criterion (Schwarz 1978), do not work well here
(selecting a K that is too small), probably because the asymptotics on which they are based do not apply. We therefore used
the following cross-validation approach to select K. For eachdata set, we masked (i.e., made missing) ∼15% of the genotypes at random (note that, in our tests of imputation accuracy
below, the data set to be analyzed would already have hadsome genotypes masked; we did not use these masked genotypes when selecting K). Then, for a range of values of K (weconsidered K p 4, 6, 8, 10, and 12; the upper limit reflects
our desire to keep the computational burden low), we used
our model to estimate the masked genotypes, as describedabove, comparing these estimates with the true genotypes. We
selected the K that maximized the number of correctly estimated genotypes.
This procedure is relatively computationally intensive, since
it involves fitting the model for several values of K, and thecomputational cost increases with K 2. To reduce the computational burden, we used only a small number of starts (as few
as three) for the EM algorithm when performing the crossvalidation. In addition, in the results we present here, we sometimes chose a single K for several data sets on a common setof individuals. For example, for the SeattleSNPs data considered later, when analyzing data on a particular set of individuals, we selected a single K for all genes, on the basis of applying the cross-validation approach to a subset of the genes.
Preliminary comparisons suggested that this approach gave
similar average accuracy to the more computationally intensive
approach of choosing K separately for each gene (data notshown). For missing-genotype estimation results from the634

HapMap project, we selected K by using only a small portionof the data, since the data sets were so large.
Because we found performance to be relatively robust to arange of values of K, in cases where the computational costof this strategy becomes inconvenient an alternative approach
would be to simply select a fixed value of K (K p 8 seemed
to perform reasonably well across a range of scenarios in our
tests). It would also be possible, and perhaps fruitful, to combine results across parameter estimates obtained using different
values of K rather than selecting a single K value.

Results
The methods described above for imputing missinggenotype data and estimating haplotypes are implemented in a software package called “fastPHASE,”
which is available for download from the Stephens Lab
Web site. Here, we compare performance of these methods with several other available methods, including
PHASE version 2.1.1 (Stephens et al. 2001; Stephens andDonnelly 2003; Stephens and Scheet 2005), GERBIL
versions 1.0 and 1.1 (Kimmel and Shamir 2005b), andHaploBlock version 1.2 (Greenspan and Geiger 2004).
The models underlying both GERBIL and HaploBlockbear some similarity to our model, being based on theidea of clusters of haplotypes, but in these models, cluster membership is allowed to change only at certain
points in the genome (“block-boundaries”), which are
estimated from the data. The model underlying PHASE
is based on the PAC model of Li and Stephens (2003),
which shares the flexibility of the model we present here
but is considerably more costly to compute. In comparisons elsewhere (Stephens and Scheet 2005), we foundthat PHASE outperformed several other methods in accuracy of both missing-data imputation and haplotype
estimation, but GERBIL and HaploBlock were not included in those comparisons. For haplotype estimation,
Kimmel and Shamir (2005b) found that GERBIL performed better than HaploBlock but slightly less well than
PHASE. However, they did not examine accuracy in
missing-data imputation, and, as far as we are aware,
our comparisons provide the first published assessmentof GERBIL and HaploBlock for this task. PHASE andGERBIL were run with their default settings, and HaploBlock was run with the ⫺W option (which produces
estimates by combining over multiple solutions, similar
to our strategy of combining estimates from multiple
runs of the EM algorithm; this option seemed to give
more accurate results but took longer to run than thealternative ⫺F option).
Missing-Data Imputation
We examined accuracy of methods for missing-data
imputation with both complete sequence data in multiple candidate genes (from the SeattleSNPs Variation

The American Journal of Human Genetics Volume 78 April 2006

www.ajhg.org

Discovery Resource) and genomewide dense SNP data
(1 SNP every 2–3 kb across entire chromosomes) fromthe International HapMap Project (International HapMap Consortium 2005).
SeattleSNPs.—We analyzed the same data from the
SeattleSNPs resequencing project as were analyzed elsewhere (Stephens and Scheet 2005); we considered sequence data on 50 autosomal genes, containing 15–230
SNPs, from 24 African Americans and 23 individuals ofEuropean descent. For each gene, we masked ∼5% of
the individual genotypes at random. We then used various methods to estimate the masked genotypes and assessed performance by computing the error rate as the
proportion of masked genotypes that were not estimatedcorrectly.
As elsewhere (Stephens and Scheet 2005), we usedeach method to analyze the data in two ways: (1) analyzing the African American and European-descentsamples separately and (2) analyzing the combined sample of 47 individuals together. For both of these approaches, we computed both the overall (total) error
rates and error rates stratified by subpopulation (table1). All four methods performed well—better than mostof the methods considered elsewhere [Stephens andScheet 2005])—with fastPHASE yielding the lowest
error rate (although differences among fastPHASE,
PHASE, and HaploBlock were not statistically significant). As in a previous study (Stephens and Scheet 2005),
our analysis of data from the African American andEuropean-descent samples combined seemed to give asmall decrease in error rate compared with analysis of
the samples separately, which illustrates the relative robustness of the methods to deviations from HWE.

We also assessed robustness of the results from fastPHASE to the number of clusters, K (table 1). For these
data, results are relatively robust across the range of K
considered here, at least provided that K is sufficiently
large (say, at least 6). Error rates generally tended todecrease as K increased from 4 to 12, although, forsufficiently large K, we would expect the error rates
to increase, and the rather small differences between
K p 8, 10, and 12 suggest that larger values of K would
not produce a substantial improvement in performance.
To examine whether accuracy could be improved bytaking account of the subpopulation of each sample, we
also analyzed the combined-sample data by using the“separate a, r” model from (12). This model produced
a very small improvement in error rate and also appeared
to make results more robust to the choice of K (lastcolumn of table 1).
The results for fastPHASE in table 1 were all obtainedby averaging over parameter estimates from T p 20
starts of the EM algorithm. We found that estimates
based on only the parameter value (among these 20 estimates) that maximized the likelihood were consistently
less accurate. For example, in the combined analysis with
K p 6, this latter (maximization) strategy gave an error
rate of 0.051, compared with 0.045 for the averaging
strategy.
CEPH HapMap data.—We obtained HapMap data forchromosomes 7 (41,018 SNPs across 159 Mb) and 22(15,532 SNPs across 35 Mb) from parents in 30 CEPH
trios (60 unrelated individuals) from a phase I data freeze(March 2005; International HapMap Project). We produced data sets with missing genotypes by masking 10%
and 25% of genotypes at random and computed error

Table 1
Error Rates for Estimation of Missing Genotypes for SeattleSNPs
ANALYZED SEPARATELYMETHOD
fastPHASE
PHASE version 2.1.1
GERBIL version 1.1
HaploBlock
fastPHASE with fixed K:
Kp4
Kp6
Kp8
K p 10K p 12ANALYZED COMBINED

AA
Error

ED
Error

Totala

AA
Error

ED
Error

TotalaTotal(Separate a,r)a

.053
.058
.067
.053

.024
.030
.030
.026

.039
.044
.048
.040

.051
.052
.063
.051

.022
.024
.026
.026

.037
.038
.045
.039

.035
NA
NA
NA

.066
.059
.058
.054
.053

.029
.024
.026
.027
.026

.048
.042
.042
.041
.039

.072
.061
.054
.051
.051

.032
.027
.024
.023
.022

.052
.045
.039
.037
.037

.046
.039
.036
.036
.035

NOTE.—Error rates are based on estimation of 9,479 missing genotypes. The best-performing
method in each column is in bold italics. The differences among fastPHASE, PHASE, andHaploBlock are not statistically significant; the differences between these three methods and GERBIL
are significant (P ! .007) on the basis of bootstrap resampling of the 50 genes. AA p African
American sample; ED p European-descent sample.
aTotal error rate for combined sample. The “Separate a,r ” total gives results from use of themodel in section 2.1.4.
www.ajhg.org

The American Journal of Human Genetics Volume 78 April 2006

635

Table 2
Error Rates for Estimation of Missing Genotypes for CEPH HapMap
Data
ERROR RATE FOR CHROMOSOME
7DATA ANALYZEDAND METHOD
Whole chromosome:
fastPHASE
fastPHASE, 1 start
Separate 150-SNP data sets:
fastPHASE
PHASE version 2
GERBIL version 1.1

22

10%
25%
10%
25%
Masked Masked Masked Masked
.034
.046

.041
.057

.033
.045

.039
.056

.036
…
.056

.044
…
.077

.035
.038
.054

.042
.049
.073

NOTE.—For each chromosome, 10% and 25% of the data weremasked, resulting in 242,481 and 606,985 missing genotypes for chromosome 7 and 93,476 and 232,731 missing genotypes for chromosome
22. Results for fastPHASE were obtained from 20 random starts ofthe EM algorithm, except for the “1 start” case, for which results wereobtained from a single random start.

rates for different methods of estimating these genotypes.
Using fastPHASE, we were able to analyze the complete
data for each chromosome. However, other methods
struggled computationally; thus, to allow comparisonswith these methods, we also split the data sets into nonoverlapping segments, each containing 150 consecutive
SNPs, and analyzed each segment separately. Even so,
HaploBlock failed to finish computing results for several
data sets after weeks of running on multiple machines
and thus was omitted from the comparisons. Because of
the amount of computation required, we applied PHASE
to only the chromosome 22 data sets.
Results are given in table 2. For these data, fastPHASE
again produced a lower error rate than those of the othermethods, with very slightly better accuracy when all data
were analyzed simultaneously instead of in 150-SNP segments. PHASE performs about as well as fastPHASE,
and the improvement of these methods over GERBIL
was more substantial for these data than for the SeattleSNPs data. The main differences between the data sets
are that (i) the HapMap SNPs tend to have a higher
minor-allele frequency, because of the way in which these
SNPs are ascertained, and (ii) the HapMap SNPs are
more widely spaced and thus would be expected to exhibit less LD. Both these factors presumably contributeto the slightly worse performance of all methods for the
HapMap data compared with the European-descentsample of the SeattleSNPs data.
The high accuracy with which genotypes could be estimated with even 25% missing data prompted us to
examine in more detail the relationship between accuracy and rates of missingness (table 3). Although thismissing-at-random pattern is not a realistic assumption
for missing data observed in real studies, the fact that
accuracy remains high (193%) even with 50% of the636

genotypes deleted illustrates both the effectiveness of themethodology and the strong correlations that exist
among SNPs at this density.
Haplotype Inference
X-chromosome data.—We assessed the accuracy of
haplotype estimates by using the X-chromosome datafrom Lin et al. (2002), also analyzed elsewhere by Stephens and Donnelly (2003) and by us (Stephens andScheet 2005). The data consist of X-chromosome haplotypes derived from 40 unrelated males. The haplotypes
comprise eight regions, which range in length from 87
to 327 kb and contain 45–165 SNPs. For each of theeight genes, we used the same 100 data sets as elsewhere(Stephens and Scheet 2005), each consisting of 20
pseudo-individuals, created by randomly pairing the 40
chromosomes.
As in previous comparisons of this type, haplotype
estimates obtained using different methods were scoredusing two error rates: the individual error (the proportion of ambiguous individuals whose haplotypes are notcompletely correct) and the switch error (the proportionof heterozygote genotypes that are not correctly phased
relative to the previous heterozygote genotype). In computing these proportions, we summed both numerator
and denominator over all 8 # 100 data sets (which isslightly different from computing an error rate for eachof the eight genes and then averaging these rates, as intable 2 in an earlier publication [Stephens and Scheet2005]). In computing these scores for each individual,
we ignored sites where one or both alleles were missing.
Results are given in table 4. For data sets like these,
which contain at least a moderate number of SNPs, estimating any individual’s haplotypes completely corTable 3
Error Rates for Estimationof Missing Genotypes with
fastPHASE for CEPH HapMap
Data, Chromosome 22
Missing Data
(%)
10
20
30
40
50
60
70
80
90

fastPHASE Error
.033
.037
.042
.051
.064
.089
.137
.227
.358

NOTE.—We masked 10%–90%of the data at random, which resulted in 93,476–837,853 missing
genotypes, and applied fastPHASE
to the entire chromosome 22.

The American Journal of Human Genetics Volume 78 April 2006

www.ajhg.org

Table 4
Individual and Switch Error Rates for Haplotype
Estimates Produced by Different Methods
for the X-Chromosome Data
MethodfastPHASE (dˆ SW)
fastPHASE (dˆ IN)
PHASE version 2.1.1
GERBIL version 1.0
HaploBlock
fastPHASE (dˆ SW) with fixed K:
Kp4
Kp6
Kp8
K p 10K p 12

Individual
Error

Switch
Error

.654
.641
.624
.660
.702

.111
.116
.113
.118
.122

.654
.642
.645
.650
.657

.111
.109
.110
.111
.113

NOTE.—The best-performing method in each errorrate column is in bold italics.

rectly is difficult. As a result, individual error rates ofall methods are high, and it could be argued that theswitch error is a more meaningful measure of performance. However, the qualitative conclusions are thesame based on either error rate. Consistent with the results of Kimmel and Shamir (2005b), PHASE slightly
outperformed GERBIL, which slightly outperformed
HaploBlock; fastPHASE produced an individual error
rate between that of PHASE and that of GERBIL andproduced the lowest switch error rate (although the difference from the switch error rate of PHASE is small
and not statistically significant). As one might hope, thepoint estimate dˆ SW from fastPHASE, which aims to minimize the switch error, produces a lower switch error
than does dˆ IN, which aims to minimize the individual
error. Conversely, dˆ IN produces a lower individual error
rate. Results from fastPHASE are again robust to a rangeof values of K.
For these data, averaging results over multiple parameter estimates obtained from multiple starts of the EM
algorithm turned out to be particularly important for
obtaining good performance. For example, for K p 4,
the point estimate dˆ SW obtained using only the single setof parameter estimates that give the largest likelihood
yielded error rates that were worse than those of any of
the other methods considered here (individual error
0.716; switch error 0.134).
It is also notable that, for these data, even the simplecluster model (which can be obtained as a special caseof our model by setting rm p 0 for all m) performs similarly to GERBIL, with an individual error rate of 0.648
and a switch error rate of 0.119. This is perhaps becauseof lower levels of historical recombination for these Xchromosome genes; thus, these data may not provide thebest guide to performance of haplotype-inference methods on autosomal data.
www.ajhg.org

HapMap data.—Marchini et al. (2006) compared several methods for haplotype inference on samples of both
unrelated and related individuals (trios of parents andchild). We consider here the data they used from unrelated individuals, which consist of three simulated data
sets and one real data set. The simulated data consist of
three “trials” (SU1, SU2, and SU4 of Marchini et al.
[2006]), which were simulated using coalescent methods
with the following conditions: for trial 1, constant-sized
population and constant recombination; for trial 2, constant-sized population with variable recombination; andfor trial 4, population demography approximating thatof a European population and variable recombination,with 2% of genotypes masked. The real data consist of
unrelated CEPH samples (60 parents from 30 trios) fromthe HapMap project (International HapMap Consortium 2005), for which the real haplotypes were determined using the trio data under the assumption of norecombination from parents to offspring. (Under thisassumption, the trio data determine phase at a large
proportion of sites; the remaining ambiguous sites wereignored in scoring methods.)
We applied fastPHASE to the unphased genotypes
from these data and sent the estimated haplotypes(dˆ SW) to J. Marchini, who independently scored the results. Table 5 compares the results from fastPHASE withthose of other methods in the original comparison. The
results from fastPHASE were consistently worse thanthose of PHASE (and those of wphase, for the simulated data) and were consistently better than those of
the other methods, HAP and HAP2. Encouragingly forour model, the performance difference between fastPHASE and PHASE was the smallest for the real data,with switch errors of 0.055 and 0.051, respectively.
Calibration of Probability Calculations
All the comparisons above are based on assessing theaccuracy of point estimates of estimated genotypes orinferred haplotypes. However, by use of our model, it isalso quick and easy to compute probabilities for eachmissing genotype and to produce samples from the conditional distribution of haplotypes, given the unphased
genotypes. These can be used to take account of uncertainty in estimated genotypes and/or haplotypes indownstream analysis—for example, by performing thesubsequent analysis on multiple sampled haplotype reconstructions to check for robustness of conclusions or,more formally, by using Bayesian statistical methods.
However, to justify this strategy, the model should ideally produce approximately calibrated predictions. For
example, of genotypes assessed to have a probability of
0.9 of being correct, ∼90% should actually be correct.
We therefore assessed the calibration of predictions
from our model, for both genotype imputation and haplotype inference. For each imputed genotype, we com-

The American Journal of Human Genetics Volume 78 April 2006

637

Table 5
Accuracy of Haplotype-Inference Methods for Simulated Data (Trials 1, 2, and 4) and Real Data from HapMap
TRIAL 1METHOD
fastPHASE (dˆ SW)
PHASE version 2wphaseHAPHAP2

TRIAL 2

TRIAL 4REAL DATA

Individual
Error

SwitchError

Individual
Error

SwitchError

Individual
Error

SwitchError

Missinga

Individual
Error

Switch
Error

.653
.355
.480
.886
.735

.045
.024
.037
.065
.069

.887
.404
.521
.971
.990

.069
.022
.037
.098
.151

.760
.620
.680
.906
.871

.058
.053
.066
.074
.087

.091
.075
.090
.116
.150

.879
.815
…
.919
.901

.055
.051
…
.066
.078

NOTE.—Results for wphase (N. Patterson, personal communication), HAP (Halperin and Eskin 2004), and HAP2 (Lin etal. 2002) were obtained by Marchini et al. (2006).
a
Genotype-imputation error rate.

puted the probability, p, under our model, that the imputed genotype was correct. We then grouped imputed
genotypes into bins, according to their value for p, and,for each bin, we compared the average value of p withthe proportion of genotypes that were actually correct.
For haplotype reconstruction, we examined the calibration of a sample of diplotype configurations, d˜ , by looking at potential switch errors—that is, by examining, ineach individual, the phase of each heterozygous site relative to the previous heterozygous site. For each suchpair of heterozygous sites, we computed the proportion,q, of the configurations in d˜ that contained the morecommon of the two possible phasings. We then grouped
these site pairs into bins, according to their value for q,
and, for each bin, we compared the average value of qwith the proportion of phasings that were actuallycorrect.
The results (fig. 2) show that our model is reasonably
well calibrated, but slightly conservative, for both tasks.
For example, of genotypes assessed a 90% chance ofbeing correct, roughly 96% were actually correct in boththe SeattleSNPs and HapMap data sets. One curious
feature of the results is the slight drop in accuracy forthe highest confidence bin (corresponding to 98% predicted probability of being correct, a value that resultsfrom the limits of 0.99 and 0.01 we imposed on elementsof v) in the SeattleSNPs data. Closer inspection reveals
that this is because of errors in imputing genotypes ofmasked heterozygotes at singleton SNPs (i.e., SNPs
where only one copy of the minor allele is present). Whensuch genotypes are masked, our model very confidently
but wrongly assesses them to be homozygous for the
only observed allele. This could be viewed as an artifact
due to the fact that the data contain only polymorphic
SNPs and that our model does not condition on themarkers being polymorphic.
Differences in Computational Requirements
The relative computation times of the different methods we consider here vary across data sets and dependon the way in which the methods were applied (e.g.,638

how many iterations were used). As we applied themethods here, fastPHASE and GERBIL require similar
amounts of computational resources, and both are considerably faster than PHASE and HaploBlock. Computation times for all results are summarized in table 6.
The increased speed of fastPHASE, compared with
PHASE, would be greater for samples with a larger number of individuals.
Discussion
We have presented a model for genetic variation among
unrelated individuals, which is computationally tractable for large-scale sets of unphased SNP genotype data.
Each data set considered here, including whole-chromosome data from phase I of the HapMap project, took
!10 h to analyze (table 6). To test the feasibility of applying our model to even larger data sets, we created adata set containing 3,150 individuals typed at 290,630
SNPs by concatenating 15 copies of the chromosome 2
genotypes of 210 unrelated individuals from phase II of
the HapMap project. Our software required 97 h (on asingle 3-GHz Xeon processor with 8 GB of RAM) to fit
the model once to these data. In addition to its computational convenience, the model is also flexible andcan capture both the sudden decline of LD that might
be expected across a recombination hotspot and themore gradual decline of LD with distance. Indeed, forthe task of imputing missing genotypes, our model performed better than any other method we considered. For
inference of haplotypes, it performed slightly less well(on the real-data comparisons) than the best of the methods we considered, PHASE version 2, but at a fractionof the computational cost. It seems slightly puzzling that,at least for CEPH HapMap data, fastPHASE appears tooutperform PHASE for missing-genotype imputation(table 2) but to perform less well than PHASE for haplotype inference (last column of table 5). We have no
good explanation for why this should be (the differences,though small, appear statistically significant because thedata sets are large).

The American Journal of Human Genetics Volume 78 April 2006

www.ajhg.org

Figure 2

Calibration of our model for predicting uncertainty in inferred genotypes and haplotypes. Points (triangles) represent probabilities
obtained by averaging over the 20 runs of the EM algorithm, as described in the text.

There are at least two factors that may contribute to
the improved performance of our model compared withthe block-based model underlying GERBIL (which isperhaps the most similar to our model among those considered here). The first is the increased flexibility of ourmodel, in that cluster membership can change continuously along the genome and not only across blockboundaries. The second is the fact that we average resultsover multiple fits of our model to improve accuracy,
whereas, as far as we are aware, GERBIL does not. Sincea similar averaging strategy also seems to improve performance of HaploBlock, it may be interesting to examine whether averaging could be used to improve performance of GERBIL and, indeed, of other methods.
Our strategy of averaging results across applicationsof the EM algorithm is slightly unusual. However, averaging results across models has often been observed
to produce improved predictive performance, datingback at least to Bates and Granger (1969) and becoming
particularly popular recently with the increased use ofBayesian methods and the development of methods suchas bagging (Breiman 1996) and boosting (Freund andSchapire 1996). Indeed, there are good theoretical reasons to expect averaging across multiple runs to produce
better performance than using results from a single runof the EM algorithm; the averaging reduces the varianceof predictions while leaving any bias unchanged. Of
course, this alone does not fully explain our empiricalfinding that averaging across multiple runs produces better performance than making predictions based on the
run with the highest likelihood. However, it is worth
noting that this latter, more standard approach is also
theoretically dubious in this setting, because the asymptotic theory that underlies maximum-likelihood estimation may not be applicable to most data sets, becauseof the large number of parameters in the model.
The averaging scheme that we use involves formationwww.ajhg.org

of equally weighted averages across multiple model fits.
One might expect that a more sophisticated averaging
scheme might further improve performance. We tried
weighting results according to the likelihood of the corresponding parameter estimates, but this typically produced a worse performance, very similar to that of selecting the single parameter values with the largest
likelihood, because one or two of the likelihoods are
typically much greater than all the others, and so theaverage is dominated by these parameter values. A
Bayesian version of our model would certainly be possible and would provide a more coherent way to averageover parameter values. However, a fully Bayesian implementation would presumably greatly increase computational cost and seems unlikely to lead to substantial
improvements in prediction accuracies.
Another novel aspect of our work is the extension of
our model to deal with samples from multiple subpopulations. For the data we considered here, involving 24
African Americans and 23 individuals of European descent, we found that a slight but consistent improvementin performance could be obtained by taking account of
the subpopulation of origin of each individual. We havealso found a similar slight improvement when analyzing
data from the four analysis panels of the HapMap project (results not shown). It seems likely that the gains ofusing this model will depend on sample size and theamount of divergence among the subpopulations, although we have not studied this dependence. Given thesimilarities between our model and those of Pritchard
et al. (2000) and Falush et al. (2003), it seems natural
to consider extending our model to the case in which
the subpopulation(s) of origin of each individual is unknown, effectively producing a method for clustering
individuals that can deal with sets of tightly linked markers. This might be especially helpful for investigating
population structure in organisms with small genomes

The American Journal of Human Genetics Volume 78 April 2006

639

Table 6
Comparison of Computation Times (in Hours) for Different Methods
MISSING-GENOTYPE ESTIMATION

SeattleSNPs(10%)

METHOD
fastPHASE
PHASE version 2
GERBIL
HaploBlock

1.3
29.3
.6
160.1

10%

25%

10%

25%

HAPLOTYPE
INFERENCE FOR
X CHROMOSOME
(7%)

9
NA
29
NA

8
NA
32
NA

5
323
10
NA

5
720
8
NA

2
151
2
265

HapMap
Chromosome 7

HapMap
Chromosome 22

NOTE.—Each method was applied on a 3-GHz Xeon processor with 1 GB of memory. The
SeattleSNPs column is for analyses of all 47 individuals together. For missing-genotype estimation,
fastPHASE was applied to the data without inference of haplotypes. Calculation times for HaploBlock with the SeattleSNPs data and for PHASE version 2 with HapMap chromosome 22 data
are based on extrapolation from a subset of the data sets. Approximate percentages of missing
data are in parentheses. Bold italics indicate the fastest method(s) for each data set.

and/or little recombination, for which finding large setsof unlinked (or loosely linked) markers will be more
problematic than in humans.
Our model has a very large number of parameters,
and, for realistic data sets, we would not expect all parameters to be well estimated. It is possible to decrease
the number of parameters by imposing constraints on
a, v, and/or r. For example, we tried constraining r anda to be constant across the genome. We also experimented with ad hoc smoothing schemes to encourage aand r to vary smoothly across the genome. However, incomparisons (not shown), these strategies did not leadto an improved performance in our applications, suggesting, perhaps surprisingly, that the large number ofparameters does not greatly diminish the predictive
power of the model.
While preparing this article, we became aware of independent work by Rastas et al. (2005) and Kimmel andShamir (2005a) on models similar to the one we present
here. In particular, these authors also pursue the strategy
of modeling cluster membership along the chromosome
by using an HMM. The main difference between theirmodels and ours is that, in our model, when a jump incluster membership occurs, the distribution of the newcluster does not depend on the current state (that is,
the right-hand side of eq. [5] does not depend on k),
whereas, in their models, it does. This results in ourmodel being less computationally complex (the procedures we describe here have computational complexity
O(K 2 ) compared with O(K 3) for the EM algorithm ofRastas et al. [2005]). A less important difference is that
we allow for the possibility of allowing the jump probability between markers in the HMM to depend on thephysical distances between markers—for example, byconstraining the rm to be equal in equation (4). However,as noted earlier, we did not find that this improved per640

formance over the unconstrained model, in which jumpprobabilities do not depend on marker spacing.
In addition to this difference between the models,there are also several further differences, both in theapplications considered (Rastas et al. [2005] apply theirmodel to haplotype inference, and Kimmel and Shamir[2005a] consider disease mapping) and in the way theseapplications are tackled. For example, Rastas et al.
(2005) estimate haplotypes for each individual by firstusing the Viterbi algorithm to find the most-probablecluster memberships and then obtaining the most probable diplotype, given these cluster memberships (see also
Kimmel and Shamir 2005b). Note that this will not, ingeneral, find the most probable diplotype for each individual, but the empirical results suggest that it is areasonable procedure. In contrast, we use Monte Carlosimulation to find haplotype estimates that attempt tominimize two different error measures. In addition, Rastas et al. (2005) and Kimmel and Shamir (2005a) use
much more complex initialization strategies for their EM
algorithms than we do here. It seems possible that ourapproach of averaging across runs of the EM algorithm
obviates the need for more-complex initialization procedures, and so it is unclear whether combining theseapproaches would improve performance.
The two applications we considered here—namely,
missing-data imputation and haplotype inference—areof direct interest in themselves, particularly as a prelude
to subsequent analyses using methods that cannot deal
with missing or unphased genotypes. We believe thatthey are also of indirect interest for the role they may
eventually play in the development of powerful association-based methods for mapping disease genes (see also
Kimmel and Shamir 2005a). Moreover, although the useof haplotypes in such methods has received more attention, we would argue that, in the longer term, methods

The American Journal of Human Genetics Volume 78 April 2006

www.ajhg.org

to accurately estimate missing and untyped genotypes
may prove to be more important. One of the primary
motivations for haplotype-based mapping methods isthat, if an untyped variant is responsible for affecting
phenotype and if this untyped variant (and therefore thephenotype) is strongly associated with a particular haplotype but not with any individual SNP, then haplotypebased tests might succeed in detecting a significant effectwhere tests based on individual SNPs fail. That is, they
are based on the idea that haplotypes may be better
predictors of untyped genotypes than are individual
SNPs. Indeed, some existing methods are explicitly basedon the idea of using haplotypes to predict genotypes at
one or more causal SNPs (e.g., Zöllner and Pritchard2005). However, in these methods, the use of haplotypes
is merely a convenient intermediate step in predicting
the untyped variants; if one had a “black box” that could
predict genotypes of untyped variants directly from unphased genotype data, then this could similarly form thebasis of methods for association mapping, bypassing theneed for haplotype estimates (see also the work of Chapman et al. [2003], for relevant discussion and the caveatthat haplotype estimation methods may still be usefulfor studying certain types of interactions among closely

linked causal SNPs). For this application, a limitationof our model is that, because it has a parameter vector(the columns of v) that must be estimated for each SNP,
it is suited to the imputation of genotypes only at SNPs
where several (and preferably many) individuals havebeen genotyped and not at sites where no individuals
have been genotyped. Although the HapMap and largescale resequencing studies, such as the SeattleSNPs project, helpfully provide reference panels of individuals
typed at large numbers of SNPs, many SNPs potentiallyinvolved in human disease are currently untyped in suchpanels, and we are working to extend our model to allow
it to impute genotypes at such SNPs. Ultimately, ourhope is that the kinds of models and methods examined
here will form the foundation for new and more-effective
statistical methods for analyzing whole-genome association studies.

Acknowledgments
We thank C. Geyer, for helpful discussions, and two anonymous referees, for comments on the submitted version. This
research was supported by National Institutes of Health grant
1RO1HG/LM02585-01.

Appendix A
Forward and Backward Algorithms
Two fundamental quantities associated with HMM computations are the so-called forward and backward probabilities:
fni (m,{k1 ,k 2 }) : p p( gi1 , … ,gim ,z•im p {k1 ,k 2 }Fn)

I{k1(k 2 }

(12)

and

bni (m,{k1 ,k 2 }) : p p( gi(m⫹1) , … ,giMFz•im p {k1 ,k 2 },n) ,where I{A} is equal to 1, if A is true, and zero otherwise. (The factor (1/2)I{k (k } is not usually included in the definitionof f, but we include it here for later notational convenience.) Although computation of these quantities for HMMs
via recursive formulas is standard, we give these formulas here because some care is needed to ensure thesecomputations have complexity O(K 2 ), rather than O(K 4 ).
The forward calculation is given by12[

fni (m ⫹ 1,{k1 ,k 2 }) p p( gi(m⫹1)Fz•i(m⫹1) p {k1 ,k 2 },n) p( Jim p 0Fn)fni (m,{k1 ,k 2 })⫹p(Jim p 1Fn)
2(

ak1(m⫹1)

冘
K

k p1

冘冘 (
K

⫹ p( Jim p 2Fn)ak1(m⫹1) ak 2 (m⫹1)

www.ajhg.org

冘 (K

fni (m,{k ,k 2 }) ⫹ ak 2 (m⫹1)
K

k1 p1 k2 p1

k p1

)

fni m,{k ,k1 })]fni m,{k1 ,k 2 }) ,The American Journal of Human Genetics Volume 78 April 2006

641

for m p 1, … ,M ⫺ 1, where fni (1,{k1 ,k 2 }) : p p( gi1Fz•i1 p {k1 ,k 2 })ak11 ak 21 and Jim is defined in appendix C. Then,KKp( giFn) may be calculated as 冘k1p1 冘k 2p1 fni (M,{k1 ,k 2 }).
The corresponding backward recursion isbni (m ⫺ 1,{k1 ,k 2 }) p p( Jim p 0Fn)p( gimFz•im p {k1 ,k 2 },n)bni (m,{k1 ,k 2 })⫹p(Jim p 1Fn)
2

(冘 (
K

kp1

p gimFz•im p {k1 ,k},n)bni (m,{k1 ,k})akm

冘(
K

⫹

kp1

)

p gimFz•im p {k,k 2 },n)bni (m,{k,k 2 })akm冘冘 (K⫹p( Jim p 2Fn)K

k1p1 k 2p1

p gimFz•im p {k1 ,k 2 })bni (m,{k1 ,k 2 })ak1 mak 2 m ,for m p 2, … ,M ⫺ 1, and with bni (M,{k1 ,k 2 }) : p 1 for all k1 ,k 2.
To obtain p(z•im p {k1 ,k 2 }Fgi ,n), we usep(z•im p {k1 ,k 2 }Fgi ,n) ∝ fni (m,{k1 ,k 2 })2 I{k (k }bni (m,{k1 ,k 2 }) ,1

冘 冘 p(z p {k ,k }Fg ,n) p 1 .
Kwith the constraint that2

K

k1p1 k 2pk1

•
im

1

2

i

Appendix B
Sampling from p (dFgii ,n)Recall that di denotes the pair of haplotypes (hia,hib) for individual i. Additionally, let wi denote the ordered pairof cluster-of-origin indicators that correspond to the haplotypes (hia,hib) . Thus, di and wi may be thought of as“phased versions” of gi and z i, respectively. To sample from p(diFgi ,n), perform the following.
1. Sample z̃•i from p(z•iFgi ,n). This involves sampling the hidden state z•i , conditional on the data gi and parametersn, which is a standard procedure for HMMs. First, sample z̃•iM ∼ p(z•iMFg,n) ∝ p( g,z•iM,n). Then, recursively form p M ⫺ 1, … ,1, sample z˜ •im fromp(z•imFz•i(m⫹1) ,gi ,n) ∝ p( gi1 , … ,gim ,z•imFn)p(z•i(m⫹1)Fz•im ,n) p fni (m,z•im )2 I{z (z } pm⫹1(z•im r z•i(m⫹1) ) ,
a
im

b
im

where pm⫹1(z•im r z•i(m⫹1)) is given by equation (9).
˜•˜ i• ,n). Since
2. Sample w̃i from p(wFgii ,zi ,n) p p(wFzi˜ i ,n) p p(wi1Fz˜ i1 ,n)p(wi2Fwi1 ,z˜ i2 ,n) … p(wiMFwi(M⫺1) ,z˜ iM,n) ,
p(wFz
i
•
•
˜ im may be sampled sequentially (for m p 2, … ,M ), given w
˜ i(m⫺1) and z˜ im
p {k1 ,k 2 }, there are,each w. Given z˜ im
at most, two possibilities for wim: (k1 ,k 2 ) and (k 2 ,k1). Thus, probabilities of these outcomes are

˜ i(m⫺1) p (k1 ,k 2 ),z•im p {k1 ,k 2 },n) ∝ pm (k1 r k1 )pm (k 2 r k 2 ) andp(wim p (k1 ,k 2 )Fw
˜ i(m⫺1) p (k1 ,k 2 ),z•im p {k1 ,k 2 },n) ∝ pm (k1 r k 2 )pm (k 2 r k1 ) .
p(wim p (k 2 ,k1)Fw
˜ i ,gi ,v) p 写 M
˜
3. Sample d˜ i from p(diFw
mp1 p(dimFwim ,v) . This is nontrivial only for heterozygous sites—that is, when
gim p 1. Then,
a
1⫺him
b
him
k2 m

a
b
a
b
˜ im
˜ im
p(dim p (him
,him
)F(w
,w
) p (k1 ,k 2 ),gim ,v) ∝ vkh1imm (1 ⫺ vk1 m )
a

642

v

1⫺hbim(1 ⫺ v )The American Journal of Human Genetics Volume 78 April 2006

k2m,

www.ajhg.org

abfor (him,him) p (0,1),(1,0).

Appendix C
EM Algorithm
Here, we describe an EM algorithm for the estimation of n p (a,v,r). To do this, we introduce latent variables
relating to “jumps” that occur in the continuous Markov jump process underlying z•i (see “Local clustering of
haplotypes” in the “Material and Methods” section). Specifically, let Jim denote the number of jumps between
markers m ⫺ 1 and m for individual i, and let Jimk denote the number of these that jump to cluster k. ThusK
Jim : p 冘kp1 Jimk ande⫺2rm d m,
jp0p( Jim p jFr) p 2 (1 ⫺ e⫺rm d m) e⫺rm d m, j p 1 .
(1 ⫺ e⫺rm d m)2 ,jp2{Now, let Q(nFn∗ ) be the expected complete-data loglikelihood, En[log p( g,z,JFn∗ )Fg] . The algorithm is first initiatedwith a random guess n (0). Then, the following is repeated for c p 1, … ,C:n (c⫹1) p argmax Q(nFn (c) ) ,nfor sufficiently large C. The maximization above is accomplished by finding solutions to [⭸Q(nFn∗)] / [⭸vkm ] p 0,[⭸Q(nFn∗)] / [⭸akm ] p 0, and [⭸Q(nFn∗)] / [⭸rm ] p 0, for all k p 1, … ,K and m p 1, … ,M (r1 is not defined). This leads
to the following estimators for n:

冘冘 I
i k

v̂km p{gim(0}(v∗km (1 ⫺ v∗k m)

I{gimp1}∗v∗km (1 ⫺ v∗k m) ⫹ v∗k m (1 ⫺ vkm))

p(z•im p {k,k  }Fgi ,n∗ )2 I{k pk}
冘冘 p(z p {k,k }Fg ,n )2
i

•
im

k



∗

i

I{k pk}

冘 E [J Fg]â p,冘冘 E [J Fg]
i

n∗

imk

n∗

imk 

(C1)

(C2)

km

i k

,

and

( 冘冘

⫺ log 1 ⫺ i kr̂m pEn∗ [JimkFg]
2n

dm

)

.

(C3)

In practice, it would be inefficient to calculate r̂ using equation (C3) only to exponentiate it as in equation (4).
In fact, when using the model of (7) or (11) in which rm is estimated separately in each marker interval, (1 ⫺e⫺rm d m ) could be replaced by a single parameter and estimated with

冘 冘 E [J Fg]
i

k

n∗

2n

imk

.

However, writing it as above (eq. [C3]) facilitates construction of EM algorithms for the constrained model inwhich all rm are equal.
www.ajhg.org

The American Journal of Human Genetics Volume 78 April 2006

643

Finally, we give expressions for terms necessary in calculating equations (C1–C3). Note that p(z•imFg,n) depends
only on the data for individual i:p(z•im p {k1 ,k 2 }Fg,n) ∝ fni (m,{k1 ,k 2 })2 {k1(k 2 }bni (m,{k1 ,k 2 }) ,with 冘k1p1 冘k 2pk1 p(z•im p {k1 ,k 2 }Fg,n) p 1 .
2
To calculate En [ JimkFg], we use 冘jp0 j # p( Jimk p jFgi ,n), which reduces to
K

K

En[ JimkFg] p

冘 (
( ) [

akm

K

p giFn k p1

p Jim p 1Fr)冘 (Kk p1

fni m ⫺ 1,{k ,k  })]⫹ 2p( Jim p 2Fr)p( gi(⭐m⫺1)Fn)ak m p( gimFz•im p {k,k  },v)bni (m,{k,k  }) .

Web Resources
The URLs for data presented herein are as follows:
GERBIL, http://www.cs.tau.ac.il/˜rshamir/gerbil/
HAP Web site, http://research.calit2.net/hap/
HaploBlock, http://bioinfo.cs.technion.ac.il/haploblock/International HapMap Project, http://www.hapmap.org/
SeattleSNPs, http://pga.gs.washington.edu
Stephens Lab Web site, http://www.stat.washington.edu/stephens/
software.html (for PHASE and fastPHASE software)

References
Akaike H (1974) A new look at the statistical model identification.
IEEE Trans Automatic Control AC 19:719–723
Bates JM, Granger CWJ (1969) The combination of forecasts. OperRes Q 20:451–468
Breiman L (1996) Bagging predictors. Mach Learn 24:123–140
Chapman J, Cooper J, Todd J, Clayton D (2003) Detecting diseaseassociations due to linkage disequilibrium using haplotype tags: aclass of tests and the determinants of statistical power. Hum Hered
56:18–31
Dempster AP, Laird NM, Rubin DB (1977) Maximum likelihood fromincomplete data via the EM algorithm. J R Stat Soc Ser B 39:1–38
Fallin D, Schork NJ (2000) Accuracy of haplotype frequency estimation for biallelic loci, via the expectation-maximization algorithm
for unphased diploid genotype data. Am J Hum Genet 67:947–959
Falush D, Stephens M, Pritchard JK (2003) Inference of population
structure using multilocus genotype data: linked loci and correlated
allele frequencies. Genetics 164:1567–1587
Freund Y, Schapire R (1996) Experiments with a new boosting algorithm.
In: Proceedings of the Thirteenth International Conference of Machine
Learning. Morgan Kaufmann, San Francisco, pp 148–156
Greenspan G, Geiger D (2004) Model-based inference of haplotype
block variation. J Comput Biol 11:493–504
Halperin E, Eskin E (2004) Haplotype reconstruction from genotype
data using imperfect phylogeny. Bioinformatics 20:1842–1849
International HapMap Consortium (2005) The International HapMap
Project. Nature 437:1299–1320
Kimmel G, Shamir R (2005a) A block-free hidden Markov model for
genotypes and its application to disease association. J Comput Biol
12:1243–1260——— (2005b) GERBIL: genotype resolution and block identification
using likelihood. Proc Natl Acad Sci USA 102:158–162

644

Koivisto M, Perola M, Varilo T, Hennah W, Ekelund J, Lukk M,
Peltonen L, Ukkonen E, Mannila H (2003) An MDL method forfinding haplotype blocks and for estimating the strength of haplotype block boundaries. In: Altman RB, Dunker AK, Hunter L, Jung
TA, Kline TE (eds) Proceedings of the Pacific Symposium on Biocomputing. Vol 8. World Scientific, Teaneck, NJ, pp 502–513
Li N, Stephens M (2003) Modeling linkage disequilibrium and identifying recombination hotspots using single-nucleotide polymorphism data. Genetics 165:2213–2233
Lin S, Cutler D, Zwick M, Chakravarti A (2002) Haplotype inference
in random population samples. Am J Hum Genet 71:1129–1137
Marchini J, Cutler D, Patterson N, Stephens M, Eskin E, Halperin E,
Lin S, Qin ZS, Munro HM, Abecasis GR, Donnelly P, for the International HapMap Consortium (2006) A comparison of phasing algorithms for trios and unrelated individuals. Am J Hum Genet 78:
437–450
Pritchard JK, Stephens M, Donnelly P (2000) Inference of population
structure using multilocus genotype data. Genetics 155:945–959
Rabiner LR (1989) A tutorial on HMM and selected applications inspeech recognition. Proc IEEE 77:257–286
Rannala B, Mountain JL (1997) Detecting immigration by using multilocus genotypes. Proc Natl Acad Sci USA 94:9197–9201
Rastas P, Koivisto M, Mannila H, Ukkonen E (2005) A hidden Markov
technique for haplotype reconstruction. Lect Notes Comput Sci
3692:140–151
Schwarz G (1978) Estimating the dimension of a model. Ann Stat 6:
461–464
Smouse PE, Waples RS, Tworek JA (1990) A genetic mixture analysis
for use with incomplete source population data. Can J Fisheries
Aquatic Sci 47:620–634
Stephens M, Donnelly P (2003) A comparison of Bayesian methods
for haplotype reconstruction from population genotype data. Am J
Hum Genet 73:1162–1169
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
Stephens M, Smith NJ, Donnelly P (2001) A new statistical method
for haplotype reconstruction from population data. Am J Hum Genet 68:978–989
Sun S, Greenwood CTM, Neal RM (2004) Haplotype inference usinga hidden Markov model with efficient Markov chain sampling [abstract 2934]. In: Proceedings and abstracts of the American Societyof Human Genetics 2004 Annual Meeting, Toronto, October 26–30
Zöllner S, Pritchard JK (2005) Coalescent-based association mappingand fine mapping of complex trait loci. Genetics 169:1071–1092

The American Journal of Human Genetics Volume 78 April 2006

www.ajhg.org