\chapter{Bayesian Inference in Mixed and Admixed Populations}


\textsl{Biologists regularly encounter populations of organisms with disparate ancestries.  Untangling the
composition of such populations is a problem for conservation biologists and wildlife
managers.  In many cases the population under question is known to consist of individuals from two
different subpopulations and their hybrids.  This occurs, for example, in hybrid zones between two species or
in regions recently colonized by exotics capable of reproducing with resident inhabitants.  This paper develops
techniques using multilocus genetic data for Bayesian clustering of individuals to purebred or
genetically-mixed categories.  The method relies on a novel application of the forward-backward recursions
in a two-component, finite mixture model.  Though developed in the context of the genetic
admixture problem, these calculations are relevant more generally to Bayesian inference in finite mixtures; they may
potentially improve mixing of the Gibbs sampler in such contexts.  
The technique is applied to genetic data on the Scottish wildcat, {\em Felis~sylvestris}, a protected species
whose distinctness from domestic housecats has been questioned.  A high proportion ($\approx .60$) of the
wild-living cats from which the sample was drawn are arguably purebred {\em F.~sylvestris}. } 

\textsl{Using the Bayes factor, we compare our new model, which allows for both purebred and admixed
individuals, to a model in which all individuals are assumed genetically admixed to some degree.  It is
difficult to accurately compute the marginal likelihood directly in these models, so we compute the Bayes
factor by reversible-jump MCMC\@. The approach follows from the original MCMC formulation of the problem,
and should help to illustrate ways in which reversible-jump methods may be implemented for comparisons
between a small set of closely-related models.}   

\section{Introduction}
Populations studied by geneticists are seldom the ideal, randomly-mating
and genetically-isolated collections of individuals for which much
genetic theory has been developed.  In particular, natural populations
may possess internal structure which prohibits random
mating, or they may be recently formed by migration
and co-mingling of individuals from two or more originally separate
populations. Such structure complicates many types of genetic
studies.  For instance, when using population-level data to map genetic
diseases, population structure, if not accounted for, may lead to
spurious associations between genetic markers and disease status
\cite{Ewe&Spi1995}.  Additionally, in the ecological study of plants
and animals there is considerable interest in population structure,
especially in regions of apparent overlap and interbreeding between
different subpopulations---so-called ``hybrid zones." For these, and
other types of problems, it is desirable to be able to infer population
structure from genetic data.  To this end, models of population genetic
{\em mixture} and {\em admixture} have been useful.  I describe the
application of such models to the inference of population structure,
focusing primarily on applications to hybrid zones of two different
groups of individuals.  Such situations are now encountered frequently as
a result of anthropogenic disturbance reducing barriers to gene flow
between formerly separate subpopulations.  Invasions of exotic species
are one pervasive example.  

As used here, a ``genetic mixture model" attributes structure in
a population to the presence of two or more
subpopulations.  Within these subpopulations individuals may mate at
random, but no interbreeding occurs between subpopulations. 
Every individual in the mixture is
considered to be a purebred descendant of only one subpopulation. 
Such models have been developed and used extensively in the field of
fisheries management \cite{Milneretal1981,Pel&Mil1987,Smouseetal1990,Millar1991,Pel&Mas}.  


On the
other hand, admixture, throughout the genetics
literature
\cite{Cav&Bod1971,Thompson1973,Wijsman1984,Long1991} refers to
interbreeding between members from different subpopulations. Accordingly,
 a ``genetic admixture model" attempts to
model a population's genetic structure by the presence of two or more
previously separate subpopulations between which there has been some
recent interbreeding. Such a population is said to be admixed. 
Additionally we will call an individual ``admixed" if it possesses
genes descended from more than one of the historically separate
populations. 
Early investigations of admixed populations sought to
estimate the relative contributions of two founding populations to the
admixed population.  These studies assumed the admixed population had
undergone enough generations of random mating to eliminate allelic
associations between loci that accompanies genetic admixture.  As a
result, the individual allele frequencies observed in the two founding
subpopulations and the admixed population were treated as sufficient
statistics.  
 It was not until recently that
statistical methods were proposed whereby the information in {\em
multilocus} data could be used to elucidate structure
in a recently admixed population 
\cite{Ran&Mou1997,Paetkauetal1995,Pritchardetal2000}.


\citeN{Pritchardetal2000} propose a versatile model for genetic
inheritance in admixed populations and use it in 
Bayesian analyses of population structure in several different
species.  A limitation of this model, however, is that it
assumes every individual is admixed to some degree.  In many
situations, such as with populations spanning hybrid zones, there is
reason to expect both purebred and admixed individuals.  A probability
model to accommodate such scenarios will include elements both of
genetic mixture models and genetic admixture models.  
In this paper I extend the methods of \citeN{Pritchardetal2000} to handle
explicitly  purebred individuals.  In sections~\ref{sec:Mix}
and~\ref{sec:Adm}, I review mixture and admixture formulations for modeling
population structure.  

In Section~3, I develop a method for making
joint, Gibbs updates of large blocks of variables in the \citeN{Pritchardetal2000} model. 
The method uses the fact that the latent allocation variables of an i.i.d.~finite mixture, with a Dirichlet
prior on mixing proportions can be shown to follow a hidden Markov chain, after integrating out the mixing 
proportions.
This computation facilitates MCMC simulation in a model, described in Section~\ref{sec:Simul}, that
allows for both purebred and admixed individuals.  
Additionally, I describe in the Discussion how such a method could help
the Gibbs sampler to escape from trapping states \cite{Robert1996}
encountered in finite mixture problems. 

I apply these techniques to data on the Scottish wildcat {\em Felis
sylvestris}. In Scotland, {\em F.~sylvestris} evolved for thousands of
years with little or no genetic exchange with cats in continental
Europe.  Within the last 2,000 years these Scottish cats have suffered
population declines due to human influences and have been exposed to
possible interbreeding with domestic cats.  It can be difficult to
distinguish {\em F.~sylvestris} from domestic cats on the basis of
morphological characters alone and conservation biologists are concerned
that the wild-living cats in Scotland may now represent an admixture of
{\em F.~sylvestris} and domestic cats.  The data were previously analyzed
by \citeN{Beaumontetal} using the method of \citeN{Pritchardetal2000}.  
However, this analysis does not address the issue of particular
interest---that of 
estimating the proportion of purebred
{\em F.~sylvestris} individuals in the population.  Nor does that analysis allow estimation
of posterior probabilities that particular individuals in the sample are purebred cats.  These questions about
the Scottish wildcat population are similar to those for many species of conservation interest to which the
present methods apply.  

Finally, using reversible-jump MCMC, it is possible to compute the Bayes factor for
comparing the new, expanded model to that of \citeN{Pritchardetal2000} given the Scottish
cat data.  While the reversible-jump sampler allows  estimation of the true Bayes factor, 
it is also possible to compute the ``pseudo-Bayes factor" \cite{Gelfandetal1992}, and assess how
accurately that estimates the Bayes factor.


\section{Genetic Mixture Models}
\label{sec:Mix}
The formulation of a genetic mixture model follows that for a general
finite mixture.  Let $N$ diploid individuals be sampled from a population
and typed at $L$ loci.  The population is assumed to consist of $J$
subpopulations indexed by $j=1,\ldots,J$.   The proportion of individuals in the mixed population
from subpopulation $j$ is the unknown parameter $\pi_j$ with
$\sum_{j=1}^J \pi_j = 1$. Assign to each individual
a latent allocation variable
$Z_i$, $i=1,\ldots,N$.  We use $Z_i=j$ to indicate that
the $i\thh$ individual is from the $j\thh$ subpopulation.

Denote the multilocus phenotype of the $i\thh$ individual by $\btY_i$. 
I use ``phenotype" as opposed to ``genotype" because the phrase
``multilocus genotype" typically implies knowledge of the gametic phase
 of the alleles present at
different loci. No
 such knowledge is available here.  The multilocus phenotype, $\btY_i$, consists of the
allelic type of each of the two alleles carried by the $i\thh$ individual
at $L$ loci.  Since we will later  identify and label specific gene
copies in an individual, we consider the two alleles carried at a locus
to be ordered.  This order may be arbitrary.  For example, it can
merely be the order in which the types of those two alleles are reported
in the data on an individual.  Thus,
$\btY_i$ can be regarded as a vector of length $2L$ with its first element
giving the allelic type of the first allele at locus~1, its second
element giving the type of the second allele at locus~1, its third
element the type of the first allele at locus~2 and so forth.  In
general, $Y_{i,t}$,
$t=1,\ldots,2L$, is the type of allele number ($t~\modtwo +1$)  at locus
number $\ceil{t/2}$ in individual $i$, where $\ceil{x}$ denotes the
smallest integer greater than or equal to
$x$.

The allele frequencies in the $j\thh$ subpopulation are denoted by 
$\btheta_j = (\bttheta_{j,1},\ldots,\bttheta_{j,L})$, where $\bttheta_{j,\ell}$,
$\ell=1,\ldots,L$, is a vector of length equal to $K_\ell$---the number of
types of alleles observed at locus~$\ell$ across all the individuals
sampled.   The frequency in
the $j\thh$ subpopulation of
the $k\thh$ allelic type of the $\ell\thh$ locus is $\theta_{j,\ell, k}$,
$k=1,\ldots,K_\ell$.  We will adopt the notation $\thetaof{j}{y_{i,t}}$ 
to mean
$\theta_{j,\ell, k}$ where $\ell = \ceil{t/2}$ and $k$ is the allelic type
of the $(t~\modtwo + 1)\thh$ allele at the $\ell\thh$ locus in the
$i\thh$ individual.


Given that an individual is from subpopulation $j$, it is assumed to have
a multilocus phenotype resulting from random mating  and linkage
equilibrium between the $L$ loci {\em within} subpopulation $j$.  Thus, 
\begin{equation}
P(\btY_i|\btheta_j,Z_i=j) = \prod_{t=1}^{2L} \thetaof{j}{\btY_{i,t}}
\end{equation}
where $P(\cdot|\cdot)$ will be used throughout to denote conditional
probability mass or density functions.  
The likelihood for $\bpi=(\pi_1,\ldots,\pi_J)$ and
$\bTheta=(\btheta_1,\ldots,\btheta_J)$, with $\bY$ denoting $(\btY_1,\ldots,\btY_N)$,
is
\begin{equation}
P(\bY|\bpi,\bTheta) = \prod_{i=1}^N P(\btY_i|\bpi,\bTheta) = \prod_{i=1}^N
\sum_{j=1}^J \pi_j P(\btY_i|\btheta_j,Z_i = j).
\label{eq:gtypProbP}
\end{equation}
Note that this formulation does not include
the familiar binomial coefficient, 2, for heterozygotes because we have
arbitrarily ordered the two alleles carried by an
individual at each locus.  This makes the likelihood in this model
comparable to that in the admixture model of \citeN{Pritchardetal2000}.

``Training" or ``learning" samples may be 
available.  They might take the form of specially tagged individuals
which, though sampled along with the rest of the mixture, may be
unambiguously assigned to a subpopulation.  Such an individual, say
$i^*$, known to come from subpopulation $j^*$ is easily accommodated by
setting $z_{i^*} = j^*$ and defining
$P(\btY_{i^*}|\btheta_j,Z_{i^*}=j)
\equiv 0$ for all $j\neq j^*$.  However, if a learning sample from the
$j\thh$ subpopulation is drawn separately (for
example,  if taken during a season when the subpopulations can be sampled
separately) it contributes a term of the form
$C\cdot\prod_{\ell=1}^L\prod_{k=1}^{K_\ell} \theta_{j,\ell, k} ^{n_{j,\ell,
k}}$ to the likelihood, where $C$ is a product of multinomial coefficients
and
$n_{j,\ell ,k}$ is the number of alleles of type $k$ observed at locus
$\ell$ in the learning sample taken separately from the $j\thh$
subpopulation.  (In the Bayesian framework, these changes are equivalent to altering the
prior for $\bTheta$ and $\bpi$ appropriately.)

Treating this mixture model from the Bayesian perspective requires prior
distributions for $\bpi$ and $\bTheta$.  The conjugate prior for
$\pi$ is the Dirichlet distribution,
$\Dir(\zeta_1,\ldots,\zeta_J)$.  Prior information could be incorporated
in the values of the $\zeta_j$, or, if no prior information is available,
the uniform distribution $\zeta_j=1$, $j=1,\ldots,J$, is a reasonable
choice.   The conjugate prior for each
$\btheta_{j,\ell}$ is $\Dir(\lambda_{j,\ell, 1}, \ldots,\lambda_{j,\ell,
K_\ell})$.  In this paper, I use uniform Dirichlet priors,
$\lambda_{j,\ell, k} = 1$ $\forall j,\ell,k$, which tend to de-emphasize
the influence of rarely-occurring allelic types.  This is a
conservative assumption, and works well when the subpopulations are
sufficiently genetically distinct.  Note, however, that
\citeN{Pritchardetal2000} discuss different Dirichlet priors and
\citeN{Pel&Mas} describe a useful Empirical Bayes approach
to assigning allele frequency priors in the application of Bayesian
methods to mixed-stock fishery analysis with closely-related
subpopulations.  

With the priors specified, the posterior distribution of $\bpi$ and
$\bTheta$, as well as other quantities of interest, may be investigated via
Gibbs sampling as described by \citeN{Die&Rob1994}.  The relevant full
conditionals are
\begin{eqnarray}
\bpi|\cdots  &\sim&  \Dir(\zeta_1 + \#\{Z=1\},\ldots,\zeta_J + \#\{Z=J\} )
\nonumber   \\ 
\rule{0pt}{5ex}
\bttheta_{j,\ell}|\cdots &\sim &\Dir(\lambda_{j,\ell, 1} + m_{j,\ell ,1} +
n_{j,\ell ,1}, \ldots, \lambda_{j,\ell, K_\ell} + m_{j,\ell, K_\ell} +
n_{j,\ell, K_\ell}), 
 \nonumber \\
& & j = 1,\ldots, J;~~\ell = 1,\ldots, L \nonumber \\
\rule{0pt}{5ex}
P(Z_i=j|\cdots) &=& \frac{\pi_j P(\btY_i|\btheta_j,Z_i=j)}
{\sum_{k=1}^J \pi_j P(\btY_i|\btheta_j,Z_i = k)}~,~
i=1,\ldots,N;~j=1,\ldots,J \nonumber
\end{eqnarray}
where $\#\{Z=j\}$ is the number of individuals currently allocated to
subpopulation $j$, $m_{j,\ell, k}$ is the number of alleles of type $k$ at
locus $\ell$ in individuals currently allocated to subpopulation $j$, and
the $n_{j,\ell, k}$ are, as before, the allele counts from the learning
samples (if any) drawn separately from the mixture sample.  


\section{A Model with Admixed Individuals}
\label{sec:Adm}
With $\bTheta$ and $\bY$ defined as in the previous section, the model of
\citeN{Pritchardetal2000} is quickly described.   Now, $j$ indexes the $J$
conceptual ``gene pools" or ``historical subpopulations" from which
individuals may be descended.  Allowing for admixed individuals requires
a different model of genetic inheritance, which, in turn, requires
different latent variables.  The $i\thh$ individual in the
sample gets  a vector of probabilities $\btQ_i = (Q_{i,1},\ldots,Q_{i,J})$, 
$\sum_{j=1}^J Q_{i,j} = 1$, which are the unobserved proportions of that
individual's genome descended from each of the $J$ gene pools. 
Also, let
$\btW_i = (W_{i,1},\ldots,W_{i,2L})$ be a vector of unobserved allocation
variables which is parallel to the the vector of allelic types
$y_i$.  Hence, $W_{i,t}=j$  indicates that the $(t~\modtwo +1)\thh$
allele at the $\ceil{t/2}\thh$ locus in the $i\thh$ individual is from
the $j\thh$ gene pool.   Given $w_{i,t} = j$ the type of allele is assumed
to be drawn randomly according to $\theta_j$.  Under this model
\begin{eqnarray}
P(y_i|\theta,w_i) = \prod_{t=1}^{2L} \thetaof{w_{i,t}}{y_{i,t}}
\end{eqnarray}
independently for each $i$.
By assigning the prior $q_i\sim\Dir(\alpha,\ldots,\alpha)$,  
$i=1,\ldots,N$, and the hyperprior $\alpha \sim \mathrm{Uniform}(0,A]$,
\citeN{Pritchardetal2000}'s model is obtained.  In effect this is a 
hierarchical model for $N$ different finite mixtures---the genes carried
by the $i\thh$ individual are a sample from a mixture with mixing proportions 
given by $q_i$, while the $q_i$ themselves $(i=1,\ldots,N)$ are drawn from
a symmetrical $\Dir(\alpha,\ldots,\alpha)$ distribution.




In this model, Gibbs sampling
proceeds  using the full conditionals
\begin{eqnarray}
q_i|\cdots  &\sim&  \Dir(\alpha_1 + \#\{w_i=1\},\ldots,\alpha_J +
\#\{w_i=J\} ),~~i=1,\ldots,N
\nonumber   \\ 
\rule{0pt}{5ex}
\theta_{j,\ell}|\cdots &\sim &\Dir(\lambda_{j,\ell, 1} + r_{j,\ell, 1},
\ldots, \lambda_{j,\ell, K_\ell} + r_{j,\ell, K_\ell}), 
 \nonumber \\
& & j = 1,\ldots, J;~~\ell = 1,\ldots, L \nonumber \\
\rule{0pt}{5ex}
P(w_{i,t}=j|\cdots) &=& \frac{q_{i,j} \thetaof{j}{y_{i,t}}}
{\sum_{k=1}^J q_{i,j} \thetaof{k}{y_{i,t}}}~,~
i=1,\ldots,N;~j=1,\ldots,J; \nonumber  \\
& & t = 1,\ldots,2L \nonumber
\end{eqnarray}
where $\#\{w_i=j\}$ is the number of alleles in the $i\thh$ individual
currently allocated to gene pool $j$ and $r_{j,\ell, k}$ denotes the number
of alleles of type $k$ at locus $\ell$
currently allocated to gene pool
$j$. 
\citeN{Pritchardetal2000} update $\alpha$ 
by a Metropolis-Hastings method as described in Section~{sec:MH}. The posterior
distribution of
$\alpha$ thus estimated provides some insight into the degree to which
admixture has occurred across individuals.  


Learning samples would be available if
there were substantial prior knowledge about the gene pools contributing
to the admixture and if known, purebred descendants from them
were separately sampled.   By assuming any
effects of genetic drift to be negligible, such samples could be treated
as learning samples in the mixture model.  The full conditional for
$\theta_{j,\ell}$ would then be modified to include the $n_{j,\ell, k}$ as
before.  

\subsection{Block-updating $w_i$ when $J=2$}
\label{sec:Polya}
In many situations involving invasions of exotic species, there is
substantial prior knowledge that the number of major subpopulations or
``gene pools" involved is two---the native population and the invading
population.  Additionally, many hybrid zones are known to be areas of
hybridization (admixture) between two species or populations.  Here I
present novel computations, feasible when only two subpopulations or gene
pools are involved, that eliminate the explicit need for the variable
$q=(q_1,\ldots,q_N)$ in implementing a Gibbs sampler.   Such a method
slightly improves mixing of the chain, but is primarily useful as it
makes possible Gibbs sampling in a simultaneous mixture and admixture
analysis as will be described in Section~\ref{sec:Simul}.

The computations themselves may be derived as follows.  Let $J=2$, so
that each allele in an individual may have originated from gene
pool~1 or gene pool~2.  Then, each $q_{i,1}$ will follow a
$\mathrm{Beta}(\alpha,\alpha)$ distribution and $q_{i,2} = 1 - q_{i,1}$.
Conditional on
$q_{i,1}$, each $w_{i,t}$ will then be independently a Bernoulli trial with
$P(w_{i,t} = 1|q_{i,1}) = q_{i,1}$.  Marginalizing over $q_{i,1}$ (not
conditioning on the data) it follows that $\#\{w_i=1\}$ follows a
beta-binomial distribution with parameters $(\alpha,\alpha)$.  Of course,
each allele in an individual is uniquely labelled so the elements of
$w_i$ may be interpreted as following a {\em labelled} beta-binomial
distribution.  Under such a distribution, the elements of $w_i$ are not
independent, but they are exchangeable \cite{deFinetti1972},
and hence their marginal distributions are invariant to permutations of
their order (and thus the arbitrary order we have imposed upon them is
acceptable).  

This labelled beta-binomial sampling mechanism is easily visualized
by a P\'{o}lya-Eggenberger urn scheme
\cite{Feller1957,Joh&Kot1977}.  Imagine an urn
initially filled with $b_1$ balls labelled ``1" and $b_2$ balls labelled
``2."  Draw a ball randomly and record
$w_{i,1} = 1$ or $2$ according to the ball's label.  Then replace the
ball to the urn, adding, at the same time, $c$ more balls of the
same type (1 or 2) as the ball just drawn.  Repeat the process, assigning
a value to
$w_{i,2}$ and so forth until $w_{i,2L}$ has also been assigned a 1 or
2\@.  If
$b_1$,
$b_2$, and $c$ were chosen to satisfy $b_1/c = b_2/c = \alpha$, then the
resulting vector $w_i$ would be a realized value from the labelled
beta-binomial distribution with parameters $(\alpha,\alpha)$.  (One should notice,
also, that this extends to a non-symmetrical beta distribution,
say $\mathrm{Beta}(\alpha_1,\alpha_2)$, by choosing $b_1/c = \alpha_1$ and $b_2/c = 
\alpha_2$.)

By such a scheme it is apparent that if $d_t$ balls of type~1 have been
drawn in the first $t$ drawings from the urn, then the
probability that the next ball drawn is a 1 is given by
\begin{equation}
\frac{b_1 + d_t c}{b_1 + b_2 + tc}.
\label{eq:transprob}
\end{equation}
And so the pairs $(w_{i,t},d_t)$, $t=1,\ldots,2L$, can be interpreted as
forming a time-inhomogeneous Markov chain in time $t$ with transition
probabilities determined by (\ref{eq:transprob}) and the obvious
fact that $d_{t+1} = d_t + 1\{w_{i,t+1}=1\}$, where $1\{x=a\}$ takes the
value one when $x=a$ and zero otherwise.  This Markov chain dependence
structure was previously noted by \citeN{Freedman1965}, who used 
it to obtain limiting distributions of quantities associated with urn models.

The foregoing has all been considered in the absence of data, $y_i$. 
However, given $\theta$, the data provide some information about the true
value of each $w_{i,t}$ by the relation $ P(y_{i,t}|w_{i,t},\theta) =
\thetaof{w_{i,t}}{y_{i,t}}$.    Therefore, conditional on $\theta$
and $y_{i,t}$, the pairs $(w_{i,t},d_t)$ participate in a {\em hidden}
Markov chain. Recognition of this fact allows application of a
``filter-forward, simulate-backward" type of algorithm which may be derived following the  
computations of \citeN{Baumetal1970} in order to
realize the elements of $w_i$ from their joint full conditional
distribution, $P(w_i|\alpha, \theta, y_i)$.  Furthermore, using the
\citeN{Baum1972} algorithm, it is possible to compute
$P(y_i|\alpha,\theta)$, effectively performing a sum over all possible
binary vectors of length $2L$ in an efficient manner.  This is described below.

Take $b_1$, $b_2$, and $c$ as defined above.  Suppressing the $i$
subscript for clarity, let $w_t \in \{1,2\}$, $t=1,\ldots,2L$, and define
$d_t = \sum_{\tau=1}^t 1\{w_\tau = 1\}$.  We adopt the notation $w_{\leq
t}$ ($w_{\geq t}$) to mean $w_1,\ldots,w_{t}$ ($w_t,\ldots,w_{2L}$) for
 components of $w$, and use the same notation with $y$ and $d$.   The
pairs
$(w_t,d_t)$ can be interpreted as following a Markov chain in $t$:
\begin{eqnarray}
P(w_{t+1},d_{t+1} | w_{\leq t},d_{\leq t})& = &P(w_{t+1},d_{t+1} | w_{
t},d_{ t})  \nonumber \\
&=& 
\frac{b_1 + d_tc}{b_1+b_2+tc}
1\{d_{t+1} = d_t + 1\{w_{t+1} = 1\} \}\nonumber .
\end{eqnarray}
The ``perturbed" or ``degraded" observations of the chain are the allelic types 
$y_1,\ldots,y_{2L}$ which
depend in hidden Markov fashion on $w$.  For notational clarity, we assume
implicit dependence on the allele frequencies $\theta$,
\[
	P(y_t | w_{\leq 2L}, d_{\leq 2L}) = P(y_t | w_t) = \thetaof{w_t}{y_t}.
\]
This dependence structure is shown in the undirected graph of
Figure~\ref{fig:undirgraph}.
\begin{figure}
\begin{center}
\includegraphics*[width=.75\textwidth]{undirgraph.eps}
\caption[Undirected graph showing the dependence between $w$, $d$]{An undirected graph showing the
dependence between $w$, $d$ and
$y$ in Section~\ref{sec:Polya}.  This graph describes hidden Markov structure
for the pairs $(w_t,d_t)$.  The dependence on $\theta$ is implicit
and not shown. }
\label{fig:undirgraph}
\end{center}
\end{figure}

In the forward step we compute and store values of
$P(w_t,d_t|y_{\leq t})$ for $w_t = 1,2$ and $d_t = 0,\ldots,t$,
recursively for $t=1,\ldots, 2L$, by the relations
\begin{equation}
P(w_{t+1}, d_{t+1} | y_{\leq t}) =
\sum_{1\leq w_t \leq 2 }
P(w_{t+1},d_{t+1}|w_t,d_t^*) P(w_t,d_t^*|y_{\leq t})
\end{equation}
where $d_t^* = d_{t+1} -
1\{w_{t+1}=1\}$,
and 
\begin{equation}
P(w_{t+1},d_{t+1} | y_{\leq t+1}) = \frac{1}{\phi_{t+1}}
P(w_{t+1},d_{t+1}|y_{\leq t}) P(y_{t+1}|w_{t+1},d_{t+1})
\end{equation}
where
\begin{equation}
\phi_{t+1} = P(y_{t+1}|y_{\leq t}) = 
\sum_{\scriptstyle 1\leq w_{t+1} \leq 2 \atop \scriptstyle 0\leq d_{t+1} \leq t+1}
P(w_{t+1},d_{t+1}|y_{\leq t}) P(y_{t+1}|w_{t+1},d_{t+1}).
\end{equation}

At the end of the forward step, notice that $\prod_{t=1}^{2L} \phi_t =
P(y_1,\ldots,y_{2L})$, which in the context of the Gibbs sampler (and
if we were to reinstate the $i$ subscript) is the desired quantity
$P(y_i|\alpha,\theta)$ for the $i\thh$ individual.  At the
end of the forward step we have also obtained the distribution
$P(w_{2L},d_{2L}|y_{\leq 2L})$.  After simulating values for $w_{2L}$
and $d_{2L}$ from that distribution, we are in a position to simulate
values for $w_t$ going backwards recursively for
$t=2L-1,2L-2,\ldots,1$, using the conditional
distributions stored during the forward step and the values just realized for
$w_{t+1}$ and $d_{t+1}$.  The backward step uses the following relations recursively to 
compute the conditional distribution from which to realize values of $(w_t,d_t)$:
\begin{eqnarray}
P(w_t,d_t | y_{\leq 2L}, w_{\geq t+1}, d_{\geq t+1}) &= & 
P(w_t,d_t|y_{\leq t},w_{t+1},d_{t+1})  \nonumber \\
&= & 
\frac{1}{\psi_t} P(w_t,d_t|y_{\leq t})P(w_{t+1},d_{t+1} | w_t,d_t),
\end{eqnarray}
where $\psi_t$ is a normalizing constant
\begin{equation}
\psi_t = 
\sum_{1\leq w_t \leq 2 }
P(w_t,d_t^*|y_{\leq t})P(w_{t+1},d_{t+1} | w_t,d_t^*)
\end{equation}
and where, again, $d_t^* = d_{t+1} -
1\{w_{t+1}=1\}$.
It is apparent that a realization of the variable $(w_1,\ldots,w_{2L})$ thus
obtained is drawn from the distribution of $w_1,\ldots,w_{2L}$
conditional on $y_{\leq 2L}$.  As such, in the context of the Gibbs
sampler, and reinstating the $i$ subscript, it is a realization from
$P(w_i|\alpha,\theta,y_i)$ for the $i\thh$ individual, as desired.  

The amount of computation required for the backward step is linear in $L$.  
The forward step at time $t$ requires a handful of elementary operations for each
of the $2t$ states that the pair $(w_{t},d_{t})$ may take.  This makes the entire
forward step  $O(L^2)$ for the case of two subpopulations.  Depending on how many loci are
available  this will typically be
computationally reasonable.  However, extending this method to $J>2$ will be computationally
difficult.  With $J>2$, $d_{t}$  becomes a vector whose elements record the number of balls of type
$1,\ldots,J$ which have been drawn up to and including time $t$.  The number of possible states for the pair
$(w_t,d_t)$ is then $ J(t+J-2)!/[(t-1)!(J-1)!]$ which gets large rapidly with $t$ and $J$.


\section{A model for simultaneous population mixture and admixture}
\label{sec:Simul}
Continuing in the case of two subpopulations ($J=2$), a common goal in applications would be to identify
purebred versus admixed individuals and to estimate
the proportion of those types in the population.  This corresponds to
partitioning one's sample into purebred and admixed groups.  The $i\thh$ individual's inclusion in
one of the two groups can be denoted by a latent variable $v_i$ taking the values 
\[
v_i = \left\{
\begin{array}{ll}
\mathrm{P}, & \mathrm{if~purebred} \\
\mathrm{A}, &
\mathrm{if~admixed~following~the~Pritchard~et~al.~(2000a)~model}
\end{array}
\right.
\] 
Using the calculation of Section~\ref{sec:Polya}, this partition problem
can be treated as a mixture problem using Gibbs sampling. The proportion
of individuals of the two types in the population are given by the new
parameter $\xi = (\xi_\mathrm{P},\xi_\mathrm{A})$, with $\xi_\mathrm{P} + \xi_\mathrm{A} =
1$.  The full conditional distribution for $v_i$ is then, for example, for
$v_i=\mathrm{P}$ 
\begin{equation}
P(v_i = \mathrm{P}|\cdots) = \frac{\xi_\mathrm{P} P(y_i|\alpha,\theta,v_i
=
\mathrm{A})}{
\xi_\mathrm{P} P(y_i|\pi,\theta,v_i =\mathrm{P}) +
\xi_\mathrm{A} P(y_i|\alpha,\theta,v_i =\mathrm{A})
}.
\label{eq:vfullcond}
\end{equation}
Calculating each of the necessary phenotype probabilities, $P(y_i|\pi,\theta,v_i =\mathrm{P})$ and 
$P(y_i|\alpha,\theta,v_i =\mathrm{A})$,  has been
described in Equation~\ref{eq:gtypProbP} and  Section~\ref{sec:Polya}. 



The conjugate prior for $\xi_P$ is 
a $\mathrm{Beta}(\delta_\mathrm{P},\delta_\mathrm{A})$ which gives the full conditional 
\begin{equation}
\xi | \cdots~\sim~\mathrm{Beta}(\delta_\mathrm{P} +
\#\{v=\mathrm{P}\},\delta_\mathrm{A}+\#\{v=\mathrm{A}\}).
\label{eq:xifullcond}
\end{equation}
I have used a uniform $(\delta_\mathrm{P} = \delta_\mathrm{A} = 1)$ prior for
$\xi_\mathrm{P}$.   This prior corresponds to each individual $i$ having prior
probability of $1/2$ of being either purebred or admixed.   

In this expanded model, which we
will call model
$\Mpa$ a sweep consists of
\begin{enumerate}
\item Gibbs update for $\pi$ using only the individuals with $v_i=\mathrm{P}$,
\item Gibbs update for $\theta$ where contributions to the full conditionals are determined
by $Z_i$ for individuals with $v_i=\mathrm{P}$ and by $w_i$ for those with $v_i =
\mathrm{A}$,
\item Gibbs updates for each individual's $Z_i$ if $v_i=\mathrm{P}$ and for $w_i$ if
$v_i=\mathrm{A}$,
\item Gibbs update for $\xi$ from  Equation~\ref{eq:xifullcond},
\item Gibbs update for each $v_i$ using Equation~\ref{eq:vfullcond},
\item Metropolis-Hastings update for $\alpha$ as described in Appendix~\ref{sec:MH}.
\end{enumerate}
The output from the resulting Markov chain can provide Rao-Blackwellized \cite{Liuetal1994}
estimates for the posterior probability that individuals in the sample are purebred or admixed as
well as esimates of the posterior distributions for $\xi$, $\pi$, $\theta$,
$\alpha$, and each $q_i$ given $v_i=A$  (though the $q_i$'s are not necessary for running
the chain, they may still be realized from their full conditional distributions and they
provide good summary statistics).

\section{Metropolis Updates for $\alpha$}
\label{sec:MH}
The method of Metropolis sampling is used to update values of $\alpha$.  A new
value for $\alpha$ denoted $\alpha^*$ is drawn from a proposal distribution.  Since $\alpha$
is constrained to the interval $(0,A]$, I use a folded normal distribution, centered at
$\alpha$.  Hence a variable $a$ is drawn from a $\mathrm{Normal}(\alpha,\sigma^2)$
distribution.  If $0<a\leq A$ then $\alpha^* = a$.  Otherwise if $-A\leq a<0$ then $\alpha^*
= -a$ and if $A<a\leq 2A$ then $\alpha^* = 2A -a$.  In all other cases ($a<-A$ or $a>2A$) the proposal is
rejected without further consideration.  The proposal density is then still symmetrical
\[
h(\alpha^*|\alpha) = \mathcal{N}(\alpha^*;\alpha,\sigma^2) 
+ \mathcal{N}(-\alpha^*;\alpha,\sigma^2) + \mathcal{N}(2A-\alpha^*;\alpha,\sigma^2)  
= h(\alpha|\alpha^*)
\]
with $\mathcal{N}$ denoting the normal density function.  The standard deviation, $\sigma$, of the
proposal distribution
 requires some tuning.  Under model $\Ma$,  $\sigma \approx .12$ seems to
work well, while when individuals may be purebred or admixed (model $\Mpa$) then $\sigma
\approx .5$ encourages better mixing with the Scottish cat data.

The proposed value $\alpha^*$ is accepted as the new value with probability given by the
minimum of 1 or the  Hastings ratio.   For \citeN{Pritchardetal2000}'s model, using, the
$q_i$'s, the acceptance probability is 
\[
\min\left\{1, 
\frac{\prod_{i=1}^N \mathcal{D}(q_i;\alpha^*,J)}
{\prod_{i=1}^N \mathcal{D}(q_i;\alpha,J)} \right\}
\]
where $\mathcal{D}(q;\alpha,J)$ denotes the density of a Dirichlet random vector $q$ of $J$
components with all $J$ parameters equal to $\alpha$.

When able to eliminate the $q_i$'s (as in Section~\ref{sec:Polya}), then with only admixed
individuals (model $\Ma$) the acceptance probability may be written as
\[
\min\left\{1,  
\frac{\prod_{i=1}^N P(y_i|\alpha^*,\theta)}
{\prod_{i=1}^N P(y_i|\alpha,\theta)} \right\}.
\]
In the model $M_\mathrm{P,A}$ which includes both purebred and admixed individuals,  the
acceptance probability is
\[
\min\left\{1, 
\frac{\prod_{i=1}^N [\xi_\mathrm{P} P(y_i|\pi,\theta) + \xi_\mathrm{A}
P(y_i|\alpha^*,\theta)]} {\prod_{i=1}^N [\xi_\mathrm{P} P(y_i|\pi,\theta) + \xi_\mathrm{A}
P(y_i|\alpha,\theta)]} \right\}.
\]

\section{MUST ADD A SECTION ABOUT MISSING DATA}
IN THE JABES ARTICLE I DIDN'T MAKE ANY MENTION OF THE FACT THAT WE CAN IMPUTE
MISSING DATA VALUES.  THERE JUST WASN'T ROOM ENOUGH IN THAT PAPER FOR ANY 
OF THAT.



\section{Data and Results}
The data from Scottish wildcats were provided by Mark Beaumont (University of Reading, UK) and are fully
described in \citeN{Beaumontetal}.  The data set is freely available at
\texttt{http://www. rubic.rdg.ac.uk/~mab/\-data.html}. Briefly, genetic samples were collected from wild-living cats
throughout Scotland by a variety of methods including trapping and tissue collection from road kills and carcasses.  Samples
were also obtained from 13 museum specimens.  In all, 230 wild-living cats were sampled and typed at eight  microsatellite loci
with numbers of alleles ranging from  nine to 17 per locus.  Additionally, 74 housecats were typed at those eight loci using
blood samples from veterinary centers in the south of England.  These 74 cats can be considered a learning sample for the
domestic cat subpopulation.  

I analyzed the data under model $\Mpa$ using runs of length 62,000 sweeps of ten different chains started
from overdispersed starting points by initializing values of all parameters $(\alpha,\theta,\xi,\pi)$ with
values simulated from their prior distributions.  All ten chains converged very quickly to the same part of the
parameter space.  The first 2,000 sweeps were discarded as burn-in, as observing
the estimated scale reduction potential factor \cite{Gelman1996} suggests this is more than
adequate burn-in.  I give the results in the next section.
I performed an analagous run under model $\Ma$ and compare the differences between the results obtained under
$\Mpa$ and $\Ma$ in Section~\ref{sec:comparison}.  For both runs I used an upper bound of $A = 3$ for the
parameter $\alpha$.  Each run took about 11 hours on a laptop computer with a 266 Mhz G3 (Macintosh) processor. 

It should be noted that the learning sample of housecats breaks the symmetry in the posterior with 
respect to permutations on the labels for the two components ({\em F.~sylvestris} and housecats) in the model.
Thus, there is not a substantial label-switching \cite{Stephens2000JRSS} problem in this case.  

\subsection{Results for model $\Mpa$}
The posterior mean estimate of $\xi_\mathrm{P}$, the proportion of purebred cats, is .65, with a 90\% credible set
spanning the range from .47 to .79.  The MCMC estimate of the posterior density of $\xi_\mathrm{P}$ is given in
Figure~\ref{fig:xipost}.
%%%%%%%%%%  FOUR PANEL FIGURE
\begin{figure}
\begin{center}
\subfigure[$P(\xi_\mathrm{P} |
y)$\label{fig:xipost}]{\includegraphics*[width=.32\textwidth]{XiPureDensOneLine.eps} }
\hspace*{.05\textwidth}
\subfigure[$P(\pi|y)$\label{fig:pipost}]{\includegraphics*[width=.32\textwidth]{PiDensOneLine.eps} } 
\\
\subfigure[$P(\alpha|y)$\label{fig:alphapost}]{\includegraphics*[width=.32\textwidth]{JumpAlphaDensOneLine.eps} }
\hspace*{.05\textwidth}
\subfigure[$\alpha$~vs.~$\xi_\mathrm{P}$\label{fig:alphaversusxi}]{\includegraphics*[width=.32\textwidth]{AlphaXiScatter.eps}
} 
\end{center}
\caption[Aggregate-level parameters for the Scottish cat 
dataset]{Graphical summaries of the aggregate-level parameters for the Scottish cat 
dataset.  (a--c) are unsmoothed posterior density estimates taken by scaling histograms 
with bin widths of 0.01 in (a) and (b) and 0.03 in (c):  (a) proportion of
purebred cats,
$\xi_\mathrm{P}$ in the population from which the cats were sampled, (b)
proportion of cats from the {\em F.~sylvestris} subpopulation, conditional on being
purebred, (c) the parameter $\alpha$.  (d) a scatterplot of 5,000 pairs
$(\alpha,\xi_\mathrm{P})$ sampled from the Markov chain.  The two parameters are clearly
correlated.  Lower values of $\alpha$ correspond to lower values of $\xi_\mathrm{P}$, as
expected. }
\label{fig:Mixfourpanel}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%
  The
distribution is long-tailed to the left.  These low values of $\xi_\mathrm{P}$ coincide with low values
of the parameter $\alpha$ (Figure~\ref{fig:alphaversusxi}).  This correlation is expected; when $\alpha$ is low,
then admixed individuals are expected to have admixture proportions near to zero or one, and hence the ability to
distinguish between admixed and purebred individuals is diminished.  The estimated posterior density for $\alpha$
itself is shown in Figure~\ref{fig:alphapost}.  It has a peak around 0.7, and tapers off with larger values,
but it is still rather high at the upper bound imposed on it of 3.  The choice of $A=3$ is clearly a choice of
prior to which the final result will be sensitive.  A larger $A$ would reduce the
posterior probability for low values of  $\xi_\mathrm{P}$, reducing the skewness of the posterior for
$\xi_\mathrm{P}$ and increasing its posterior mean estimate.  This issue will be taken up
again in the Discussion. 

Figure~\ref{fig:pipost} gives the estimated posterior density for the probability that a cat is {\em
F.~sylvestris} conditional on its being of purebred type.  The posterior mean is .83 with a 90\% credible
interval from .73 to .94.   This suggests that a
large proportion ($>60\%$) of the wild-living cats in Scotland are purebred {\em F.~sylvestris}.  On the other
hand, there is evidence that between 21\% and 53\% of the wild-living cats are admixed individuals with 
ancestry from both {\em F.~sylvestris} and domestic cats.  Further, it cannot be ruled out
that some wild-living cats are pure housecats that have gone feral. 

Figure~\ref{fig:purefork} summarizes the results for individual cats.
%%%%%%%%  PURE FORK FIGURE
\begin{figure}
\begin{center}
\includegraphics*[width=\textwidth]{PureFork.eps}
\end{center}
\caption[$P(v_i=\mathrm{P})$ versus $P(Z_i =\mathrm{``}\mathit{ F.~sylvestris}\mathrm{"}|v_i=\mathrm{P})$
for the Scottish cat data]{A plot of posterior mean estimates for 
$P(v_i=\mathrm{P})$ on the horizontal axis against
$P(Z_i =\mathrm{``}\mathit{ F.~sylvestris}\mathrm{"}|v_i=\mathrm{P})$ on the vertical
axis.  Each open circle represents one of the 230 wild-living cats in the
sample.  The cluster in the upper right includes 120 individuals all with posterior
probability of being purebred {\em F.~sylvestris} greater than $.80$.  At far left are
some eight individuals with high posterior probability of being admixed.  For cats with
intermediate estimates of $P(v_i=\mathrm{P}|y)$, the credible sets tend to be quite
wide (not shown).    }
\label{fig:purefork}
\end{figure}
%%%%%%%%%%%%
  On the horizontal axis, is the posterior probability
of being purebred.  On the vertical axis is the posterior probability of being {\em
F.~sylvestris} conditional on being purebred.  A cluster in the upper right represents 102 of the cats in the
sample, all with posterior probability of being pure greater than .80.  Given that these cats are pure, they have
posterior probability close to one of being {\em F.~sylvestris}.  Also evident is a small cluster of cats with
$P(v_i=\mathrm{P}|y) > .65$ but which, if they are purebred cats, are almost certainly domestic cats.  At the
other end of the scale are several cats with very high probability of being admixed.  


\subsection{Comparison of results for models $\Mpa$ and $\Ma$\label{sec:comparison}}
For parameters shared by $\Ma$ and $\Mpa$, the estimates
differ between models most for
$\alpha$.  Under $\Ma$,  $\alpha$ is much smaller, so as to accommodate the
purebred cats as admixed individuals with admixture proportions close to 0 or 1 (Figure~\ref{fig:AllAlphaDens}).
%%%%%  ALL ALPHA DENS
\begin{figure}
\begin{center}
\includegraphics*[width=.45\textwidth]{AllAlphaDens.eps}
\end{center}
\caption[Posterior densities for $\alpha$ from the Scottish cat data]{ Posterior densities for $\alpha$
from the Scottish cat data.  Broken line shows
$P(\alpha|y)$ under model~$\Ma$. Solid line shows
$P(\alpha|y)$ under model~$\Mpa$. 
It appears that under $\Ma$, most of the information constraining values
of $\alpha$ in the posterior comes from individuals of pure, or mostly
pure origin.   }
\label{fig:AllAlphaDens}
\end{figure}
Additionally, a significantly smaller proportion of the alleles in the sample get
allocated to the housecat population under $\Mpa$ than under $\Ma$.  Histograms of the proportion of alleles allocated to the
housecat subpopulation for $\Mpa$ and
$\Ma$ are shown in Figure~\ref{fig:PoA}.   
\begin{figure}
\begin{center}
\includegraphics*[width=.45\textwidth]{PoA100s.eps}
\end{center}
\caption[Proportion of alleles allocated to the housecat
subpopulation]{Comparison of the proportion of all the alleles in the sample allocated to the housecat
subpopulation.  Solid line and circles are a histogram from $\Mpa$; broken line and open circles are a
histogram from $\Ma$.  There is little overlap between the two.  A higher proportion of alleles is
allocated to the housecat population under $\Ma$.}
\label{fig:PoA}
\end{figure}
MAYBE OUGHT TO PUT SOMETHING IN HERE ABOUT PARTICULAR INDIVIDUAL CATS.
THERE IS A FIGURE HERE AT LEAST FOR INDIVIDUAL NUMBER 4.  MUST ADD A BIT MORE IN THE TEXT ABOUT 
THAT.



%%%%%%%  INDIV 4 Q HISTS
\begin{figure}
\begin{center}
\mbox{
\subfigure[From model $\Mpa$]{\includegraphics*[width=.4\textwidth]{Indiv4JumpQHist.eps}} 
\hspace*{.05\textwidth}
\subfigure[From model $\Ma$]{\includegraphics*[width=.4\textwidth]{Indiv4AdmixQHist.eps}}
}
\end{center}
\caption[Proportion of genome from the {\em
F.~sylvestris} subpopulation for cat \#4]{Histograms for the proportion of genome from the {\em
F.~sylvestris} subpopulation for cat \#4 from the Scottish Cat dataset.
(a) Under $\Mpa$ the proportion of genome includes the weighted average of that  which would be assigned
to {\em
F.~sylvestris} given  that it was purebred and the histogram of $q_4$ conditional on $v_4 = \mathrm{A}$.
(b)
Under $\Ma$ this is the histogram of $q_4$.    Note that $\Mpa$ allows
for a mode  in the posterior distribution of $q_4$ which is in the center of its range.      }
\end{figure}



\section{Bayesian Model Comparison} 
Once able to entertain the model $\Mpa$, it is natural
to ask whether that expanded model has gained us anything.  One way to pose the question is to ask
whether the data provide more support for $\Mpa$ than for the model we will call $\Ma$ which requires all
individuals to be admixed and governed by a single $\alpha$.  A rough estimate of the Bayes factor, 
$B = P(y|\Mpa) /
P(y|\Ma)$,   might be obtained by observing the proportion of time the Markov chain defined under $\Mpa$ spends in
states with zero or almost zero individuals allocated to the purebred group (since retricting
$\xi_\mathrm{P}$ to zero in
$\Mpa$ essentially gives
$\Ma$).  However, this is unsatisfactory as there is no prior probability mass
on the point $\xi_\mathrm{P} = 0$.  Furthermore, the chain may visit states with low
$\xi_\mathrm{P}$ so infrequently, that it is impossible to get a good estimate of $B$
that way.  

\citeN{Gelfandetal1992} suggest approximating $B$ by the ``pseudo-Bayes factor" formed as the
product over all observations of the ratio of cross-validitation predictive densities 
under the two models.  The cross-validation predictive density for the $i\thh$ individual,
may be approximated by the harmonic mean of the values 
$P(y_i|\alpha_s,\theta_s)$ under $\Ma$ and the values $P(y_i|\alpha_s,\pi_s,\theta_s,\xi_{\mathrm{P}_s})$, 
computed as the denominator in (\ref{eq:vfullcond}), under $\Mpa$, where $s$ subscripts the states visited by the
chain over which the harmonic mean is taken.  
\citeN{Raftery1992} cautions that 
the pseudo-Bayes factor, being akin to a pseudo-likelihood, may be an inaccurate approximation
to the Bayes factor and should not be used for model comparison 
if the latter is available.  However as discussed by \citeN{Pritchardetal2000}, it is 
difficult to reliably estimate the marginal likelihood $P(y|\Ma)$, and the same 
is true for $P(y|\Mpa)$, making computation of the Bayes factor by that route challenging.



As an alternative, I have developed a reversible-jump MCMC scheme  \cite{Green1995} for
computing the Bayes factor $B = P(y|\Mpa)/P(y|\Ma)$.  While reversible jump methods have recently received
widespread attention for sampling over numerous models in complex model spaces
\cite{Rue&Hur1999,Del&For1999,Giu&Gre1999}, it seems they have been used less often when a small number
of closely-related models are being considered, as in the present case.  The posterior distributions  estimated
from separate runs under
$\Mpa$ and
$\Ma$ can guide us in devising reversible-jump proposals that are easy to implement and offer a good estimate of
$B$. Details appear in Appendix~\ref{sec:RJ}. This circumvents the potential problem of instability in a direct,
sampling-based estimate of $P(y|\Mpa)$ or
$P(y|\Ma)$, and affords an opportunity to compare the pseudo-Bayes factor
to the full Bayes factor in the comparison of two complex, hierarchical models. 

\section{Reversible Jump MCMC for Model Comparison}
\label{sec:RJ}
We may compute the Bayes factor, $B$, by reversible jump MCMC \cite{Green1995}.  This method
allows for the construction of a Markov chain that may jump between state spaces of varying dimension.  In our
case we construct a chain which takes values in two spaces indexed by $m=1$ or $2$.  If $m=1$ then the
chain is currently in the space associated with model $\Ma$, and it moves to new values within that space
as described in Section~\ref{sec:Adm}.  If $m=2$, then the chain is currently in the state space associated with
model $\Mpa$, and it moves to new values within that space as described in Section~\ref{sec:Simul}.  Since $\Mpa$
includes the variables $\xi$ and $\pi$  ($\xi$ has one degree of freedom, $\xi_\mathrm{P}$, and, in the case of
$J=2$,
$\pi$ has one degree of freedom as well) which are absent in model $\Ma$, there are two extra degrees of
freedom when
$m=2$.  For this reason, reversible-jump moves are required to move from $m=1$ to $m=2$.  The formulation of
these moves is such that detailed balance is satisfied, ensuring that the proportion of time the chain spends
with $m=1$ converges to
$P(\Ma|y)$ as the chain is run for infinite time, and so, for a run of the chain of length $n$, the quantity 
\begin{equation}
\frac{\sum_{i=1}^n  1\{m_i = 1\}}
{ \sum_{i=1}^n  1\{m_i = 2\}   } 
\label{eq:posteriorodds}
\end{equation}
estimates the posterior odds, which, upon division by the prior odds, gives $B$.

For a reversible jump move from $m=2$ to $m=1$ we leave $\theta$ unchanged and  propose a new value for $\alpha$,
say $\alpha^\prime$, that is a deterministic, many-to-one, function $g$ of the current values of $\alpha$,
$\xi_\mathrm{P}$, and $\pi$.  We are at liberty to choose any appropriate and suitable $g$.  For the Scottish cat
problem, by examining the posterior distribution of $\xi_\mathrm{P}$, $\pi$, and $\alpha$ under $\Mpa$, and by
surmising that high values of
$\xi_\mathrm{P}$ in $\Mpa$ should correspond to low values of $\alpha^\prime$ in $\Ma$, I empirically chose
\begin{equation}
\alpha^\prime = g(\alpha,\xi,\pi) = 0.0925 + 0.13638 \alpha - 0.21\sin^{-1}( \xi_\mathrm{P}^2).
\end{equation}
In this case, $\sin^{-1}(\xi_\mathrm{P}^2)$ was chosen, since that transformed variable has a simpler (i.e. more linear)
relationship with 
$\alpha$ than does 
$\xi_\mathrm{P}$, itself, in the MCMC output from $\Mpa$ (see Figure~\ref{fig:alphaversusxi} in the main text). 
Notice that
$\pi$ does not actually appear in the function
$g$, since this simplifies the Jacobian, and it does not seem essential (i.e. there is not large correlation
between $\alpha$ and $\pi$).    Taking the 5000 pairs
$(\alpha,\xi_\mathrm{P})$ that were plotted in Figure~\ref{fig:alphaversusxi} and applying $g$ to them gives
values of $\alpha^\prime$ summarized by their histogram in Figure~\ref{fig:alphaprimehist}.  This histogram
resembles the posterior distribution of $\alpha$ under $\Ma$ (broken line in Figure~\ref{fig:AllAlphaDens}), as
desired.  


%%%%%%    ALPHAPRIME HISTOGRAM AND PROPOSAL DENSITY FIGURES
\begin{figure}
\begin{center}
\subfigure[
$\alpha^\prime$ histogram\label{fig:alphaprimehist}]{\includegraphics*[width=.30\textwidth]{AlphaPrimeHist.eps}}
\\
\subfigure[$f_\xi(\xi_\mathrm{P})$ \&
$f_\pi(\pi_0)$\label{fig:fdens}]{\includegraphics*[width=.30\textwidth]{XiPiProposalDens.eps}}
\end{center}
\caption[Proposal densities for reversible-jump moves]{Elements of the reversible-jump proposals. 
(a)~Histogram of values of
$\alpha^\prime$ computed as $g(\xi_\mathrm{P},\alpha,\pi)$ for 5,000 points visited by a Markov chain run
under $\Mpa$. The function
$g$ was chosen so that this histogram would be similar to the posterior density for $\alpha$ under $\Ma$, shown in
Figure~\ref{fig:AllAlphaDens}.
(b)~The proposal densities: solid line is $f_\xi(\xi_\mathrm{P})$; broken line is  $f_\pi(\pi_0)$ 
These were chosen to represent overdispersed versions of $P(\xi_P|y)$ and $P(\pi_0|y)$ under $\Mpa$, shown 
in Figures~\ref{fig:xipost} and~\ref{fig:pipost}.
 }
\end{figure}

To propose the reverse move from $m=1$ with a current value $\alpha^\prime$ to $m=2$ with proposed values for the
parameters
$\alpha$, $\xi_\mathrm{P}$ and $\pi$ requires simulating new values for $\xi_\mathrm{P}$ and $\pi$ from known
densities and then using those values and the inverse of the function $g$ to determine what value of $\alpha$
shall be proposed.  The known densities were chosen to approximate the posterior density estimates of
$\xi_\mathrm{P}$ and $\pi$ under $\Mpa$.  Letting
$\pi_0$ denote the proportion of purebred cats that are {\em F.~sylvestris}, the densities used were
$f_\xi(\xi_\mathrm{P})\equiv \mathrm{Beta}(8,4)$ and $f_\pi(\pi_0) \equiv .8\mathrm{Beta}(30,9) +
.2\mathrm{Beta}(15,2)$.  These densities are shown in Figure~\ref{fig:fdens}.  Comparison to Figures~\ref{fig:xipost}
and~\ref{fig:pipost} shows that they resemble overdispersed versions of $P(\xi_\mathrm{P}|y)$ and $P(\pi_0|y)$. With values of
$\xi_\mathrm{P}$ and
$\pi_0$ drawn from these densities,
$\alpha$ is determined by 
\[
\alpha = g^{-1}(\alpha^\prime,\xi_\mathrm{P},\pi_0) = \frac{0.21\sin^{-1}( \xi_\mathrm{P}^2) + \alpha^\prime -
.0925}{.13638}.
\]

We may propose a reversible jump move at the end of each sweep.  Thus, if $m=1$, after a sweep updating all the
variables associated with $\Ma$, we propose a jump up to $m=2$.  If $m=2$, then after a sweep updating all the
variables associated with $\Mpa$ we propose a jump down to
$m=1$.  Under such a scheme, the acceptance probability for a proposed move from $m=1$ with
$\alpha=\alpha^\prime$ to
$m=2$ and parameter values $(\alpha,\xi_\mathrm{P},\pi_0)$, is $\min\{1,\rjaccept\}$, with  $\rjaccept$ reducing
to 
\begin{eqnarray}
\rjaccept &=& \frac{P(\Mpa)}{P(\Ma)} \times
\frac{P(\alpha|\Mpa)}{P(\alpha^\prime|\Ma)} \times
\frac{P(\xi_\mathrm{P}|\Mpa) P(\pi_0|\Mpa)}{f_\xi(\xi_\mathrm{P}) f_\pi(\pi_0)}  \nonumber \\
\rule{0pt}{5ex}
 & \times &
\frac{\prod_{i=1}^N [\xi_\mathrm{P} P(y_i|\pi,\theta) + \xi_\mathrm{A}
P(y_i|\alpha,\theta)]}
{\prod_{i=1}^N P(y_i|\alpha^\prime,\theta)}  \times
\frac{1}{0.13638}
\end{eqnarray}
where $P(\cdot|M)$ denotes prior densities for parameters under different models $M$.  If proposing a move down
from $m=2$ with current values $(\alpha,\xi_\mathrm{P},\pi_0)$ to $m=1$ with $\alpha=\alpha^\prime$, the
acceptance probability is $\min\{1,\rjaccept^{-1}$\}.  The factor of $(0.13638)^{-1}$ is the Jacobian from the transformation
$g$.

Figure~\ref{fig:aprobMa} shows a trace of $\log\rjaccept$ from a chain forced to stay in $m=1$ (i.e. it makes
proposals to $m=2$ but is not allowed to accept them) using the Scottish cat data with learning samples,
and assuming prior odds for the models 
$P(\Mpa)/P(\Ma) = 1$.  Figure~\ref{fig:aprobMpa} shows a similar trace of
$\log\rjaccept^{-1}$ for a chain restricted to $m=2$.   
\begin{figure}
\begin{center}
\mbox{
\subfigure[Trace of $\log\rjaccept$ under $\Ma$\label{fig:aprobMa}]{\includegraphics*[width=.35\textwidth]{aprobsMa.eps} } 
\hspace*{.05\textwidth}
\subfigure[Trace of $-\log\rjaccept$ under $\Mpa$\label{fig:aprobMpa}]{\includegraphics*[width=.35\textwidth]{aprobsMpa.eps} }
}
\end{center}
\caption[Values of $\log\rjaccept$ under $\Ma$ and $\Mpa$] {With prior odds of 1, (a)~a trace of values
of
$\log\rjaccept$ plotted as unconnected points for 10,000 sweeps of a chain with
$m$ fixed at 1 (model $\Ma$).  The majority of points lie above 5, indicating that most proposals to move
to $\Mpa$ from $\Ma$ by the proposed reversible-jump move would be accepted.  Note that some values of
$\log\rjaccept$ are greater than 20.  So, even with prior log-odds of $\log [P(\Mpa)/P(\Ma)] \approx -20$,
proposals from $m=1$ to $m=2$ will occasionally be accepted.  (b)~a trace of $-\log\rjaccept$ for 10,000 sweeps
of a chain restricted to $m=2$.  Many values are less than  $-20$.  However, again, with $\log [P(\Mpa)/P(\Ma)]
\approx -20$ proposals from $m=2$ to $m=1$ will be occasionally accepted. } 
\end{figure}

It is apparent from these traces that, without imposing strong
prior support for $\Ma$, it is unlikely that a chain in $m=2$ would ever move to $m=1$.   Thus, I made three
different runs with prior log-odds,  $\log[P(\Mpa)/P(\Ma)]$, equal to $-19$, $-20$, and $-21$.  From each of
these runs, I estimated the posterior log adds by taking the log of (\ref{eq:posteriorodds}).  The value of the
posterior log-odds calculated as the average over ten chains started from overdispersed starting points as a
function of sweep number is shown for the three different prior odds in Figure~\ref{fig:PostLogOddsTraces}. 
Though the chains may not have been run long enough for an extremely precise estimate of the posterior log-odds,
an order-of-magnitude estimate can clearly be made.  Subtracting the prior log-odds from the estimated posterior
log-odds gives, for each of the three different prior odds used, an estimate of $\approx 20.3$ for the log
of the Bayes factor.  Hence, $2\log B > 40$, indicating overwhelming support in the data for model $\Mpa$ over
$\Ma$.
The log of the pseudo-Bayes factor is 12.3 which, quite notably, differs by eight 
from the true $\log B$.  








%%% MDL CMP TRACES FIGURE
\begin{figure}
\begin{center}
\includegraphics*[width = .55\textwidth]{PostLogOddsTraces.eps}
\caption[Traces of estimated posterior log-odds,  $\log[P(\Mpa|y)/P(\Ma|y)$]{Estimated posterior
log-odds, 
$\log[P(\Mpa|y)/P(\Ma|y)]$ for three different prior log-odds.  Each line shows the estimate as a
function of number of sweeps.  Ten parallel chains started from overdispersed points were used for each
estimate.  Hence the estimates shown are the $\log$ of the expression in (\ref{eq:posteriorodds}) for
$n=10$ times the number of sweeps.     The numbers at the right of each line are the prior log-odds
assumed.  $\Mpa$ is so highly favored, that in order to get the reversible-jump sampler to mix, huge
prior weight must be given to $\Ma$.  The log of the Bayes factor,
$\log B$, may be calculated by subtracting the prior log-odds from the estimated posterior log-odds.  For
all three values of the prior odds this gives about 20.3. }
\label{fig:PostLogOddsTraces}
\end{center}
\end{figure}



\section{Discussion}
In applications to conservation biology and ecology, populations of interest may be pure mixtures of
two subpopulations, or they may contain admixed individuals from two originally separate
subpopulations.  Genetic data, in conjunction with statistical models of genetic mixture and 
admixture have been useful for clustering individuals and genes from different subpopulations.
This paper presents a novel application of the ``filter-forward-simulate-backward"
algorithm akin to the computations presented in \citeN{Baumetal1970} to the population-admixture model of
\citeN{Pritchardetal2000}.  This computation makes it possible to expand that model to one
that allows for individuals to be either purebred or admixed. Such an expanded model ($\Mpa$) is vastly more
supported by the Scottish cat data than one including admixture only.  It is likely that $\Mpa$ will 
fit other datasets much better as well, because samples from recently admixed populations
will typically include some purebred individuals.  

While the dramatic improvement of model fit is encouraging, it also raises some issues that
bear further investigation.  The first of these is that $\Mpa$ may fit the data better not
simply because it allows separate classes of purebred and admixed individuals.  It may be
that a great deal of improvement comes from including the parameter $\pi$ which allows different
contributions of pure cats from the two subpopulations.  This contrasts to the formulation in 
$\Ma$ where, due to the symmetry of the $\mathrm{Beta}(\alpha,\alpha)$ prior for the $q_i$'s, the
marginal probabilities are equal that any gene copy is from the housecat or the {\em F.~sylvestris} subpopulation.
That is to say, under $\Ma$, $P(w_{it} = 1|\alpha) = P(w_{it} = 2|\alpha) = .5$ for all $i,t$, and $\alpha$.  
By contrast, under $\Mpa$, for different values of $\xi,\pi,$ and $\alpha$, the marginal probability that any gene
copy is from the housecat population is not constrained to equal the marginal probability that it
is from the {\em F.~sylvestris} population.    
The symmetry imposed by $\Ma$ might explain some systematic biases for estimates of $q_i$
that \citeN{Pritchardetal2000AJHG}
report for simulated data with unequal admixture proportions.  Furthermore, the
issue may have implications for model $\Mpa$.  If the population
admixture proportions depart from $.5$, then $P(y_i|\alpha,\theta)$ might be inflated for
individuals with large amounts of ancestry from the lesser-represented subpopulation, and
deflated for individuals with more ancestry from the greater-represented subpopulation.  For
this reason, in the Scottish cat problem, one might expect that the posterior probability of
being a purebred individual will be overestimated for cats that resemble {\em F.~sylvestris} and
underestimated for cats that appear to be housecats.  This may also induce some bias in the posterior estimate of
$\xi_\mathrm{P}$.  All this suggests that a fruitful extension to the model $\Ma$
of \citeN{Pritchardetal2000} would be to allow population-specific $\alpha$'s.  For example, 
in the case of $J$ subpopulations, $q_i\sim\Dir(\alpha_1,\ldots,\alpha_J)$.


The results also suggest that estimation in $\Mpa$ may be sensitive to the upper bound, $A$, chosen for $\alpha$.
Had $A$ been chosen greater than three, then values of $\alpha > 3$ would surely have been visited in the 
MCMC simulation, and the resulting estimate for $\xi_\mathrm{P}$ would have been somewhat larger, since 
$\alpha$ and $\xi_\mathrm{P}$ are positively correlated.  This is observed in a separate run made with $A=10$---the chain
visits values of $\alpha$ between 3 and 10 quite frequently.  In fact, the
estimated posterior density 
for $\alpha$ decreases only slightly between 3 and 10.  However, the effect on the other parameters is not overwhelming.
For example, with $A=10$, the posterior mean (90\% credible interval) for $\xi_\mathrm{P}$ was 
.71 (.53,~.82), as opposed to .65 (.47,~.79) with $A=3$. 
It is interesting that  the choice of
$A$ has little effect in the poorer-fitting model $\Ma$, because that model tries to fit purebred cats as admixed individuals. 
This  keeps $\alpha$ low regardless of $A$.  Under $\Mpa$, however, once the purebred
individuals are removed from the admixed class there is little information left for estimating
$\alpha$.  So, paradoxically, to use the better-fitting model $\Mpa$ requires imposing more
prior information.  In the case of $A$, however, biological knowledge can guide the choice.

I chose $A=3$ because, with only two
subpopulations, large values of
$\alpha$ indicate that admixed individuals carry close to half of their ancestry from one subpopulation and half from 
the other.  In a population like this, the most plausible explanation for such a pattern would be that the
admixed individuals were all first-generation (F1) hybrids between individuals from the two
subpopulations.  If this is the case, then, at each locus, an admixed individual will carry exactly
one allele from one subpopulation, with the other allele coming from the other subpopulation. 
This condition can be used to compute the posterior probability that an individual in the sample is an
F1 hybrid.  I leave the description of such a procedure to a separate
paper, however, I have found that  none of the individuals in the sample had  posterior probability 
greater than .5 of being an F1 hybrid.  In fact, for all but seven of the individuals, the
posterior probability of being an F1 hybrid was below .10.   For this reason, it seemed
implausible that
$\alpha$ should be allowed to range past about three.

The Bayesian model comparison revealed that $\Mpa$ is a much better model for the Scottish cat data. 
Computing Bayes factors in these models is often difficult because calculating the marginal likelihood
can require a difficult computation of an unkown normalizing constant.  Rather than directly computing the 
marginal likelihood, Appendix~\ref{sec:RJ} gives an example of how approximations to posterior densities of several
parameters in  different models may be used to formulate reversible-jump moves between a small set of
closely-related models.  This gives us a good approximation to the Bayes factor.  In turn, that allows us to compare the
true Bayes factor to the pseudo-Bayes factor (a product of ratios of cross-validation predictive densities). The pseudo-Bayes
factor has been advocated as a computationally manageable  approximation to the Bayes factor.   While cross-validation and
predictive densities offer  a fine level of detail for exploring which observations, in particular, are poorly fit by a model,
their use in overall model comparison via the pseudo-Bayes factor should be done with reservation.
In the current problem, the pseudo-Bayes factor provided a poor approximation of the Bayes factor.  

Quite apart from genetic mixtures, the forward-backward computation here may be useful in more
general mixture problems.  Sometimes,  the Gibbs sampler mixes poorly in the Bayesian analysis of
mixtures. \citeN{Robert1996} describes this in terms of {\em trapping states} in finite normal mixtures:
when only one or a few observations are allocated to a component, the parameters for that component fit
the few observations so tightly that few if any of the other observations would likely get allocated to
the component.  Reparametrizing the normal mixture model, as done by \citeN{Men&Rob1995},  corrects the
problem by keeping the component-specific parameters from fitting the observations in a near-empty component
too tightly.  However, this does not address trapping states that may occur simply because the mixing proportion
for a component becomes small.  If the mixing proportion of a component happens to reach a value near zero, then 
the probability of allocating any observations to that component will also be small, and the component may remain
empty through many iterations of the chain.


The block-updating scheme of Section~\ref{sec:Polya} can provide a useful Gibbs move that could
be executed to restore empty components, by the following rationale:  in a $J$-component finite mixture
with a $\Dir(\zeta_1,\ldots,\zeta_J)$ prior on the mixing proportions, the latent allocation variables, $Z_i$,
marginally follow a labelled compound multinomial-Dirichlet distribution. Consequently, conditional
on current values of all the $Z_i$'s, the subset of those having any two values  will follow a labelled beta-binomial 
distribution \cite{Johnsonetal1997}.  That is, the marginal distribution of $\{Z_i : Z_i~=~j_a~\cup~Z_i~=~j_b,j_a\neq j_b, i =
1,\ldots,N\}$ follow a labelled beta-binomial distribution with parameters
$(\zeta_{j_a},\zeta_{j_b})$.  Thus, the methods of Section~\ref{sec:Polya} could be applied to redistribute elements
amongst the two  components $j_a$ and $j_b$, having marginalized over the mixing
proportions $\pi_{j_a}$ and $\pi_{j_b}$.  And so, observations may be reallocated to component $j_a$ (or $j_b$), according
to their full conditional distributions, even if $\pi_{j_a}$ (or $\pi_{j_b}$) is close to zero.  

\section*{Acknowledgments}
This work developed out of conversations during a brief visit
to University of Oxford, where the author was hosted by the lab of Peter Donnelly in the 
Department of Statistics. Partial support came from
National Science Foundation Grant BIR--9807747.  The author thanks Jonathan Pritchard, Elizabeth Thompson, and
Matthew Stephens for helpful discussion and comments on earlier drafts. Special thanks to Mark Beaumont for providing
the data on Scottish wildcats.



