\chapter{$\lambda$ and a P\'{o}lya Urn Model for Genetic Inheritance}
\label{chap:Urns}

\section{Introduction}
In the previous chapter, we saw the Wright-Fisher model for genetic
inheritance in a randomly-mating population.  This is one of the
simplest stochastic models possible for the transmission of genes between
generations without a specified pedigree structure. Under the
Wright-Fisher model, the marginal distribution of the number of offspring of an
individual  is binomial.  Its use in determining the effective size of a
natural population involves a considerable degree of abstraction---the
many complicated interactions, stochastic events, and fitness
considerations that lead to extra-binomial variance in family size or a
certain rate of increase in inbreeding are presumed to be 
adequately accounted for by considering that the natural population can be
modeled as a Wright-Fisher population of a particular size. 

Explicitly modeling the many factors  influencing inter-generation
genetic dynamics in a natural population would be very challenging. 
Fortunately, for many natural populations, the Wright-Fisher population
of size $N_e$ provides an admirable approximation.  However, when data
are available on the census size, $C$, of a population, biologists are
often interested in estimating the ratio of effective size to census
size, $N_e/C$, which we shall denote $\lambda$. Pursuing likelihood or
Bayesian inference for $\lambda$ while strictly adhering to the
Wright-Fisher model requires a somewhat inelegant interpretation of
population size and leads to difficulties in applying Markov chain Monte Carlo
methods to the problem.  Namely, the probability of the unobserved allele counts in
the population  is not well-defined for different values of $\lambda$ under the
Wright-Fisher model.  For this reason,   I propose a new model for genetic inheritance based
on an urn sampling scheme with stochastic replacements.  Given the census size of a breeding
population, this model has a single parameter, $s$, which can change the effective size of
the population, without altering the number of individuals in that population.  This
provides  several advantages in the inference problems pursued here, simplifying the
interpretation of the modeled process and the computation necessary for estimating
$\lambda$.  Results from this model, however, may be easily translated back to an
interpretation from within the familiar Wright-Fisher perspective.
 
In this chapter I will briefly describe the Wright-Fisher based model for
estimating $\lambda$, and a single-site-updating scheme for a Markov chain Monte Carlo
(MCMC) algorithm suitable for approximating likelihood ratios or posterior
probabilities for $\lambda$. I will then illustrate the shortcomings of
the Wright-Fisher model in this context.  This motivates the development
in Section~\ref{sec:UrnDev} of the urn model.  I devote
Section~\ref{sec:TwoStage} to alternative interpretations and developments of
the urn model which provide a wider range of biological interpretations for it. 
In Section~\ref{sec:UrnCompare}, I derive expressions for the inbreeding and
variance effective sizes of the urn-model population, investigate
probabilities of allele fixation in the urn model, and show the relationship
between census sizes
$C_t$ through time, the parameter of the urn model, $s$, and the parameter
$\lambda$.  In Section~\ref{sec:Fixation}, I argue that in terms of allele
fixation probabililities the Wright-Fisher model corresponds to possibly unrealistic
assumptions about the distribution of family sizes.  Finally, I briefly describe
the use of the urn model to develop a simple MCMC scheme for Bayesian inference
of
$\lambda$.  I apply this to the data of \citeN{Begonetal1980} from the preceding
chapter.  Elaborations on this MCMC method involving different sampling times
and schemes, and different life-histories, are presented in the following
chapter.

\section{Estimating the Ratio of Effective to Census Population Size}
In this section, I consider the scenario in which census sizes of breeding
adults have been recorded each generation over some time period in a
semelparous population with discrete generations and genetic samples are
taken with replacement at some intervals of time either from those breeding
adults or from the juveniles descended from them (other sampling models
will be entertained in the following chapter).  Much of this treatment is adapted
from
\citeN{And&Tho1999}.


At generation $t$,  $C_t$ diploid
individuals reproduce, giving rise to $C_{t+1}$ individuals
at generation $t+1$.  We take genetic samples of size
$S_{0},\ldots,S_{T}$ (assume $S_0>0,\
S_T>0$) individuals, and find counts of the $K$ different
allelic types at a locus,
$\bY = (\btY_{0},\ldots,\btY_{T})$ where
$\btY_{t} = (Y_{t,1},\ldots,Y_{t,K})$.  For estimating $\lambda$ the
population at
$t$ is modelled as
$\floor{\lambda C_t}$ ideally-reproducing (\ie via a Wright-Fisher
inheritance model) adults, where
$\floor{x}$ is the largest integer $\leq x$.  Underlying the data are
latent allele counts $\bX = (\btX_{0},\ldots,\btX_{T})$ with
$\btX_{t} = (X_{t,1},\ldots,X_{t,K})$. $\{\btX_t,\;t\geq 0\}$ is a
first-order Markov chain with transition probabilities
$P_\lambda(\btX_{t+1}|\btX_{t})$ being multinomial with cell
probabilities $\btX_{t}/\floor{2\lambda C_{t}}$ and number of
trials $\floor{2\lambda C_{t+1}}$.  The genetic data at time $t$ are
 assumed to be samples from the gamete pool produced by the
$\floor{\lambda C_{t}}$ adults.  This assumption is equivalent to sampling
adults with replacement, and closely approximates the sampling of
juveniles, when the number of juveniles is large. Hence, for
$S_t>0$,
$P_\lambda(\btY_{t}|\bX) = P_\lambda(\btY_{t}|\btX_{t})$ is
multinomial with parameters $\btX_{t}/\floor{2\lambda C_{t}}$ and
$2S_{t}$.  For $S_t=0$, $P_\lambda(\btY_{t}|\btX_{t})\equiv 1$.  Summing out the
nuisance parameters
$\btX_0$ over a uniform prior $P_\lambda(\btX_0)$ gives the likelihood
\begin{eqnarray}
L(\lambda) = P_\lambda(\bY) & = & \sum_\bX P_\lambda(\bY,\bX) 
\label{eq:ch3Like} \\
& = &
\sum_{\btX_0,\ldots,\btX_T} P_\lambda(\btX_0)P_\lambda(\btY_{0}|\btX_{0})
\prod_{t=1}^T P_\lambda(\btX_t|\btX_{t-1})
P_\lambda(\btY_{t}|\btX_{t}). \nonumber
\label{eq:lik}
\end{eqnarray}
With
$K=2$, the sum over $\bX$  may be evaluated exactly.  With larger
$K$, however, the huge space of possible
$\btX_t$'s makes this infeasible.

\subsection{MCMC likelihood for $\lambda$\label{sec:MCMCLikeLambda}}
To obtain an efficient Monte Carlo estimate of $L(\lambda)$, one may
consider the likelihood ratios $L(\lambda)/L(\lambda_0)$, 
\cite{Tho&Guo1991,Gey&Tho1992}
\begin{equation}
\frac{L(\lambda)}{L(\lambda_0)} =
\frac{P_\lambda(\bY)}{P_{\lambda_0}(\bY)} = \sum_\bX 
\frac{P_\lambda(\bY,\bX)}{P_{\lambda_0}(\bY,\bX)}
P_{\lambda_0}(\bY|\bX)=
E_{\lambda_0}\left(\left.\frac{P_\lambda(\bY,\bX)}
{P_{\lambda_0}(\bY,\bX)}\right|\bY\right)
\label{eq:ThoandGuo}
\end{equation}
which may be estimated by $\frac{1}{m}\sum_{i=1}^m 
P_\lambda(\bY,\bX^{(i)})/P_{\lambda_0}(\bY,\bX^{(i)})$ where each
$\bX^{(i)}$ is  realized from $P_{\lambda_0}(\bX|\bY)$. This is an
efficient Monte Carlo estimator of the likelihood ratio provided
$\lambda$ is near to
$\lambda_0$. Independent
samples of
$\bX$ are not available because
$P_{\lambda_0}(\bX|\bY)$ is known only up to scale.  Instead, values of 
$\bX^{(i)}$ can be  realized from a Markov chain with limit distribution
$P_{\lambda_0}(\bX|\bY)$ using 
 a component-wise Metropolis-Hastings
algorithm (Hastings 1970):  Start with initial values of $\bX$;
Select a pair
$(X_{t,k},X_{t,\ell})$,
$k\neq
\ell$  at random from $\bX$;  Propose updating the pair to
$(X^*_{t,k},X^*_{t,\ell})=(X_{t,k}-w,X_{t,\ell}+w)$, where $w$ is a random integer
drawn with probability $q(w|X_{t,k},X_{t,\ell})$;   
accept the proposal with probability 
$\min\{1,[q(-w|X^*_{t.k},X^*_{t,\ell})P_{\lambda_0}(\bY,\bX^*)]/
[q(w|X_{t,k},X_{t,\ell})  P_{\lambda_0}(\bY,\bX)] \}$.  After initial 
updates for burn-in, $\bX^{(i)}$'s are sampled as the state of $\bX$
at a spacing of $u$ updates. 


When estimating a curve for $L(\lambda)$, the range of $\lambda$'s of
interest may be large.  In such a case it does not suffice to realize
$\bX^{(i)}$'s under a single $\lambda_0$. Instead, one must sample from
several chains, each indexed by a different $\lambda_0$,
$\lambda_0\in\Lambda$.  \citeN{GeyerTech} describes a reverse logistic
regression method for reweighting the samples from each chain and
estimating the whole likelihood surface.

\subsection{A Bayesian approach\label{sec:ABayesianApproach}}
The MCMC sampler described above is useful for Bayesian inference of
$\lambda$ as well.  Suppose that we are interested in computing the
posterior probabilities for each $\lambda\in\Lambda$.  Let $P(\lambda)$
be the prior distribution for $\lambda$ assigning probability mass to each
of the discrete points $\lambda\in\Lambda$.  We shall assume for brevity
and simplicity that $\sum_{\lambda\in\Lambda} P(\lambda) = 1$.  It is
possible to also propose changes to $\lambda$ in the MCMC sampler, and
thereby sample from the posterior distribution of $\lambda$ by the
following scheme:  given the current state of the chain, ($\bX,\lambda$),
propose a new value $\lambda^*$ with probability $h(\lambda^*|\lambda)$
from the proposal distribution $h(\cdot|\lambda)$.  Then, accept the proposal with
probability 
\[
\min\left\{1,\frac{h(\lambda|\lambda^*)}{h(\lambda^*|\lambda)}
\frac{P(\lambda^*)}{P(\lambda)}
\frac{P_{\lambda^*}(\bY,\bX)}{P_{\lambda}(\bY,\bX)} \right\}
\]
The successive values $\lambda^{(i)}$ generated in this way are a
dependent sample from $P(\lambda|\bY)$ and may accordingly
be used to estimate that posterior distribution. 

\subsection{Shortcomings of the Wright-Fisher model in this case}
The first difficulty in using the Wright-Fisher model to formulate the
likelihood given in (\ref{eq:ch3Like}) is the dependence of the sampling
terms $P_\lambda(\btY_t|\btX_t)$ on $\lambda$.  The parameter $\lambda$ is
supposed to affect the transitions $P_\lambda(\btX_{t+1}|\btX_{t})$, so it
is inelegant to also have the dependence appear in
$P_\lambda(\btY_t|\btX_t)$.  This dependence results from the
fact that a given allele count $X_{t,k}$ corresponds to a different allele
{\em frequency} depending on the value of $\lambda$.  Since the likelihood
is the sum over all $\bX$, this dependence has little real effect,
but its presence is enough to cause people to take pause\footnote{The
audience at the 2000 PMMB short course in Berkeley apparently found this
unwanted dependence to be worrying when Elizabeth
presented this work.}. 

The undesireable effect of this dependence becomes even more
clear when one considers the possibility of sampling adults without
replacement {\em before} reproduction---Sampling Plan~II of
\citeN{Nei&Taj1981}, also described in \citeN{Waples1989}.  Under such
sampling, in the context of a population with very high variance in family
size, it is possible that the effective size estimated for a population
could be smaller than the number of adults sampled without replacement
from the population modeled as a Wright-Fisher population.  While this presents
no problem conceptually---it is quite possible that the effective size
of a population may be smaller than the number of individuals sampled from it---it exposes
how difficult it would be to correctly model Sampling Scheme~II within the MCMC scheme of
the previous chapter.  Trying to explicitly model sampling of adults without replacement in
the above model becomes difficult because $\btX_t$ is affected by
the sampling process, and the degree to which it is affected depends not on $\floor{\lambda
C_t}$ but on $C_t$ itself.  Suffice it to say that pursuing the explicit modeling of Sampling
Scheme~II within the context of the Wright-Fisher model of genetic inheritance, while
possible, would not be straightforward.  

A further shortcoming of the Wright-Fisher model is even more
practically problematic---different values of $\lambda$ imply different sizes of the
state space of the latent variable $\bX$, and this makes the computation
of $L(\lambda)/L(\lambda_0)$ or of $P(\lambda|\bY)$ more difficult and
less efficient.   Consider first the computation of
$L(\lambda)/L(\lambda_0)$ using the identity (\ref{eq:ThoandGuo}).  This
identity only holds when the support of $\bX$ under $\lambda_0$ is
equal to or contains the support of $\bX$ under $\lambda$.  This only
transpires when $\lambda_0 > \lambda$, since if $\lambda> \lambda_0$, it is
possible that a value of $\bX$ for which $P_\lambda(\bY,\bX) > 0$  might
include a $t$ for which 
$\sum_{k=1}^{K-1}X_{t,k} > \floor{2\lambda_0C_t} $, and hence
$P_{\lambda_0}(\bY,\bX)$ would be zero.  As a consequence, the importance
sampling algorithm implied by (\ref{eq:ThoandGuo}) can be employed only
for $\lambda_0>\lambda$.  This is inefficient because it is not possible
to use the realizations generated at values of $\lambda_0<\lambda$ to
estimate $L(\lambda)/L(\lambda_0)$.  It should be pointed out that a
similar problem would occur with the importance sampling method of the
previous chapter if one tried to combine the realizations from
$P^*_{N_e}(\bX)$ with realizations from another distribution, say
$P^*_{N_{e0}}(\bX)$ ($N_{e0}< N_e$) to estimate the likelihood at $N_e$. 

Additionally, it is not immediately clear how best to compute  $P_\lambda(\bY,\bX^{(i)}) /
P_{\lambda_0}(\bY,\bX^{(i)})$ even
when
$\lambda_0>\lambda$.  It is easy to compute $P_{\lambda_0}(\bY,\bX^{(i)})$, but not
so for $P_\lambda(\bY,\bX^{(i)})$ when $\lambda \neq \lambda_0$.  One
approach we have used is to ``let the last allele at each locus take up the slack."  In
other words,  when $\bX^{(i)}$ has been simulated from a chain under $\lambda_0$, then 
$P_\lambda(\bY,\bX^{(i)})$ is computed by replacing each component, $\btX^{(i)}_t$,
of $\bX^{(i)}$ with a new vector differing in the number of alleles subscripted by
$K$.  That is, the new vector is
$(X^{(i)}_{t,1},\ldots,X^{(i)}_{t,K-1},X^{(i)}_{t,K}-\floor{2\lambda_0 C_t}
+\floor{2\lambda C_t} )$.  While this sort of scheme seems to work
reasonably well in practice,  the quantity
$P_\lambda(\bY,\bX^{(i)})/P_{\lambda_0}(\bY,\bX^{(i)})$ computed in this way is not
invariant to permutations of the order of alleles, and so it is
unattractive.  Since the ratio
$P_{\lambda^*}(\bY,\bX)/P_{\lambda}(\bY,\bX)$ must also be computed for the
Bayesian approach to estimating $\lambda$ described above, this problem
also plagues the Bayesian approach using the Wright-Fisher model.

In short, to compute the likelihood or the posterior
distribution for
$\lambda$ using MCMC, it is desirable to be able to propose a reasonable
model under which changes in the value of $\lambda$ do not change the size
of the state space of
$\bX$.  Such a model is derived in the next section.


\section{The Urn Model\label{sec:UrnDev}}
From the preceding section, it should be clear that we desire a model of genetic
inheritance that allows a population of effective size $N_e$, {\em but census
size} $C$, to give rise to a new generation with genetic dynamics (increase in
allele frequency variance, inbreeding, \etc) characteristic of a Wright-Fisher
population of size $N_e$.  The important, and more difficult part, is that we
require doing this {\em without modeling the population as having any size other
than
$C$}. If our primary concern was that of matching the increase in allele
frequency, we would seek a genetic inheritance model with the same number of
individuals that we observe in the census, but with an increased variance in
progeny allele frequency.   An obvious way to do this would be to allow
allele frequency in the following generation to follow a scaled beta-binomial
distribution, rather than a scaled binomial distribution, because the
beta-binomial is a classical discrete distribution with extra-binomial
variance.  While such a distribution is convenient from the perspective of
allele frequency, it is also the case that it arises from a system of mating
that we can specify in terms of a sampling process on individual gene copies.
In this section I will describe an urn sampling scheme from which it arises. In
Section~\ref{sec:TwoStage} I will describe another genesis of the same
distribution.  

Let us consider intergenerational sampling from generation~1, with $C_1$
reproductive adults, to their offspring who become $C_2$ reproductive
adults.  Let there be $K$ alleles, indexed by $i$, with numbers of each at
time~1 given by $n_i$ $(\sum_{i=1}^K n_i = 2C_1$) and at time~2 by $x_i$
$(\sum_{i=1}^K x_i = 2C_2$).  By our urn scheme, the alleles in
generation~2 are drawn by the following scheme:  a gene copy is drawn at
random from those present in generation~1, and a copy of it is placed in
generation~2.  Then the original gene copy is returned to generation~1,
along with $s$ new copies of it.  This defines a multivariate
P\'{o}lya-Eggenberger urn scheme \cite{Joh&Kot1977}, where the genes of different
allelic type may be regarded as balls of different color.  By this scheme,
$\btX=(X_1,\ldots,X_K)$, conditional on the $n_i$, follows the compound multinomial
Dirichlet distribution, which is a multivariate generalization of the
well-known beta-binomial distribution \cite{Johnsonetal1997}.

The compound multinomial distribution with $2C_2$ trials arises as the
marginal distribution of $\btX$ from the hierarchy
\begin{eqnarray}
\bm{Q} & \sim & \mathrm{Dirichlet}(\alpha_1,\ldots,\alpha_K)  \nonumber \\
\btX|\bm{Q} & \sim & \mathrm{Mult}_K(2C_2,\bm{Q}), \nonumber
\end{eqnarray}
and has the probability mass function
\begin{equation}
P(\btX|2C_2;\alpha_1,\ldots,\alpha_K) =
\frac{(2C_2)!\;\Gamma(\alpha_\bigdot)}{\Gamma(2C_2 +
\alpha_\bigdot)}
\prod_{i=1}^K \biggl( \frac{\Gamma(X_i + \alpha_i)}{X_i!\;
\Gamma(\alpha_i)}
\biggr)
\label{eq:CMD}
\end{equation}
where $\Gamma(y) = \int_0^\infty t^{y-1}e^{-t}dt$ is the gamma function, $\alpha_\bigdot =
\sum_{i=1}^K
\alpha_i$, and, as before
$\sum_{i=1}^K x_i =   2C_2$.  For the urn sampling scheme described
above, $\btX$ has p.m.f (\ref{eq:CMD}) with $\alpha_i = n_i/s$,
$i=1,\ldots,K$.

Conceptually, it is apparent that $s$ plays an important role in
determining the variance in family size in this urn-model for genetic
inheritance. If $s=0$ we have sampling with replacement---the
fundamental feature of the Wright-Fisher model.  However, if $s>0$ then
each time a gene copy in generation~1 is selected to produce a daughter
in generation~2, the probability is increased that that particular gene
copy will have more daughters (and the probability is increased that
other individuals will have fewer daughters).  So, when
$s>0$ the variance in offspring number is increased (and the effective
size of the population is thus decreased).  The situation of $s<0$
corresponds to an effective size which is larger than the census size of
the population, as is sometimes achieved by a prescribed system of
breeding in captive populations.  $s=-1$ is hypergeometic sampling of gene
copies into the next generation. In an urn-model population of
constant size, this would correspond to every gene copy being copied
exactly one time into the next generation.  In the remainder of this
chapter, however, we will concentrate on cases of $s\geq 0$.  Though this 
model was developed, conceptually, in terms of stochastic replacements that were in
integer units of gene copies, there is no mathematical restriction that 
$s$ be an integer; non-integer values for $s$ are permissible.

It is possible to derive standard results about the
increase of allele frequency variance and the probability of identity by
descent from a generation of genetic sampling via this urn model.  This will be done in
Section~\ref{sec:UrnCompare}. These calculations are simplified by the fact that the marginal
distribution of the
$i\thh$ component in the compound multinomial distribution has a
beta-binomial distribution with parameters $2C_2$, $\alpha_1$ and
$\alpha_\bigdot - \alpha_i$.  Thus it has  p.m.f.
\begin{eqnarray}
P(X_i |2C_2;\alpha_1,\ldots,\alpha_K)& =&  P(X_i
|2C_2;\alpha_i,\alpha_\bigdot - \alpha_i) \nonumber \\ 
& = &
\frac{(2C_2)!\;\Gamma(\alpha_\bigdot)}{\Gamma(2C_2 + \alpha_\bigdot)}\times
\frac{\Gamma(X_i + \alpha_i)}{X_i!\; \Gamma(\alpha_i)}\times  \mbox{}
\nonumber \\
& &
\frac{\Gamma(2C_2-X_i + \alpha_\bigdot -
\alpha_i)}{(2C_2-X_i)!\;\Gamma(\alpha_\bigdot-\alpha_i)}
\label{eq:betabomb}
\end{eqnarray}
and first two central moments
\begin{eqnarray}
\Exp X_i& = &2C_2\frac{\alpha_i}{\alpha_\bigdot}  \label{eq:ExpBB} \\
\Var(X_i) &=& 2C_2\biggl(\frac{2C_2 + \alpha_\bigdot}{1 + \alpha_\bigdot}
\biggr)\biggl(\frac{\alpha_i}{\alpha_\bigdot}\biggr)
\biggl(\frac{\alpha_\bigdot-\alpha_i}{\alpha_\bigdot}\biggr).
\label{eq:VarBB}
\end{eqnarray}
Any two components of a compound multinomial distribution are correlated, having
covariance 
\[
\Cov(X_i,X_j) = 2C_2 \biggl(\frac{2C_2 + \alpha_\bigdot}{1 +
\alpha_\bigdot}\biggr)
\biggl(\frac{\alpha_i}{\alpha_\bigdot}\biggr)
\biggl(\frac{\alpha_j}{\alpha_\bigdot}\biggr).
\]
Derivation of the preceding results may be found in
\citeANP{Johnsonetal1997}~(\citeyearNP{Johnsonetal1997}, pg.~81,82), where the
authors also point out that the variance-covariance matrix for a compound
multinomial random variable can be written as the product 
\label{page:FormOfCMDVar} of $(2C_2 +
\alpha_\bigdot)/(1 + \alpha_\bigdot)$ and the variance covariance matrix for a
multinomial random variable with $2C_2$ trials and cell probabilities of
$\alpha_i/\alpha_\bigdot$, $i=1,\ldots,K$.

\section{Other Interpretations of the Urn Model\label{sec:TwoStage}}
The urn model described above is a special case of more general classes of
genetic inheritance models.  Recognizing this provides us with a better idea of
why this model might be suitable from a biological perspective, and will help in the
analysis of fixation probabilities later in the chapter.

First, this urn model is a special case of a conditional branching
process model.  In such a model, each gene in a population independently produces a random
number $k$ of offspring, with $k$ following the same distribution for each gene.  The final
result, in the following generation, however, is made conditional upon the population size
at that time being  $C_2$ diploids.  By this conditioning, the numbers
of offspring of each gene are no longer independent, but they are still exchangeable.  This
type of model was introduced by \cite{Mor&Wat1958}, and studied in great detail by
\citeN{Kar&McG1965}.  The urn model corresponds to a conditional branching process model in
which offspring number has the negative binomial distribution, a versatile distribution
which has been previously employed to model the distribution of family sizes
\cite{Raoetal1973}.  

The negative binomial probability mass function may be parameterized in
terms of
$\alpha$ and $\beta$,  shape and scale parameters, respectively, analogous to
the parameters of a gamma distribution.  The pmf of a 
$\mathrm{NegBin}(\alpha,\beta)$ rv is
\[
P(X|\alpha,\beta) = \biggl(\frac{\Gamma(X+\alpha)}{\Gamma(X+1)\Gamma(\alpha)}
\biggr)
\biggl( \frac{1}{\beta+1}
\biggr) ^\alpha
\biggl( \frac{\beta}{\beta+1}
\biggr) ^X.
\]
Further, if $X_1$ and $X_2$ have negative binomial distributions with common scale
$\beta$ and shape parameters $\alpha_1$ and $\alpha_2$, then their sum has a
$\mathrm{NegBin}(\alpha_1+\alpha_2,\beta)$ distribution.  This, and the fact that the ratio
of a negative binomial variable with a sum of itself and other independent negative binomial
variables with common scale is a beta-binomial random variable,  makes it easy to verify
that the conditional branching process model with negative binomial offspring distribution
is the same as the urn model.


Another  general class of models, which I call ``two-stage" models
may be formulated as follows:
\begin{enumerate}
\item Let there be $C_1$ diploids in the current generation
\item Form a gamete (or juvenile) pool in which  the $i\thh$ gene copy
$(i=1,\ldots,2C_1)$ is represented by $H_i$ copies of itself, where $H_i$ is a
random number drawn from some distribution $P_H$, and where $P_H$ is the same
for all $i$.  
\item Sample $2C_2$ gametes with replacement from the gamete pool to form the
next generation.
\end{enumerate}
This model will behave differently
according to the type of distribution and the parameters specified for $P_H$. 

When $P_H$ is chosen to be a gamma distribution, $\mathrm{Gamma}(\alpha,\beta)$,
the two-stage model is identical to the urn model described in the previous
section. This is so because the ratio  of a gamma random variable to a sum of
itself and other independent gamma random variables with the same scale parameter is a beta
random variable. This relation generalizes to the multivariate case:
if
$Y_1,\ldots,Y_K$ follow independent gamma distributions, each with their own shape parameter
$\alpha_i$, but all with the same scale parameter, then
\[
\bm{Q} = \biggl(\frac{Y_1}{Y_1 + \cdots +Y_K}~,~
\ldots~,~
\frac{Y_K}{Y_1 + \cdots +Y_K}
\biggr)
\]
is a Dirichlet random vector with parameters $(\alpha_1,\ldots,\alpha_K)$. 
Knowing this, it is straightforward to verify that the urn model with parameter
$s$ is identical to the two-stage model with $P_H$ being a gamma distribution
with shape parameter $\alpha = 1/s$, and arbitrary scale.

The gamma distribution might seem a poor choice for the number of
gametes produced by a gene, since this will yield non-integer values for the
number of gametes. However, the mathematical tractability achieved by using a
continuous distribution for gamete number is substantial.  Further, the gamma
distribution is a continuous analogue which closely resembles the negative binomial
distribution.  Obviously, the continuous
approximation will be better when the scale is larger, as will be the case with
organisms that produce many offspring.  Nonetheless, even with fairly small
offspring number, the continuous approximation is good.

If we allow that the urn model is a good approximation to the two-stage model
with negative binomial gamete numbers, then we have yet another interpretation
of the model owing to the genesis of the negative binomial as a gamma-weighted
mixture of Poisson random variables.  That is, if
$\phi\sim\mathrm{Gamma}(\alpha,\beta)$ and $H\sim\mathrm{Poisson}(\phi)$, then,
marginally $H\sim\mathrm{NegBin}(\alpha,\beta)$.  Therefore, another biologically
reasonable interpretation of the urn model is that the number of offspring copies
of a gene  surviving  to childhood follows a Poisson distribution; however, there is
heterogeneity in the population, so that the expected number of offspring
surviving to childhood varies randomly across the pool of parental gene copies,
following a gamma distribution.  Children are then
chosen randomly (and with replacement) to survive to adulthood. So long as the number
of juveniles is large, sampling with replacement is a
reasonable approximation to the actual hypergeometric sampling that would
actually transpire.

\section{Comparison to the Wright-Fisher Model\label{sec:UrnCompare}}
For different values of $C_1$, $C_2$, and $s$, how does this model compare
to the Wright Fisher model?  In this section we investigate some fundamental
quantities in the urn model: variance in offspring number, probability of
identity by descent, and allele frequency variance.  From these calculations, it
is possible to derive the inbreeding and variance effective sizes of ideal populations
reproducing via the urn model.  I also investigate the probability of allele fixation in
these models, and finally describe how one may relate the stochastic replacement
quantity $s$ in a population of changing census size over time, to the quantity
$\lambda$.  

\subsection{Offspring number}
The expected number of offspring copies of a parental gene may be computed
simply from (\ref{eq:ExpBB}).  We consider an urn model in which one ball
(gene copy) is white and the other $2C_1-1$ balls are black.  The expected
number of offspring genes of a single gene copy is then the expected
number of white balls obtained when $2C_2$ total balls are drawn in our
urn sampling scheme.  In this case $n_i = 1$, so $\alpha_i = 1/s$ and
$\alpha_\bigdot = 2C_1/s$.  So the expected number of offspring
genes is
\begin{equation}
\Exp(\mathrm{\#~of~offspring~of~a~single~gene}) = 2C_2 \cdot
\frac{1/s}{2C_1/s} =
\frac{C_2}{C_1}.
\label{eq:OffExp}
\end{equation}
Note that this is the same expectation that one would
get for Wright-Fisher sampling (s=0) in a population of changing size.

The variance in offspring number is obtained by substituting the proper
expressions into (\ref{eq:VarBB}), giving
\begin{eqnarray}
\Var(\mathrm{\#~of~offspring~of~a~single~gene}) & = &
\frac{2C_2\biggl(2C_2 + \frac{2C_1}{s}\biggr)
\biggl(\frac{1}{s}\biggr)
\biggl(\frac{2C_1-1}{s}\biggr)}
{\biggl(1+ \frac{2C_1}{s}\biggr)
\biggl(\frac{2C_1}{s}\biggr)^2}.  \nonumber \\
& = & \frac{C_2}{C_1}\cdot
\frac{2C_1 + 2C_2 s - 1 - sC_2/C_1}{2C_1 + s}
\label{eq:OffVar}
\end{eqnarray}
It is instructive to notice that in the case of $C_1=C_2=C$ this reduces
to
\[
\mathrm{Variance~of~offspring~number}
= \biggl(1-\frac{1}{2C}\biggr)\biggl(\frac{s+1}{1+\frac{s}{2C}}\biggr)
\]
which is precisely the variance in offspring number in a Wright-Fisher
population, inflated by the factor $2C(s+1)/(2C+s)$.  Thus, if
$s=0$, we obtain the Wright-Fisher variance in offspring
number, as we ought to.

\subsection{Identity by descent}
We may similarly determine the probability that a randomly chosen pair of
gene copies drawn from the $2C_2$ gene copies of generation~2 are
identical by descent (IBD), \ie they are both copies of a single gene copy in generation~1. 
Let the random variable $X_j$ denote the number of offspring of a single, specific, gene
copy, say the $j\thh$ gene copy, in generation~1.  Conditional on
$X_j=x_j$, the probability that a pair drawn from generation~2 is IBD
{\em and} that both members of the pair are copies of the $j\thh$
gene in generation~1 is
$(\frac{x_j}{2C_2})(\frac{x_j-1}{2C_2-1})$.  Therefore, the probability
that a pair is IBD and are copies of the $j\thh$ gene may be written as
\[
P(\mathrm{IBD}~\mathit{and}~\mathrm{copies~of~gene}~j) = 
\sum_{x_j=1}^{2C_2}
\biggl(\frac{x_j}{2C_2}\biggr)\biggl(\frac{x_j-1}{2C_2-1}\biggr)
P(X_j=x_j),
\]
where $P(X_j=x_j)$ is, as given in (\ref{eq:betabomb}), the beta-binomial
probability of gene offspring number.  This may be rewritten as 
\begin{eqnarray}
\mbox{} & = &
\biggl(\frac{1}{2C_2(2C_2-1)}\biggr)\sum_{x_j=1}^{2C_2}(x_j^2 -
x_j)P(X_j=x_j)  \nonumber \\
\mbox{} & = & \biggl(\frac{1}{2C_2(2C_2-1)}\biggr)\biggl(
\Var(X_j) + (\Exp X_j)^2 - \Exp X_j \biggr)
\end{eqnarray}
which, upon substituting expressions (\ref{eq:OffExp}) and
(\ref{eq:OffVar}) for the mean and the variance, and factoring out a
factor of $C_2/C_1$ simplifies to 
\begin{eqnarray}
& = & \frac{1}{2C_1(2C_2-1)}
\biggl( \frac{2C_1 + 2C_2 s - 1 - s\frac{C_2}{C_1}}{2C_1 + s}
+ \frac{C_2}{C_1} - 1 \biggr)  \nonumber \\
& = & \frac{1}{2C_1(2C_2-1)}
\biggl( \frac{(2C_2 -\frac{C_2}{C_1} - 1)s - 1 }{2C_1 + s} +
\frac{C_2}{C_1}
\biggr).
\label{eq:IBDj}
\end{eqnarray}

Equation~\ref{eq:IBDj} gives the probability that a randomly drawn pair
is IBD {\em and} both  are copies of a particular gene $j$.  Since the
events of being copies of a particular gene $j$ and copies of a gene $k$
are disjoint for $j \neq k$, the probability of IBD may be found by
summing  (\ref{eq:IBDj}) over $j$ from 1 to $2C_1$, and so
\begin{equation}
P(\mathrm{IBD}) = 
\frac{1}{2C_2-1}
\biggl( \frac{(2C_2 -\frac{C_2}{C_1} - 1)s - 1 }{2C_1 + s} +
\frac{C_2}{C_1}
\biggr).
\label{eq:IBD}
\end{equation}
It is once again instructive to consider the case $s = 0$, in which the
above expression simplifies to $1/(2C_1)$.  This does not depend at all on
$C_2$, which is true with Wright-Fisher sampling.  With respect to inbreeding in the
Wright-Fisher model, it does not matter how much a population has grown
over the past generation; what counts is how small the population was in
the parental generation.  The same is not true in the urn
model with
$s\neq 0$.  This results from the fact that the relative effect of $s$
depends on how many gene copies are initially in the urn {\em and} how many are
drawn from it.   

If the urn-model population were of constant size $C$ ($C_1=C_2=C$), then
(\ref{eq:IBD}) reduces to 
\[
P(\mathrm{IBD}) = \frac{s+1}{2C+s}.
\]
Thus, for an urn-model population of constant size $C$ to have an inbreeding
effective size of $N_e$, $s$ would have to be chosen so that
\begin{equation}
s = \frac{2C - 2N_e}{2N_e - 1}.
\label{eq:sForInbNe}
\end{equation}


\subsection{Allele frequency variance}
We can consider an allele of type ``$i$" which is present in $n_i$ copies
out of the
$2C_1$ gene copies of generation~1.  The number of copies of allelic type
$i$ in the $2C_2$ gene copies of the following generation is then the
random variable
$X_i$.  Substituting the appropriate quantities for $\alpha_\bigdot$ and
$\alpha_i$ into (\ref{eq:VarBB}) gives
\begin{eqnarray}
\Var(X_i) & = & \frac{ 2C_2 
\biggl(2C_2+\frac{2C_1}{s}\biggr)
\biggl(\frac{n_i}{s}\biggr)
\biggl(\frac{2C_1-n_i}{s}\biggr)
}{\biggl(\frac{2C_1}{s}\biggr)^2
\biggl(1+\frac{2C_1}{s}\biggr)}  \\
& = & \frac{2C_2}{4C_1^2} \cdot n_i(2C_1-n_i) 
\frac{2C_1/s + 2C_2}{2C_1/s +1}. 
\label{eq:CountVar}
\end{eqnarray}
Hence the corresponding allele {\em frequency}, $X_i/(2C_2)$ will have
variance which is
$1/(4C_2^2)$ of the expression in (\ref{eq:CountVar}).  Denoting $n_i/(2C_1)$ by $p$, we may
write
\begin{equation}
\Var\biggl(\frac{X_i}{2C_2}\biggr) =
\frac{p(1-p)}{2C_2}\biggl(\frac{2C_1 + 2C_2s}{2C_1 + s}\biggr).
\label{eq:FreqVar}
\end{equation}
Once again, with $s=0$ we have the binomial variance of the Wright-Fisher
model.  Also, it is clear from (\ref{eq:FreqVar}) that for $s>0$ the allele
frequency variance is inflated over that of a Wright-Fisher population of
size $C_2$.  The magnitude of the effect of $s$ depends on the relative
and absolute sizes of $C_1$ and $C_2$.  With $C_1=C_2=C$
(\ref{eq:FreqVar}) reduces to 
\[
\Var\biggl(\frac{X_i}{2C}\biggr) = p(1-p)
\frac{1+s}{2C+s}.
\]
This means that for an urn-model population of constant size $C$ to have a
variance effective size of $N_e$, $s$ would be 
\begin{equation}
s = \frac{2C - 2N_e}{2N_e - 1},
\label{eq:sForVarNe}
\end{equation}
the same as (\ref{eq:sForInbNe}).

\subsection{Probability of allele fixation in one generation\label{sec:proboffix}}
So far, we have seen that both the inbreeding and variance effective sizes
of an urn-model population of constant size $C$ are $N_e$ when $s$ is
chosen to be $(2C-2N_e)/(2N_e-1)$. 
Given this, it is tempting to imagine that the eigenvalue effective size
\cite{Ewens1979}  in
an urn model population of size
$C$ will be the same as that in a Wright-Fisher population of size $N_e$
when
$s=(2C-2N_e)/(2N_e-1)$. 
The eigenvalue  effective size is the size of a Wright-Fisher
population with the same largest non-unit eigenvalue of its transition
probability matrix, and it determines the rate of loss of rare alleles over time
in the population. I have not
yet computed the eigenvalue effective size of this urn model. Such an analysis should not be
difficult, however, as this urn model is an exchangeable model in the sense of
\citeN[,~p.~77]{Ewens1979}, and so the general theory of
\citeN{Cannings1974} for computing eigenvalues as the expected value of products
of offspring numbers should hold.  
I have not yet pursued this, as I am primarily interested in
the application of this urn model to inference problems.  Nonetheless, future
work along this line might permit some analytical results of the model to be 
obtained directly from the earlier work of
\citeN{Kar&McG1965}.



However, I have performed some related numerical investigations.  Rather than determining
the eigenvalue effective size of the population, directly, I have  calculated, for the
Wright-Fisher model and the urn model, the probability that an allele is lost or becomes
fixed {\em in a single generation}.  This quantity will be related to the eigenvalue
effective size, and is simpler to compute.
These
numerical investigations (using the probabilities of the zero class
computed from Equation~\ref{eq:betabomb} and from the binomial
distribution) show that the probability of allele fixation in a single generation of
reproduction in an urn model population of census size $C$ and variance effective size
$N_e<C$ is always smaller than the probability of allele fixation in a single generation in
Wright-Fisher population of size $N_e$.  An example of this is given in
Figure~\ref{fig:UrnFixCompare}  which  shows the 
probability of allele loss for different values of the variance effective size.  Dots are
plotted for the probability that an allele found in five copies in an urn-model population
of size 100 diploids is lost in the next generation. 
\begin{figure}
\begin{center}
\includegraphics*[width=.7\textwidth]{UrnFixCompare.eps}
\end{center}
\caption[Comparison of one-generation fixation probabilities]{Probability that an
allele at frequency .025 is lost from a population in a single generation of reproduction. 
Crosses show the single-generation fixation probabilitities for a Wright-Fisher
population of  size given on the horizontal axis. Filled circles show
the single-generation fixation probabilities for an urn model population with census size
$C=100$ and variance
 effective size given on the horizontal axis.   As is apparent, rare alleles are more
likely to be lost in a single generation from a Wright-Fisher population of size $N_e$ than
from an urn model population of variance effective size $N_e$.  }
\label{fig:UrnFixCompare}
\end{figure}
Plotted on the same graph are crosses
showing the probability that an allele at frequency $5/200 = .025$ is lost in one
generation of Wright-Fisher sampling in a population of the same variance
effective sizes.
It is apparent from the figure that the single-generation fixation probabilities are much
higher in the Wright-Fisher model.  For example, with $\lambda=.3$ ($N_e=30$ in the
graph), it is three times more probable that zero copies of an allele will
appear in the following generation in the Wright-Fisher population than in the
urn-model population.

It might seem unfortunate that the urn-model population does not correspond
closely to the Wright-Fisher population in terms of single-generation allele fixation
probabilities.  I argue, however, that this results from a deficiency of the
Wright-Fisher model, rather than an ``inaccuracy" in the urn model I have
proposed.  To some extent the discrepancy must be due to the fact that in a
Wright-Fisher population of size
$N_e$, there are only
$2N_e+1$ possible values for the number of copies of any one allele, as opposed
to $2C+1$ possible values in an urn model population.  The implication of this
is that in an urn-model population, alleles may persist at lower frequencies than
would be allowed in a Wright-Fisher population.  Since almost all natural
populations have census sizes which are larger than their effective sizes,
alleles in those natural populations may survive at lower frequencies than would
be allowed in a Wright-Fisher population of comparable variance or inbreeding
effective size.  In other words, modeling a population
by a Wright-Fisher model of comparable variance effective size
will, in almost all cases, lead one to overestimate the probability that an allele is lost
from the population.  I devote Section~\ref{sec:Fixation} to this topic.
 

\subsection{Comparable $\lambda$}
If census sizes are known, $C_0,\ldots,C_T$, and we want
to find an urn model that is comparable to one in which the reproduction is like
a Wright-Fisher population with effective sizes $\floor{\lambda
C_0},\ldots,\floor{\lambda C_T}$, then we can define stochastic
replacements
$s_t$ for each pool of $C_t$ adults, $t=1,\ldots,T$.  In this, we shall make
our main concern that of matching the increase in allele frequency
variance.  In a generation of Wright-Fisher sampling from a gamete pool
with allele frequency $p$ and drawing $\floor{2\lambda C_{t}}$ gene
copies, the variance will be $p(1-p)/\floor{2\lambda C_t}$.  Setting this
equal to the allele frequency variance derived for the urn model
(\ref{eq:FreqVar}), and changing $C_1$ and $C_2$ to $C_{t-1}$ and $C_t$,
respectively, gives
\[
\frac{p(1-p)}{\floor{2\lambda C_t}} = 
\frac{p(1-p)}{2C_{t}}\biggl(\frac{2C_{t-1} + 2C_{t}s}{2C_{t-1} + s}\biggr).
\]
Then, disregarding the floor function, which resulted from a discretization
imposed by adherence to the Wright-Fisher model, this may be solved for
$s$, giving
\begin{equation}
s_t = \frac{2C_{t-1}(1-\lambda)}{2\lambda C_{t}-1}~~~~~,~~~~~t =
1,\ldots,T.
\label{eq:sfromlambda}
\end{equation}
This will form the basis for making inference about $\lambda$ using the urn
model.


\section{Allele Fixation in Population-Genetic Models\label{sec:Fixation}}
In \ref{sec:proboffix} we saw that an urn model with a given variance effective
size has a much smaller probability of allele loss in one generation than the corresponding
Wright-Fisher model.   The purpose of this section is to demonstrate that the
probability of fixation in a single generation of reproduction in a Wright-Fisher model of
size
$N_e$ is probably too high for a natural population of variance effective size $N_e$.  This
is done indirectly, by showing that the two-stage model (Section~\ref{sec:TwoStage}),
even under an extreme assumption about $P_H$, exhibits single-generation fixation
probabilities for a given variance effective size that are lower than those in the
Wright-Fisher model.

The argument proceeds as follows: let $P_H$ belong to $\mathcal{C}$, the class of
distributions with compact support on non-negative values.  This
corresponds to individuals having some maximum number $O_\mathrm{max}$ of
offspring that may survive to the gamete stage in the two-stage model.  The urn
model can approximately fit into this scheme if we define the scale parameter of
the corresponding gamma distribution to be small enough that the probability that
offspring number is greater than
$O_\mathrm{max}$ is very small.  Within this class, one may have discrete or
continuous, unimodal or multimodal distributions, \etc\ \ Through computer
simulations, I have found that, for a given variance effective size, the
$P_H\in\mathcal{C}$  that gives rise to the largest probability of allele
fixation in a single generation of reproduction seems to be the scaled Bernoulli
distribution---with probability
$p$ an individual has $O_\mathrm{max}$ offspring, and with probability $1-p$ it has
zero offspring.  (Of course, the value of $O_\mathrm{max}$ is irrelevant in this
case and may be set to unity without loss of generality.)  I have not 
yet tried seriously to prove that $P_H\sim\mathrm{Bernoulli}$ gives the maximum
single-generation fixation probability for a two-stage model of given variance effective
size, but a proof should be possible.  The probability of fixation in a Wright-Fisher model
of size
$N_e$ however, is still always greater than that in the two-stage model with
$P_H\sim\mathrm{Bernoulli}(p)$ for given variance effective size (which depends
on $p$).  An example of this appears in Figure~\ref{fig:Bernoulli}, which is
similar to Figure~\ref{fig:UrnFixCompare}, except that it includes a series of
dots for the two-stage model with Bernoulli offspring distribution.  
\begin{figure}
\begin{center}
\includegraphics*[width=.7\textwidth]{TwoStageBernoulli.eps}
\end{center}
\caption[Fixation probabilites in a two-stage model with
$P_H\sim\mathrm{Bernoulli}$]{Probabilities that an allele at frequency .025 is lost in one
generation of reproduction as a function of variance effective size under three different
models of genetic inheritance.  The dashed line is for a Wright-Fisher model of size $N_e$. 
The solid line is for an urn model of size 100, with $s$ set to obtain the corresponding
variance effective size.  The dots are simulation results using a two-stage model in a
population of size 100 adults with
$P_H\sim\mathrm{Bernoulli}(p)$.  The parameter $p$ increases from $.04$ to 1.0, from left to
right in the figure---higher values of
$p$ give rise to higher variance effective sizes.}
\label{fig:Bernoulli}
\end{figure}


Figure~\ref{fig:UrnFixCompare} shows us that the probability of allele loss in a single
generation in a Wright-Fisher model {\em exceeds} that of a two stage model with Bernoulli
offspring number distribution for all variance effective sizes.  This suggests that  using
a  Wright-Fisher model of size $N_{e_v}$ to calculate the probability of allele extinction
in one generation in a natural population of census size $C$  and estimated variance
effective size
$N_{e_v}$ will likely overestimate that probability.  This is apparent, because the
single-generaton fixation probability calculated from a Wright-Fisher model is greater than
that for any two-stage model of census size $C$.
The degree to which the single-generation probability of allele fixation computed from  a
Wright-Fisher model of size $N_{e_v}$ overestimates the true fixation probability depends on
the actual distribution of offspring numbers in the natural population.  If the natural
population is one in which a few individuals have many offspring, and the rest have zero,
then the Wright-Fisher model will not be grossly inaccurate. If, however, the population
is one in which the distribution of offspring number is not bimodal (for
example if the distribution of offspring number is shaped like the negative binomial
distribution), then the Wright-Fisher model of size $N_{e_v}$ will poorly represent that
natural population in terms of the rate of loss of rare alleles.  In such a case,
the urn model will provide a much more faithful representation.
While the former scenario may be appropriate
for organisms in which families all live or die together on the basis of some environmental
factors, it seems unlikely that it would apply well to many species.  

Finally, these
concerns are important from an inference perspective.  In making likelihood inference about
effective size, the information in fixation probability is used (since one makes full use of
the entire distribution of allele counts---not just their variance).  As a result, for most
natural populations, using a likelihood derived from the Wright-Fisher model will yield an
estimate of effective size which is larger than the estimate that would be obtained under a
likelihood model based on the urn-sampling scheme using an estimate of the census size of
the population (as described in the following section).  Such an effect, however, is quite
small, as we see below.

\section{MCMC for Bayesian Estimation Under this Urn Model }
The urn model allows easier implementation of the sampler described in
Section~\ref{sec:ABayesianApproach} for Bayesian inference of $\lambda$.  The formulation
still follows directly from Sections~\ref{sec:MCMCLikeLambda}
and~\ref{sec:ABayesianApproach} as presented, but with the definition of
$P_\lambda(\btX_t,\btX_{t-1})$ made in terms of transition probabilities dictated by the urn
model with a stochastic replacement corresponding to a particular value of $\lambda$ and the
observed census sizes.  Thus, for a particular value of $\lambda$ and the census sizes
$C_{t-1}$ and $C_t$ at times $t-1$ and
$t$,
$P_\lambda(\btX_{t}|\btX_{t-1})$ is written in terms of $s_t$, following 
\ref{eq:CMD}:
\begin{eqnarray}
P_\lambda(\btX_{t}|\btX_{t-1}) &=& P(\btX_{t}|\btX_{t-1},s_t,C_t,C_{t-1})\nonumber \\
& = &
P(\btX_t|2C_t;\alpha_1,\ldots,\alpha_K) \nonumber \\
& = & 
\frac{(2C_t)!\;\Gamma(\alpha_\bigdot)}{\Gamma(2C_t +
\alpha_\bigdot)}
\prod_{i=1}^K \biggl( \frac{\Gamma(X_{t,i} + \alpha_i)}{X_{t,i}!\;
\Gamma(\alpha_i)}
\biggr)
\end{eqnarray}
where $\alpha_i = X_{t-1,i}/s_t$, $\alpha_\bigdot = \sum_{i=1}^K\alpha_i$, and $s_t$ is
computed from Equation~\ref{eq:sfromlambda}.  I will briefly describe a simple
implementation of an MCMC sampler here.  This is a special case of the sampler described in
the following chapter where a more comprehensive treatment is given.  

I implemented such a sampler to compute the posterior probability for values of $\lambda \in
\Lambda =
\{\lambda_{\mathrm{min}},\lambda_1,\ldots,\lambda_n,\lambda_\mathrm{max}\}$ using a simple
random walk proposal distribution for
$\lambda$ (\eg $h(\lambda^*|\lambda) = 1/2$ for $\lambda^* = \lambda_{i-1}$ or
$\lambda^* = \lambda_{i+1}$, $\lambda \not\in\{\lambda_\mathrm{min},\lambda_\mathrm{max}\}$;
$h(\lambda^* = \lambda_1|\lambda_\mathrm{min}) = 1$; and $h(\lambda^* =
\lambda_n|\lambda_\mathrm{max})=1$), and a discrete uniform proposal distribution
$q(w|X_{t,k},X_{t,\ell})$ between the integers $-L$ and $L$, inclusive, where 
$L = 3 + \sqrt{\min\{X_{t,k},X_{t,\ell}\}}$.  I proposed $E$ updates to random components of
$\bX$ for each update proposed to $\lambda$, where $E$ was chosen so that each component of
$\bX$ was expected to receive one proposal for each proposal to change $\lambda$.
The proportion of time during the MCMC run in which $\lambda= \lambda_i$ estimates the
posterior probability that $\lambda=\lambda_i$.\footnote{This uses the Monte Carlo estimator
for probabilities  (\ref{eq:ProbMC}).  In the following chapter I present a more
sophisticated method for updating $\lambda$ that will also permit a Rao-Blackwellized
estimator (\ie Equation~\ref{eq:RB}) for the posterior distribution of $\lambda$.}

I used this sampler to compute a posterior
distribution for $\lambda$ given the data of \citeN{Begonetal1980}.  They reported census
sizes of tens of thousands of flies.  I thus set $C_t = 10,000$ and computed a posterior for
$\lambda\in\Lambda$ assuming a uniform prior on the values in $\Lambda$.  The points in
$\Lambda$ were chosen to correspond to the values of $N_e$ for which
log-likelihood values were computed in Chapter~\ref{chap:ForBack}.  Starting values for the
$\bX$ were obtained using a realization from the forward-backward sampler of
Chapter~\ref{chap:ForBack}.  This made burn-in essentially unnecessary. The sampler was run
for 100,000 updates of
$\lambda$, with the collection interval $u=1$.  This took 12 hours on a laptop computer with
a 266 Mhz G3 (Macintosh) processor.    Since this is a case with a one-dimensional parameter
with uniform prior, I was able to convert the posterior distribution for $\lambda$ into a
log-likelihood curve for
$\lambda$ and hence for
$N_e$.  Thus, I was able to compare it to the log-likelihood curve computed in
Chapter~\ref{chap:ForBack}. The result appears in Figure~\ref{fig:BegonBayes}.

\begin{figure}
\begin{center}
\includegraphics*[width=.8\textwidth]{BayesChainCompare.eps}
\end{center}
\caption[Log-likelihood curve for $N_e$ using the urn model]{Log-likelihood (scaled so that
the maximum is zero) for
$N_e$ given the data of \citeN{Begonetal1980}.  Open circles show the result 
from the importance sampling scheme of Chapter~\ref{chap:ForBack}.  Filled circles show the
log-likelihood curve computed by an MCMC scheme using the urn model of the present
chapter and an assumed census size of 10,000 flies each generation.  The two methods give
comparable results, but the MCMC scheme requires less computer time. The curve for the MCMC
scheme is shifted slightly to the left relative to the curve from importance sampling.  This
may be due to the differences in fixation probabilities between the Wright-Fisher model
used in the importance sampling scheme, and the urn model used for the MCMC approach.}
\label{fig:BegonBayes}
\end{figure}

The two curves are nearly identical.  However, under the urn model the log-likelihood is not
as low for small values of $N_e$ as for the importance sampling method under the
Wright-Fisher model.  This slight difference may
arise from the different probabilities that the two models assign to the event of allele
fixation as described in Section~\ref{sec:Fixation}.

The curve obtained by MCMC required less computational time than that obtained
by importance sampling. It was also significantly simpler to implement the MCMC scheme than
the importance sampling scheme of Chapter~\ref{chap:ForBack}.  The MCMC scheme and sampling
models here are as simple as possible. In the following chapter they  will both  be extended
to handle a wider range of scenarios.  



\section{Discussion}

In this chapter, I have described and characterized an alternative to the Wright-Fisher
model for genetic inheritance.  This alternative, which I call the ``urn model" carries
particular benefits for estimating $\lambda$, the per-generation ratio of effective breeders
to the census number of breeders, when the census number is known.  The urn model allows the
formulation of a sensible likelihood for $\lambda$ and it makes MCMC calculation of that
likelihood possible.  The likelihood was derived under three assumptions that merit
particular attention here.  The first is the assumption that $\lambda$ is constant over
time.  The remaining two concern questions of the census size: ``What ages or
stages of individuals should be counted in the census size of the population?" and ``What if
census sizes are not known without error, but are, rather, themselves estimates with some
uncertainty?"


\subsection{Constancy of $\lambda$}
In the likelihood for $\lambda$ developed in this chapter, the assumption is made that
$\lambda$ is constant over time.  In practice, this assumption will likely be violated.  In
natural populations, one would expect that  $\lambda$
could be influenced by time-varying factors like population size or climatic conditions.  In
managed or controlled populations it would  be natural for $\lambda$ to change
over time due to changes in harvest regimes or breeding practices or pesticide use.  The
methods I have developed could be extended to allow the value of $\lambda$ to vary over the
different inter-sample intervals.  For example, $\lambda_{[0,t_1]}$ could represent the
value of
$\lambda$  that applied to the generations between the sample at time 0 and the sample at
time $t_1$ when the next sample is drawn from the population, and $\lambda_{[t_1,t_2]}$
would apply to the  interval between the sample at time $t_1$ and the next sample in time,
and so forth.  Then, these specific $\lambda_{[\cdot,\cdot]}$'s could all be estimated as
separate parameters.  Some precision will be
lost because, though all the $\lambda_{[\cdot,\cdot]}$'s would be estimated jointly, only a
fraction of all the data will apply directly to each  $\lambda_{[\cdot,\cdot]}$.
Since $\lambda$ may also be affected by the census size of the population (for example, if 
there are limited nesting sites, then a generation descended from a large group of parents
might be expected to have a smaller $\lambda$ than a generation descended from only a few
parents who did not have to compete for nest sites), one could also formulate a model in
which $\lambda$ in any generation, $t$, is a function of $C_t$ or $C_{t-1}$.  In the
Bayesian setting it would be possible to design a reversible jump MCMC \cite{Green1995}
sampler for comparing different models which correspond to different assumptions about how
$\lambda$ varies over time or as a function of census size. 

It is also not unlikely that $\lambda$ would vary  each generation {\em within} an
inter-sample interval; however, it would be folly to try to model this with a separate
$\lambda$ for each generation because, within a single inter-sample interval, separate
$\lambda$'s for each generation would not be practically identifiable.   For the
case of
$\lambda$ varying from generation to generation within the inter-sample interval, it does
not seem straightforward to derive an approximate expression for what the maximum likelihood
estimator for $\lambda_{[\cdot,\cdot]}$ estimates.  It is not, for example, the harmonic
mean of the $\lambda$ values each generation, except in the case of constant census size.

When $\lambda_{[\cdot,\cdot]}$'s for all the inter-sample intervals are constrained in the
model to be equal to an overall $\lambda$ as I have presented the model earlier in the
chapter, another question arises: ``How is the information from different intervals weighted
and combined to arrive at  a single overall estimate of $\lambda$?" The likelihood-based
or Bayesian  estimate of
$\lambda$  will be influenced most by the intervals in which census sizes are
small.  This occurs because the amount of information about $\lambda$ in the data, relative
to the noise from the random process of drawing genetic samples, increases when the census
size descreases.  This observation does suggest that in future analyses estimating $\lambda$
in populations with census sizes that fluctuate greatly, it would be prudent to carry out an
alternative analysis under a model in which intervals containing very small census sizes had
a
$\lambda_{[\cdot,\cdot]}$ which was a separate parameter from the
$\lambda_{[\cdot,\cdot]}$'s in intervals with much larger census sizes.

The observation also exposes a generic difficulty in estimating $N_e$ from multiple samples
in time.  The first formal method described for doing so, that of \citeN{Pollak1983}, assumes
that $N_e$ is constant over the entire period from the first to the last sample.  This is
likely to be untrue, especially if the census size fluctuates greatly
over time.  The likelihood analyses of the \citeN{Begonetal1980} dataset presented in this
chapter and the preceding one are similar in that they assume a constant $N_e$ or $C_t$ from
$t=0$ to $T$.  In all of the above cases, asking the question of how to weight information
from different intervals shows us that esimating a single $N_e$ over a period of time which
includes more than two sampling episodes is not a particularly well-defined problem.  When
census sizes are available, then estimating $\lambda$, instead of $N_e$ provides an
alternative.  

\subsection{How do we define census sizes?}
In populations with discrete generations and in which the reproductively active adults may
be readily counted and sampled separately from the rest of the population, it is quite clear
what the census size $C_t$ should be---the number of actively reproducing adults at 
generation $t$.  For example, with pink salmon ({\em Oncorhynchus gorbuscha}), $C_t$ should
count the number of spawning adults at time $t$. In other species, in which the distinction
betwen mature adult and developing youngster is not so clear, it is also less clear what
quantity (\ie total number of animals, number of females, \etc) should be chosen to be
represented by
$C_t$.  Still, a number of sensible choices could be used, depending on the life-history
features of the species under study.  For example, in a population with discrete generations
in which only the individuals above a certain age are reproductively active, $C_t$ should
count the number of such individuals at time $t$---it should not be the whole population
size---adults and youngsters together.  Likewise, the genetic samples should be drawn from
the reproductively active individuals or their immediate offspring.  

Many organisms do not have discrete generations, of course.  In such cases, the choice of
which population census quantities to define as $C_t$ and which segments of the population
to sample should be made within the context of a stochastic model that is faithful to the
life history of the organisms under study.  Only within the context of such a life-history
based model, will the meaning of $\lambda$ be clear.  Chapter~\ref{chap:overlap}
describes the elaboration of the urn model of this chapter to the life history pattern of
Pacific salmon that mature at different ages.  

There is a somewhat lively literature and debate involving the estimation of the ratio
$N_e/N$, where $N_e$ in the numerator represents a ``long term" effective size of a
population over time and
$N$ in the  denominator represents some sort of ``long term" census population size
\cite{Nunney1995,Hus&Bar1995,Vucetichetal1997}.   The goal in this series of papers is quite
different from mine here. Under the perspective of the above authors, fluctuations in
population size and overlapping generations are lumped together with variance in offspring
number as inseparable factors which influence the effective size of a population.  In
contrast, the P\'{o}lya urn model here allows one to define the parameter $\lambda$ to
measure the degree to which genetic change is influenced by variation in offspring number
{\em separately} from fluctuations in population size (and, in the next chapter, separately
from the effect of overlapping generations in the population).  

It should be clear why being able to estimate $\lambda$ is advantageous: it is often
possible to directly observe, and therefore account for, fluctuating population size and
overlapping year classes.  It is harder to observe $\lambda$, and in fact, in organisms with
high juvenile mortality, it must be estimated using genetic data.  Once an estimate of
$\lambda$ has been made, however, it can be used to predict genetic change in a population
given patterns of fluctuating population size or overlapping generations observed in the
future.  This sort of ``predictive analysis", applied to census data collected in the
future, but using an estimated $\lambda$ from genetic and census data collected in the
past, is not available if one computes ``long-term" $N_e/N$ as other authors have pursued.  

\subsection{Census sizes estimated with error}
Thoughout this chapter (and the next) I assume that the census sizes $C_t$ are known without
error.  While census sizes of some organisms can be determined quite accurately, they are
seldom known with certainty. For some populations, in fact, census estimates may be very
imprecise.  It would be worthwhile to include this uncertainty in census sizes into the
estimation procedure for $\lambda$. I have not pursued that here, but leave it as an open
problem.  

An {\em ad hoc} approach to propagating the uncertainty in census size estimates to the
estimates of $\lambda$ would probably be worth investigating, since treating the problem
fully from the likelihood or Bayesian perspective would be very challenging,
computationally.  For example, if the true, unknown census
size values were modeled as latent variables, then  in an MCMC scheme, the Hasting's ratio
involved in proposing changes to those latent census sizes would depend not only on the
census size data, but also on the genetic data at all the loci.  Not only that, but changes
to the latent census size would also change the size of the space of the latent $\bX$
variables.  Designing an MCMC  sampler that mixed well in this context would, I believe, be
exceedingly difficult.  