\chapter{$\lambda$ and Overlapping Generations}
\label{chap:overlap}

\section{Introduction}
The reproduction of natural populations is not always well-characterized
by a model with discrete generations.  In particular, of the species of
Pacific salmon, only the pink salmon, {\em Oncorhynchus gorbuscha}, has a
discrete-generation life history; all pink salmon, within their native range, mature
at two years of age.  The other species of Pacific salmon mature,
reproduce, and senesce at a variety of ages.  For
example, a spawning collection of chinook salmon might consist of three-,
four-, five-, and six-year old fish.  Each different year class has
descended from a different collection of reproducing parents.  

These factors complicate the estimation of effective size in salmon
populations---the Wright-Fisher model simply does not describe their life history
very well. In simulation studies, however, 
\citeN{Wap&Tee1990} and  
\citeN{Waples1990(2)} show that many quantities of interest,
such as allele frequency variance, rate of loss of heterozygosity, and the rate of
loss of rare alleles, in a salmon population all depend on the average
generation length and the effective number of breeders per year, $N_b$.   
\citeN{Waples1990} demonstrates that there is an approximately linear relationship
between Wright's $F$-statistic and $1/N_b$ in salmon populations.  He then shows
how that relationship may be used to estimate the harmonic mean $N_b$ from
genetic samples of juveniles descended from temporally-spaced brood years.   

The goals of this chapter are different.  Rather than estimating
an overall effective number of breeders for the population, the interest here is in
estimating a $\lambda$-like quantity---a ratio of
effective spawners to  the census number of spawners---given data on the census sizes
of fish of different age-groups and genetic data either from adults or juveniles or
both.  This goal is pursued within the context of a long time series of demographic and
genetic data of the sort that should become increasingly available due to the falling
costs of genotyping and the ability to amplify DNA from archived fish materials \cite{Nielsenetal1999}. 
For example,
\citeN{Ardren1999} describes extensive fish scale collections from two
intensively-studied steelhead ({\em Oncorhynchus mykiss}) populations on the West
Coast.  These fish scales, taken from both spawning adults and outmigrating juveniles
allow the age of each fish to be determined.   Also, as
\citeN{Mil&Kap1997} and 
\citeN{Ardren1999} have shown, microsatellite loci may be reliably amplified
from these fish scales. Furthermore, Canadian fisheries agencies together with other
scientists have proposed launching a program of close genetic monitoring of a
``reference" stream on the coast of Vancouver Island, in which spawners are carefully
counted and samples from the population are genotyped on a regular basis (William
Ardren, pers.\ comm.).  The data-analysis framework described in this chapter would be
very appropriate for such monitoring programs.   

While we will conceptually think in terms of a $\lambda$ for each age group of
adults, we will rely heavily on the urn model for
genetic inheritance, described in the last chapter, in order to derive a probability model and develop
Markov chain Monte Carlo (MCMC) methods for computing the  posterior probabilities of the parameters. 
Having such  a model in which the census number of breeders is considered known, and is used in the
probabilistic model for the population,  but in which  the corresponding effective size may be altered
by changing a simple parameter  which does not alter the census sizes, is crucial to formulating a
reasonable probability model.  In the following section I develop the probability model and several
extensions to accommodate different sampling strategies and the occurrence of null alleles.  In 
Section~\ref{sec:MCMCoverlap}, I exploit the simple neighborhood structure in the model to develop
single-site Metropolis-Hastings updates for the latent variables in the model.  These updates form the
basis of a Markov chain from which we may sample from the posterior distribution of the parameters of
interest.  I represent the dependence structures using the intuitively appealing ``language" of
graphical models.  Since I use only the simplest results from the theory of graphical models, it should
be self-explanatory to most.  However, the reader interested in learning more about graphical models in
statistics is referred to the comprehensive text by
\citeN{Lauritzen1996}.
Finally, in Section~\ref{sec:OverlapSim}, I demonstrate the potential of the method in
several small trials on genetic data simulated using census size estimates of chinook
salmon from a Snake River tributary.  The results suggest that the method 
works under such conditions.  However, future work assessing  the robustness of the
method to departures from the assumed model and  characterizing the mixing properties of the
sampler under different data scenarios is warranted.

Earlier work that I did on this topic involved an extension of the methods of
Chapter~\ref{chap:ForBack} to a case with a Pacific-salmon-like life history.  I did
not pursue that approach any further, but  I include a brief description of it in
Appendix~\ref{app:MVN}.  The urn model provides a superior approach.


\section{Overlapping Generations via an Urn Model}
The urn model for genetic inheritance described in the previous chapter provides
a good  mechanism for modeling genetic drift in populations with
complex life histories, like those of Pacific salmon.  This section describes how
it may be applied in such a context.  First we shall examine a model for the
conditional dependence structure of the variables in such a population, without
reference to specific probability distributions.  We then ``clothe that
backbone" with the specific probability distributions
chosen to represent the  population-genetic sampling, as well as the taking of
genetic samples from juveniles and adults.

\subsection{Dependence structure with the Pacific salmon life history}
  We consider a
population of dioecious, diploid, semelparous organisms, in which adults may
mature and mate between the ages of
$a^-$ and
$a^+$, inclusive, and from which it is straightforward to sample and count the
reproductive adults separately from the rest of the population (as is the case
with Pacific salmon).    For example, a pink salmon population would have
$a^-=a^+=2$, while for a species like chinook salmon in some rivers $a^-$ might
be 3 and
$a^+$ might be 5 or 6.  Assume that accurate estimates of the census sizes
of adults of different age classes are available over a specific time period
beginning at
$t=0$ and ending at $t=T$.  The census of $a$-year-old adults breeding at time $t$ is
denoted $\CA_{t,a}$.  We shall regard these
estimates as known without error.  Additionally, we shall assume that the number
of juveniles each year has been estimated, or can be specified (to within a rough
approximation, at least) based on the number of adults giving rise to them. We  denote
the estimated juvenile population size at time $t$ by  $\CG_t$.  We will assume that
this population behaves in genetic terms as if it were an ideal population governed
by a parameter $\blambda$ which can be construed as a vector having several
components---one for each age group $ (\lambda_{a^-},\ldots,\lambda_{a^+})$ and one
for the sampling from adult  into juvenile or gamete stages, say  $\lambda^{(w)}$. 
This will become more clear when we actually start assigning probability distributions
in this model. 



Furthermore, assume that genetic samples are
available from the adult and juvenile populations at
$t=0$,
$t=T$, and at least some (and preferably many) time points in between.  It must be 
possible to determine the age of adults, so that the genetic samples can be regarded as
drawn from adults of known ages.  Adult ages can be determined from scales or otoliths taken from
individuals.  Likewise, when sampling juveniles we shall assume that it is possible to sample reliably
from a single age class of juveniles, so that they are known to have descended from a particular brood
year of adults.   This is possible, in practice, because juveniles  of many species of salmon will
migrate to the ocean  at a single, early age; thus, the juveniles in a stream in a given season will
all be of a known age class.  For species, like steelhead, in which the freshwater residence time of
juveniles   may vary widely from individual to individual,
juvenile age, like adult age, can be determined from scales or otoliths as well.  

The
genetic samples involve typing individuals at $L$ loci  assumed to be independently
segregating.  In such a case, it is easy to combine data from the multiple loci, so I
will describe the methodology in detail for a single locus only, and then later
describe how to combine data from multiple loci.  From this single locus, let $K$
alleles be observed in the genetic samples from adults and juveniles.   $\SA_{t,a}$
denotes the sample size of adults of age $a$ taken at time $t$, and $\btY_{t,a} =
(Y_{t,a,1},\ldots,Y_{t,a,K})$ is a vector of allele counts for the $K$ different
alleles observed in the sample of
$a$-year-olds at time $t$.  Likewise, we denote sample sizes from juveniles at time
$t$ by
$\SG_t$, and the observed numbers of alleles from a sample of juveniles at time
$t$ by the
$K$-vector,
$\btZ_t =  (Z_{t,1},\ldots,Z_{t,K})$.

The unobservable, or {\em latent} variables in this model are the allele counts
in the adults of different ages at each of the times $t$, $\btX_{t,a} = 
(X_{t,a,1},\ldots,X_{t,a,K})$, and the allele counts amongst the juveniles at the
different times $t$, $\btW_t =  (W_{t,1},\ldots,W_{t,K})$. Note that the sum of the
$K$ components of $\btX_{t,a}$ is
$2\CA_{t,a}$, and the sum of the components of $\btW_t$ is $2\CG_t$.

Before specifying probability distributions for the observed genetic samples and
for the transitions between the latent variables, it is helpful to simply
consider the conditional dependence structure between the variables, given the
overlapping year-class nature of the population's reproduction. 
We will first investigate this dependence structure under the assumption that
individuals sampled from amongst the adults are not then precluded from reproducing
themselves.  This corresponds to Sampling Scheme~I of \citeN{Nei&Taj1981}  (without the
restriction that the census size is equal to the effective size of the population). 
This sort of sampling would be realized if non-invasive genetic sampling (\eg fin
clips) was used, or if adults were sampled destructively {\em after} spawning. I will
consider Sampling Scheme~II in Section~\ref{sec:SS2}.


Figure~\ref{fig:overlapdigraph} shows an acyclic directed graph for a hypothetical
population in which $T=7$, $a^-=2$, and
$a^+ = 4$.  
%%%%%%%%%% OVERLAP DIGRAPH
\begin{sidewaysfigure}
\includegraphics*[width=\textheight]{OverlapGraphBold.eps}
\caption[Dependence structure with overlapping generations]{Acyclic directed graph
describing the conditional dependence structure in the urn model with overlapping
generations.  Adults reproduce at ages 2, 3, and 4.  The vertical, dotted lines separate 
time into different spawning years.  The dotted circles represent latent variables in the set
$\bX_\mathcal{P}$ (see text).}
\label{fig:overlapdigraph}
\end{sidewaysfigure}
%%%%%%%%%%%%%%%%
In this graph,
the arrows may be taken to represent a temporally-defined dependence. 
That is, $c\longrightarrow d$ may be read to mean ``$c$ is a variable that
`occurs' before $d$ in time, and upon which the distribution of $d$ depends."\footnote{The variable
$c$ is said to be a ``parent" of $d$, and variable $d$ is called a ``child" of variable $c$. This
terminology will be used later in the context of moralizing directed graphs to find neighborhoods of
variables. } The form of the graph thus follows exactly from what we know about reproduction in a
population of Pacific salmon from which we sample both juveniles and adults. The shape of the graph also
admits a simple factorization of the joint probability of all the variables involved.  To express this
succinctly, the following notation will be useful: let the set of relevant times and ages be
denoted
$\mathcal{T} =
\{(t,a):0\le t\le T, a^- \le a \le a^+\}$.  The set of times and ages which are
``initial points" are those for which we must posit a prior distrubution for
adult allele counts over which we will integrate.  This set is $\mathcal{P} =
\{(t,a)\in\mathcal{T}: t-a < 0 \}$, and we will use the shorthand
$\XA_\mathcal{P}$ to refer to the latent allele counts in adults of those ages
and times.  In the graph of Figure~\ref{fig:overlapdigraph}, the elements of
$\bX_\mathcal{P}$ are surrounded by dotted circles.   We will refer to the set
of pairs,
$(t,a)$ which are not in
$\mathcal{P}$ as being in the set $\mathcal{P}^c =
\{(t,a)\in\mathcal{T}:t-a\ge 0\}$.  We shall denote by
$\syset =
\{(t,a)\in\mathcal{T}:\SA_{t,a} > 0\}$ the set of times and ages
for which we have drawn genetic samples from the adults.  Similarly, the set of
all times for which a genetic sample from the juveniles has been taken will be
denoted by
$\szset =
\{t:0\leq t\leq T,\SG_t > 0\}$.  And finally let the bold roman versions of each
variable refer to sets of variables as follows:
$\bY=\{\btY_{t,a}:(t,a)\in\syset\}$, $\bZ=\{\btZ_t:t\in\szset\}$,
$\bX=\{\btX_{t,a}:(t,a)\in\mathcal{T}\}$, and $\bW=\{\btW_t:0\le t\le T\}$.

The joint
probability of the observed and latent variables may then be written as
\begin{eqnarray}
P_\blambda(\bY,\bZ,\bX,\bW) &  = & P_\blambda(\bX_\mathcal{P})  
\label{eq:jointoverlap} \\
 \mbox{} &\times & \prod_{(t,a)\in\syset} P(\btY_{t,a}|\btX_{t,a})
\times \prod_{t\in\szset} P(\btZ_t|\btW_t)  \nonumber \\
\mbox{} & \times &
\prod_{(t,a)\in\mathcal{P}^c} P_\blambda(\btX_{t,a}|\btW_{t-a}) \times
\prod_{0\le t\le T} P_\blambda(\btW_t|\btX_{t,a^-},\ldots,\btX_{t,a^+})  \nonumber
\end{eqnarray}
where $P(\cdot|\cdot)$ denotes a conditional probability distribution function
not depending on $\blambda$, $P_\blambda(\cdot|\cdot)$ a conditional distribution
depending on
$\blambda$ and 
$P_\blambda(\bX_\mathcal{P})$ is the prior probability of
$\bX_\mathcal{P}$, which
also depends on $\blambda$.  This prior distribution,
$P_\blambda(\bX_\mathcal{P})$,  must necessarily be a joint distribution on the
components of $\bX_\mathcal{P}$, since we expect that those components will be
dependent.  I will treat this in more detail in Section~\ref{sec:OverPriors}, but for now we take
the joint prior distribution as given.  The two terms on the second line of
(\ref{eq:jointoverlap}) are the probabilities of the observed allele counts in all the
samples of adults and juveniles, respectively.  The two terms on the third line of the
equation are 1) the probabilities due to population-genetic sampling of the latent
allele counts in the adult groups given the juvenile cohorts to which they belonged,
and 2) the probability of the latent allele counts amongst a juvenile cohort given all
the adult age classes contributing to it. 

\subsection{Dependence structure under Sampling Scheme~II and with null alleles
\label{sec:SS2}}
The dependence structure described in the previous section applies to many
situations, but one may encounter other cases which require extensions to that
basic dependence structure.  Here I will deal with two such cases: 1) that when
the genetic sampling is destructive and occurs before reproduction, so that individuals
which are sampled do not have the opportunity to contribute offspring to the following
years, and 2) the case of alleles that are not codominantly expressed. 

\citeN{Nei&Taj1981} used the name ``Sampling scheme~II" for the case when the census
size is larger than the effective size, and the genetic samples are
destructively obtained before the organism is able to  reproduce.  \citeN{Waples1989}
showed that the two sampling schemes could be handled within the same general
$F$-statistic framework, with only a slight difference in the formulae for converting
estimates of $F$ to estimates of
$N_e$.  In our case, using a probability model derived from the urn
model of the previous chapter,  if the census size of the population is known, then the
two different sampling plans can be treated using the different probability
distributions that they give rise to.  

The dependence structure of Figure~\ref{fig:overlapdigraph} applies to Sampling
Scheme~I.  For Scheme~II, the dependence structure is different.  Because
the sampling is destructive, the gene copies sampled are not available
to contribute gametes to the gamete pool.  Hence, $W_t$ will depend upon both
$X_{t,a}$ {\em and} 
$Y_{t,a}$ for $a^- \leq a \leq a^+$.   Figure~\ref{fig:SP2graph} shows the
dependence structure between the variables in a year $t$ under Sampling
Scheme~II.  
%%%%%%%%%% SAMPLING PLAN II DIGRAPH
\begin{figure}
\begin{center}
\includegraphics*[width=.15\textwidth]{SP2digraph.eps}
\end{center}
\caption[Dependence structure under Sampling Plan~II]{Acyclic directed graph describing
the conditional dependence structure in the probability model for overlapping
generations at year $t$, with three age classes of adults (2,3,4), under Sampling
Scheme~II---sampling adults destructively and before reproduction.  The arrows
connecting these variables to other times are omitted in this figure.}
\label{fig:SP2graph}
\end{figure}
%%%%%%%%%%%%%%%%
The arrows between years are not shown, though they occur in the same
places and directions as in Figure~\ref{fig:overlapdigraph}.  Note the inclusion
of the arrows (shown with finely dotted lines) from the samples to the gamete pool. 
This implies a modification of (\ref{eq:jointoverlap}), changing the last factor
to be 
\begin{equation}
\prod_{0\le t\le T}
P_\blambda(\btW_t|\btX_{t,a^-},\ldots,\btX_{t,a^+},\btY_{t,a^-},\ldots,\btY_{t,a^+}).
\end{equation}

It would not be difficult to modify the dependence structure further to account for
the destructive genetic sampling of juveniles. 
However, I do not pursue that here, assuming instead that the gamete/juvenile pool
from which samples are drawn is large enough that the effect
of removing a sample of juveniles has little impact on the allele frequencies which
will occur in the spawning populations of mature organisms. 

Another complication which is frequently encountered is the occurrence of alleles that
are not codominantly expressed.  In this case, it is often possible to detect
homozygotes of a particular allele, but the heterozygotes appear to be homozygotes of
an alternate allele.  This adds another layer
of complexity to the model.  The reason for this is that, when some alleles cannot
be reliably detected in heterozygote form, it is not possible to actually observe
allele counts $\btY_{t,a}$ in samples taken from the adults. Instead one
observes only the counts of phenotypes (heterozygotes and apparent homozygotes) of
different types, which I shall denote by
$\bG^{(Y,o)}_{t,a}$ for $(t,a)\in\syset$.  The superscript $^{(Y,o)}$ refers to the
fact that these are the observed phenotypes in the sample from adults.  Similarly, the
samples from juveniles permit only the observation of phenotypes which will be
denoted by
$\bG^{(Z,o)}_{t}$, $t\in\szset$.  Part of $\bG^{(Y,o)}_{t,a}$ and $\bG^{(Z,o)}_{t}$ should be
thought of as symmetrical matrices with $(i,j)^\thh$ element equal to $(j,i)\thh$
element and giving the number of observed phenotypes with a copy of allele $i$ and
a copy of allele
$j$ ($i$, $j$  codominant).   One additional category of phenotypes must be  included in both
$\bG^{(Y,o)}_{t,a}$ and $\bG^{(Z,o)}_{t}$.  For this we use
$G^{(Y,o)}_{t,a,-}$ and
$G^{(Z,o)}_{t,-}$, to denote the number of individuals in the samples from adults
and juveniles, respectively, in which no bands on a gel were detected.  For example, 
if only allele
$i$ at a locus  was  undetectable, then for the sample from juveniles at time $t$,
$G^{(Z,o)}_{t,i,j}$ would  be zero,
$G^{(Z,o)}_{t,j,j}$ would be the sum of the number of $(j,j)$  genotypes and the
number of heterozygotes of $i$ and $j$, and $G^{(Z,o)}_{t,-}$ would be the number of
$(i,i)$ homozygotes.  

Computing the probability of $\bG^{(Y,o)}_{t,a}$ given
$\btX_{t,a}$ or $\bG^{(Z,o)}_{t}$ given $\btW_t$ would require a sum over all possible
unobserved genotypes consistent with the observed phenotypes.  To
avoid having to do that sum directly, we will introduce more latent variables and
effectively sum over them using MCMC.  This also greatly simplifies
the joint probability function in the case of Sampling Scheme~II in the presence of
null alleles.  The new latent variables are  $\bG^{(Y,\ell)}_{t,a}$ and
$\bG^{(Z,\ell)}_{t}$, which are analogous to the symmetrical matrix portions of 
$\bG^{(Y,o)}_{t,a}$ and
$\bG^{(Z,o)}_{t}$, except that they count the number of different types of genotypes
that would be observed if all the alleles were fully penetrant.  Note that there is a
many-to-one map from the space of $\bG^{(Y,\ell)}_{t,a}$ to that of
$\bG^{(Y,o)}_{t,a}$, and similarly from the space of $\bG^{(Z,\ell)}_{t}$ to 
$\bG^{(Y,o)}_{t}$.


So long as the genetic transmission processes we consider are exchangeable, the
dependence structure between these new variables within a year is given by the graph of
Figure~\ref{fig:GraphsWithGs}(a) for Sampling Scheme~I and
Figure~\ref{fig:GraphsWithGs}(b) for Sampling Scheme~II.  
%%%%%%%%%% SAMPLING PLAN I DIGRAPH
\begin{figure}
\begin{center}
\mbox{\subfigure[]{
\includegraphics*[width=.2\textwidth]{SP1_WithGs_digraph.eps}}
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
\subfigure[]{
\includegraphics*[width=.2\textwidth]{SP2_WithGs_digraph.eps}}
}
\end{center}
\caption[Dependence structure with null alleles]{Acyclic directed graphs describing
the conditional dependence structure in the probability model for overlapping
generations at year $t$ with three age classes (2,3,4) and with some alleles not
codominantly expressed.  The arrows connecting these graphs to other times are omitted
in this figure. (a)  Sampling Scheme~I, sampled adults still contribute offspring to
future generations  (b) Sampling Scheme~II, adults sampled destructively.}
\label{fig:GraphsWithGs}
\end{figure}
%%%%%%%%%%%%%%%%
The joint distribution of all the variables involved
 can then be written similarly to (\ref{eq:jointoverlap}).  Using the notation
$\bG^{(Y,o)} =
\{\bG_{t,a}^{(Y,o)}: (t,a)\in\syset\}$ and $\bG^{(Z,o)} =\{\bG_{t}^{(Z,o)}:
t\in\szset\}$ along with
$\bG^{(Y,\ell)}=\{\bG_{t,a}^{(Y,\ell)}: (t,a)\in\syset\}$ and 
$\bG^{(Z,\ell)}=\{\bG_{t}^{(Z,\ell)}: t\in\szset\}$, we have
\begin{eqnarray}
\lefteqn{P_\blambda(\bG^{(Y,o)},\bG^{(Z,o)},\bG^{(Y,\ell)},\bG^{(Z,\ell)},\bX,\bW)  
=} 
\label{eq:NullOverlap} \\
 P_\blambda(\bX_\mathcal{P})  &\times & \prod_{(t,a)\in\syset}
P(\bG^{(Y,o)}_{t,a}|\bG^{(Y,\ell)}_{t,a})~
\times \prod_{t\in\szset} P(\bG^{(Z,o)}_t|\bG^{(Z,\ell)}_t)  \nonumber \\
 \mbox{} &\times & \prod_{(t,a)\in\syset} P(\bG^{(Y,\ell)}_{t,a}|\btX_{t,a})~
\times \prod_{t\in\szset} P(\bG^{(Z,\ell)}_t|\btW_t)  \nonumber \\
& \times &\prod_{(t,a)\in\mathcal{P}^c} P_\blambda(\btX_{t,a}|\btW_{t-a})~~ \times
\;\prod_{0\le t\le T}
P_\blambda(\btW_t|\btX_{t,a^-},\ldots,\btX_{t,a^+}) 
\nonumber 
\end{eqnarray}
for Sampling Scheme~1.  For Sampling Scheme~II, the final term in the product must be
replaced by 
\[
\prod_{0\le t\le T}
P_\blambda(\btW_t|\btX_{t,a^-},\ldots,\btX_{t,a^+},\bG^{(Y,\ell)}_{t,a^-},\ldots,
\bG^{(Y,\ell)}_{t,a^+}).
\]
Specification of the probability functions specific to sampling with recessive
alleles is deferred until Section~\ref{sec:NullProbs}.


\subsection{Specifying probability distributions\label{sec:SpecProbDist}}
The graph of Figure~\ref{fig:overlapdigraph} and the corresponding
factorization of Equation~\ref{eq:jointoverlap} (as well as their extensions
for the special cases described above) indicate that the probability model here may be
fully defined by assigning distributions to
$P_\blambda(\bX_\mathcal{P})$ and the different $P_\blambda(\cdot|\cdot)$ and $P(\cdot|\cdot)$
distributions.  Specifying these distributions  requires several assumptions
to be made about how reproduction occurs.  In general, I shall model
population-genetic sampling by P\'{o}lya urn models, and the drawing of
genetic samples by sampling without replacement from the populations.  In this context,
sampling ``without replacement" is {\em not} referring to whether or not sampled
individuals are able to reproduce; it is referring to how the genetic samples are
obtained.  Certainly, destructive genetic sampling will occur without replacement,
but even non-invasive sampling could occur without replacement since any previously-sampled fish
will bear marks (for example the loss of a fin clipped for genetic sampling) that should prevent
it from being sampled twice.  For
the samples taken from a large pool of juveniles, there will be little difference
between the multivariate hypergeometric sampling implied by sampling without
replacement and the multinomial sampling implied by sampling with replacement.

The simplest transition to model is that described by
$P_\blambda(\btX_{t,a}|\btW_{t-a})$.  This is the population-genetic sampling
which occurs when juveniles from time $t-a$ are ``selected" or ``sampled" to
survive to be reproducing adults of age $a$ at time $t$.  This depends on
$\lambda_a$, and, following Equation~\ref{eq:sfromlambda} on
Page~\pageref{eq:sfromlambda},  may be parametrized   in terms of a stochastic
replacement quantity $\SRA_{t,a}$ which depends on the juvenile census size,
$\CG_{t-a}$, the adult census size,
$\CA_{t,a}$, and $\lambda_a$:
\begin{equation}
\SRA_{t,a} = \frac{2\CG_{t-a}(1-\lambda_a)}{2\lambda_a\CA_{t,a} - 1}.
\end{equation}
Thus, given $\btW_{t-a}$,
$\btX_{t,a}$ follows the compound multinomial distribution (\ref{eq:CMD}). 
The probability mass function may be
expressed, similarly to (\ref{eq:CMD}), as a normalizing constant times a product of
$K$ terms corresponding to the $K$ different alleles:
\begin{eqnarray}
P_\blambda(\btX_{t,a}|\btW_{t-a}) & = & 
P(
\btX_{t,a}|\btW_{t-a},\lambda_a,\CA_{t,a},\CG_{t-a})  \nonumber \\
& = & P(\btX_{t,a} |\btW_{t-a},\varphi_{t,a},\CA_{t,a}) 
\nonumber \\ & = &
\frac{(2\CA_{t,a})!\;\Gamma(\alpha_\bigdot)}{\Gamma(2C_{t,a} +
\alpha_\bigdot)}
\prod_{i=1}^K \biggl( \frac{\Gamma(X_{t,a,i} + \alpha_i)}{X_{t,a,i}!\;
\Gamma(\alpha_i)}
\biggr)
\end{eqnarray}
where $\alpha_i = W_{t-a,i}/\SRA_{t,a}$ and $\alpha_\bigdot = \sum_{i=1}^K
\alpha_i$.

Modeling the stochastic process and distribution for
$P_\blambda(\btW_t|\btX_{t,a^-},\ldots,\btX_{t,a^+})$ is more difficult,
and requires that more assumptions be made about reproduction and survival in the
population.  The particular problem that arises is that the distribution of
$\btW_{t}$ depends not only on the vagaries of sampling alleles from within each
age class of adults (\ie non-multinomial sampling of gene copies from amongst
the $C_{t,a}$ $a$ year-olds), but also on the fact that adults of different age classes
may produce different mean numbers of juvenile offspring, either by producing more gametes
or by producing individuals with higher survival to the juvenile stage.  This second
effect is akin to that discussed in
\citeN{Rym&Lai1991}, in which the inbreeding effective size of a population is
decreased due to the higher survivorship of a segment of the population included in a
supportive breeding program.  Since it is impossible to determine the
age of the parent of any gene copy sampled amongst juveniles, these two sources of variation in
$\btW_t$ are confounded and may not be separated.  Rather than include these two
confounded processes in a model which is not identifiable, I assume an ideal model for
the production of juveniles from adults of different age classes, and then account for
both of the above-mentioned processes by a single parameter in an urn model scheme.

This ideal model assumes that each adult at time $t$ produces an age-specific number of
gametes, and the survivors to the juvenile stage are sampled from those gametes by an
urn scheme with stochastic  replacement parameter $\SRG_t$.  More specifically, each
diploid adult of age $a$ contributes 
$\fitness_a$ copies of each of its two gene copies to the gamete pool.  Thus, the
counts of the different alleles in the gamete pool at time $t$ are given by the
$K$-vector  $\btB_t = (B_{t,1},\ldots,B_{t,K}) \label{page:Bdef}
=\sum_{a=a^-}^{a^+}\fitness_a\btX_{t,a}$.  Then, the
$2\CG_t$  gene copies in the juveniles are sampled from this gamete pool via a
P\'{o}lya urn scheme in which the stochastic replacement quantity depends on the
parameter
$\lambda^{(w)}$---the conceptual ratio of ``effective juveniles" to the census number
of juveniles.  Letting $B_{t,\bigdot}$ denote the total number of gametes in the gamete
pool at time $t$ ($B_{t,\bigdot} =
\sum_{a=a^-}^{a^+} 2\fitness_a\CA_{t,a}=\sum_{i=1}^K B_{t,i}$), then, once again by
Equation~\ref{eq:sfromlambda} on Page~\pageref{eq:sfromlambda}, we have the stochastic
replacement 
\begin{equation}
\SRG_t = \frac{B_{t,\bigdot}(1-\lambda^{(w)})}{2\lambda^{(w)}\CG_t - 1}.
\end{equation}
And so, the conditional probability
$P_\blambda(\btW_t|\btX_{t,a^-},\ldots,\btX_{t,a^+})$  may now be
expressed as
\begin{eqnarray}
P_\blambda(\btW_t|\btX_{t,a^-},\ldots,\btX_{t,a^+}) & = &
P(\btW_t|\btX_{t,a^-},\ldots,\btX_{t,a^+},\CA_{t,a^-},\ldots,\CA_{t,a^+},
\bfitness,\CG_t)  \nonumber \\
& = &
P(\btW_t|\btB_t,\SRG_t,\CG_t)   \nonumber \\
& = & 
\frac{(2\CG_{t})!\;\Gamma(\alpha_\bigdot)}{\Gamma(2\CG_{t} +
\alpha_\bigdot)}
\prod_{i=1}^K \biggl( \frac{\Gamma(W_{t,i} + \alpha_i)}{W_{t,i}!\;
\Gamma(\alpha_i)}
\biggr)
\end{eqnarray}
where $\alpha_i = B_{t,i}/\SRG_{t}$ and $\alpha_\bigdot = \sum_{i=1}^K
\alpha_i$.

The quantities $\bfitness = (\fitness_{a^-},\ldots,\fitness_{a^+})$ may be interpreted
as fitness measures for different age classes expressing how successful they are at
producing juveniles of sampling age.   In practice,
$\fitness_a$ can be chosen to reflect the biology of the situation.  For
example, a reasonable choice for salmon would be one half the fecundity of age $a$
females.  It should be clear from the above expression, that the absolute magnitudes
of the $\fitness_a$'s are actually irrelevant; the parametrization of $\SRG_t$ in terms
of $\CG_t$ and the relationship between $\alpha_i$, $B_{t,i}$, and $\SRG_t$ ensure
that the relative sizes of the $\fitness_a$'s are all that matter.  Nonetheless, it is
computationally convenient to think of the $\fitness_a$'s in terms of the number of
gametes produced.

Another, and a possibly more elegant, interpretation of this population-genetic sampling
scheme for juveniles is provided by the conditional branching process model of
\citeN{Kar&McG1965} with negative binomial distributions of offspring number (see
Section~\ref{sec:TwoStage}).  In this interpretation, the total number of juvenile gene
copies is fixed to be
$2\CG_t$, however the distribution of the number of
copies of each gene within an age $a$ adult appearing among the juveniles is
exchangeably negative binomial with arbitrary (but equal for all genes) scale
parameter
$\beta$, and shape parameter
$\fitness_a/\SRG_t$.  By such an interpretation it is perhaps even more clear that
$\SRG_t$,  the stochastic replacement quantity for reproduction into juveniles at time
$t$, represents both non-Wright-Fisher sampling within age classes, but also a
departure from our best guess as biologists as to the fitnesses/fecundities of adults
of different age classes.  Since both the
non-Wright-Fisher sampling within age classes, and the unkown differential
survival between age classes reduce effective size of
populations, and will therefore affect $\lambda$, it seems quite reasonable that both
are accounted for in  the parameter $\lambda_t$.  

In the case of Sampling Scheme~II, in which adults are destructively sampled before
reproduction, defining the probability function
$P_\blambda(\btW_t|\btX_{t,a^-},\ldots,\btX_{t,a^+},\btY_{t,a^-},\ldots,\btY_{t,a^+})$
requires only a simple modification to the above mechanism for transmission of genes
to juveniles.  Since sampled adults do not contribute to future generations, we
need merely define $B_{t,i}$ so as to reflect that.  Namely, $B_{t,i} =
\sum_{a=a^-}^{a^+} \fitness_a(X_{t,a,i} - Y_{t,a,i})$, \label{page:BSS2} and
$B_{t,\bigdot}$ must be modified accordingly: ($B_{t,\bigdot} =
\sum_{a=a^-}^{a^+} 2\fitness_a(\CA_{t,a}-\SA_{t,a})=\sum_{i=1}^K B_{t,i}$).  For
Sampling Scheme~II in the presence of null alleles, $Y_{t,a,i}$ in the immediately
preceding sentence may be replaced by the quantity $Y_{t,a,i}^{(\ell)}$
described in the next section.

Finally, we only have to specify probability distributions for the genetic samples
drawn from adults and juveniles,  $P(\btY_{t,a}|\btX_{t,a})$
and $P(\btZ_t|\btW_t)$.  As stated at the beginning of  this section, sampling without
replacement is a good model for the acquisition of genetic samples.  With multiple
alleles, this leads to the multivariate hypergeometric distribution (see
\citeNP{Johnsonetal1997}, Chapter~39).  This distribution may also be written as a
normalizing constant multiplied by a product of $K$ terms, one for each of the $K$
alleles.  So, for the genetic samples taken from adults we have
\begin{equation}
P(\btY_{t,a}|\btX_{t,a}) = \frac{(2\CA_{t,a}-2\SA_{t,a})!
(2\SA_{t,a})!}{(2\CA_{t,a})!}
\prod_{i=1}^K \frac{X_{t,a,i}!}{(X_{t,a,i}-Y_{t,a,i})!Y_{t,a,i}!}.
\end{equation}

For genetic samples taken from the juveniles, we can also use the multivariate
hypergeometric distribution
\begin{equation}
P(\btZ_{t}|\btW_{t}) = \frac{(2\CG_{t}-2\SG_{t})!
(2\SG_{t})!}{(2\CG_{t})!}
\prod_{i=1}^K \frac{W_{t,i}!}{(W_{t,i}-Z_{t,i})!W_{t,i}!},
\label{eq:JuveHyperGeom}
\end{equation}
or, since the number of juveniles is typically large, modeling the process as sampling
with replacement will yield essentially the same result, and so the multinomial
probability distribution is appropriate:
\begin{equation}
P(\btZ_{t}|\btW_{t}) = (2\SG_{t})!
\prod_{i=1}^K\frac{[W_{t,i}/(2\CG_t)]^{Z_{t,i}}} {Z_{t,i}!}.
\label{eq:JuveMultinom}
\end{equation}
Notice that (\ref{eq:JuveMultinom}) also includes a simple product of terms over
alleles.


\subsection{Probabilities with recessive alleles\label{sec:NullProbs}}
For recessive or null alleles at a locus with $K$ alleles, I assume Hardy-Weinberg equilibrium and  a
simple penetrance model which may be summarized by the matrix $\bA$ having elements
$a_{i,j}$, $1\leq i,j \leq K$.  $a_{i,j} = 0$ implies that an allele of type $i$ is
detectable (\ie leaves a band on a gel) when it occurs in the same individual as an
allele of type $j$.  If $i$ subscripts a null allele, than $a_{i,i}=0$ and also
$a_{i,j}=0$ for all other $j$.  This penetrance model can also account for other simple
dominance relationships between alleles (\eg $a_{i,j} = 0$ but $a_{i,i} = 1$).  

Given the latent genotypes of sampled juveniles, $\bG^{(Z,\ell)}_t$, the observed
phenotypes can be found by $G^{(Z,o)}_{t,i,i} = a_{i,i}G^{(Z,\ell)}_{t,i,i} +
\sum_{j\neq i} a_{i,j}|a_{j,i}-1|G^{(Z,\ell)}_{t,i,j} $ and, for $j\neq
i$, by $G^{(Z,o)}_{t,i,j} = a_{i,j}a_{j,i}G^{(Z,\ell)}_{t,i,j} $.  The number of
individuals showing no bands, $G^{(Z,o)}_{t,-}$, is found by subtraction, being
half the number of gene copies not otherwise accounted for.  
Since there is a deterministic map from $\bG_{t,a}^{(Y,\ell)}$ to
$\bG_{t,a}^{(Y,o)}$, 
$P(\bG_{t,a}^{(Y,o)}|\bG_{t,a}^{(Y,\ell)})$ will take the value one whenever
$\bG_{t,a}^{(Y,\ell)}$ is consistent with $\bG_{t,a}^{(Y,o)}$ and zero otherwise.  
The map from
$\bG^{(Y,\ell)}_{t,a}$ to $\bG^{(Y,o)}_{t,a}$ works similarly.  Notice also that the
allele counts in the sample may be easily obtained from
$\bG^{(Y,\ell)}_{t,a}$ or
$\bG^{(Z,\ell)}_{t}$.  We will denote these as $\btY^{(\ell)}_{t,a}$ and
$\btZ^{(\ell)}_t$, respectively.

Deriving the probability distribution $P(\bG^{(Y,\ell)}_{t,a}|\btX_{t,a})$
under the assumption of sampling without replacement requires some combinatoric
calculations. Since the fates of gene copies in the genetic transmission and sampling
models adopted here are exchangeable, the probability of every ordering of gene copies into the adults
in the population, and therefore into the sampled adults from the population, is the
same.  Therefore
$P(\bG^{(Y,\ell)}_{t,a}|\btX_{t,a})$ may be found by counting the ways of drawing
particular combinations of pairs of genes,
$\bG^{(Y,\ell)}_{t,a}$, from the allele counts in the adults $\bX_{t,a}$, and
dividing by the total number of ways of drawing any $\SA_{t,a}$ pairs from the
population.  I show this below, suppressing the $_{t,a}$ subscript and
$^{(Y,\ell)}$ superscript on elements of $\bG^{(Y,\ell)}_{t,a}$, the  $_{t,a}$
subscript on elements of $\btX_{t,a}$ and on the population and sample sizes
$\CA_{t,a}$ and $\SA_{t,a}$, and the
$^{(\ell)}$ superscript and
$_{t,a}$ subscript on elements of $\btY_{t,a}^{(\ell)}$.

First, the denominator of the probability $P(\bG^{(Y,\ell)}_{t,a}|\btX_{t,a})$ is the number of ways of
drawing an unordered collection of
$\SA$ unordered pairs from a population of
$2\CA$ gene copies, which is 
\begin{equation}
\frac{1}{\SA!} \prod_{i=0}^{\SA-1} \smallbico{2\CA-2i}{2} =
\frac{(2\CA)!}{2^{\SA}\SA!(2C-2S)!}.
\label{eq:Pairs}
\end{equation}
The product of binomial coefficients arises from sequentially choosing unordered pairs
without replacement, and the $1/(\SA!)$ accounts for the different orders
in which those pairs may be drawn. 

The numerator of $P(\bG^{(Y,\ell)}_{t,a}|\btX_{t,a})$ may be written as
\begin{equation}
\prod_{i=1}^K \left(\smallbico{X_i}{Y_i} \cdot
\frac{Y_i!}{(2G_{i,i})!\prod_{j\neq i} G_{i,j}!} \cdot
\frac{(2G_{i,i})!}{2^{G_{i,i}}G_{i,i}!} \cdot
\prod_{j<i} G_{i,j}!
\right)
\label{eq:PairsNumerator}
\end{equation}
and explained as follows: we have a product over alleles of four factors; the first
factor is a binomial coefficient that counts the number of ways of choosing $Y_i$ gene
copies of type $i$ from a population having $X_i$ such gene copies.  The second factor
is a multinomial coefficient which counts the ways of partitioning those $Y_i$ gene
copies into the groups of genes participating in
the different categories of genotypes.  The third factor counts the number of ways
$2G_{i,i}$ gene copies of allelic type $i$ can be paired up into $G_{i,i}$
unordered homozygous genotypes (this is a special case of (\ref{eq:Pairs})). And,
finally,
$G_{i,j}!$ is the number of ways of making $G_{i,j}$ heterozygote genotypes from
$G_{i,j}$ copies of alleles of type
$i$ and $G_{i,j}$ copies of alleles of type $j$.  The product of $G_{i,j}!$ is taken
over
$j<i$ since, in combination with the other product from $i=1$ to $K$, this leads to the product
over all heterozygote classes. 

Equation~\ref{eq:PairsNumerator} simplifies modestly so we may write our desired probability
as 
\begin{equation}
P(\bG^{(Y,\ell)}_{t,a}|\btX_{t,a}) = 
\frac{
\prod_{i=1}^K \left(
\frac{X_i!}{(X_i-Y_i)!} \cdot
\frac{1}{2^{G_{i,i}}\prod_{j=1}^K G_{i,j}!} \cdot
\prod_{j<i} G_{i,j}!
\right)
}
{
\frac{(2\CA)!}{2^{\SA}\SA!(2C-2S)!}
}.
\label{eq:DesiredProbability}
\end{equation}
This is 
(\ref{eq:PairsNumerator}) divided by (\ref{eq:Pairs}).  The same is true for
$P(\bG^{(Z,\ell)}_{t}|\btW_{t})$ using $\CG$ and $\SG$, and with $Z$'s replacing $Y$'s
and
$W$'s replacing $X$'s, under the assumption that sampling from juveniles is done
without replacement.  Under the assumption that sampling from juveniles is done with
replacement, $P(\bG^{(Z,\ell)}_{t}|\btW_{t})$ is a simple expression given by a
multinomial distribution with cell probabilities being the genotype frequencies
expected under Hardy-Weinberg equilibrium.

\subsection{The prior distribution for allele counts\label{sec:OverPriors}}
The prior distribution $P(\bX_\mathcal{P})$ presents some interesting difficulties.
Ideally, we would like to use some sort of stationary distribution of allele
counts $\bX_\mathcal{P}$ for the salmon population under study. However, this is
difficult, first, because with fluctuating sizes, the population allele frequencies
won't strictly have a stationary distribution, and second, because even if we knew the
historical sizes of the population, it would not be straightforward to determine the
distribution of $\bX_\mathcal{P}$.  Below, I present, in series, several different ways
of handling the prior,  $P_\blambda(\bX_\mathcal{P})$, starting with the most naive. 
In practice, some combination of the methods described below will probably work
best.   The choice of which to use is a matter of balance between reflecting the
reality of the situation and imposing too much (and possibly incorrect) structure on
the latent variables, which will affect the inferences made. 


The most naive approach would be to use independent priors for the components of
$\bX_\mathcal{P}$.  Independent, uniform priors, for example, would assert very little
{\em a priori} structure on the model.  This is naive because some of the components
of $\bX_\mathcal{P}$ reflect fish that have matured from the same pool of juveniles. 
Clearly, allele frequencies amongst adults matured from the same cohort of juveniles will be
correlated.  Fortunately, if a large sample has been taken from every time and age in $\mathcal{P}$,
then the choice of prior may have little effect since those data (let us call them
$\bY_\mathcal{P}$) will constrain $\bX_\mathcal{P}$ considerably.

An improvement on the above can be made by adding to the model the allele counts in
juvenile pools contributing to the adult populations of $\mathcal{P}$.  For example,
extending the graph in Figure~\ref{fig:overlapdigraph} to include juvenile pools in
``negative time," we could have the variable $\btW_{-1}$, from which $\btX_{1,2}$,
$\btX_{2,3}$, and
$\btX_{3,4}$ are drawn; $\btW_{-2}$ parental to $\btX_{1,3}$ and $\btX_{2,4}$ in the graph; and
$\btW_{-3}$ parental to $\btX_{1,4}$.  Then, even with independent prior distribution on
$\btW_{-3}$, $\btW_{-2}$, and $\btW_{-1}$, the correlation between $\btX_{1,2}$,
$\btX_{2,3}$, and
$\btX_{3,4}$ would be modeled, as well as that between the other elements of
$\btX_\mathcal{P}$.  In general, this approach requires specifying $a^+$ new
variables,  $\btW_{-a^+},\ldots,\btW_{-1}$, and giving them independent priors.  While
this is a great improvement over the first approach, it still does not account for the
correlation that is bound to exist between the allele counts in the juvenile pools in
``negative time."  An {\em ad hoc} approach to doing so is described in the following
method.

An approximate relationship between the variables $(\btW_{-a^+},\ldots,\btW_{-1})$ can
be derived using the work of
\citeN{Waples1990} which explores the expected $F$-statistics
between the allele frequencies corresponding to 
$\btW_{-a^+}$ and the remaining components, as a function of the effective number of
breeders $N_b$ and the proportion of fish maturing at different ages.  Consider
a salmon population progressing through time with effective numbers of spawners $N_b$,
possibly changing each year, and with
$\btf=(f_{a^-},\ldots,f_{a^+})$ being a vector of proportions giving the probability
that a fish matures at a particular age.  Through computer simulation of such a
population,
\citeN{Waples1990} found a linear relationship between the expected value of $F$
computed from allele frequencies in the gamete pools separated by $t$ years and the
inverse of twice the harmonic mean effective number of breeders ($\tilde{N}_b$) in the
$t$ years between the gamete pools considered.  He also showed that  the slope of this
linear relationship depends on $t$ and the proportions $\btf$.  We will denote this
slope by the function
$\Delta_t(\btf)$.  \citeN{Tajima1992} gives a convenient recursive algorithm for
computing $\Delta_t(\btf)$.  In our case,  denoting the allele {\em frequency} in a
juvenile or gamete pool at time $i$ by $p_i$, \citeN{Waples1990} empirically shows
that
\begin{equation}
\frac{\Exp[(p_{-a^+} - p_i)^2]}{p_{-a^+}(1-p_{-a^+})} \approx
\frac{\Delta_{i+a^+}(\btf)}{2\tilde{N}_b}
\label{eq:WaplesEmpirical}
\end{equation}
for $i = -a^+ + 1, \ldots, -2,-1$.  If we assume that the juvenile/gamete pools in
these years are all of the same size, $\CG^{(-)}$ diploids (the superscript $^{(-)}$
refers to these being in ``negative time"), and the expected value of each of
$\btW_{-a^++1},\ldots,\btW_{-1}$ is $\btW_{-a^+}$, then, by (\ref{eq:WaplesEmpirical}), we have for
allele $j$ at time
$i$
\begin{eqnarray}
\Var(W_{i,j}|W_{a^+,j}) &=& \Exp[(2\CG^{(-)})^2(p_{-a^+} - p_i)^2]  \nonumber \\
 &\approx &  (2\CG^{(-)}) p_{-a^+}(1-p_{-a^+})
\left(\frac{(2\CG^{(-)}) \Delta_{i+a^+}(\btf)}{2\tilde{N}_b}\right).
\label{eq:VarFromWaples}
\end{eqnarray}

This is the variance of a binomial random variable with $2\CG^{(-)}$ trials
and success probability $p_{a^+}$, multiplied by the term in the large parentheses. 
That, in turn, is the form of the variance of a beta binomial random variable.
From the discussion in Section~\ref{page:FormOfCMDVar} (Page~\pageref{page:FormOfCMDVar}) of the
relationship between the variance of beta-binomial and binomial random variables, it may be seen that a
distribution satisfying the variance relationship in (\ref{eq:VarFromWaples}) is the
beta-binomial distribution with $2\CG^{(-)}$ trials and  parameters $\alpha_j$ and
$\alpha_\bigdot-\alpha_j$ such that $a_j/\alpha_\bigdot = p_{-a^+}$ and 
\begin{equation}
\frac{2\CG^{(-)} + \alpha_\bigdot}{1+\alpha_\bigdot} = 
\frac{(2\CG^{(-)}) \Delta_{i+a^+}(\btf)}{2\tilde{N}_b}.
\label{eq:WaplesOverdisperse}
\end{equation}
This suggests that the following would be a reasonable way to construct a prior for
the vector
$(\btW_{-a^+},\ldots,\btW_{-1})$:
\begin{enumerate}
\item Assume reasonable values for the proportions of individuals maturing at
different ages, $\btf = (f_{a^-},\ldots,f_{a^+}$).
\item Let $\btW_{-a^+}$ follow a discrete uniform prior (since $\CG^{(-)}$ does not
change in the MCMC  simulations, this term in the distribution will conveniently never
change, either).
\item Given $\btW_{-a^+}$, assume that  $\btW_{i}$ ($i=-a^++1,\ldots,-1$) are drawn
from the gene copies in the juvenile pool at time $-a^+$ via independent  P\'{o}lya urn
schemes with stochastic replacement quantities
$\SRG_i^{(-)} = 2\CG^{(-)} / \alpha_{\bigdot_i}$, where $\alpha_{\bigdot_i}$ is
calculated according to (\ref{eq:WaplesOverdisperse}):
\begin{equation}
\alpha_{\bigdot_i} = \frac{2\tilde{N}_b - \Delta_{i+a^+}(\btf)}
{\Delta_{i+a^+}(\btf) - \tilde{N}_b/\CG^{(-)}}.
\end{equation}
\end{enumerate}

It can be shown that conditional on $\btW_{-a^+}$ such a distribution will have the
variance of (\ref{eq:VarFromWaples}).  Although it will not properly reflect the
covariance between the gamete pools, it should be a very reasonable
approximation.  In practice, of course, the harmonic mean effective number of breeders
will be unknown, but one should be able to make a reasonable estimate at the
harmonic mean census number of spawners of all age groups in the $a^+$ years before
data started being recorded for the population.  Denoting that quantity as
$\tilde{\CA}^{(-)}$, a simple way of estimating $\tilde{N}_b$ given 
$\tilde{\CA}^{(-)}$ and $\blambda$ is $\tilde{N}_b = \lambda^{(-)}\tilde{\CA}^{(-)}$
where $\lambda^{(-)} = \sum_{a=a^-}^{a^+} f_a\lambda_a$.  This is the way in which the
prior distribution depends on $\blambda$.  

Finally, if census sizes (or estimates thereof) of the different aged fish in the
population are known in the years before  the genetic data started being collected,
that information can similarly be used to help define a prior distribution
for
$\bX_\mathcal{P}$.  Doing so is simple---one merely defines time 0 to be the time at
which the first census size data are available.  Then, everything from the previous
several paragraphs still applies for constructing a prior on initial gamete pools, but
one also has several years of census data over which $\btX_{t,a}$'s and $\btW_t$'s may
be sampled in an MCMC sampler, helping to more accurately reflect the joint
distribution of
$\btX_{t,a}$'s when genetic samples are finally taken.

In concluding this section, I point out that while the {\em ad hoc} approach described above
is reasonable and practical, it is not deeply satisfying.  The derivation of an elegant prior
$P_\blambda(\bX_\mathcal{P})$ remains an interesting, open problem.    
 
\section{A Bayesian Formulation and MCMC Simulation from $P(\blambda|\bX,\bW)$\label{sec:MCMCoverlap}}
A Bayesian formulation of this problem is obtained by assigning a prior distribution
$P(\blambda)$ for
$\blambda$. This leads to the posterior distribution
\begin{equation}
P(\blambda|\bY,\bZ) = 
\frac{P(\blambda)\sum_{\bX,\bW}P_\blambda(\bY,\bZ,\bX,\bW)}
{\int_\blambda P(\blambda)\sum_{\bX,\bW}P_\blambda(\bY,\bZ,\bX,\bW)d\blambda}
\end{equation}
where the integral in the denominator is over all values of $\blambda$ and the sum is over
all possible values of $\bX$ and
$\bW$.  This sum and integral are intractable.  However, it is possible to simulate values of
$\blambda$ from this posterior distribution using the Metropolis-Hastings algorithm, and thus, the
posterior distribution may be evaluated by Markov chain Monte Carlo.  This is presented in overview in
the following two paragraphs, and in detail in the remainder of the section.  

Given current values of $\bX$, $\bW$, and $\blambda$, a Metropolis-Hastings update step for $\bX$
involves simulating a new value $\bX'$ from a proposal distribution $q_\bX(\bX'|\bX,\cdots)$ that depends
on
$\bX$, and possibly on the current values of other variables in the model (denoted by ``$\cdots$").  A
uniform random variable on the unit interval $U$ is then drawn. If $U<H_\bX$ then the proposal is
accepted and the value of $\bX$ is changed to $\bX'$.  If $U>H_\bX$, then the value of $\bX$ remains
unchanged.  If
\begin{equation}
H_\bX = \frac{ q_\bX(\bX|\bX',\cdots)P_\blambda(\bY,\bZ,\bX',\bW)} 
 {q_\bX(\bX'|\bX,\cdots)P_\blambda(\bY,\bZ,\bX,\bW)},
\label{eq:XHR}
\end{equation}
then, if $q_\bX$ is such that successively applying the updates using $U$ and $H_\bX$ above leads to an
irreducible Markov chain of $\bX$ values,  that Markov chain will have limit
distribution
$P_\blambda(\bX|\bY,\bZ,\bW)$.  Similarly, updates to $\bW$ can be made by proposing new values $\bW'$
from $q_\bW(\bW'|\bW,\cdots)$, drawing $U$ and accepting the proposal if $U$ is less than 
\begin{equation}
H_\bW = \frac{ q_\bW(\bW|\bW',\cdots)P_\blambda(\bY,\bZ,\bX,\bW')} 
 {q_\bW(\bW'|\bW,\cdots)P_\blambda(\bY,\bZ,\bX,\bW)}.
\label{eq:WHR}
\end{equation}
In the same way, updates to $\blambda$ are made with a proposal distribution
$q_\blambda(\blambda'|\blambda,\cdots)$ and accepted according to the Hastings ratio
\begin{equation}
H_\blambda = \frac{ q_\blambda(\blambda|\blambda',\cdots) P(\blambda')P_{\blambda'}(\bY,\bZ,\bX,\bW)} 
 {q_\blambda(\blambda'|\blambda,\cdots)P(\blambda)P_\blambda(\bY,\bZ,\bX,\bW)}.
\label{eq:LHR}
\end{equation}

Applying these updates in series (update $\bX$, update $\bW$, update $\blambda$, update
$\bX$, update $\bW$, and so on\ldots) leads to a Markov chain with limit distribution
$P(\blambda,\bX,\bW|\bY,\bZ)$.  Sampling $n$ values of
$\blambda$ visited by this chain gives a sequence $\blambda^{(1)},\ldots,\blambda^{(n)}$ which may be
used to estimate $P(\blambda|\bY,\bZ)$ by Monte Carlo.  The following three sections provide
greater  detail on the calculations involved.  Section~\ref{sec:Nghds} shows how to exploit the
conditional dependence structure of the graph in Figure~\ref{fig:overlapdigraph} to simplify the
calculation of Hastings ratios for $\bX'$ and $\bW'$.  Then Section~\ref{sec:propXW} gives a
prescription for the proposal distributions $q_\bX$ and $q_\bW$.  Finally, in
Section~\ref{sec:propLambda}, proposal distributions for $\blambda$ are considered, and a
Rao-Blackwellized estimator for $P(\blambda|\bY,\bZ)$ is given.  


\subsection{Neighborhood structures and joint probability ratios\label{sec:Nghds}} 
Making updates to $\bX$ and $\bW$
requires repeated calculation of the Hastings ratios (\ref{eq:XHR}) and (\ref{eq:WHR}).  This task is
made easy by proposing changes only to small parts (two components, for example) of either $\bX$ or $\bW$
at any one time.  The neighborhood structure inherent in the graph of Figure~\ref{fig:overlapdigraph} and
the fact that the probabilities described above can all be written in terms of a product over the $K$
alleles make this particularly attractive, as described below. 

Let the data from the genetic samples be considered fixed at $\bY$ and
$\bZ$,  and suppose the current values for $\bX$ and $\bW$ are denoted by  $\bX$
and
$\bW$, respectively.  Let $\bX^\prime$ and $\bW^\prime$ differ from $\bX$ and $\bW$
only at an arbitrary, single component subscripted by
$(t^\prime,a^\prime)\in\mathcal{T}$ for $\bX'$ and by $t'\in\{0,\ldots,T\}$ for $\bW'$.
In doing MCMC  we will make frequent use of the ratios 
\begin{equation}
\frac{P_\blambda(\bY,\bZ,\bX^\prime,\bW)}{P_\blambda(\bY,\bZ,\bX,\bW)}~~~~\mathrm{and}
~~~~
\frac{P_\blambda(\bY,\bZ,\bX,\bW^\prime)}{P_\blambda(\bY,\bZ,\bX,\bW)}.
\label{eq:TheRatios}
\end{equation} 
Calculating such ratios is done quickly by noting that they are functions only of a
small collection of variables adjacent in the graph to the altered component.   The
variables adjacent to the altered component in the graph are members of its
neighborhood, and the factors in the joint density including those neighbors are the
only ones that are changed by the alteration in that component.  Hence, the other
factors cancel out in the ratio.  The neighborhoods can be graphically
found and represented via the moralized, undirected  graph associated with the directed
graph \cite{Lauritzen1996}. The moralized subgraph around $\btX_{t',a'}$ (or $\btW_{t'}$) is formed  by
starting with the subgraph containing all variables which are either connected to  $\btX_{t',a'}$ (or
$\btW_{t'}$) by arrows in either direction or which 
are parents of any children of $\btX_{t',a'}$ (or
$\btW_{t'}$), and  then converting all the arrows between those variables to undirected edges and
moralizing the subgraph.  Moralizing is done by including edges between any unconnected parents in the
directed graph.  Directed and moralized versions of the subgraphs around $\btX_{6,3}$ and $\btW_3$ from
Figure~\ref{fig:undirgraph} are shown in Figure~\ref{fig:Moralized}.
\begin{figure}
\begin{center}
\subfigure[]{\includegraphics*[width = .49\textwidth]{AdultNghd.eps}}
\\
\subfigure[]{\includegraphics*[width = .49\textwidth]{AdultNghdUndirected.eps}}
\\
\subfigure[]{\includegraphics*[width = .49\textwidth]{GameteNeighborhood.eps}}
\\
\subfigure[]{\includegraphics*[width = .49\textwidth]{GamNghdsUndirected.eps}}
\end{center}
\caption[Neighborhoods for the allele count amongst the juveniles and adults]{Neighborhoods for the
allele count amongst the juveniles and adults. (a) and (b) are respectively the directed and the
moralized, undirected subgraphs for the relevant neighborhood in 
$\bX'$ with
$t'=6$ and
$a'=3$. (c) and (d) are the same for the neighborhood around $\btW_{3}$ (\ie $t'=3$). }
\label{fig:Moralized}
\end{figure}
The corresponding distribution associated with each undirected graph in the figure may
be factorized by their cliques (maximally connected subgraphs).  Therefore
the ratios in (\ref{eq:TheRatios}) may be written as ratios of terms corresponding to
the cliques.  Using the notation $\bX_{\{t',\backslash a'\}}$ to refer to the set
$\{\btX_{t',a^-},\ldots,\btX_{t',a^+}\}$, {\em excluding} $\btX_{t',a'}$, we have
\begin{eqnarray}
\frac{P_\blambda(\bY,\bZ,\bX^\prime,\bW)}{P_\blambda(\bY,\bZ,\bX,\bW)} 
& = &
\frac{ P_\blambda(\btW_{t'}|\btX'_{t',a'}, \bX_{\{t',\backslash a'\}}  )   }
     { 
      P_\blambda(\btW_{t'}|\btX_{t',a^-}, \ldots, \btX_{t',a^+}  )}
		   \times \frac{P(\btY_{t',a'}|\btX'_{t',a'})}{ P(\btY_{t',a'}|\btX_{t',a'})} 
\label{eq:XRatios}\\
 &  \times &\frac{P_\blambda(\btX'_{t',a'}|\btW_{t'-a'}) } 
{P_\blambda(\btX_{t',a'}|\btW_{t'-a'})  }    \nonumber  
\end{eqnarray}
for the ratio involving an altered version of $\bX$.  Note that if
$(t',a')\in\mathcal{P}$ then a term corresponding to the prior $P(\bX_\mathcal{P})$
would also appear in the ratio.  For the ratio involving the altered version of $\bW$
we have
\begin{eqnarray}
\frac{P_\blambda(\bY,\bZ,\bX,\bW')}{P_\blambda(\bY,\bZ,\bX,\bW)} 
&=& \frac{ P_\blambda(\btW'_{t'}|\btX_{t',a^-}, \ldots, \btX_{t',a^+}  )   }
     { 
      P_\blambda(\btW_{t'}|\btX_{t',a^-}, \ldots, \btX_{t',a^+}  )}  \times
		    \frac{P(\btZ_{t'}|\btW'_{t'})}{ P(\btZ_{t'}|\btW_{t'})} 
\label{eq:WRatios} \\
 &  \times & \prod_{(t,a)\in\mathcal{W}_{t'}}\frac{P_\blambda(\btX_{t,a}|\btW'_{t'})
}  {P_\blambda(\btX_{t',a'}|\btW_{t'})  }    \nonumber  
\end{eqnarray}
where $\mathcal{W}_{t'}$ represents the times and ages of adults descended from the
juvenile pool at time $t'$.  That is, $\mathcal{W}_{t'} = \{(t,a)\in\mathcal{T}: t-a
= t'\}$.

Let us now make the further restriction that $\bX'$ differs from $\bX$ only in two
components of $\bX'_{t',a'}$.  That is to say $X'_{t',a',i}$ and $X'_{t',a',j}$ can
take any non-negative values so long as $X'_{t',a',i}+X'_{t',a',j} =
X_{t',a',i}+X_{t',a',j}$.  We shall make a similar restriction on $\bW'$. In such a
case, the probability ratios in (\ref{eq:XRatios}) and (\ref{eq:WRatios}) simplify
further still, with the normalizing constants and the terms for unaltered allele counts
cancelling out. Hence we have
\begin{equation}
\frac{P(\btY_{t',a'}|\btX'_{t',a'})}{P(\btY_{t',a'}|\btX_{t',a'})} = 
\frac{X'_{t',a',i}!(X_{t',a',i}-Y_{t',a',i})!}
     {X_{t',a',i}!(X'_{t',a',i}-Y_{t',a',i})!} \cdot
\frac{X'_{t',a',j}!(X_{t',a',j}-Y_{t',a',j})!}
     {X_{t',a',j}!(X'_{t',a',j}-Y_{t',a',j})!} 
\end{equation}
when $(t',a')\in\syset$ and 1 otherwise.
A similar expression applies to 
$\frac{P(\btZ_{t}|\btW'_{t})}{P(\btZ_{t}|\btW_{t})}$ for sampling without
replacement from juveniles.   For sampling with replacement from juveniles we have
\begin{equation}
\frac{P(\btZ_{t'}|\btW'_{t'})}{P(\btZ_{t'}|\btW_{t'})} = \biggl(
\frac{W'_{t',i}}{W_{t',i}}\biggr)^{Z_{t',i}}
\biggl(
\frac{W'_{t',j}}{W_{t',j}}\biggr)^{Z_{t',j}}
\end{equation}
for $t'\in\szset$, and 1 otherwise.
For the terms having to do with population genetic sampling into the adult stage, we
have
\begin{equation}
\frac{P_\blambda(\btX'_{t',a'}|\btW_{t'-a'}) } 
{P_\blambda(\btX_{t',a'}|\btW_{t'-a'})  } 
=
\frac{\Gamma(X'_{t',a',i}+W_{t'-a',i}/\SRA_{t',a'})X_{t',a',i}!}
{\Gamma(X_{t',a',i}+W_{t'-a',i}/\SRA_{t',a'})X'_{t',a',i}!}
\cdot
\frac{\Gamma(X'_{t',a',j}+W_{t'-a',j}/\SRA_{t',a'})X_{t',a',j}!}
{\Gamma(X_{t',a',j}+W_{t'-a',i}/\SRA_{t',a'})X'_{t',a',j}!}
\end{equation}
and
\begin{eqnarray}
\frac{P_\blambda(\btX_{t'+a,a}|\btW'_{t'}) } 
{P_\blambda(\btX_{t'+a,a}|\btW_{t'})  } 
& = &
\frac{\Gamma(X_{t'+a,a,i}+W'_{t',i}/\SRA_{t'+a,a})\Gamma(W_{t',i}/\SRA_{t'+a,a})}
{\Gamma(X_{t'+a,a,i}+W_{t',i}/\SRA_{t'+a,a})\Gamma(W'_{t',i}/\SRA_{t'+a,a})}
\\
& \times &
\frac{\Gamma(X_{t'+a,a,j}+W'_{t',j}/\SRA_{t'+a,a})\Gamma(W_{t',j}/\SRA_{t'+a,a})}
{\Gamma(X_{t'+a,a,j}+W_{t',j}/\SRA_{t'+a,a})\Gamma(W'_{t',j}/\SRA_{t'+a,a})}.
\nonumber
\end{eqnarray}
For the terms having to do with population sampling into the juvenile stage we
compute the two relevant ratios using the quantity $\btB_t$ (defined on
Page~\pageref{page:Bdef} in Section~\ref{sec:SpecProbDist}) and its altered version $\btB'_t$ when
necessary. Thus we have
\begin{eqnarray}
\frac{ P_\blambda(\btW'_{t'}|\btX_{t',a^-}, \ldots, \btX_{t',a^+}  )   }
     { 
      P_\blambda(\btW_{t'}|\btX_{t',a^-}, \ldots, \btX_{t',a^+}  )} &=&
\frac{ P_\blambda(\btW'_{t'}|\btB_{t'}  )   }
     { 
      P_\blambda(\btW_{t'}|\btB_{t'}  )}  \\
& = & 
\frac{\Gamma(W'_{t',i}+B_{t',i}/\SRG_{t'}) W_{t',i}!}
      {\Gamma(W_{t',i}+B_{t',i}/\SRG_{t'}) W'_{t',i}!} \cdot
\frac{\Gamma(W'_{t',j}+B_{t',j}/\SRG_{t'}) W_{t',j}!}
      {\Gamma(W_{t',j}+B_{t',j}/\SRG_{t'}) W'_{t',j}!}  
\nonumber
\end{eqnarray}
and 
\begin{eqnarray}
\frac{ P_\blambda(\btW_{t'}|\btX'_{t',a'}, \bX_{\{t',\backslash a'\}}  )   }
     { 
      P_\blambda(\btW_{t'}|\btX_{t',a^-}, \ldots, \btX_{t',a^+}  )} &=&
\frac{ P_\blambda(\btW_{t'}|\btB'_{t'}  )   }
     { 
      P_\blambda(\btW_{t'}|\btB_{t'}  )}  \\
& = & 
\frac{\Gamma(W_{t',i}+B'_{t',i}/\SRG_{t'})\Gamma(B_{t',i}/\SRG_{t'})}
     {\Gamma(W_{t',i}+B_{t',i}/\SRG_{t'})\Gamma(B'_{t',i}/\SRG_{t'})}
\nonumber \\
& \times &
\frac{\Gamma(W_{t',j}+B'_{t',j}/\SRG_{t'})\Gamma(B_{t',j}/\SRG_{t'})}
     {\Gamma(W_{t',j}+B_{t',j}/\SRG_{t'})\Gamma(B'_{t',j}/\SRG_{t'})}.
\nonumber
\end{eqnarray}
Both of the above extend immediately to the case of Sampling Scheme~II with
$\btB_{t}$ defined appropriately (\ie with the $\btY_{t,a}$'s subtracted out as on
Page~\pageref{page:BSS2} in Section~\ref{sec:SpecProbDist}).

The derivation of the ratios of joint probabilities when non-penetrant alleles are present proceeds
in similar fashion to the above treatment, but is omitted for brevity.

\subsection{Proposal distributions for $\bX'$ and $\bW'$\label{sec:propXW}}
In the preceding, we saw that it is advantageous to consider changes to pairs of alleles
at a single time and age for $\bX$ and a single time for $\bW$.  Consequently the 
proposal distribution $q_\bX$ can be a function just of those components, and can be written 
$q_\bX(X'_{t,a,i},X'_{t,a,j}|X_{t,a,i},X_{t,a,j},\ldots)$.  Since $X'_{t,a,i}+X'_{t,a,j}$ must equal
$X_{t,a,i} + X_{t,a,j}$, the proposal distribution is simply a distribution on $X'_{t,a,i}$, imposing,
for uniqueness of reverse moves in this sampler, the condition $i<j$.


The proposal distribution $q_\bX(X'_{t,a,i}|X_{t,a,i},X_{t,a,j},\cdots)$  should reflect a compromise
between  statistical efficiency and computational efficiency.  From a statistical perspective, it is
most efficient to simulate $X'_{t,a,i},X'_{t,a,j}$ from their full conditional distribution.  However,
not much efficiency is gained this way, and calculating the
full conditional distribution incurs a heavy computational cost.  Instead, I define $q_\bX$ to be a
uniform distribution with width determined by the current values $X_{t,a,i}$ and $X_{t,a,j}$. That is,
$X'_{t,a,i}$ is drawn from a uniform distribution on the integers (excluding the current value,
$X_{t,a,i}$) between $X_\mathrm{lo}$ and $X_\mathrm{hi}$, inclusive, where the values of $X_\mathrm{lo}$
and
$X_\mathrm{hi}$ are chosen as a linear function of the approximate standard deviation of $X_{t,a,i}$
conditional only upon its parents in the graph.  Namely $X_\mathrm{lo} = X_{t,a,i} - w$ and
$X_\mathrm{hi} = X_{t,a,i} + w$ where $w$ is the greatest integer less than or equal to
\begin{equation}
 2 \beta \biggl(\frac{L + \alpha}{1 + \alpha} 
X_{t,a,i}( 1 - X_{t,a,i}/L) \biggr)^{1/2}
\end{equation}
where $L= X_{t,a,i} + X_{t,a,j}$, $\alpha = L/\SRA_{t,a}$, and $\beta$ is a scaling factor that may be
adjusted to achieve a desired acceptance proportion.  It can be tuned automatically during run time if
desired.  The width of $q_\bW$ may be tuned similarly.

It is also desirable to include some checking in the computer code to ensure that $q_\bX$ does not give
positive probability to any values of
$X'_{t,a,i}$ which would be incompatible with the descendants of $X_{t,a,i}$ and $X_{t,a,j}$ in the
graph.  


\subsection{Proposal distributions for $\blambda$\label{sec:propLambda}}
To make updates to $\blambda$, we consider changes to just one of its components at a time, $\lambda_a$
for the discussion here. A naive, computationally simple proposal distribution for $\lambda_a$ is less
desirable than a full conditional update for $\lambda_a$, because the latter allows for a
Rao-Blackwellized (see Section~\ref{sec:RBIntro} on Page~\pageref{sec:RBIntro}) Monte Carlo estimator of
$P(\blambda_a|\bY,\bZ)$.  This does require that the parameter space for $\blambda$ be discretized. 
This has little effect on the final inferences one can make if the discretization is fine enough.  For
example one could choose to consider $n$ values for
$\lambda_a$ say,
$\lambda_{a,0},\lambda_{a,1}\ldots,\lambda_{a,n}$ where $\lambda_{a,i} = .02 * i$.  For most
situtations, this will be a fine enough discretization.  Writing $\Lambda_a$ for the set
$\{\lambda_{a,0},\ldots,\lambda_{a,n}\}$, and $\blambda'$ for $\blambda$ with its $a\thh$ component set
to $\lambda'_a$,  I use for
$q_\blambda(\lambda'_a)$ the full conditional distribution
\begin{equation}
q_\blambda(\lambda'_a|\cdots) = P(\lambda_a|\bY,\bZ,\bX,\bW) = P(\lambda'_a|\bX,\bW) 
=
\frac{P(\lambda') P_{\blambda'}(\bX,\bW)}{\sum_{\lambda'_a\in\Lambda_a}P(\lambda')
P_{\blambda'}(\bX,\bW)}
\label{eq:lambdaFullCond}
\end{equation}
For each update to $\lambda_a$, (\ref{eq:lambdaFullCond}) must be computed for all
$\lambda'_a\in\Lambda_a$.  This is computationally expensive, but that is more than offset by the
fact that at the $i\thh$ update, one is realizing the values $P(\lambda'_a|\bX^{(i)},\bW^{(i)})$ for each
$\lambda'_a\in\Lambda_a$  with
$(\bX^{(i)},\bW^{(i)})$ being simulated from their posterior distribution given $\bY$ and $\bZ$. 
Therefore the successive values $P(\lambda'_a|\bX^{(i)},\bW^{(i)})$ may be averaged over the course of a
run of the Markov chain to yield an efficient, Rao-Blackwellized estimate
of
$P(\lambda_a|\bY,\bZ)$.  Furthermore, since
$q_\blambda(\lambda_a)$ is the full conditional distribution for 
$\lambda_a$, the above scheme defines a  
Gibbs sampling proposal for $\lambda_a$, and  the
Hastings ratio
$H_\blambda$ reduces to unity, always.  

 In empirical
tests, this method of updating $\blambda$ takes more computational time, but yields far superior
estimates of
$P(\lambda_a|\bY,\bZ)$ than a naive (\ie uniform) proposal distribution for each $\lambda_a$ in fewer
updates of the chain.  


\section{Special Cases}
There are some situations in which it might be advantageous or imperative to consider a model which is
simpler and has fewer parameters than the one just described.  One obvious simplification would be to
restrict the $\lambda_a$'s of each age group to be equal. This would be appropriate if data were only
available on juveniles, since in that case the different $\lambda_a$'s would be
unidentifiable.  It might also be desirable if data are
relatively sparse, and/or if one has prior reason to believe that $\lambda_a$'s would not differ
greatly over different age classes.

Similarly, it is possible to restrict $\bW$ to match perfectly the allele frequencies implied 
by $\bX$.  This corresponds to the assumption that juvenile populations are very large and all of the
population-genetic sampling that is {\em not} random with respect to an individual's family of origin
occurs in the mortality between juvenile and adult stages.  While this restriction would not allow the
independent estimate of a
$\lambda$ for juveniles, it would be prudent in the case when data are available only on adults or
only on juveniles.  Any non-random (with respect to family) sampling that occurred before the juvenile
stage would then be estimated as part of the population-genetic sampling from juvenile to adult.  


\section{Simulated Data\label{sec:OverlapSim}}
It is always the case that the process of formulating and  describing this type of MCMC method is far
simpler and takes less time than the cycle of implementation, testing, and debugging that is required to
produce software to actually carry out the Markov chain simulations described.  I have not yet been able
to test all the parts of the current version of the software implementing the
method described in the chapter.  However, I am satisfied that some of its modules are functioning
properly, and the results are sufficiently promising to present the method's performance on simulated
data.  For this purpose I have used the census  data  in
Table~\ref{tab:Inmaha}, from the Inmaha Creek chinook salmon population.  These data appear in
\citeN{Beamsderferetal1998}, and were kindly provided to me in electronic format by Robin Waples at the
National Marine Fisheries Service.   As in Chapter~\ref{chap:ForBack}, the purpose of this demonstration
is not to assess the bias and variance of the Bayesian estimator for $\lambda$ derived in this chapter. 
Doing so would require a prohibitive amount of computing in order to average the results over a large
number of simulated datasets.  Rather, this section demonstrates that the MCMC method itself is able to
provide a good approximation to the posterior probability for $\lambda$ given a single dataset.  

Inmaha Creek is a tributary of the Snake River in which the chinook salmon population has declined
dramatically in the last four decades.  Fish return at ages 3, 4, and 5.  The 3-year-olds are almost all
small males called ``jacks."   Census size estimates, broken
down by age class, are available for the years 1954 to 1999.  While jacks certainly contribute some to
future generations, it is unlikely that the contribution, on a per-fish basis, is nearly as great as that
of four and five year-olds.  Further, since there are  so few of them, and because their occasional zero
census size estimates cause conflicts with the current version of my software, they were
excluded from the dataset.  

\begin{table}
\caption[Census sizes of Inmaha Creek chinook salmon]{Estimates of the number of
chinook spawners returning to Inmaha Creek (a tributary on the Snake River
drainage) in years 1954 to 1999.  Age 3 fish are young males known as ``jacks."  
(Data source: \citeN{Beamsderferetal1998})
}
\label{tab:Inmaha}
\begin{center}
\input{InmahaTab.tex}
\end{center}
\end{table}

Genetic data were simulated for a single locus given these census sizes by initializing a juvenile pool
in year 1949 with 5 alleles having counts in the proportions $(.4,.2,.2,.1,.1)$.  Allele
counts in the juvenile pool in years 1950 to 1953 were then considered to be $\btW_{-4}$ to
$\btW_{-1}$ and were drawn according to the urn scheme describing the prior distribution
(Section~\ref{sec:OverPriors}) with
$\tilde{C}$ being 900.  In other words, $\btW_{-4}$ to
$\btW_{-1}$ were simulated by independent, random draws from an urn containing  alleles
in the initial  frequencies, $(.4,.2,.2,.1,.1)$.   The remaining
latent variables
$\bX$ and
$\bW$ were  simulated throughout the graph via the urn scheme described in this chapter, with
$\lambda_4=\lambda_5 = 0.4$, and using the age-specific fitnesses of $\gamma_4 = 450$ and
$\gamma_5=650$. These values were obtained by using rough fecundity/length, length/age, and juvenile
survivorship relationships for chinook salmon (both stream- and ocean-type combined) given in
\citeN{Healey1991}.   The latent data were simulated under the assumption that no
genetic drift occurs between the adult and the juvenile stage.     Genetic data were not simulated  from
1954 to 1963.  However, genetic drift was simulated in the population during that interval.  This
allowed the allele frequencies between different years to settle closer to their joint stationary
distribution before starting the simulated sample collection.   Other simulations (Robin Waples,
National Marine Fisheries Service, unpublished result) show that 20 years is  sufficient to allow the
alleles frequencies to ``warm-up."   For the purposes of the present demonstration, ten years should be
sufficient.

   From  1963 to 1988 I simulated datasets with samples of varying sizes drawn every year from the same
simulated set of latent variables.  The three different sample sizes considered were:
\begin{enumerate}
\item $\SA_{4}=\SA_5 = 10$  and 
$\SG=30$
\item $\SA_{4}=\SA_5 =25$ and
$\SG=60$
\item
$\SA_{4}=\SA_5 =60$ and
$\SG=125$
\end{enumerate}
In years when the sample size would have been larger than one half the census size of the population of
a particular age (4 or 5), the sample size for that age group was decreased to be one half of the census
size of the population. Data were not simulated and used for the last eleven years (1989--1999) of the
census data because the small population sizes in those years meant that even with very small samples
from the adult populations, a good estimate of
$\blambda$ was possible, and I wanted a more challenging scenario for demonstrating the method.


The simulated data were analyzed under the assumption that $\lambda_4 = \lambda_5$ (which shall
hereafter be referred to as $\lambda$) and that no drift occurs in the transmission of genes to the
juvenile stages; hence $\btW_t$ is completely determined by $(\btX_{t,a^-},\ldots,\btX_{t,a^+})$. 
$\lambda$ values in the set
$\{.02,.04,\ldots,.98\}$ were considered. To reduce burn-in time,
$\bX$ was initialized to the value that was realized in the original simulation.  I have subsequently
verified that the burn-in time required for other reasonable starting configurations (like all allele
frequencies of $\bX$ initialized to the average frequency of the alleles observed in the samples) is
short.  A sweep of the algorithm consisted of:
\begin{enumerate}
\item $E$ updates in series, first with a random pair $(X_{t,a,i},X_{t,a,j})$ between the years 1963
and 1988 and then with  a random pair  $(W_{t,i},W_{t,j})$ from the juvenile pools used to construct
the prior distribution in the years 1958 to 1962.
\item An update of $\lambda$.
\end{enumerate}
$E$ was chosen so that each component of $\bX$ was updated twice on average during a sweep.  For the
different scenarios I simulated, I performed 70,000 sweeps of the algorithm.  Inspection of the
estimated posterior for $\lambda$ suggests that the estimate changed imperceptibly over the last
50,000 sweeps of the algorithm.  70,000 sweeps required 2 hours on a laptop computer with a 266 Mhz
G3 (Macintosh)
processor. 

At each of the different sampling intensities I analyzed the data under the assumption that all the
samples were available (Figure~\ref{fig:4HistsTogether}), and also under the assumption that only the
adult samples were available (Figure~\ref{fig:NoZY}(a)).  I also did one simulation in which
samples from the adults were not available, but samples of size
$\SG=125$ from juveniles at all years were available (Figure~\ref{fig:NoZY}(b)).  For comparison, I have
plotted each of these posterior distributions next to the posterior distribution that one
would obtain if $\bX$ and $\bW$ were known without error. 

\begin{figure}
\begin{center}
\includegraphics*[width=.9\textwidth]{4HistsTogether.eps}
\end{center}
\caption[ Posterior distributions for $\lambda$ with samples from all groups available]{Plot of the
posterior probabilities for
$\lambda =\lambda_4 =\lambda_5$ from genetic data simulated on the Inmaha Creek chinook population
(Table~\ref{tab:Inmaha}), when data were available on both adults and juveniles.  The first line in the
graph corresponds to the estimate with the latent variables known without error.  The other three lines
correspond to the different sample sizes of adults and juveniles.  The true $\lambda_4 = \lambda_5 =
0.4$.}
\label{fig:4HistsTogether}
\end{figure}


%\begin{figure}
%\begin{center}
%\subfigure[$\SA = 60$, $\SG = 125$]{\includegraphics*[width=.8\textwidth]{ProbSample60_125.eps} } \\
%\subfigure[$\SA = 25$, $\SG = 60$]{\includegraphics*[width=.8\textwidth]{ProbSample25_60.eps} } \\
%\subfigure[$\SA = 10$, $\SG = 25$]{\includegraphics*[width=.8\textwidth]{ProbSample10_30.eps} }
%\end{center}
%\caption{300 distributions drawn from the sampler\ldots }
%\end{figure}


\begin{figure}
\begin{center}
\subfigure[Samples taken only from adults]{\includegraphics*[width=.8\textwidth]{4Hists_NoZ.eps} } \\
\subfigure[Samples taken only from juveniles]{\includegraphics*[width=.8\textwidth]{1HistNoXs.eps} }
\end{center}
\caption[Posterior distributions for $\lambda$ under different sampling scenarios]{Graphs as described in
Figure~\protect\ref{fig:4HistsTogether}, but under different sampling scenarios.  (a) Juvenile sample
sizes all zero, and adult sample sizes as shown.  (b) Adult sample sizes not all zero, but samples of
125 taken each year from the juveniles. }
\label{fig:NoZY}
\end{figure}


The results, as shown in Figures~\ref{fig:4HistsTogether} and~\ref{fig:NoZY}, suggest that
the MCMC sampler is functioning appropriately and computing the posterior distribution for $\lambda$. 
The curvature generally decreases with sample size, reflecting the loss of information, as it
should.\footnote{The posterior distribution for sample sizes $\SA = 25$, $\SG=60$, being more peaked
than the posterior distribution for   $\SA = 60$, $\SG=125$, is an exception to the trend.  This results
from the fact that for the particular set of data simulated for $\SA = 25$, $\SG=60$, the estimated
$\lambda$ happens to be smaller than for the data simulated with $\SA = 60$, $\SG=125$.  The credible
set will be smaller for a lower estimated value of $\lambda$ because, when $\lambda$ is smaller, then
the amount of genetic drift expected will be larger, relative to the amount of error due to random
sampling of genes.  In other simulations (not shown) in which the maximum {\em a posteriori} estimate
of $\lambda$ for the simulated data with
$\SA = 25$,
$\SG=60$ was closer to that for the simulated data with $\SA = 60$, $\SG=125$, the credible set is wider
for the dataset with smaller samples. }  Furthermore, in all cases except one, the 90 percent credible
interval for
$\lambda$ overlaps the true value of 0.4.  While the narrowest posterior distributions occur with
samples  from both adults and juveniles, there still seems to be a substantial signal in the data, even
when samples are taken either only from adults or only from juveniles. 


\section{Discussion}
These results are encouraging.  They
demonstrate that the MCMC sampler devised here is able to compute the posterior probability distribution
for $\lambda$, suggesting  that the method presented in this chapter permits use of data over
multiple years from a salmon population with known census sizes to estimate the ratio
$\lambda$ with good precision.   It should be kept in mind that
these simulations exploit the data from only a single locus with five alleles.  Narrower credible sets
would be obtained with data on multiple loci.  The posterior distribution for $\blambda$ given data on
multiple, independently segregating loci is proportional to the product of the posterior probabilities
for $\blambda$ from each of the loci treated separately, as described here. Therefore, the extension to
multiple loci is simple.


The method developed herein would be particularly appropriate for estimating $\blambda$ in hatchery
populations of salmon where the census sizes of spawning adults are well known.  As in the previous
chapter, the method thus far developed in the current chapter assumes that $\blambda$ remains
constant over time.    Future work is required to assess how robust this estimate is to departures from
the underlying model.   However, like the method of Chapter~\ref{chap:Urns}, it would
also be possible here to propose new models in which $\blambda$ varied over time, and to compare those
models within a Bayesian framework using reversible jump MCMC \cite{Green1995}.  Such a method would be
well-suited to using genetic data to detect the impact of supportive breeding programs
\cite{Rym&Lai1991,Hansenetal2000} on
$\blambda$ in salmon populations.  