\chapter{Explicit Modeling of the Hybridization Process}
\label{chap:Hybrids}

The previous chapter described an extension to the model proposed by 
\citeN{Pritchardetal2000} for individuals
from structured populations.  The \citeANP{Pritchardetal2000} model is a
mathematically convenient model which seems to apply well to many
situations of population structure.  However, when more information is
available about the nature of the admixing process, then a more detailed
analysis of the situation is possible by using a more detailed model
of the process.  This chapter describes such a detailed model for a
commonly encountered sitatuation in conservation biology---the case of
two species which are known to have been hybridizing for a limited number
of generations or in which the hybrids have fitness which is reduced to
the extent that all of the admixed individuals in regions of overlap
between the species are the result of recent hybridization.   The model
applied is one in which individuals belong to one of many different
hybrid categories (\eg $\F_1$, $\F_2$, and various backcrosses), and
inference will be done in a similar way to that in the previous
chapter.   Using Markov chain Monte Carlo (MCMC) we compute
for each individual the posterior
probability  of inclusion in a particular hybrid
category, and simultaneously compute the posterior distribution of the
allele frequencies in the populations under study of the two different
species. 

 I present these methods in the context of
hybridization between sympatric populations of a species~\A and a
species~\B, and develop the probability model  that arises from
explicitly modeling the hybridization process and using data on multiple,
unlinked markers.  This
probability is the product of probabilities for single-locus genotypes
conditional on the hybrid category.  There is a simple calculus for these
conditional probabilities which is described below.     With the model
thus specified, I describe how MCMC proceeds in similar fashion to the
previous chapter. 

The main goals of this chapter are to derive the method, apply 
it to real and simulated data, and then describe several
modifications and extensions suggested by the results of these data analyses.  
I apply  the method to real genetic data 
(a sample of juvenile salmonids containing rainbow trout ({\em Oncorhynchus mykss}), cutthroat
trout ({\em O.~clarki}), and their hybrids) and then to simulated data with many relatively
uninformative markers and again with fewer very informative markers.
The results of these analyses, in Section~\ref{sec:CuttResults}, suggest that
distinguishing between different hybrid categories beyond the $\F_1$ stage
 requires a substantial
amount of data, and clear genetic differentiation between the species.
They also underscore the importance of having an established sampling scheme
when trying to detect hybrids.  These issues suggest further extensions to the
method which are beyond the current scope of the thesis.  

\section{Population and Probability Model} 
We consider a group  of individuals in the wild which consists of
sympatric populations of two species~\A and~\B, and hybrids of the two
species which have occurred from $n$ potential generations of
interbreeding.  We take $n$ to be known or assumed.  We also assume that
we have a sample of $M$ individuals drawn from this group for genetic
analysis.  For now, we shall assume that individuals are sampled randomly
and independently of whether they are of species~\A, or  species~\B, or
are a hybrid of the two species.  This sort of sampling would arise,
if, for example, juveniles were sampled and were very difficult to
distinguish on the basis of morphology between the two species or hybrids
thereof. In Section~\ref{sec:CuttDiscussion},  I briefly discuss a modeling/sampling approach for
relaxing that assumption.  We have genotype information from
$L$ unlinked loci on the individuals sampled.  Let the
$\ell\thh$ locus possess
$K_\ell$ alleles detected in the sample.  We denote the allele
frequencies in species~\A and~\B, $n$ generations ago, by $\bTheta_\A$
and
$\bTheta_\B$, respectively.  Each of these $\bTheta$'s is a collection of
vectors, with each vector giving the allele frequencies at a particular
locus.  For example, for species~\A, $\bTheta_\A =
(\bttheta_{\A,1},\ldots,\bttheta_{\A,L})$, where $\bttheta_{\A,\ell} =
(\theta_{\A,\ell,1},\ldots,\theta_{\A,\ell,K_\ell})$ are allele
frequencies at the $\ell\thh$ locus.  The alleles found in individuals
from species~\A and species~\B are assumed to be drawn randomly from the
allele frequencies
$\bTheta_\A$ and
$\bTheta_\B$ respectively,
$n$ generations ago.  Likewise, individuals $n$ generations before sampling are
assumed to be in Hardy-Weinberg and linkage equilibrium with reference
to their contemporaneous conspecifics.  

Each individual in the sample is genotyped at $L$ loci.  The
gene copies carried at any locus are considered to be ordered,
though that order may be arbitrary; for example it may merely be
the order in which the genetic data on that locus in that individual
happened to be recorded.  This ordering arises because we will shortly
introduce latent data that applies to each particular gene copy.  The
allelic types of the two gene copies at locus $\ell$ in individual $i$
are denoted by $\btY_{i,\ell} = (Y_{i,\ell,1},Y_{i,\ell,2})$, with each
of $Y_{i,\ell,1}$ and $Y_{i,\ell,2}$ taking an integer value between 1
and $K_\ell$, inclusive, corresponding to the possible allelic types at
the $\ell\thh$ locus.  The
$L$ single-locus genotypes in the $i\thh$ individual are denoted by
$\bY_i = (\btY_{i,1},\ldots,\btY_{i,L})$, and all of the genetic data,
over all
$M$ individuals in the sample is $\bY = (\bY_1,\ldots,\bY_M)$.  We do not
know from which species each of an individual's gene copies descended,
but we denote that unknown information by the latent variable
$\btW_{i,\ell} = (W_{i,\ell,1}, W_{i,\ell,2})$.  $W_{i,\ell,1}$ takes
the value 0 if the first gene copy at the $\ell\thh$ locus of the
$i\thh$ individual originated from the species~\A population, and it takes
the value 1 if that gene copy originated from species~\B.  $W_{i,\ell,2}$
takes values analogously, depending on the origin of the second gene
copy.  We use $\bW_{i} = (\btW_{i,1},\ldots,\btW_{i,L})$ to denote all the
latent gene origin indicators in the $i\thh$ individual, and $\bW$
denotes the latent gene origin indicators in all the individuals. 


 To
develop the probability for the observed data, $\bY$, we must
explicitly state what is meant by a  ``hybrid category." For unlinked
markers such a definition is given in the following subsection in terms of
``genotype frequency classes." 

\subsection{Hybrid categories}
When hybridization between two species has been potentially occurring for
$n$ generations, the possible hybrid categories into which  an individual
may fall can be enumerated and described by considering the possible
arrangements of different species amongst the founders in an
$n$-generational pedigree.  The individual of
interest is taken to be the member at the bottom of the pedigree, and is
assumed to be non-inbred over the last $n$ generations; hence we assume
there are no loops in its $n$-generational pedigree.
Figure~\ref{fig:n=2ped} illustrates this for the case of
$n=2$.
\begin{figure}
\begin{center}
\includegraphics*[width=\textwidth]{2peds.eps}
\caption[Hybrid classes with $n=2$]{Six arrangements of founders on
pedigrees of $n=2$ two generations each.  Each box represents a locus. 
The circles within each box represent the two gene copies possessed by the
diploid organism at the locus. The founders are the individuals in the
top row of each pedigree.  Black gene copies amongst founders are those
originating from the Species~$A$ population, and the white gene copies
are from Species~$B$.   The individual at the bottom of each pedigree
belongs to a different hybrid category, determined by the arrangement of
species amongst the founders.  (a) through (f) represent six distinct
{\em hybrid classes}.  (a) through (f) also represent six distinct {\em
genotype frequency classes}.   There are, however, only five distinct
{\em gene frequency classes}; the individuals at the bottoms of pedigrees
(c) and (f) are both in the same gene frequency class. }
\label{fig:n=2ped}
\end{center}
\end{figure}
Depending on the type of genetic data available (\eg linked versus
unlinked loci), it may not be possible to resolve all the possible
arrangements of founders on the pedigree.  I describe three ways of
classifying these different arrangements (and hence the individuals at
the bottom of the pedigree) the last two of which correspond to the level
of resolution that is available in unlinked versus linked marker data.   

At the simplest level, we can classify
individuals into different {\em gene frequency classes} which correspond
to the expected proportion of gene copies 
originating from one species or another within that individual.  This
is determined by the number of founders from  each species on the
pedigree.  Since there are $2^n$ founders for a pedigree with $n$
generations, there are $2^n + 1$ different gene frequency classes.  Each
class is determined by the number,  $a$,  of founders originating from
the species~$A$ population ($a = 0,1,\ldots,2^n$).  In
Figure~\ref{fig:n=2ped}, both (c) and (f) belong to the gene frequency
class with $a=2$. The individuals at the bottoms of the other pedigrees
belong to the remaining four distinct gene frequency classes.   

Thinking in terms of gene frequency classes provides some perspective on
the  latent variable $Q_i$ described in the previous chapter.  
\citeN{Pritchardetal2000} introduced $Q_i$ as the proportion of genome of
the
$i\thh$ individual originating from Population~0.   $Q_i$
reflects membership of the $i\thh$ individual in a gene frequency class as
$n\rightarrow\infty$.  As $n\rightarrow\infty$ the number of possible
gene frequency classes becomes infinite---hence the continuous nature of
$Q_i$ in the previous chapter.  In the present chapter, we will
use $Q$ to denote the proportion of an individual's genome derived from
the Species~A population; however, here, $Q$ will be discrete in nature,
rather than continuous.  We refer to $Q$ as the ``genetic heritage
proportion."

In recently hybridized populations, there is more
information available than that used by considering only gene frequency
classes.  We may instead
consider the different {\em genotype frequency classes} into which an
individual may fall.
The members of a genotype frequency class all share the same expected
proportion of the three possible single-locus genotypes with respect to
origin of the gene copies.  These three genotypes are: 1) both gene copies
originating from the Species~\A population; 2) one gene copy from
Species~\A, the other from Species~\B; and 3) both gene copies from
Species~\B\@.  Enumerating these genotype frequency classes, and computing
the expected proportions of the genotypes follows from Mendel's laws. 
Since each individual receives one gene copy randomly selected from the
two in its mother, and another randomly selected from the two in its
father, the expected proportions of the genotypes in an individual are
determined by the gene frequency classes to which its parents belong. 
We let $\btG_g = (G_{g,1},G_{g,2},G_{g,3})$ denote the expected
proportions of the three different genotypes in an individual in genotype
frequency class $g$; these proportions are 
\begin{eqnarray}
G_{g,1} &= &Q_mQ_f  \\
G_{g,2} &=& 2Q_m (1-Q_f) \nonumber \\
G_{g,3} &= & (1-Q_m)(1-Q_f) \nonumber
\end{eqnarray}
where $Q_m$ and $Q_f$ are the genetic heritage proportions of
the individual's mother and father, respectively.  Straightforward
algebra verifies that two individuals $i$ and $j$ will belong to the same
genotype frequency class if and only if the parents of $j$ belong to the
same gene frequency classes as the parents of $i$.  Consequently, the
number of distinct genotype frequency classes after $n$ generations of
possible interbreeding between the species is the number of unordered
pairs of different gene frequency classes after $n-1$ generations: 
$(2^{n-1}+1)(2^{n-1}+2)/2$. For $n\geq 2$ there are always more genotype
frequency classes than there are gene frequency classes. With data on
multiple unlinked loci, it is possible to distinguish between individuals
in different genotype frequency classes.  This is one of our inference
goals, and in the Bayesian context will be pursued by computing the
posterior distribution of a latent variable
$Z_i$ which takes the value $g$ if the $i\thh$ individual belongs to
genotype frequency class $g$. 

Another classification of hybrids may be made into what I refer to as
different {\em hybrid classes}.  Members of the same  hybrid class share
the same number of founders from each species {\em and} the same
arrangement of those founders, up to changes of branching order at any
node on the binary tree of the pedigree.  For $n>2$ there are always more
hybrid classes than there are genotype frequency classes.  As an example,
with $n=3$,
$\F_2$ and $\F_3$ hybrids are in different hybrid classes as shown by the
different arrangements on the three generational pedigrees of
Figure~\ref{fig:F2andF3}.  
%%%%%%%%%% F2 and F3 Peds with n=3
\begin{sidewaysfigure}
\mbox{
\subfigure[$\F_2$ hybrid class]{
\includegraphics*[width=.4\textheight]{F2_3ped.eps}}
\hspace*{.15\textheight}
\subfigure[$\F_3$ hybrid class]{
\includegraphics*[width=.4\textheight]{F3_3ped.eps}}
}
\caption[Pedigrees of the $\F_2$ and $F_3$
hybrid classes with $n=3$]{Pedigrees of the $\F_2$ and $\F_3$
hybrid classes with $n=3$.  (a) is the pedigree of an $\F_2$ individual
following three generations of potential interbreeding, and (b) is the
pedigree of an $\F_3$ hybrid.  Note that the gene frequency classes of
the parents of the individual at the bottom of each pedigree are the
same.  Thus, the bottom individuals on each pedigree belong to the same
genotype frequency class, even though they belong to different
hybrid classes.}
\label{fig:F2andF3}
\end{sidewaysfigure}
%%%%%%%%%%%%%%%%
Nonetheless, they are in the same genotype
frequency class.  With only unlinked markers, it is not possible to
distinguish individuals that are in different hybrid classes, but in the
same genotype frequency class.   For the formulation with unlinked
markers we will deal  exclusively with the genotype
frequency classes.  A probability model with more complex dependence structure 
would be required for dealing with linked markers.



\subsection{Probability of the data}
 Similar to the previous chapter we introduce notation here
to avoid awkward subscripting: let $\thetaAof{\A}{i}{\ell}{1}$ denote
the frequency in the species~\A population of the allele possessed by the
$i\thh$ individual at the first gene copy of its $\ell\thh$ locus.  This
is a shorthand for the doubly-subscripted
$\theta_{\A,\ell,Y_{i,\ell,1}}$. Similarly, for the second gene copy we
will write $\thetaAof{\A}{i}{\ell}{2}$  and for the frequencies in the
population from species~\B we have $\thetaAof{\B}{i}{\ell}{1}$ and
$\thetaAof{\B}{i}{\ell}{2}$. 

Given the population allele frequencies, the gene origin indicators, and 
the genotype frequency class to which an individual belongs, it is straightforward to
compute the probability of that individual's single-locus genotype at the
$\ell\thh$ locus. For our purposes later, it is more useful to have an
expression for the joint probability of the genotype and the gene origin
indicators. In the
$i\thh$ individual in the
$g\thh$ genotype frequency class, this is 
\begin{equation}
P(\btY_{i,\ell},\btW_{i,\ell}|Z_i =
g,\bttheta_{\A,\ell},\bttheta_{\B,\ell}) =
\left\{
\begin{array}{ll}
    \thetaAof{\A}{i}{\ell}{1}  \thetaAof{\A}{i}{\ell}{2}  
    G_{g,1} , &
    ~~\mathrm{if}~W_{i,\ell,1}= W_{i,\ell,2}= 0 
\\ 
     \thetaAof{\A}{i}{\ell}{1}  \thetaAof{\B}{i}{\ell}{2}  
     G_{g,2}/2 , &
     ~~\mathrm{if}~W_{i,\ell,1}=0,~W_{i,\ell,2}= 1   
\\ 
    \thetaAof{\B}{i}{\ell}{1}  \thetaAof{\A}{i}{\ell}{2}  
     G_{g,2}/2 , &
     ~~\mathrm{if}~W_{i,\ell,1}=1,~W_{i,\ell,2}= 0
\\
    \thetaAof{\B}{i}{\ell}{1}  \thetaAof{\B}{i}{\ell}{2}  
     G_{g,3}, &
     ~~\mathrm{if}~W_{i,\ell,1}=W_{i,\ell,2}= 1.
\end{array}
\right.
\label{eq:JointWandY}
\end{equation}
The product of
the two allele frequencies in the above expressions follows from the
assumption that each gene copy in the founders ($n$ generations ago) of
the 
$i\thh$ individual is sampled randomly from the
alleles present in its population of origin.  Then, $G_{g,1}$ is
the probability that an individual in genotype frequency class $g$ has
both gene copies originating from species~\A, $G_{g,2}/2$ is the
probability that the first (second) gene copy originates from species~\A
and the second (first) originates from species~\B, and $G_{g,3}$ is the
probability that both gene copies originated from speces~\B.


For a given genotype frequency class, the marginal probability of the
$i\thh$ individual's genotype at locus
$\ell$ is computed by summing (\ref{eq:JointWandY}) over the latent gene
origin indicators:
\begin{equation}
P(\btY_{i,\ell}|Z_i =
g,\bttheta_{\A,\ell},\bttheta_{\B,\ell})
= 
\sum_{\scriptstyle 0\leq W_{i,\ell,1} \leq 1 \atop \scriptstyle 0\leq
W_{i,\ell,2} 
\leq 1}
  P(\btY_{i,\ell},\btW_{i,\ell}|Z_i =
g,\bttheta_{\A,\ell},\bttheta_{\B,\ell}).  
\label{eq:YprobWsOut}
\end{equation}
And finally, under the assumption of unlinked markers in Hardy-Weinberg
and linkage equilibrium among conspecifics $n$ generations ago, the
probability of the $i\thh$ individual's multilocus genotype is just the
product over the
$L$ single-locus genotype probabilities:
\begin{equation}
P(\bY_i|Z_i = g, \bTheta_\A, \bTheta_\B)
= \prod_{\ell = 1} ^ L
P(\btY_{i,\ell}|Z_i =
g,\bttheta_{\A,\ell},\bttheta_{\B,\ell}).
\end{equation}
This gives us an expression for the probability of the data on a single
individual. We now must derive the probability for the data on all
$M$ individuals in the sample. 

Given $n$ generations of potential interbreeding between the species,
there are
$\Gn = (2^{n-1}+1)(2^{n-1}+2)/2$  genotype frequency classes that members
of our genetic sample may fall into.  We model the individuals
in the sample as being randomly and independently drawn from a mixture of
individuals, each belonging to one of the $\Gn$ genotype frequency
classes with probability
$\pi_g$, $g = 1,\ldots,\Gn$, $\sum_{g=1}^\Gn
\pi_g = 1$.      Using
$\bpi$ to denote the vector of mixing proportions,
$(\pi_1,\ldots,\pi_\Gn)$, we may now write the probability of all the
observed data,
$\bY$, conditional on $n$, $\bTheta_\A$, $\bTheta_\B$ and $\bpi$ as  the
product over the members of the sample of the probability of each of
their multilocus genotypes:
\begin{equation}
P(\bY| \bTheta_\A, \bTheta_\B, \bpi) = \prod_{i=1}^M
\biggl(\,
\sum_{g=1}^\Gn \pi_g P(\bY_i|Z_i = g, \bTheta_\A, \bTheta_\B)
\biggr).
\label{eq:HybLike}
\end{equation}


\section{A Bayesian Specification}
Equation~\ref{eq:HybLike} is the likelihood for $\bTheta_\A$,
$\bTheta_\B$, and $\bpi$.   To pursue Bayesian inference in this problem
requires prior distributions $P(\bTheta_\A)$,
$P(\bTheta_\B)$, and $P(\bpi)$, so that
the posterior distribution may be computed.   We wish to make inferences
not only about $\bTheta_\A$,
$\bTheta_\B$, and $\bpi$, but also the latent variables $\bW$ and
$\btZ=(Z_1,\ldots,Z_M)$, so we are concerned with the posterior
distribution
$P( \bTheta_\A, \bTheta_\B, \bpi,\btZ,\bW|\bY)$ and the marginalizations
thereof.  That posterior distribution is proportional to the joint
probability of all those variables, $P(\bY,\bTheta_\A, \bTheta_\B,
\bpi,\btZ,\bW)$.  With the latent variables present, this joint density
factorizes neatly:
\begin{eqnarray}
P(\bY, \bTheta_\A, \bTheta_\B, \bpi,\btZ,\bW) & =  & 
P(\bTheta_\A)
  P(\bTheta_\B)  P(\bpi)  \\
\mbox{} & \times & \prod_{i=1}^M P(\bY_i|\bW_i,\bTheta_\A,\bTheta_\B)
P(\btW_i|Z_i)P(Z_i|\bpi),
\nonumber
\end{eqnarray}
which, as we will see in Section~\ref{sec:HybMCMC}, allows
straightforward MCMC sampling.


It is computationally convenient and biologically reasonable to take the
specific form of the prior distributions $\bTheta_\A$ and $\bTheta_\B$ to
be Dirichlet distributions, independent over the $L$ unlinked
loci.  That is,
$P(\bttheta_{\A,\ell})$ is 
$\mathrm{Dirichlet}(\lambda_{\A,\ell,1},\ldots,\lambda_{\A,\ell,K_\ell})$.  
Since the Dirichlet distribution is the conjugate prior for the
multinomial distribution, this choice facilitates simulation from the
full conditional distributions for
$\bTheta_\A$ and
$\bTheta_\B$.  The Dirichlet distribution is also 
the multivariate generalization of the beta distribution which arises
theoretically as the equilibrium distribution for gene frequencies in the
presence of genetic drift and linear pressure from migration or mutation
\cite{Wright1938,Wright1952}.  Specification of the parameters
$\blambda_{A,\ell} =
(\lambda_{\A,\ell,1},\ldots,\lambda_{\A,\ell,K_\ell}) $ and
$\blambda_{\B,\ell} =
(\lambda_{\B,\ell,1},\ldots,\lambda_{\B,\ell,K_\ell})$ provides a way to
incorporate prior information about the allele frequencies among the two
species at the $\ell\thh$ locus.  For example, if at locus $\ell$,
previous studies have indicated that species~\B has a very low frequency of
allele $j$ while species~\A has a high frequency, then
$\lambda_{\B,\ell,j}$ should be chosen small, relative to the other
components of $\blambda_{\B,\ell}$, while $\lambda_{\A,\ell,j}$ should be
chosen large.
If one had very strong prior evidence that different species were fixed for
different alleles, then this could also be incorporated in the prior in a similar
manner.
 If, on the other hand, very little prior knowledge is
available about allele frequencies in the two species, then a sensible
choice of prior may be the Jeffreys prior (see \citeN{Gelmanetal1995}). 
For species~\A this would be $\lambda_{A,\ell,j} = 1/K_\ell$ for
$j=1,\ldots,K_\ell$.  Another ``not-so-informative" prior density  is the
uniform Dirichlet distribution with $\lambda_{A,\ell,j}=1$,
$j=1,\ldots,K_\ell$.  This prior de-emphasizes the importance given to those
alleles that only appear in
several copies in the whole sample.  
 
The conjugate  prior for $\bpi$ is also a Dirichlet distribution, so we
shall let  $P(\bpi ) \equiv
\mathrm{Dirichlet}(\zeta_1,\ldots,\zeta_\Gn)$.  Here too, prior knowledge
of the biology of the situation could be incorporated into the prior. For
example, if it was well known that backcrosses between $\F_1$ hybrids and
species~\A had low fitness, then that could be reflected in the prior for
$\bpi$.  Additionally if hybridization and backcrossing was known to be
fairly rare, then this prior knowledge could be reflected by having
smaller $\zeta_g$'s for those genotype frequency classes that
required more episodes of interbreeding
within the last $n$ generations.  In the absence of prior
information on hybridization rates, the Jeffreys prior, $\zeta_{g} =
1/\Gn$, $g=1,\ldots,\Gn$, is once again a suitable choice.

\section{MCMC Simulation from the Posterior
Distribution\label{sec:HybMCMC}}
It is not possible to compute directly the posterior distributions for
the variables that we are interested in. However, simulating from the
joint posterior distribution of all the variables by MCMC can be done via
Gibbs sampling in a manner similar to that for normal finite mixture
models \cite{Die&Rob1994}.
After a sufficient period of burn-in, a sample of variables drawn from this
joint posterior distribution allows Monte Carlo estimation of 
the posterior distribution of any subset of variables of interest, either
marginally, or conditional on the values taken by another subset of
variables.   Given initial starting values for all the variables in the
model, Gibbs sampling proceeds by successively simulating new values for
particular variables in the model from their full conditional
distributions
\cite{Gem&Gem1984}.  I shall denote full conditional distributions by
$P(\cdot|\cdots)$. 

We shall refer to a standard iteration of our MCMC algorithm as a
``sweep."  A sweep consists of a series of steps in which each of the
variables in the probability model (except for the data, $\bY$, which are
fixed) is updated once.   Here, the steps in a single sweep are
\begin{enumerate}
\item For $\ell=1,\ldots,L$, simulate new values for $\btheta_{\A,\ell}$
and
$\bTheta_{\B,\ell}$ from their full conditional distributions,
$P(\bTheta_\A|\cdots)$ and 
$P(\bTheta_\B|\cdots)$, respectively, 
\item Simulate a new value of $\bpi$ from $P(\bpi|\cdots)$,
\item For $i=1,\ldots,M$ and $\ell=1, \ldots,L$, simulate a new value
of
$\btW_{i,\ell}$ from $P(\btW_{i,\ell}|\cdots)$,
\item For $i=1,\ldots,M$, simulate a new value of $Z_i$ from
$P(Z_i|\bY_i,\bTheta_\A,\bTheta_\B)$.
\end{enumerate}
By sampling the current states of all the variables after each sweep, one
aquires a dependent sample suitable for Monte Carlo estimation of most
quantities of interest.  In particular, a Rao-Blackwellized  Monte Carlo
estimate of the posterior probability that individual  $i$ is of the
$g\thh$ genotype frequency class is obtained by averaging the values of 
$P(Z_i=g|\bY_i,\bTheta_\A,\bTheta_\B)$ computed during each sweep.

The full conditional distributions are easily derived.  By conjugacy, 
\begin{equation}
P(\bttheta_{\A,\ell}|\cdots)  \equiv
\mathrm{Dirichlet}(\lambda_{\A,\ell,1} + r_{\A,\ell,1},\ldots
,\lambda_{\A,\ell,K_\ell}+ r_{\A,\ell,K_\ell})
\end{equation} where $r_{\A,\ell,j}$ is the number of gene
copies of allelic type $j$ at the
$\ell\thh$ locus currently allocated to species~\A.  (\ie gene copies
of allelic type $j$ for which the corresponding $W_{\cdot,\ell,\cdot} =
0$).  An analogous expression exists for
$P(\bttheta_{\B,\ell}|\cdots)$.  Again by conjugacy 
\begin{equation}
P(\bpi|\cdots) \equiv \mathrm{Dirichlet}(\zeta_1 + s_1,\ldots,\zeta_\Gn +
s_\Gn)
\end{equation}
where $s_g$ is the number of individuals in the sample currently allocated
to genotype frequency class $g$.  

The full conditional distribution for the pair of gene origin indicators
$\btW_{i,\ell}$ in the $i\thh$ fish currently included in the $g\thh$
genotype frequency class is obtained by Bayes Law:
\begin{equation}
P(\btW_{i,\ell}|\cdots)  = 
\frac{P(\btY_{i,\ell},\btW_{i,\ell}|Z_i =
g,\bttheta_{\A,\ell},\bttheta_{\B,\ell})}
{P(\btY_{i,\ell}|Z_i =
g,\bttheta_{\A,\ell},\bttheta_{\B,\ell})}
\end{equation}
where the numerator and denominator are  given in (\ref{eq:JointWandY})
and (\ref{eq:YprobWsOut}), respectively.   
The full conditional distribution for $Z_i$ would be 
$P(Z_i|\btW_{i,\ell},
\bttheta_{\A,\ell},\bttheta_{\B,\ell})$, however, it is more efficient to
simulate new values of $Z_i$ from $P(Z_i|\bY_i,\bTheta_\A,\bTheta_\B)$,
since this represents a marginalization over the four
possible values of $\btW_{i,\ell}$.   
$P(Z_i|\bY_i,\bTheta_\A,\bTheta_\B)$ is easy to compute, because the
quantities needed to calculate it by Bayes Law  have already been computed
in Step~3 of the sweep.  By Bayes Law
\begin{equation}
P(Z_i = j|\bY_i,\bTheta_\A,\bTheta_\B) = 
\frac{\pi_g P(\btY_{i,\ell}|Z_i =
j,\bttheta_{\A,\ell},\bttheta_{\B,\ell})}{
\sum_{g=1}^\Gn \pi_g P(\btY_{i,\ell}|Z_i =
g,\bttheta_{\A,\ell},\bttheta_{\B,\ell})}.
\end{equation}

It is a good idea to run multiple chains from different starting
values to diagnose mixing problems.  In this case, it is easy to assign
overdispersed starting values by simulating values of $\bTheta_\A$,
$\bTheta_\B$, and $\bpi$ from their prior distributions rather than their
full conditional distributions in Steps~1 and~2 of the first sweep.  I
have used \citeN{Gelman1996}'s estimated scale reduction potential factor
to assess how quickly chains in this problem converge to the target
distribution.  In the case of data from cutthroat trout and steelhead trout
described in the following section this occurs very
rapidly.  Burn in then requires little time.   

\section{Steelhead $\times$ Cutthroat Trout Hybrids in Whiskey Creek}
Hybrids of steelhead trout, {\em O.~mykiss}, and coastal cutthroat trout, {\em O.~clarki}. 
in various streams on the West Coast have been detected using genetic data in several studies
\cite{Cam&Utt1985,Neillands1990}.  In a large genetic survey conducted by the
National Marine Fisheries Service and the Washington Department of Fish and Wildlife (as 
part of the preparation of a status report to determine extinction risk for 
coastal cutthroat trout) widespread evidence for hybridization between the two species was found.
\citeN{Johnsonetal1999} summarize this survey, as well as the available literature from the field and the
laboratory  on hybridization between {\em O.~clarki} and  {\em O.~mykiss}.  They report that 
no severe developmental abnormalities occur in hybrids of the two
species; hybrid offspring are clearly viable.  However, hybrids may possess morphological
and behavioral traits that reduce their fitness in natural environments.  This
accords well with the observation that hybrid individuals are typically detected among
juvenile trout, but adult hybrids are seldom observed, and with the observation that
although hybrization may occur each year (in cases where it has been monitored over time it has been found
to be ongoing)  the two species still remain distinct.  Nonetheless, at some locales,
some fish sampled and analyzed possess genotypes suggesting they belong to
a hybrid class involving more than just one generation of hybridization.  Here, I apply the
methods of this chapter to the survey data from Whiskey Creek in Western Washington to see if 
it is possible to distinguish between $\F_1$ and later hybrids.  

As described in \citeN{Johnsonetal1999}, 74 juvenile trout, believed to be cutthroat, were
sampled from Whiskey Creek.  The sample was not a random sample of the
juvenile fish inhabiting Whiskey Creek because the biologists were expressly trying to sample cutthroat
trout only.  It is quite evident, however, that they were unsuccessful at capturing a sample only of cutthroat
trout.   In the following, the sample is treated as if it were a random sample, and, given the difficulty of
distinguishing the two species as juveniles, this may not be very inaccurate. It should be kept in mind,
however, that careful estimates of the frequency of hybrid individuals would require data from a carefully
planned study conducted for the purpose of identifying hybrids (rather than, as with the Whiskey
Creek data, a study intended to estimate the allele frequencies in cutthroat trout alone).  I briefly discuss
the statistical modeling of different sampling scenarios in the discussion.  

  The 74 fish were genotyped at 50
enzyme loci.  The occurrence at some loci of several allelic types typically found only at low frequency in
cutthroat populations but at high frequency in steelhead populations led
\citeN{Johnsonetal1999} to separate the 74 fish into a group of 48 putative cutthroat, 21 putative steelhead,
and 5 putative hybrids.  These classifications were based on two rules:
1) fish homozygous for the steelhead common allele at \locus{ADA-2}, \locus{mAH-2}, and \locus{CKA-2} were
classified as steelhead; and
2) fish possessing a steelhead
common allele in four or more of the eight loci, \locus{sAAT-4}, \locus{ADA-2}, \locus{mAH-2}, 
\locus{mAH-3}, \locus{CKA-2},
\locus{IDDH-1}, \locus{sIDHP-2} and \locus{PEPA}, were classified as hybrids.
David Teel of the National Marine Fisheries Service kindly provided me with the data
on those 74 fish, typed at 50 loci from Whiskey Creek. Somehow, this dataset did not include
the four supposedly informative loci  \locus{mAH-2}, \locus{CKA-2}, \locus{sAAT-4}, or \locus{PEPA}.  Even
without those loci, however, the dataset is still suitable for demonstrating some important points 
about the methods developed within this chapter.  In the remainder of this section I describe
four different analyses, two with real data and two with simulated data, that I carried out.  The results of
these analyses appear in Section~\ref{sec:CuttResults}.

\subsection{Four analyses for demonstration}
I performed four analyses:
\begin{enumerate}
\item I first analyzed the data from Whiskey Creek using the same method as that used for the
Scottish cat data
in Chapter~\ref{chap:MixAdmix}.  The probability model underlying this analysis allows
fish to belong either to a pure cutthroat or pure steelhead category, or to a generic
``admixed" category. The analysis also provides
a posterior distribution for  the putatively admixed individuals' genetic heritage
proportions.  Of the 50 loci available, 30 were polymorphic in the sample from Whiskey Creek
and were used for the analysis.  For the MCMC I used 5 different chains, each with 2000 sweeps for 
burn in and then 40,000 sweeps during which samples were collected every sweep from each chain.  This
required 12 hours on a laptop computer with a 266 MHz G3 (Macintosh) processor.  
\item I then re-analyzed the data from Whiskey Creek using the methods developed in the present
chapter.  For this analysis, I assumed there were 8 genotype frequency classes to which individuals might
belong---the six classes arising from $n=2$ generations of potential interbreeding, as well as two more
corresponding to offspring generated by breeding between two backcrossed individuals. 
Table~\ref{tab:GtypFreqClasses} lists the expected proportions of the different single locus genotypes in
these eight classes, and also gives the names that I use to refer to them.  I use these names only as
representative ``type"  names---for example, it should be kept in mind that the ``$\F_2$" genotype frequency
class contains other indistinguishable hybrid classes as well, such as $\F_3$.   
\begin{table}
\caption[Genotype frequency classes assumed for the analyses]{Genotype frequency classes assumed for the
analyses.  1--5 and 7 are the six  genotype frequency classes that arise given $n=2$ generations of potential
interbreeding.
$G_{g,1}$, $G_{g,2}$, and $G_{g,3}$ are as described in the text.  The final column gives names that I use to
refer to these genotype frequency classes.  Classes 6 and 8 require $n=3$ generations of potentical
interbreeding, as they would be formed by Cutt Bx $\times$ Cutt Bx and St Bx $\times$ St Bx
pairings, respectively.}
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}llcccc}
\hline\hline
$g$ & $Q$  &  \underline{$G_{g,1}$ (\A,\A)}  & \underline {$G_{g,2}$ (\A,\B) or (\B,\A)} & 
\underline{$G_{g,2}$ (\B,\B)}  & 
Name 
\\  \hline
1 & 1.00 &  1.0000  &  0.0000 & 0.0000  & Pure Cutt   \\
2 & 0.00 &  0.0000  &  0.0000 & 1.0000  & Pure St   \\
3 & 0.50 &  0.0000  &  1.0000 & 0.0000  & $\F_1$   \\
4 & 0.50 &  0.2500  &  0.5000 & 0.2500  & $\F_2$   \\
5 & 0.75 &  0.5000  &  0.5000 & 0.0000  & Cutt Bx   \\
6 & 0.75 &  0.5625  &  0.3750 & 0.0625  & Cutt $(\mathrm{Bx})^2$   \\
7 & 0.25 &  0.0000  &  0.5000 & 0.5000  & St Bx   \\
8 & 0.25 &  0.0625  &  0.3750 & 0.5625  & St $(\mathrm{Bx})^2$   \\
\hline \hline
\end{tabular*}
\label{tab:GtypFreqClasses}
\end{table}
Though prior information is available on allele frequencies in other steelhead and cutthroat populations, I
chose to do this analysis using the Jeffreys prior ($\lambda_{\ell,j} = 1/K_\ell$, $j=1,\ldots,K_\ell$) for
allele frequencies.  I also used the Jeffreys prior ($\zeta_{g} = 1/8$, $g = 1,\ldots,8$) for the mixing
proportions ($\bpi$) of the different genotype frequency classes.  100,000 sweeps of five chains started from
overdispersed starting values were run.  This required 5 hours on the same laptop computer with the 266 Mhz G3
(Macintosh) processor.
\item From the posterior mean estimates in Analysis~2 of the allele frequencies in the cutthroat and the
steelhead populations of Whiskey Creek, I simulated a new sample of 391 fish.  The purpose of this was to see
whether inferences about genotype frequency classes could be made any more sharply with many more individuals
in the sample, but still with loci that were roughly as informative as those in the Whiskey Creek dataset. It
is also, of course, instructive to analyze simulated data in which the ``truth" is known.  In this case I
included in the sample 200 Pure Cutt, 150 Pure St, 20 $\F_1$, 10 $\F_2$, 5 St Bx, and 2 each of Cutt Bx, Cutt
$(\mathrm{Bx})^2$, and St $(\mathrm{Bx})^2$ individuals.  I once again used a Jeffreys prior for allele
frequency, but used the uniform Dirichlet prior ($\zeta_g = 1$, $g=1,\ldots,8$) for the mixing proportions
$\bpi$.  35,000 sweeps of five chains with different starting values required 8 hours.
\item The last dataset analyzed is another simulated dataset with 391 individuals and the same numbers of
individuals from each of the 8 genotype frequency classes as in Analysis~3.  However, genetic data were
simulated for 12 diallelic loci with very large frequency differences between the two species---at each locus,
species~\A has an allele at frequency .995, which is at frequency .005 in species~\B.  This was designed to
mimic the situation in which a researcher has a moderate number of nearly diagnostic loci.  The
analysis was once again carried out using the Jeffreys prior for allele frequencies and a uniform Dirichlet
prior on the mixing proportions. 45,000 sweeps of five chains, each with overdispersed starting
values, required 5.1 hours.    
\end{enumerate}



\section{Results\label{sec:CuttResults}}
\subsection*{Analysis 1} 
The results of the simultaneous mixture/admixture analysis of Chapter~\ref{chap:MixAdmix} on the Whiskey
Creek dataset motivated, to some extent, the development of the techniques in this  chapter. 
Figure~\ref{fig:Anal1}(a)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\begin{center}
\subfigure[$P(\mathrm{PURE})$~by~Chapter~\ref{chap:MixAdmix}~Methods]{\includegraphics*[width=.9\textwidth]
{CuttsProbPure.eps} } \\
\subfigure[Histograms showing $P(Q_i|\bY)$ for 7 fish]{\includegraphics*[width=.97\textwidth]
{CuttQHists.eps}}
\end{center}
\caption[Results from analysis 1]{Results from analysis 1.  (a) shows the posterior probability that an
individual is purebred (Pure Cutt or Pure St) for all 74 fish ordered by that posterior probability.  The
labels of the fish (as given in the dataset) are on the horizontal axis.  Four fish have high posterior
probability of being admixed---19, 48, 16, and 72, while three more have probability of being admixed
greater than .45---52, 21, and 49.  (b)
$Q_i$ for the  seven fish that appear most admixed: 48,19, and  16 have genetic heritage proportions around
that expected for $\F_1$ or $\F_2$ hybrids while 49, 21, and 52 have $Q_i$'s closer to those expected for
backcrossed hybrids.  }
\label{fig:Anal1}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 shows the posterior probability that each of the 74 fish in the sample is
Pure Cutt or Pure St.  Seven of those fish have posterior probabilities of being purebred less than .55. 
In Figure~\ref{fig:Anal1} are posterior distributions of the genetic heritage proportions $Q_i$ of these
seven fish, suggesting that perhaps they belong to different hybrid classes.  Fish 48, 19, and 16, have
posterior means for
$Q_i$ that are near .5.  This would be the genetic heritage proportion of either an $\F_1$ or $\F_2$ hybrid;
the analysis of Chapter~\ref{chap:MixAdmix} is not able to distinguish those two different hybrid classes. 
Fish 49, has a $Q_i$ near $.25$, which would correspond to the St Bx or St $(\mathrm{Bx})^2$ categories.  Fish
21 and 52, on the other hand, appear to have roughly 3/4 of their ancestry from the cutthroat species
suggesting that they may belong to the Cutt Bx or Cutt $(\mathrm{Bx})^2$ categories.  Once again, it is not
possible to distinguish between those possibilities without explicitly modeling the hybridization process as
done in the next analysis.  




\subsection*{Analysis 2}
In this analysis, each fish is assigned a posterior probability of
belonging to each one of the 8 different genotype frequency classes.  Two of those classes are purebred
categories (Pure Cutt and Pure St).  All the individuals with a posterior probability of being Pure Cutt
greater than .5 have a negligible posterior probability of being Pure St, and vice versa.  For these
individuals, Figure~\ref{fig:Anal2}(a)
%%%%%%%%%%%%%%%%%%%
\begin{figure}
\begin{center}
\subfigure[$P(\mathrm{PURE})$~with 8 genotype frequency classes]{\includegraphics*[width=.9\textwidth]
{HYC_PureOrdered.eps}}  \\
\subfigure[Posterior probabilities of genotype frequency class for 3
fish]{\includegraphics*[width=.9\textwidth] {3RealFishHYCProbs.eps} } 
\end{center}
\caption[Results from analysis 2]{Results from Analysis~2.  Description in text.}
\label{fig:Anal2}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
shows their posterior probabilities of being Pure Cutts (open circles) or Pure St (closed cirles),
respectively.  Only three fish have posterior probability of being purebred less than .55.  These three are
48, 19, and 16.  Fish 72, 49, 21 , and 52 are no  longer among them---they all have posterior probabilities
greater than .7 of being purebred.  

Figure~\ref{fig:Anal2} shows the posterior probabilities that fish 48, 19, and 16 belong to each of the
eight different genotype frequency classes.  Fish~19 has posterior probability of being an $\F_1$ hybrid that
is more than three times that of being in the $\F_2$ category.  Nonetheless, the posterior probability that
it belongs to a genotype frequency class other than $\F_1$ is .36; there is not sufficient evidence to assign
it to a single genotype frequency class with any sort of confidence.  This is even more true for fish~48
and~16.  



\subsection*{Analysis 3}
The posterior probability that each of the 200 simulated cutthroat trout belong to the Pure Cutt class is
shown in Figure~\ref{fig:Anal3}(a).  
%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\begin{center}
\subfigure[$P(\mathrm{Pure~Cutt})$ for 200 simulated cutthroat trout]{\includegraphics*[width=.7\textwidth]
{200PureCutt_Post.eps} } \\
\subfigure[$P(\mathrm{Pure~St})$ for 150 simulated steelhead trout]{\includegraphics*[width=.7\textwidth]
{150PureSt_Post.eps}} \\
\subfigure[Probability of pure descent of 41 simulated, hybrid juvenile
trout]{\includegraphics*[width=.7\textwidth] {41Hybrids.eps}}
\end{center}
\caption[Probability of pure descent of fish in analysis 3]{Results from Analysis 3.  Description in
text.}
\label{fig:Anal3}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%
Most of these fish have high posterior probability of being Pure Cutts.
The same is true for the 150 simulated steelhead---all but a few of them have posterior probability greater
than .9 of being in the Pure St category.  In other words, with data similar to those from Whiskey Creek, and
with many purebred individuals of each species sampled, it is unlikely that any purebred individual will
receive high posterior probability of being in  a non-purebred genotype frequency category. 
Figure~\ref{fig:Anal3}(c) shows the probability of being in either of the Pure Cutt or Pure St categories
for the 41 simulated hybrid fish.  The different symbols denote the true genotype frequency class of each
fish.  The $\F_1$ and $\F_2$ individuals are most readily indentified as being hybrids of some sort.  The 
various backcrossed categories, however, contain some members with high posterior probabilities of being Pure
Cutt or Pure St---they are not easily detected as hybrids.    

It appears difficult to make clear distinctions between the different hybrid genotype frequency classes with
the genetic  data used.  Figure~\ref{fig:MoreAnal3}(a)
%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\begin{center}
\subfigure[20 simulated $\F_1$'s]{\includegraphics*[width=.9\textwidth]
{20SimmedF1_stack.eps} } \\
\subfigure[21 simulated hybrids as indicated]{\includegraphics*[width=.9\textwidth]
{21SimmedHYC_Probs.eps}}
\end{center}
\caption[Posterior probabilities of genotype frequency class for simulated hybrids]{More results from
Analysis~3.  Description in text.}
\label{fig:MoreAnal3}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
shows the posterior probabilities of inclusion in the 8 genotype frequency classes for the 20 $\F_1$
hybrids in descending order of the posterior probability that they are $\F_1$'s.  While many of them have
high posterior probability of being $\F_1$, there is also one with posterior probability greater than .5 of
being in the Pure St category.  
For the non-$\F_1$ hybrids, the situation is even
less promising.  Of the 10
$\F_2$ individuals (the first 10 individuals in Figure~\ref{fig:MoreAnal3}(b)), not one of them has posterior
probability greater than .5 of being in the $\F_2$ category.  Finally, for the backcrossed genotype frequency
classes, as the final 11 columns in Figure~\ref{fig:MoreAnal3}(b) reveal, there is little relationship
between a fish's true genotype frequency class and the estimated posterior probability of being in that
class.  


\subsection*{Analysis 4}
This analysis reveals that, even if all the loci available have great allele frequency differences between
the species, many loci are still required to distinguish individuals from genotype frequency classes beyond 
the two pure categories and $\F_1$'s.  All the simulated Pure Cutt and Pure St fish had very high posterior
probabilities of being from their respectively correct categories.  The same is true of the $\F_1$'s as
shown in Figure~\ref{fig:Anal4}(a).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\begin{center}
\subfigure[20 simulated $\F_1$'s]{\includegraphics*[width=.9\textwidth]
{20_GoodData_F1s.eps} } \\
\subfigure[21 simulated hybrids as indicated]{\includegraphics*[width=.9\textwidth]
{21_GoodData_OtherHybs.eps}}
\end{center}
\caption[Posterior probabilities for hybrids with 12 very informative markers]{Analysis 4 results.
Description in text.}
\label{fig:Anal4}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%
Each of them has posterior probability near 1 of being from the $\F_1$ category.  However, for the $F_2$'s
and beyond, the story is mixed.  Four of the $\F_2$'s have posterior probability greater than .9 of  being in
the $\F_2$ category.  However, four of them also have posterior probabilities lower than .5.  Within the backcrossed
categories, distinctions are even less reliable.  It is not difficult to see why this might be the case with
as few as 12 loci, even if all the loci were perfectly diagnostic (\ie there were only private ``cutthroat"
and ``steelhead" alleles at each).  For example, an individual in the  St
$(\mathrm{Bx})^2$ category, with probability $(1-.0625)^{12} = .46$, will have no single locus genotypes
homozygous for cutthroat origin, and will thus look very much like a St Bx individual.  
For all of these 21 second-generation (or beyond) hybrids the posterior probability of being in a genotype
frequency class ``beyond"
$\F_1$ is quite high ($>.99$ for all but one).  In that sense, they may be reliably categorized as ``$\F_2$ or
backcrossed" by the method.  

\section{Discussion \label{sec:CuttDiscussion}}
I have described a method for Bayesian inference in populations of recently hybridized species.  The four
short exercises in data analysis using the method bring three main issues to attention.  The first is that
it requires many markers, and substantial genetic differentiation between the species at those markers to be
able to reliably distinguish hybrids in classes beyond the $F_1$ class.  For many problems, it may be
preferable to include only four categories analogous to Pure Cutt, Pure St, $\F_1$, and ``$F_2$ and beyond." 
Alternatively, one could plan to collect data on many loci.  Doing so, however, increases the chance that
the assumption of unlinked loci will be violated.  It would be possible to model the case of linked loci if
one knew the recombination fractions between loci.  The genotype frequency class at one marker would then
depend on the genotype frequency classes at the other markers on the chromosome.  Just like the
identity-by-descent process along homologous chromosomes, this ``genotype frequency class" process  would be
non-Markovian.  It is from this non-Markovian process that one would gain the information
(though probably very little) to distinguish different hybrid classes within the same genotype frequency
class using linked markers.  In a strict analysis, linked markers and the possibility of linkage
disequilibrium among conspecifics at those markers would also necessitate the use of haplotype frequency
parameters   rather than the allele frequency parameters $\bTheta_\A$ and $\bTheta_\B$.  While the
necessary computations for this sort of model could all be done in terms of the underying vector of
segregation patterns at each locus,  the use of some sort of approximation, for example a Markov
approximation to the ``genotype frequency class" process along the chromosome as taken by
\citeN{McKeigue1998} in a related problem, is also a possibility. 

Another issue is the importance of the sampling model used.  In the case of the model as described in this
chapter, individuals are assumed drawn at random from a population which is a mixture in the proportions
$\bpi$ of individuals from the different genotype frequency classes.  This is violated in
the case of the cutthroat trout data, because when the biologists collected the specimens, they were
explicitly trying to obtain pure cutthroat, and hence throwing back those individuals that looked like
steelhead or hybrids. Clearly, it must not be easy to distinguish cutthroat juveniles from steelhead
or hybrid juveniles on the basis of morphological characters.  In order to estimate accurately the proportion
of hybrids in a locale, or even to estimate accurately the posterior probability that an individual is a
hybrid, it is imperative to design the study with those goals in mind.  Having an explicit model, like the one
described in this chapter, that includes the sampling of the organisms is an asset, since the model may be
tailored to particular sampling schemes.  For example, it would be possible to model stratified sampling in
which sampled organisms were first put into ``possibly hybrid" and ``probably purebred"  categories on the
basis of their morphological traits, and  then a random subset of individuals from each of those
categories was genetically typed.  Extending the present model to such a scenario  requires only that
two different mixing proportions be used for the respective samples, \ie a vector $\bpi_h$ for the ``possibly
hybrid" sample and a $\bpi_p$ for the ``probably pure" sample.  The specification of the latent data and
allele frequency parameters would remain the same, and the MCMC would proceed with very little modification. 
Such a sampling scheme might be very useful for species in which hybridization is rare, and leads to some
morphological distinction of the hybrid individuals that could be used to concentrate sampling upon them.  

Finally, the analysis of the real data exposes a difficulty in the Bayesian analysis of finite mixtures. In
this case, there are only three individuals with high posterior probability of being hybrids, yet there are 6
different hybrid genotype frequency classes.  Hence, much of the time during the running of the
Markov chain for MCMC, some of the components corresponding to those genotype frequency classes will be
empty, and they will never have many individuals allocated to them.  This makes the results more
sensitive to the prior chosen for $\bpi$ than it would be if all the genotype frequency classes were well
represented in a large sample, and provides another argument for reducing the number of components
used in the model to two pure categories, an $\F_1$ category, and then the ``$\F_2$ and beyond" category.  
Fortunately, however, in the Gibbs sampling, individuals still mix easily between the empty or nearly-empty
genotype frequency classes.  This contrasts with the Bayesian analysis of normal mixtures in which trapping
states are often encountered with nearly-empty components \cite{Die&Rob1994}.  The distinction here occurs
because the component-specific parameters are the fixed quantities $\btG_g$, which are not affected by the
number of fish allocated to each component.  This feature would also make it straightforward to implement a
reversible jump MCMC \cite{Green1995} sampler to allow the number of genotype frequency classes to be modeled
as an unknown random variable whose posterior distribution was to be determined.  
