Bayesian inference in Y-linked two-sex branching processes with mutations: ABC approach
Miguel González Rodrigo Martínez Cristina Gutiérrez
stat.CO q-bio.PE








González, M. Gutiérrez, C. Martínez, R.Department of Mathematics and ICCAEx,  University of Extremadura,
              Avda. Elvas s/n, 06006-Badajoz, SPAIN.

              Tel.: +34-924-289-300 (Ext. 57571)

              Fax: +34-924-272-911

              cgutierrez@unex.es          

              González, M.

              mvelasco@unex.es
 

              Martínez, R.

              rmartinez@unex.es





Bayesian inference in Y-linked two-sex branching processes with mutations: ABC approach
    Miguel González Cristina Gutiérrez Rodrigo Martínez
    Received: date / Accepted: date
=======================================================================================





A Y-linked two-sex branching process with mutations and blind
choice of males is a suitable model for analyzing the evolution of
the number of carriers of an allele and its mutations of a
Y-linked gene. Considering a two-sex monogamous population, in
this model each female chooses her partner from among the male
population without caring about his type (i.e., the allele he
carries).

In this work, we deal with the problem of estimating the main
parameters of such model developing the Bayesian inference in a
parametric framework. Firstly, we consider, as sample scheme, the
observation of the total number of females and males
up to some generation as well as the number of males of each genotype at last generation. Later, we introduce the information of the mutated males only in the last generation obtaining in this way a second sample scheme. For both samples, we apply the Approximate Bayesian Computation (ABC) methodology to approximate the posterior
distributions of the main parameters of this model. The accuracy
of the procedure based on these samples is illustrated and discussed
by way of simulated examples.








§ INTRODUCTION


In <cit.>, a stochastic model in the field of branching
processes was introduced with the aim of describing the evolution
of the number of carriers of a Y-linked gene and its mutations in a
two-sex monogamic population. This model allows to study the interesting and important problem of how mutations of Y-linked genes evolve in a population. In a general sense, we use the term mutation for any change in the genetic material
which gives rise to the transmission of a different trait. We consider a population where two types of alleles could coexist. We denote them as  and 𝗋. The -allele is considered a marker allele or an allele which transmits a trait of interest (not expressed in the phenotype of the male) and the 𝗋-allele is considered an allele which transmits any other trait different of that transmitted by . Moreover, we assume that
-allele could mutate transmitting a different trait of  and therefore, we also denote this mutated allele as 𝗋. That is, in our context, 𝗋-allele means all alleles which transmit a trait different of that transmitted by , stemming or not from mutations.
We also assume that backmutation is not allowed, i.e. the 𝗋-allele never can return to the -form. Therefore, there could exist a flow from  to 𝗋 but not vice versa.  Notice that, if in the population there would only be -alleles, it could appear later 𝗋-alleles which would stem from mutations.

This model, called Y-linked two-sex branching process (Y-BBP) with
mutations, considers a population formed by females and males who mate with blind choi­ce to produce offspring, i.e. each female chooses her partner from among the male
population without caring about his genotype (because the trait is not expressed in the phenotype of the male or it is not decisive at mating time). Applying the genetic inheritance rules, every couple
gives birth to females and males, with every male progeny
inheriting the genetic material corresponding to the Y-chromosome
from his father. But, during reproduction, there could occur a
mutation in the transmitted allele by a father with -allele, altering
the characteristic of the son with respect
to his progenitor. Hence, under these assumptions, a male with -allele could
give birth either a male offspring who is a clone of his genetic
material (the same allele) or a mutant with a new type of allele (𝗋).

As important example of such mutations, one could suppose that an
alteration in the allele might impair the individuals reproductive
capacity. In this way, the process could be applied to model
problems of fertility. In particular, it would allow one to study
the case of mutations which end in different levels of fertility
including total infertility (aspermia). A particular case of this
situation is presented in <cit.>, in which it is suggested
that a mutation in the USP9Y Y-chromosomal gene causes the absence
of sperm in semen. Another possibility is that the mutation may
represent the beginning of a new paternal lineage, as for example
the one that gave rise to the haplogroup I which is related to
risk of suffering coronary disease, see <cit.>.

The
aforementioned work  <cit.> should be consulted for
further background motivation and information about conditions
guaranteeing a positive probability of survival of the alleles in the population. Such conditions depend on several parameters of the model:
the reproduction mean of each genotype, the probability of being
female and the probability of mutation. Therefore, from a practical point of view, it is necessary to develop estimation procedures for these parameters.

The aim of this paper is to develop the Bayesian inferential
theory for a Y-BBP with mutations considering an enough informative and realistic
sample scheme (in the sense of the minimum amount of information that it is necessary to be observed in order to obtain accurate estimates). The branching process theory has
usually assumed that the entire family tree is needed to be
observed in order to make accurate inferences. However, to observe
such quantity of information is hard in practice. In this sense, the authors published in a previous work a study about the inference of the parameters of a Y-BBP model without mutations (see <cit.>),
based on a more realistic sampling scheme where the
total number of females and males
up to some generation as well as the number of males of each genotype
in the last generation is observed.
Carrying on with these ideas, in this paper and for the Y-BBP with mutations, we consider firstly the same sample.
However, in contrast with the model without mutation, the son's genotype is not determined directly from the father's one.
As consequence, this sample could determine the global behavior of alleles in the population, but might not provide enough information in order to make inference on the parameters of the model with mutations.
Therefore, to overcome this lack of information, some knowledge about the number of
mutated males in the last generation should be added. This will be considered as the second sampling scheme.

Moreover, in the Bayesian framework, a Markov cha­in Monte Carlo (MCMC) me­thod was used for the model without mutations in <cit.>, with very good results. However, although in general MCMC me­thod works well in many substantive problems, it can perform poorly when is applied to large data sets or complex models, as the model presented in this paper. In fact, the approximation to this problem using the MCMC methodology has provided poor results, failing to provide accurate posterior approximations in a reasonable computational time. Besides, at least in our context, such methodology often needs to make use
of the conjugate family theory representing a lack of generality.

Due to these limitations, we are interested in applying a
different statistical tool to solve this incomplete data problem,
the Approximate Bayesian Computation (ABC) methodology (see, for example, <cit.>, <cit.> or <cit.> for a recent survey). This
method is being developed during last decades as an alternative to
such more traditional MCMC methods. These likelihood-free
techniques are very well-sui­ted to models for which the likelihood
of the data are either mathematically or computationally
intractable but it is easy to simulate from them, so that they
look very appropriate, a priori, for studying the inference of the Y-BBP with mutations.

Besides this Introduction, the paper is organized in 8 sections as follows. In Section 2, it is described
in detail the Y-BBP with mutations as well as the asymptotic
behavior of the different types of alleles in the population. Section 3 is devoted to introduce the Tolerance Rejection-ABC Algorithm. We apply it in Section 4 to a simulated example based on the sample des­cribed in this Introduction. In Section 5 we set out a more informative sampling scheme.  We apply again the algorithm, but now with this new sample, in Section 6, developing a series of simulated examples which cover the different situations that can be observed in the sample. After that, in Section 7, we examine the robustness of the methodology, and in Section 8 we use the approximation of the posterior distributions of the parameters to infer the predictive posterior distribution of the size of future generations. Finally, in Section 9, we provide some concluding remarks.



§ DESCRIPTION OF THE MODEL


The genetic frame we model is given by a Y-linked gene
which presents two allelic forms, denoted as  and 𝗋, where  can
mutate giving rise to new (different) alleles, all denoted also as
𝗋. This allele represents the
transmission of any trait different from the characteristic
transmitted by the 𝖱-allele (stemming or not from its mutations).

Since the Y-chromosome is specific to males, we deal with a
two-sex population formed by females, by males which carry the
𝖱-allele (called 𝖱-males), and by 𝗋-males which carry the 𝗋-allele. It is assumed that each individual
mates with only one individual of the opposite sex if available
(perfect fidelity or monogamous mating), forming a couple.
Therefore, in the population one could find two types of couples,
denoted by 𝖱- and 𝗋-couples, depending on
whether its male is of type 𝖱 or of type 𝗋,
respectively.

According to the rules of genetic inheritance, and taking into
account the possibility of mutation, an 𝖱-couple can
give birth to females, 𝖱-males, and 𝗋-males,
whereas, given the assumption of no backmutation and that mutations of 𝗋-allele are also named as 𝗋, an
𝗋-couple gives birth to females and 𝗋-males.

Assuming non-overlapping generations and given the number of
𝖱- and 𝗋-couples in generation 𝗇,
denoted by  and , respectively, the number of females,
males, and couples of each genotype in the (𝗇+1)st
generation is determined by considering a two-stage structure,
reproduction and mating, similarly as it was described in <cit.> and <cit.> for others Y-BBP without mutations.

In the reproduction phase, couples of the 𝗇th
generation produce offspring independently of each other and
according to certain reproduction law which is the same for a
given genotype but may be different for different genotypes since
the mutation could affect the reproductive capacity. Moreover,
these reproduction laws are independent of the generation the
couples belong to. Mathematically, the number of females and males of each genotype
stemming from each type of couple is identified with the following
independent sequences of independent, identically distributed,
non-negative, and integer-valued random vectors:

    {(,,), 𝗂=1,2,...;𝗇=0,1,...}

and 
    {(,𝖬_𝗇𝗃^𝗋→ 𝗋),
    𝗃=1,2,...;𝗇=0,1,...}.

Here,  and  are,
respectively, the number of females stemming from the
𝗂th 𝖱-couple and the 𝗃th
𝗋-couple of generation 𝗇;  is the
number of males stemming from the 𝗂th
𝖱-couple of the 𝗇th generation which have
preserved the 𝖱-allele, and  is the
number of males stemming from the 𝗂th
𝖱-couple of the 𝗇th generation, whose
alleles have mutated and now are of type 𝗋; and
finally, 𝖬_𝗇𝗃^𝗋→ 𝗋 is the number of males
stemming from the 𝗃th 𝗋-couple of the
𝗇th generation, and which therefore carry also the
𝗋-allele.

We assume that the distributions of ++ and
+𝖬_𝗇𝗃^𝗋→ 𝗋 have finite means,
𝗆_𝖱 and 𝗆_𝗋, respectively, and variances.


Moreover, the conditional distribution of the vector
(,,) given ++=k is multinomial
with parameters (k, α,
 (1-α)(1-β), (1-α)β), for
k≥ 0, and 0<α<1, 0≤β <1 with
α representing the probability for an offspring to
be female and β the probability of mutation. Then,
in accordance with this multinomial scheme, the average numbers of
females, -males, and 𝗋-males generated by an
-couple are, respectively, α 𝗆_𝖱,
(1-α)(1-β) 𝗆_𝖱 and (1-α)β
𝗆_𝖱. Notice that, if β=0, then mutations do not happen a.s. so, if in the population both alleles coexist, 𝗋-allele stems from the 𝗋-couples in the initial generation and one has the Y-BBP without mutation studied in <cit.>. The case β=1 is not considered in this paper because in such case, from the first generation on, only the 𝗋-allele would survive in the population a.s. and then one has the classical bisexual branching process introduced by <cit.> describing the evolution of this allele.

With respect to the mutant-allele, the conditional distribution of
(,𝖬_𝗇𝗃^𝗋→ 𝗋) given
+𝖬_𝗇𝗃^𝗋→ 𝗋=l is also multinomial with
parameters (l,α, (1-α)), for
l≥ 0, and 0<α<1, with
α the same for both genotypes, i.e., the gene has
no influence on sex designation. Then, the average numbers of
females and 𝗋-males are, respectively, α
𝗆_𝗋 and (1-α)𝗆_𝗋.

At the end of the reproduction phase, one has the total number of
females, 𝖱-males, and 𝗋-males, denoted by
, , and , respectively, which together constitute
the (𝗇+1)th generation. Specifically, one obtains such
variables by means of the following expressions:

    =∑_𝗂=1^+∑_𝗃=1^,


    =∑_𝗂=1^ =+𝖬^𝗋→ 𝗋_𝗇+1,
 where

    =∑_𝗂=1^    𝖬^𝗋→
    𝗋_𝗇+1=∑_𝗃=1^𝖬_𝗇𝗃^𝗋→ 𝗋,

with the empty sum defined as 0, and  and 𝖬^𝗋→
𝗋_𝗇+1 denoting the total number of males with
𝗋-genotype in generation 𝗇+1 which stemming
from 𝖱- and 𝗋-couples, respectively.

Given the total numbers of females, 𝖱-males, and
𝗋-males in the (𝗇+1)st generation, the
number of couples of each type ( or 𝗋) in this generation is
determined in the mating phase as follows: perfect fidelity mating is assumed, hence
if the total number of females is greater than or equal to the
total number of males then every male finds a mate in the female
population resulting in 𝖹^𝖱_𝗇+1= couples of type
𝖱 and 𝖹^𝗋_𝗇+1= couples of type
𝗋. On the other hand, every female mate when  the total number of males
exceeds the total number of females. Moreover, since it is assumed that the genotype has no impact on the mating mechanism, females choose its mate in a blind way. Hence, the total number of 𝖱-couples in the (𝗇+1)th generation, 𝖹^𝖱_𝗇+1, follows a hypergeometric distribution with parameters ,
=+, and , while the total number of
𝗋-couples in this generation equals the number of
remaining females, i.e., 𝖹^𝗋_𝗇+1=-𝖹^𝖱_𝗇+1, whose distribution is also hypergeometric with parameters , , and .

The bivariate sequence {(,)}_n≥0, describing the
evolution of the number of couples of each type over
generations, is called Y-linked two-sex bran­ching process with
mutations and blind choice of males. It is
shown in <cit.> that the process above is a homogeneous
multitype Mar­kov chain and that each genotype shows the dual
behavior typical for bran­ching processes known as the
extinction-explosion dicho­tomy. However, the behavior of the 𝗋-allele dependents on the behavior of the -allele. In concrete, if the 𝖱-allele becomes extinct, the survival or not of the 𝗋-allele depends on its own reproductive capacity. Whereas, considering β>0, if the 𝖱-allele explodes, the 𝗋-allele also explodes due to the mutations, independently of the 𝗆_𝗋 value, so that the coexistence set is a.s. {𝖹^𝖱_𝗇→∞,𝖹^𝗋_𝗇→∞}={𝖹^𝖱_𝗇→∞}. Moreover, this set has a positive probability if min{α,(1-α)}(1-β)𝗆_𝖱>1 (see <cit.> for details).

In <cit.>, a simulation-based study was developed
to determine the behavior of the different types of alleles
in the population on the coexistence set. So, we established that the asymptotic behavior of the 𝗋-allele depends on the relation between 𝗆_𝗋, and (1-β)𝗆_𝖱. In particular, when 𝗆_𝗋≥ (1-β)𝗆_𝖱, the 𝗋-genotype is the dominant one in the sense that,
a.s. on {𝖹^𝖱_𝗇→∞}, the sequence
{/}_𝗇≥ 0, converges to infinity.
In the case 𝗆_𝗋<(1-β)𝗆_𝖱, there is no
dominant genotype because the previous sequence converges, a.s. on {𝖹^𝖱_𝗇→∞}, to a positive and finite value.

Specifically, when 𝗆_𝗋>(1-β)𝗆_𝖱, for 𝗇 large enough, it can be stated that

    /≃(𝗆_𝗋/(1-β)𝗆_𝖱)^n 𝖶,
 with 𝖶 a certain non-degenerate random variable.
When 𝗆_𝗋=(1-β)𝗆_𝖱, the sequence  {/}_𝗇≥ 0 also grows a.s. to infinity, however now it does so linearly, that is, for 𝗇 large enough, it is satisfied that

    /≃𝗇β/1-β + 𝖶^*,

where 𝖶^* is a non-degenerate random variable. Finally, for 𝗇 large enough,

    /≃β 𝗆_𝖱/(1-β)𝗆_𝖱-𝗆_𝗋

in the case 𝗆_𝗋<(1-β)𝗆_𝖱, that is,  {/}_𝗇≥ 0 converges a.s. to the constant β 𝗆_𝖱((1-β)𝗆_𝖱-𝗆_𝗋)^-1 which had been determined empirically.

Moreover, we have determined computationally the asymptotic ratio of the quotient between the total number of 𝗋-couples in consecutive generations,  𝖹^𝗋_𝗇+1/, and we have concluded that such ratio is, a.s. on {𝖹^𝖱_𝗇→∞},  min{α,(1-α)}max{𝗆_𝗋,(1-β)𝗆_𝖱}. Finally, it was proved in <cit.> that the asymptotic ratio of 𝖹^𝖱_𝗇+1/ is, a.s. on {𝖹^𝖱_𝗇→∞}, min{α,(1-α)}(1-β)𝗆_𝖱.

Based on these previous results and a deeper study of the simulations, it is easy to deduce the rates of growth of every type of couple in every case, on the set where both genotypes survives. The knowledge of such ratios is important for the development of the results of this paper. So, when 𝗆_𝗋>(1-β)𝗆_𝖱, the sequence {}_𝗇≥ 0 grows geometrically at a rate τ_𝖱=min{α,(1-α)}(1-β)𝗆_𝖱 while
{}_𝗇≥ 0 grows, also geometrically, at a rate τ_𝗋=min{α,(1-α)}𝗆_𝗋, i.e., each type of couple have a different rate of growth, being the 𝗋-allele the dominant one.

On the other hand, when 𝗆_𝗋<(1-β)𝗆_𝖱, {}_𝗇≥ 0 and {}_𝗇≥ 0 have the same rate of geometric growth given by τ_𝖱, and moreover it is verified that, as 𝗇 tends to infinity, the limit of /τ_𝖱^𝗇 is, a.s. on {𝖹^𝖱_𝗇→∞}, proportional to the limit of /τ_𝖱^𝗇 with proportionality constant β 𝗆_𝖱((1-β)𝗆_𝖱-𝗆_𝗋)^-1.

Finally, when 𝗆_𝗋=(1-β)𝗆_𝖱,  {}_𝗇≥ 0 grows at a geometric rate of τ_𝖱 while the sequence that normalizes {}_𝗇≥ 0 is {𝗇τ_𝖱^𝗇}_𝗇≥ 0. Moreover, as 𝗇 tends to infinity, the limit of /𝗇τ_𝖱^𝗇 is, a.s. on {𝖹^𝖱_𝗇→∞}, proportional to the limit of /τ_𝖱^𝗇 with proportionality constant β/(1-β).

As we indicated at the Introduction, our aim in this paper is to apply the ABC methodology to obtain accurate approximations to the posterior distributions of the parameters of the model, that is, of α, β, 𝗆_𝖱 and 𝗆_𝗋 and to verify that this methodology works adequately in all the possible situations given by the explained above relations between 𝗆_𝗋 and (1-β)𝗆_𝖱, always on the coexistence set. To do that, previously, we must select the sample we are going to observe. We are interested in finding a sufficiently informative sampling scheme observing the minimum amount of information that leads us to obtain good estimates. Related to this question, as we also indicated at the Introduction, the authors published (see <cit.>) a study about the estimation of the main parameters of a Y-BBP (without considering mutations) based on a sample where only the total number of females and males (without knowing the genotype of the males) up to some generation 𝖭 as well as
the different types of males only in the last generation 𝖭 were assumed to be observed. Following these ideas, initially we set out in this paper the Bayesian estimation of the parameters of the Y-BBP with mutations based on that same sample.



§ APPROXIMATE BAYESIAN COMPUTATION


Let ℱℳ_𝖭 denote the observed data until
generation 𝖭 which is assumed that has been generated from a
model with parameter vector
θ=(α,β,𝗆_𝖱,𝗆_𝗋). In particular

    ℱℳ_𝖭={𝖥𝖬_0,𝖥𝖬_1,...,𝖥𝖬_𝖭-1,𝖥𝖬𝖱𝗋_𝖭},
 where
𝖥𝖬_𝗇=(,), 𝗇=0,...,𝖭-1, is the vector given by the total
number of females and males in generation
𝗇 and 𝖥𝖬𝖱𝗋_𝖭=(𝖥_𝖭,𝖬^𝖱_𝖭,𝖬^𝗋_𝖭) is the vector given by the total
number of females and males of each genotype at last generation. Note that 𝖥𝖬_0 could be fixed -initial generation
at an experiment- or random -representing the first generation one
observes, non necessarily the initial fixed generation.
Henceforward, we shall focuss on the first interpretation. Moreover, we shall assume that 𝖥_𝖭>0, 𝖬^𝖱_𝖭>0 and 𝖬^𝗋_𝖭>0. Notice that this assumption implies that >0 and >0, for all 𝗇=1,...,𝖭-1 and also implies that both genotypes have coexisted at least in the last generation.

The aim of Bayesian approach is to derive the posterior
distribution of the parameter vector,
θ|ℱℳ_𝖭. ABC methodology offers good approximations to the posterior
distributions of parameters for models which ha­ve intractable
likelihoods but are easy to simulate.

The use of ABC ideas initially comes from the field of population
genetics (see  <cit.>,  <cit.> and <cit.>), although  these were quickly extended to a great variety
of scientific applications areas. The basic ideas are to simulate
a large number of data from a model depending on a parameter
vector that is drawn from a prior distribution and compare the
simulated data with the values from the observed sample. The aim
of the ABC methodology is to provide samples from a posterior-type
distribution (in the sense that it includes the sample
information) which is a good (enough) approximation of the
posterior distributions of the parameters of the model. Several
algorithms have been proposed in the literature to solve the
problem of how to choose this approximation, surveys on ABC
algorithms can be read in <cit.>, <cit.> and <cit.>.

These general ideas can be properly adapted to our model which is very easy to simulate given the parameter vector, some information about the initial
generation, as for example, the total number of females and males
of each type, and the family of probability distributions the reproduction laws belongs to. In our case, as we have a complete absence of knowledge on the reproduction laws of the model that has generated the observed data, we will assume, for simplicity, a parametric setting with Poisson distributions as reproduction laws. This distribution is frequently used as offspring distribution, see for example <cit.>, <cit.>, <cit.>, <cit.>, <cit.> or <cit.>. Another parametric reproduction law could also be considered without substantial changes in the estimates (see the sensitivity analysis showed in Section <ref>).

Moreover, in our case, it is not possible to calculate explicitly the
likelihood function, f(ℱℳ_𝖭|θ), because the complete branching structure cannot be derived due to the fact that the total number of males of each genotype, the total number of 𝗋-males stemming from 𝖱-couples and the total number of each type of couple are not observed in each generation.



 §.§ Description of the algorithm


In our particular case, the proposed algorithm is the Tolerance
Rejection-ABC Algorithm which is an adaptation of that proposed in
<cit.> which works as follows. For a
Y-BBP with mutations, assuming observed the sample in
(<ref>), it is easy to simulate for each specific
vector of parameters θ (sampled from a prior
distribution π(θ)) the entire
family tree up to the current 𝖭th generation
and to obtain the random vectors (, , , ,
𝖬^𝗋→ 𝗋_𝗇, , ), 𝗇=0,…,𝖭. Then,  using Equations
(<ref>) and (<ref>), can be obtained a simulated sample of
(,), 𝗇=0,…,𝖭-1 and (𝖥_𝖭,𝖬^𝖱_𝖭,𝖬^𝗋_𝖭),
renamed as

    ℱℳ_𝖭^={𝖥𝖬_0^,𝖥𝖬_1^
    ,...,𝖥𝖬_𝖭-1^,𝖥𝖬𝖱𝗋_𝖭^}.


Notice, ℱℳ_𝖭^ depends on  ℱℳ_𝖭 only through 𝖥𝖬_0.
Actually, 𝖥_0^=𝖥_0 and the vector (𝖬^𝖱_0^,𝖬^𝗋_0^)
is simulated from the uniform distribution, subject to the constraint 𝖬^𝖱_0^+𝖬^𝗋_0^=𝖬_0.
Moreover, we consider only paths simulated by the algorithm where both alleles have coexisted in the last generation, i.e., where 𝖥_𝖭^>0,  𝖬^𝖱_𝖭^>0 and 𝖬^𝗋_𝖭^>0, as it occurred in the observed sample ℱℳ_𝖭.


Now, for a
given ε>0, known as tolerance level, and a distance,
ρ(·,·), the algorithm compares (in terms of metric)
the simulated paths, ℱℳ_𝖭^, with the observed sample ℱℳ_𝖭.
This allows us to obtain an approximation of θ|ℱℳ_𝖭 by the distribution

    θ|ρ(
    ℱℳ_𝖭^,ℱℳ_𝖭)≤ε,
 using a small enough ε. In our case, we shall use a small enough quantile of the sample of the distances as it is usual in ABC studies (see, for example, <cit.>).


To quantify the distance between ℱℳ_𝖭^ and
ℱℳ_𝖭 we use

    ρ(ℱℳ_𝖭^,ℱℳ_𝖭)=      
          (∑_𝗇=1^𝖭(^/-
    /^)^2+∑_𝗇=1^𝖭-1(^/- /^)^2.
          +.(^/- /^)^2+(^/-
    /^)^2)^1/2


Notice that we have re-scaled each coordinate of the vectors since their magnitudes can be extremely different, depending on generation, sex and genotype (see <cit.> and <cit.>).

Then, the Tolerance Rejection-ABC Algorithm is formulated as,




Tolerance Rejection-ABC Algorithm




10cm
For i=1 to m do

repeat

generate (α^𝗌𝗂𝗆,γ^𝗌𝗂𝗆,ϕ^𝗌𝗂𝗆)∼ U(0,1)×U(0,1)× U(0,1)

generate β^𝗌𝗂𝗆=0 with probability γ^𝗌𝗂𝗆 and

    β^𝗌𝗂𝗆∼π(β) with probability 1-γ^𝗌𝗂𝗆


generate 𝗆_𝗋^𝗌𝗂𝗆=0 with probability ϕ^𝗌𝗂𝗆 and

    𝗆_𝗋^𝗌𝗂𝗆∼π(𝗆_𝗋) with probability 1-ϕ^𝗌𝗂𝗆

generate 𝗆_𝖱^𝗌𝗂𝗆∼π(𝗆_𝖱)

let
θ=(α^𝗌𝗂𝗆,β^𝗌𝗂𝗆,𝗆_𝖱^𝗌𝗂𝗆,𝗆_𝗋^𝗌𝗂𝗆)


simulate ℱℳ_𝖭^𝗌𝗂𝗆 from the likelihood
f(ℱℳ_𝖭|θ)

until
ρ(ℱℳ_𝖭^𝗌𝗂𝗆,ℱℳ_𝖭) ≤ε,

set θ^(i)=θ

end for 



Note that, we generate the parameter α^𝗌𝗂𝗆 from a uniform distribution on (0,1) and the parameter 𝗆_𝖱^𝗌𝗂𝗆 from a generic prior distribution π(𝗆_𝖱) on (0,∞). This is consistent with the fact that >0 and >0.
On the other hand, taking into account that in the model, 𝗋-allele can mean an allele different from , β^𝗌𝗂𝗆 could be null. Moreover,  𝗆_𝗋^𝗌𝗂𝗆 could also be null even being >0 (in this case β^𝗌𝗂𝗆>0, see (<ref>)). Therefore, we generate the parameters β^𝗌𝗂𝗆 and 𝗆_𝗋^𝗌𝗂𝗆 from prior distributions which are mixture of distributions: one degenerated at 0 (in order to consider the possibility that β^𝗌𝗂𝗆 and 𝗆_𝗋^𝗌𝗂𝗆 takes exactly the value 0)
and the other one π(β) on (0,1) and π(𝗆_𝗋) on (0,∞), with weights given by (γ^𝗌𝗂𝗆, 1-γ^𝗌𝗂𝗆) and (ϕ^𝗌𝗂𝗆,1-ϕ^𝗌𝗂𝗆), respectively.
Since we do not have information about the possible value of these parameters, we consider γ^𝗌𝗂𝗆 and ϕ^𝗌𝗂𝗆 following a uniform distribution on (0,1).



§ A SIMULATED EXAMPLE BASED ON THE OBSERVED SAMPLE ℱℳ_𝖭


Now, the previous algorithm is implemented using as observed data a sample which has been obtained by simulation. We analyze first the case where the relation between the parameters is 𝗆_𝗋≥(1-β)𝗆_𝖱 although it is worth to remind here that we are searching for a general method which works independently of the relation between the parameters.



 §.§ Case 𝗆_𝗋≥(1-β)𝗆_𝖱


Our objective is to approximate the posterior distribution
θ|ℱℳ_𝖭, where ℱℳ_𝖭 is
an observed sample which has been
simulated from a Y-BBP with mutation with parameter vector
θ=(α,β,𝗆_𝖱,𝗆_𝗋)=(0.46,0.005,3.2,4)
(notice that with those values the relation
𝗆_𝗋≥(1-β)𝗆_𝖱 is satisfied)
and initial vector (𝖥_0, 𝖬_0^𝖱,
𝖬_0^𝗋)= (10,5,5). For such a model with
this set of parameters and initial values, we proved in
<cit.> that there exists a positive probability of
survival of both genotypes.



We simulate 15 generations of this Y-BBP with mutations assuming
that reproduction laws of both genotypes follow the non-parametric offspring distributions with finite support given in Table <ref>, with means 𝗆_𝖱=3.2 and
𝗆_𝗋=4, respectively. The observed data can be seen in Table
<ref> and are denoted by ℱℳ_15.



We apply the Tolerance Rejection-ABC Algorithm generating the parameter vector
assuming independent non-informative prior distributions. In
particular, for α^𝗌𝗂𝗆 and β^𝗌𝗂𝗆 (when it is positive) uniform distributions in the interval (0,1) and for 𝗆_𝖱^𝗌𝗂𝗆 and 𝗆_𝗋^𝗌𝗂𝗆 (when it is positive), uniform distributions in the interval (0,10). We have chosen, obviously, 0 as the minimum value for the support of the latter uniform distributions and 10 as the maximum value because we consider that number high enough for the number of offspring of many animal species although this
number could be adapted to any specific situation. After that, we simulate Y-BBPs with mutations until generation 15 using, as  and 𝗋 reproduction laws, Poisson distributions with parameters, respectively, 𝗆_𝖱^𝗌𝗂𝗆 and 𝗆_𝗋^𝗌𝗂𝗆 (recall we use this generic type of distribution for the offspring laws because we know nothing about the true reproduction laws). We generate a pool of 50 millions of simulated paths. To compare the observed sample and the simulated ones we consider a tolerance level equal to the 0.00002 quantile of the sample
of the distances, so that the size of ABC samples
to approximate the posterior distribution θ|ℱℳ_15 is 1000.



In Figure <ref>, we present the approximate posterior
distribution of every parameter, that
is,

    δ|ρ(
    ℱℳ_15^,ℱℳ_15)≤ε,
 with
δ equal to α, β (in this case, only paths where β^𝗌𝗂𝗆>0 have been considered), 𝗆_𝖱 and
𝗆_𝗋, the corresponding in every case, together with the
true value of the parameters (vertical solid line) and 95%
HPD sets (vertical dotted lines). We can appreciate first that the approximate posterior distribution for α is very accurate. Actually, this happens in
every example we present in the paper and it is due to the
fact that the quotient between the total number of females and
the total number of individuals, which are observed, converges precisely to α when the number of generations tends to infinity (see <cit.>), therefore, the similitude between the chosen simulated paths and the observed sample
makes the estimation for α good enough.

However, we can observe in Figure <ref> that the estimation of
the posterior distribution of the parameter β is not very
accurate because P(β=0|ℱℳ_15) is very high, estimated by 0.716, despite the real value of β is strictly greater than 0. Note at this point that, to estimate β is a difficult task. First, because in general, its value, although positive, is very small in real situations (0.005 in
our example), close to zero (which represents the non-mutation). Secondly, due to the fact that from the total number of males with 𝗋-allele, 𝖬_𝗇^𝗋, is not possible to know, without some additional information, how many of those come from
mutations, 𝖬_𝗇^𝖱→ 𝗋. On the other hand, the corresponding estimates of 𝗆_𝖱|ℱℳ_15 and 𝗆_𝗋|ℱℳ_15  are enough accurate, with the last one better than the first. This is due to the fact that we are in the case in which the 𝗋-allele is the dominant one
(𝗆_𝗋≥(1-β)𝗆_𝖱) and therefore one has more
information about males with the mutant allele, in spite of the noise produced by the non-observed variable 𝖬_𝗇^𝖱→ 𝗋.
Anyway, if one goes more in deep analyzing these data, one has to consider that both approximations of the posterior distributions (of 𝗆_𝖱 and 𝗆_𝗋) are related with β, as it is shown in Figure <ref>, and therefore the estimates of these parameters inherit in some sense the inaccuracy of β.
Actually, if we consider only the simulated paths where β^𝗌𝗂𝗆=0, the kernel density estimate of 𝗆_𝖱 is really accurate (of course, since the true value of β is close to zero).




To give a measure of the accuracy of the method for the different
parameters, we consider the relative mean square error (RMSE), which
was also proposed in <cit.> and <cit.>,
calculated by

    1/n∑_k=1^n(δ^𝗌𝗂𝗆_k-δ)^2/δ^2,
 with
n=1000, δ the true value of α, β,
𝗆_𝖱 or 𝗆_𝗋, the corresponding in each case,
and δ^𝗌𝗂𝗆_k the corresponding value of δ
on the kth simulated path chosen by the method.




In particular, Table <ref> shows the RMSE of the estimates of α,
β, 𝗆_𝖱 and 𝗆_𝗋 given by the Tolerance
Rejection-ABC Algorithm when the sample ℱℳ_15 is observed and considering
all chosen simulated paths (i.e. those simulated paths such that ρ(ℱℳ_15^𝗌𝗂𝗆,ℱℳ_15)≤ε), all chosen simulated paths where β^𝗌𝗂𝗆=0 and all chosen simulated paths where β^𝗌𝗂𝗆>0. One can appreciate that, in general, the RMSE for α, 𝗆_𝖱 and 𝗆_𝗋 are very similar in all cases and very close to 0.

However, the RMSE for β when only simulated paths where β^𝗌𝗂𝗆>0 are considered, takes a high value, considerably greater than the value when all simulated paths are considered, even being the first one the real situation. This is due to the fact that P(β=0|ℱℳ_15) is very high and that the true value of β is very close to 0.

Anyway, note that although the methodology cannot
provide an adequate approximate posterior distribution of the
parameter β and consequently of the parameters 𝗆_𝖱 and 𝗆_𝗋 either, it can
provide very accurate approximate posterior distributions of the
rates of growth of both alleles, see Figure <ref>, since the total number of males of each genotype is observed. In particular, as it was indicated in Section 2, when α<0.5 and 𝗆_𝗋>(1-β)𝗆_𝖱 (as it is the
case of our example), on the set of coexistence of both alleles, the
rate of growth of the mutant allele is equal to
α𝗆_𝗋, while the rate of growth of the
𝖱-allele is equal to α(1-β)𝗆_𝖱.

Notice here that, although we do not know, a priori, whether ℱℳ_𝖭 will provide or not accurate estimates of the parameters of the model in other cases different from 𝗆_𝗋≥ (1-β)𝗆_𝖱, we are looking for a unified estimation procedure whose behavior does not depend on the parameters relation. This is why in the next section we will modify the previous sample scheme including additional information. This new sample scheme will be use in the rest of the paper.







§ INTRODUCING ADDITIONAL INFORMATION: A NEW SAMPLE SCHEME


Up to now, we have used the ABC algorithm to estimate the main
parameters of the model given ℱℳ_𝖭. However, we have seen in the example given in Section <ref> that
the estimate of the parameter β is not very accurate,
having the approximate posterior distribution huge variability. In the example we have observed an atom at zero of size 0.716 when actually the true parameter is really small but not null.

Therefore, it seems reasonable to think that it is necessary to get some information about the number of mutant alleles stemming from 𝖱-fathers, at least in
some generation, in order to obtain a more accurate approximation of the posterior distribution of β. In particular, we introduce this kind of information for the last generation. Hence, we assume from now on that the available
sample consists of the sample given by (<ref>) as well as the total number of 𝗋-males stemming from 𝖱-fathers in the
last generation, that is, 𝖬^𝖱→ 𝗋_𝖭. Moreover, as 𝖬^𝗋_𝖭 is known, the total number of 𝗋-males stemming from 𝗋-fathers in the last generation, 𝖬^𝗋→ 𝗋_𝖭, is also derived (see (<ref>)). Notice that, to obtain 𝖬^𝖱→ 𝗋_𝖭 and 𝖬^𝗋→ 𝗋_𝖭, it would be necessary to know who is every 𝗋-male's father of the generation 𝖭.
Therefore, it is plausible to assume that 𝖬^𝖱_𝖭-1 and 𝖬^𝗋_𝖭-1 are also observed, including males whom do not produce descendants. From now on, we denote this sample as ℱℳ_𝖭. Therefore,  
    ℱℳ_𝖭={ℱℳ_𝖭,𝖬^𝖱_𝖭-1,𝖬^𝗋_𝖭-1,𝖬^𝖱→ 𝗋_𝖭,𝖬^𝗋→ 𝗋_𝖭}.




§ A SERIES OF SIMULATED EXAMPLES BASED ON THE OBSERVED SAMPLE ℱℳ_𝖭


In the following subsections, we will illustrate, by means of simulated examples, how the Tolerance Rejection-ABC Algorithm works to approximate the posterior distribution θ|ℱℳ_𝖭. We will consider different situations depending on whether some variables of the sample ℱℳ_𝖭 are positive or null.



 §.§ Observing 𝖬^𝖱_𝖭>0, 𝖬^𝖱→ 𝗋_𝖭>0 and 𝖬^𝗋→ 𝗋_𝖭>0


We first consider the situation in which 𝖬^𝖱_𝖭>0, 𝖬^𝖱→ 𝗋_𝖭>0 and 𝖬^𝗋→ 𝗋_𝖭>0. This implies that 𝖬^𝖱_𝖭-1>0 and 𝖬^𝗋_𝖭-1>0.
Moreover, this assumption also implies that β and 𝗆_𝗋 are strictly positive and then their posterior distributions are not concentrated at zero value, which simplifies the Tolerance Rejection-ABC Algorithm described in Subsection <ref> because only simulated paths where β^𝗌𝗂𝗆>0 and 𝗆_𝗋^𝗌𝗂𝗆>0 will be considered. On the other hand, the metric is slight more complex, including the new observed variables 𝖬^𝖱_𝖭-1, 𝖬^𝗋_𝖭-1, 𝖬^𝖱_𝖭, 𝖬^𝖱→ 𝗋_𝖭 and 𝖬^𝗋→ 𝗋_𝖭 in the same way as previously. In particular, the distance between the simulated path, ℱℳ_𝖭^, and the observed data, ℱℳ_𝖭, is defined as

    ρ^*(ℱℳ_𝖭^,ℱℳ_𝖭)=      
          (∑_𝗇=1^𝖭(^/-
    /^)^2+∑_𝗇=1^𝖭-2(^/- /^)^2.
          +(^/- /^)^2+(^/-
    /^)^2
          +(^/- /^)^2+(^/-
    /^)^2
          +.(^/- /^)^2)^1/2




  §.§.§ Case 𝗆_𝗋≥ (1-β)𝗆_𝖱


To illustrate how to approximate the posterior distribution
θ|ℱℳ_𝖭, first we study again the case 𝗆_𝗋≥(1-β)𝗆_𝖱,
considering the same observed sample given in Table <ref> and also
assuming that now it is observed that 𝖬^𝖱→ 𝗋_15=6 (i.e. from 45850 males with 𝗋-allele in generation 15, 6 of them come from mutations), that 𝖬^𝗋→ 𝗋_15=45844, that 𝖬^𝖱_14=754 and that 𝖬^𝗋_14=24687.
With this new
information, we apply the Tolerance Rejection-ABC Algorithm using the
metric ρ^*(· ,·).

In Figure <ref>, we present the approximate posterior
distributions of all parameters together with the true value of
the parameters (solid line). One can appreciate
how the approximate posterior distribution of β has improved compared with the corresponding approximation given in Figure <ref>. Now, the true value of all parameters are into the 95% HPD sets and the corresponding RMSE for β and 𝗆_𝖱 are, respectively, 18.392 and 0.035, considerably smaller than that given in Table <ref> (1557.8 and 0.1141, respectively) where only simulated paths with β^𝗌𝗂𝗆>0 were considered, as it is now our case. For the rest of the parameters, the approximate posterior distributions in Figure <ref> are very similar to that given in Figure <ref> being the corresponding RMSE for α and 𝗆_𝗋, 0.037 and 0.053, respectively, similar values to that given in Table <ref> (0.0077 and 0.0076, respectively) where only simulated paths with β^𝗌𝗂𝗆>0 were considered.
Moreover, since the range of the posterior
distribution of β is very small, its estimation does not affect to the estimation of neither 𝗆_𝖱 nor 𝗆_𝗋, which are positively correlated (see contour plots showed in Figure <ref>).








  §.§.§ Case 0<𝗆_𝗋<(1-β)𝗆_𝖱


Next we illustrate how the algorithm works to approximate the posterior distribution
θ|ℱℳ_𝖭, in the case
0<𝗆_𝗋<(1-β)𝗆_𝖱. To this end, we consider a second
simulated example with initial
vector (𝖥_0, 𝖬_0^𝖱, 𝖬_0^𝗋)=
(10,5,5) as in the previous case, and parameter vector θ=(α, β,
𝗆_𝖱, 𝗆_𝗋)=(0.45, 0.01, 3.5, 2.6). For a Y-BBP with mutations with this set of parameters
and initial values, we also proved in <cit.> that there
exists a positive probability of survival of both genotypes.



We simulate 15 generations of this Y-BBP with mutations assuming
that reproduction laws of both genotypes follow the non-parametric offspring distributions with finite support given in Table <ref>, with means 𝗆_𝖱=3.5 and
𝗆_𝗋=2.6. The simulated data can be seen in Table
<ref> and they are denoted by
ℱℳ_15.



We now plot (see Figure <ref>) the approximate posterior
distributions of the parameters, once the algorithm have been applied. Again, we can appreciate that the methodology provides accurate approximations to the posterior distributions of all parameters in this new context, being the RMSE for α, β, 𝗆_𝖱 and 𝗆_𝗋, respectively, 0.0285, 1.7918, 0.0131, 0.0719. Notice that, in this case, the RMSE for β and 𝗆_𝖱 are smaller than that given in the case 𝗆_𝗋≥ (1-β)𝗆_𝖱 in the previous subsection, since in the case 𝗆_𝗋<(1-β)𝗆_𝖱 we have more information on these parameters because the rate that define the growth is essentially (1-β)𝗆_𝖱 that allows to obtain more accurate approximations of β|ℱℳ_15 and 𝗆_𝖱|ℱℳ_15.

Therefore, as final conclusion of this Subsection <ref>, we can establish that the proposed Tolerance Rejection-ABC Algorithm works adequately to estimate the parameters of a Y-BBP with mutations, given the information provided by the sample ℱℳ_𝖭 with 𝖬_𝖭^𝖱>0, 𝖬_𝖭^𝖱→ 𝗋>0 and 𝖬_𝖭^𝗋→ 𝗋>0 whichever the relation between the parameters.





 §.§ Observing 𝖬^𝖱_𝖭>0, >0 and =0


In the previous subsection it was considered that 𝗆_𝗋>0 since  was assumed to be non-null. Now, we research the situation where =0. Obviously, this event occurs in models with 𝗆_𝗋=0 but it can be also observed in models with 𝗆_𝗋>0. Due to this fact, the estimation of 𝗆_𝗋 in this case can be a difficult task. Most probably the approximate posterior distribution of 𝗆_𝗋|ℱℳ_𝖭 will present an atom at zero with non-null probability. This kind of problems are usual in branching process theory. For example, based on the observation of a Galton-Watson process it is difficult to make inference on whether the extinction or explosion of such process will occur (see <cit.> or <cit.>).

The algorithm works in the same way as it was described in Subsection <ref>, now
using the metric ρ^*(·,·). In this case we only consider simulated paths ℱℳ_𝖭^ such that  𝖬^𝗋→ 𝗋_𝖭^=0. Therefore, the last sum term of ρ^*(·,·) is deleted.

To illustrate this particular case, we fix the para­meter vector
θ=(α, β, 𝗆_𝖱, 𝗆_𝗋)=(0.45, 0.10,
3, 0) and initial vector (𝖥_0, 𝖬_0^𝖱,
𝖬_0^𝗋)= (10,5,5). For a Y-BBP with mutations with
this set of parameters and initial values, we proved in
<cit.> that there exists a positive probability of
survival of the 𝖱-genotype and therefore also of the
𝗋-genotype.




We simulate 15 generations of this Y-BBP with mutations assuming
that the reproduction law of 𝖱-geno­type follows non-parametric offspring distribution with finite support given in Table <ref>, with mean 𝗆_𝖱=3. The simulated data can be seen in Table <ref> and they are denoted, as in the
previous cases, by ℱℳ_15.





In Figure <ref>, we present the approximate posterior densities of every parameter. The approximations of α|ℱℳ_15, β|ℱℳ_15 and 𝗆_𝖱|ℱℳ_15 have been calculated considering all chosen simulated paths. In those cases, the algorithm provides accurate approximations to the posterior densities of all parameters with the 95% HPD sets containing their true values and with small values of their RMSE (see Table <ref>). On the other hand, the approximation of 𝗆_𝗋|ℱℳ_15 has been obtained considering only chosen simulated paths where 𝗆_𝗋^𝗌𝗂𝗆>0. In this case, we also represent in such figure the P(𝗆_𝗋=0|ℱℳ_15) (area of the vertical solid bar) which is estimated by 0.504.

At this point, from the estimates and the observed sample, one can wonder about the following hypothesis test:

    H_0:𝗆_𝗋=0   H_1: 𝗆_𝗋>0.
 Considering that we have assumed in the implementation of the algorithm that 𝗆_𝗋^𝗌𝗂𝗆 could take the value 0 with probability ϕ^𝗌𝗂𝗆, being ϕ^𝗌𝗂𝗆∼ U(0,1), we consider its expected value at calculating the Bayes factor, K, and therefore, it is verified that P(𝗆_𝗋>0)=P(𝗆_𝗋=0) and then 
    K=P(𝗆_𝗋=0|ℱℳ_15)P(𝗆_𝗋>0)/P(𝗆_𝗋>0|ℱℳ_15)P(𝗆_𝗋=0)=0.504/0.496=1.06.

Although the Bayes factor is greater than 1 and this leads us to conclude that 𝗆_𝗋=0 is supported by the observed sample, it is also true that the value of K is very close to 1 and then, the acceptance of H_0 is not strongly supported.
For that reason, in Figure <ref> we present a comparison of the approximate posterior densities of the parameters α, β and 𝗆_𝖱 considering simulated paths where 𝗆_𝗋^𝗌𝗂𝗆=0 (dotted line) and simulated paths where 𝗆_𝗋^𝗌𝗂𝗆>0 (solid line). Notice that the true values of the three parameters are into 95% HPD sets in both cases. Moreover, in Table <ref> are presented the RMSE for the estimates of all parameters for these cases. One can appreciate that in all cases the RMSE for α is very similar and  close to 0. Moreover, the RMSE for β and 𝗆_𝖱 take their smaller values when only simulated paths where 𝗆_𝗋^𝗌𝗂𝗆>0 are considered. This is due to the close relation of these parameters so, when 𝗆_𝗋^𝗌𝗂𝗆>0 the values of β^𝗌𝗂𝗆 and 𝗆_𝖱^𝗌𝗂𝗆 are smaller than in the case 𝗆_𝗋^𝗌𝗂𝗆=0 since the 𝗋-males do not stem only from mutations.





We finally estimate the difference in means of the approximate posterior densities of each parameter between these two groups (we name Group A1 to the set of all chosen simulated paths where 𝗆_𝗋^𝗌𝗂𝗆=0 and name Group A2 to the set of all chosen simulated paths where 𝗆_𝗋^𝗌𝗂𝗆>0) using the Bayesian alternative to the t test (see <cit.>). We obtain that the 95% HPD for β and 𝗆_𝖱 are, respectively, (0.02236, 0.0536) and (0.0472,0.2286) which do not include zero. The Bayes factors are, respectively, 4704.56 and 7.15, and the probabilities that the true values of the differences are greater than zero are, respectively, 100% and 99.8% which leads us to conclude that there exist significant differences in the means of the approximate posterior densities of β|ℱℳ_15 and 𝗆_𝖱|ℱℳ_15 between Groups A1 and A2. However, the 95% HPD for α is (-0.0122, 0.0162) which includes zero. The Bayes factor in this case is 0.104 and the probability that the true value of the differences is greater than zero is 60%. This leads us to conclude that there are no significant differences in the means of the approximate posterior density of α|ℱℳ_15 between Groups A1 and A2. Visually, one can appreciate such differences in Figure <ref>.



 §.§ Observing 𝖬^𝖱_𝖭>0, =0 and >0


In a similar way than in previous subsection, next we describe the algorithm when it has been observed that 𝖬^𝖱_𝖭>0, =0 and 𝖬^𝗋→ 𝗋_𝖭>0. Obviously, this event occurs in models with β=0, but it can be also observed in models with β>0. Due to this fact, as we pointed out previously, the estimation of β in this case can be a difficult task. Now, the approximate posterior distribution of β|ℱℳ_𝖭 will present an atom at zero with non-null probability. The algorithm works in the same way as it was described in Subsection <ref>, now using the metric ρ^*(·,·). In this case we only consider simulated paths ℱℳ_𝖭^ such that  𝖬^𝖱→ 𝗋_𝖭^=0. Therefore, the next-to-last sum term of ρ^*(·,·) is deleted.

To illustrate this particular case, we fix the para­meter vector
θ=(α, β, 𝗆_𝖱, 𝗆_𝗋)=(0.65, 0.01,
3, 3.5) and initial vector (𝖥_0, 𝖬_0^𝖱,
𝖬_0^𝗋)= (10,5,5). For a Y-BBP with mutations with
this set of parameters and initial values, we proved in
<cit.> that there exists a positive probability of
survival of the 𝖱-genotype and therefore also of the
𝗋-genotype.




We simulate 15 generations of this Y-BBP with mutations assuming
that reproduction laws of  𝖱 and 𝗋-genotypes follow non-parametric offspring distributions with finite support given in Table <ref>, with means 𝗆_𝖱=3 and 𝗆_𝗋=3.5. The simulated data can be seen in Table
<ref> and they are denoted, as in the previous cases, by ℱℳ_15.









In Figure <ref>, we show the approximate posterior densities of every parameter. The approximations of α| ℱℳ_15, 𝗆_𝖱|ℱℳ_15 and 𝗆_𝗋|ℱℳ_15 have been obtained considering all chosen simulated paths. In all those cases, the algorithm provides accurate approximations to the posterior densities of all parameters with the 95% HPD sets containing their true values and with small values of their RMSE (see Table <ref>). On the other hand, the approximation of β|ℱℳ_15 has been obtained considering only chosen simulated paths where β^𝗌𝗂𝗆>0. In this case the 95% HPD set also contains the true value of the parameter so the approximation to the posterior density is also considered enough accurate. Moreover, we also represent the P(β=0|ℱℳ_15) (area of the vertical solid bar) which is estimated by 0.152.

As in the previous subsection, from the estimates and the observed sample, one can wonder about the following hypothesis test:

    H_0:β=0   H_1: β>0.
 Considering that we have assume in the implementation of the algorithm that β^𝗌𝗂𝗆 could take the value 0 with probability γ^𝗌𝗂𝗆, being γ^𝗌𝗂𝗆∼ U(0,1), we consider its expected value at calculating the Bayes factor, K, and then, it is verified that P(β>0)=P(β=0), thus 
    K=P(β=0|ℱℳ_15)P(β>0)/P(β>0|ℱℳ_15)P(β=0)=0.152/0.848=0.18.

That value of the Bayes factor leads us to conclude (see <cit.>) that there are substantial evidences against the null hypothesis, and then β>0 is more supported by the observed sample, which is the real situation.

In Figure <ref> we present the approximate posterior densities of α|ℱℳ_15, 𝗆_𝖱|ℱℳ_15 and 𝗆_𝗋|ℱℳ_15 calculated considering only paths where β^𝗌𝗂𝗆>0. Notice that, the true values of the all three parameters are into 95% HPD sets. Moreover, in Table <ref> are presented the RMSE for the estimates of those parameters. Notice that the RMSE are very similar to those calculated when all chosen simulated paths are considered and they are very close to 0.



As final conclusion of Subsections <ref> and <ref> we establish that, if in the observed sample, one of the random variables 𝖬^𝖱→ 𝖱_𝖭 or 𝖬^𝖱→ 𝗋_𝖭 is equal to 0 then we apply the Tolerance Rejection-ABC Algorithm and solve the corresponding hypothesis test considering the approximate posterior distributions conditioned to the decision given by this test.

In order not to extend the paper, we have not considered explicitly in subsection <ref> an example where =0 and 𝗆_𝗋>0, however an example of this kind of situation is considered in subsection <ref> where =0 and β>0. Analogously, it has not been considered explicitly in subsection <ref> an example where =0 and β=0, however an example of this kind of situation is consider in subsection <ref> where =0 and 𝗆_𝗋=0. In both cases, the results are analogous to those shown in the paper.

The case =0, i.e. =0 and 𝖬^𝗋→ 𝗋_𝖭=0, is not illustrated in the paper. A sample where =0 is observed could belong to a coexistence path although it would not be guaranteed. Anyway, to make inference about the parameters in this situation, the Tolerance Rejection-ABC Algorithm would be applied and then both hypothesis test, in (<ref>) and (<ref>), should be solved.

 The case =0 is not illustrated either in the paper, since this case represents the extinction of the 𝖱-allele and then the behavior of 𝗋-allele is described by a two-sex Galton-Watson process (see <cit.>).



§ SENSITIVITY ANALYSIS


In this section we examine the sensitivity of inferences depending on the probability distribution used to generate the simulated paths. We apply  the Tolerance Rejec­tion-ABC Algorithm to the examples in Subsections 6.1.1 and 6.1.2. but now generating the pool of simulated paths, instead from the Poisson distribution, from negative binomial distribution laws with different value of size (𝗄) since these kind of distributions have been also used in practical cases (see <cit.>, <cit.> or <cit.>).

In Tables <ref> and <ref> we present the point estimates of  α, β, 𝗆_𝖱 and 𝗆_𝗋, under squared error loss as
well as their 95% HPD sets for the two examples, respectively. It can be
seen that in all cases the HPD sets contain the true values of the parameters, being very similar among them for different distributions.
This allow us to conclude that this is a robust methodology against the probability distribution used to simulate the processes.








§ PREDICTION OF THE FUTURE POPULATION SIZE


Finally, once that the algorithm has been proved to be a useful tool to obtain accurate approximations of the posterior distributions of the parameters, from them, we can also estimate others random variables related to the process. For instance, from a practical standpoint, it is  of interest to infer the size of future generations. Next, we apply a Monte Carlo procedure, proposed in <cit.>, to approximate the predictive distributions.
In particular, for each θ^(i), i=1,…,m sampled from θ|ℱℳ_𝖭, one can simulate s process until the 𝗅th generation, which started with  females,
 𝖱-males and  𝗋-males, obtaining values to approximate the predictive posterior distributions (𝖥_𝖭+𝗅, 𝖬^𝖱_𝖭+𝗅, 𝖬^𝖱→ 𝗋_𝖭+𝗅, 𝖬^𝗋→ 𝗋_𝖭+𝗅)|ℱℳ_𝖭
and (𝖹^𝖱_𝖭+𝗅, 𝖹^𝗋_𝖭+𝗅)|ℱℳ_𝖭 by Gaussian kernel estimators.

To illustrate this procedure, we consider the example given en subsection 6.1.2 considering the observed sample ℱℳ_15 given in Table <ref>,
m=1000, s=2000 and 𝗅=1. Concretely, we simulate a generation of 2000 processes started with (𝖥_15,𝖬^𝖱_15,𝖬^𝗋_15)=(5437,6351,258), for each parameter θ^(i), i=1,…,1000.

Figure <ref> shows the approximated predictive posterior distributions for 𝖥_16, 𝖬^𝖱_16, 𝖬^𝖱→ 𝗋_16 and 𝖬^𝗋→ 𝗋_16, given ℱℳ_15 in Table <ref>. Notice that these estimates are in accordance with the relation between the parameters and with the observed sample, where the -allele is the dominant one.



The software environment  for statistical   computing and graphics 𝐑 (“GNU  S”, see <cit.>) has been used to perform the ABC methodology and the simulation study. To calculate the kernel density estimation the GenKern  package (see <cit.>) and sm package (see <cit.>) have been used. To implement the Bayesian t-test, the BayesFactor and BEST packages (see <cit.> and <cit.>, respectively) have been applied.




§ CONCLUDING REMARKS


The aim of this work has been to develop Bayesian
inference theory for a Y-linked two-sex branching process with blind
choice which is useful to model the evolution of the number of carriers of two alleles (named as  and 𝗋) of a Y-linked gene considering the possibility of mutations from -allele to the 𝗋-allele.

We have focussed mainly on approximating the posterior distributions of the main parameters of such mo­del considering for that, at the beginning, a realistic sampling scheme where the observation of the total number of females and males in each generation is assumed as well as the observation of the total number of each type of males (males with -allele and males with 𝗋-allele) in the last generation. Then, we have described the development of a method based on the Approximate Bayesian Computation (ABC) methodology (Tolerance Rejection-ABC Algorithm) to approximate the posterior distributions of the model parameters based on such sample scheme.

We have shown throughout a simulated example that
the methodology presents difficulties to estimate the posterior
distribution of the probability of mutation, β, due to the fact that with the observed sample it is not possible to know how many of the observed 𝗋-alleles stem from mutations. For that reason, we consider another sampling scheme where also is observed, in the last generation, the total number of 𝗋-males stemming from -fathers as well as the total number of each type of males in the penultimate generation.

We have illustrated how the Tolerance Rejection-ABC Algorithm works based on this sampling scheme and considering different situations which can be observed in the sample in the case of coexistence of both alleles. In this sense, we have considered special situations which can be observed in the last generation of the sample: when there are the two types of males (i.e. >0 and >0), when there are not 𝗋-males stemming from 𝗋-fathers (i.e. =0) and when there are not 𝗋-males stemming from -fathers (i.e. =0). In all cases, we have obtained accurate approximations to the posterior densities of all parameters with the 95% HPD sets containing the true values of the parameters.

The case where =0 is the special interest because it is no possible to know whether the mean number of individuals stemming from 𝗋-couples, 𝗆_𝗋, is equal to 0 or strictly positive. Analogously, the case where =0 is interesting because it is no possible to know whether the probability of mutation is equal to 0 or strictly positive. In both cases, after applying the ABC methodology, we have proposed a hypothesis test to decide the more plausible option (see (<ref>) and (<ref>)). In the two considered examples, the Bayes factor has lead us to conclude that the true situation was the supported one by the observed sample.

Notice that we have taken 15 generations in the sample schemes of all simulated examples considered in the paper. We considere that this is a balanced number in the sense that it is big enough to observe whether one of the alleles is the dominant and also it is a feasible number to be observed in many animal populations with sex reproduction. Moreover, in the examples, we have covered all possible situations between the parameters in the coexistence set taking into account the different magnitudes of the rates of growth.

We have also studied the robustness of the methodology by mean of a general simulated experiment where we have applied the methodology for different base distributions concluding that this is a robust methodology against the probability distribution used to simulate the processes.

Finally, we have been able to predict the future population size approximating the predictive distributions of the random variables related to the total number of females and the total number of each type of males in the following generation to the last one observed.

Note that the Approximate Bayesian Computation
is a proved statistical tool very useful for inference in
parameters of complex models in population genetics as is our case.
It is easy and fast to simulate from our model, and therefore,
in this case has been more convenient that the Gibbs sampler.


This research was supported by Grant MTM2015-70522-P (MINECO/FEDER, UE) and Grant IB16­103 (Junta de Extremadura / Fondo Europeo de Desarrollo Regional, UE).










31


urlstyle

[Beaumont et al.(2002)Beaumont, Zhang, and Balding]beaumont
M.A. Beaumont, W. Zhang, and D.J. Balding.
Approximate Bayesian Computation in population genetics.
Genetics, 162:0 2025–2035, 2002.

[Bertoin et al.(2008)Bertoin, Fontbona, and Martínez]Bertoin08
J. Bertoin, J. Fontbona, and S. Martínez.
On prolific individuals in a supercritical continuous-state branching
  process.
J. Appl. Probab., 3:0 714–726, 2008.

[Blumberg and Lloyd-Smith(2013)]Bl
S. Blumberg and J.O. Lloyd-Smith.
Inference of r0 and transmission heterogeneity from the size
  distribution of stuttering chains.
PLOS Comput. Biol., 9(5):0 1–17, 2013.

[Bowman and Azzalini(2014)]sm
A.W. Bowman and A. Azzalini.
R package : nonparametric smoothing methods,
  2014.
URL <URL http://www.stats.gla.ac.uk/ adrian/sm,
  http://azzalini.stat.unipd.it/Book_sm>.
R package version 2.2-5.4.

[Charchar et al.(2012)Charchar, Bloomer, Barnes, Cowley, Nelson, Wang,
  Denniff, Debiec, Christofidou, Nankervis, Dominiczak, Bani-Mustafa,
  Balmforth, Hall, Erdmann, Cambien, Deloukas, Hengstenberg, Packard,
  Schunkert, Ouwehand, Ford, Goodall, Jobling, Samani, and
  Tomaszewski]Charchar
F.J. Charchar, L. Bloomer, T.A. Barnes, M.J. Cowley, C.P. Nelson, Y. Wang,
  M. Denniff, R. Debiec, P. Christofidou, S. Nankervis, A.F. Dominiczak,
  A. Bani-Mustafa, A.J. Balmforth, A.S. Hall, J. Erdmann, F. Cambien,
  P. Deloukas, C. Hengstenberg, C. Packard, H. Schunkert, W. H. Ouwehand,
  I. Ford, A.H. Goodall, M.A. Jobling, N.J. Samani, and M. Tomaszewski.
Inheritance of coronary artery disease in men: an analysis of the
  role of the Y chromosome.
The Lancet, 379:0 915 – 922, 2012.
ISSN 0140-6736.
10.1016/S0140-6736(11)61453-0.

[Daley(1968)]daleyA
D.J. Daley.
Extinction conditions for certain bisexual Galton-Watson branching
  processes.
Z. Wahrscheinlichkeitsth., 9:0 315–322, 1968.

[Farrington and Grant(1999)]fg
C.P. Farrington and A.D. Grant.
The distribution of time to extinction in subcritical branching
  processes: applications to outbreaks of infectious disease.
J. Appl. Probab., 36:0 771–779, 1999.

[Farrington et al.(2003)Farrington, Kanaan, and Gay]fkg
C.P. Farrington, M.N. Kanaan, and N.J. Gay.
Branching process models for surveillance of infectious diseases
  controlled by mass vaccination.
Biostatistics, 4(2):0 279–295, 2003.

[González et al.(2006)González, Hull, Martínez, and
  Mota]ghmm06
M. González, D.M. Hull, R. Martínez, and M. Mota.
Bisexual branching processes in a genetic context: The extinction
  problem for Y-linked genes.
Math. Biosci., 202:0 227–247, 2006.

[González et al.(2009)González, Martínez, and Mota]gmm09
M. González, R. Martínez, and M. Mota.
Bisexual branching processes to model extinction conditions for
  Y-linked genes.
J. Theor. Biol., 258:0 478–488, 2009.

[González et al.(2012)González, Gutiérrez, and
  Martínez]ggm2012a
M. González, C. Gutiérrez, and R. Martínez.
Extinction conditions for Y-linked mutant-alleles through two-sex
  branching processes with blind-mating structure.
J. Theor. Biol., 307:0 104–116, 2012.

[González et al.(2013a)González, Gutiérrez, and
  Martínez]ggm2013a
M. González, C. Gutiérrez, and R. Martínez.
Parametric Bayesian inference for Y-linked two-sex branching
  models.
Stat. Comput., 23(6):0 727–741, 2013a.

[González et al.(2013b)González, Gutiérrez,
  Martínez, and del Puerto]ggmp2013b
M. González, C. Gutiérrez, R. Martínez, and I.M. del Puerto.
Bayesian inference for controlled branching processes through MCMC
  and ABC methodologies.
RACSAM, 107:0 459–473, 2013b.

[González et al.(2016)González, Gutiérrez, Martínez,
  Minuesa, and del Puerto]ggmmp2016
M. González, C. Gutiérrez, R. Martínez, C. Minuesa, and I. del Puerto.
Bayesian analysis for controlled branching process.
Branching Processes and Their Applications (del Puerto, I.,
  González, M., Gutiérrez, C., Martínez, R., Minuesa, C. Molina, M.,
  Mota, M. and Ramos, A., eds.) Lecture Notes in Statistics, 219:0
  185–205, Springer–Verlag, 2016.

[Gutiérrez(2012)]cgutierrez
C. Gutiérrez.
Two-Sex Branching Processes in a Genetic Context.
PhD thesis, Department of Mathematics. University of Extremadura,
  2012.

[Guttorp and Perlman(2013)]Guttorp2013
P. Guttorp and M.D. Perlman.
Predicting extinction or explosion in a Galton-Watson branching
  process.
J. Stat. Infer. Stoch. Proc., 16:0 113–125, 2013.

[Guttorp and Perlman(2015)]Guttorp2015
P. Guttorp and M.D. Perlman.
Predicting extinction or explosion in a Galton-Watson branching
  process with power series offspring distribution.
J. Stat. Plan. Infer., 167:0 193–215, 2015.

[Jeffreys(1961)]jef
H. Jeffreys.
Theory of probability (3rd ed.).
Oxford University Press, 1961.

[Knuth(2005)]knuth
D. Knuth.
The art of the computing programing. Vol. 1.
Pearson Education, Inc., 2005.

[Kruschke and Meredith(2017)]BEST
J.K. Kruschke and M. Meredith.
BEST: Bayesian Estimation Supersedes the t-Test, 2017.
URL <https://CRAN.R-project.org/package=BEST>.
R package version 0.5.0.

[Lintusaari et al.(2017)Lintusaari, Gutmann, Dutta, Kaski, and
  Corander]lintusaari
J. Lintusaari, M.U. Gutmann, R. Dutta, S. Kaski, and J. Corander.
Fundamentals and recent develpments in Approximate Bayesian
  Computation.
Syst. Biol., 66(1):0 e66–e82, 2017.

[Lucy and Aykroyd(2010)]GenKern
D. Lucy and R. Aykroyd.
GenKern: Functions for generating and manipulating binned
  kernel density estimates, 2010.
URL <http://CRAN.R-project.org/package=GenKern>.
R package version 1.1-10.

[Marin et al.(2012)Marin, Pudlo, Robert, and Ryder]marin
J.M. Marin, P. Pudlo, C.P. Robert, and R.J. Ryder.
Approximate Bayesian Computational methods.
Stat. Comput., 22:0 1167–1180, 2012.

[Mode and Sleemam(2000)]modesleemam
C.J. Mode and C.K. Sleemam.
Stochastic processes in epidemiology.
World Scientific, 2000.

[Morey and Rouder(2015)]BayesFactor
R.D. Morey and J.N. Rouder.
BayesFactor: Computation of Bayes Factors for Common Designs,
  2015.
URL <https://CRAN.R-project.org/package=BayesFac­tor>.
R package version 0.9.12-2.

[Pakes(2003)]pakes
A.G. Pakes.
Biological applications of branching processes.
Handbook of Statistic. Stochastic Processes: Modelling
  and Simulation (Shanbhag, D.N. and Rao, C.R., eds.), Vol. 21. Chapter
  18:0 693–773, 2003.

[Pritchard et al.(1999)Pritchard, Seielstad, Pere-Lezaun, and
  Feldman]pri
J. Pritchard, M. Seielstad, A. Pere-Lezaun, and M. Feldman.
Population growth of human Y chromosomes: a study of Y
  chromosome microsatellites.
Mol. Biol. Evol., 16:0 1791–1798, 1999.

[R Development Core Team(2011)]r
R Development Core Team.
R: A Language and Environment for Statistical Computing.
R Foundation for Statistical Computing, Vienna, Austria, 2011.
URL <http://www.R-project.org/>.
ISBN 3-900051-07-0.

[Sun and Heitman(2012)]sun
S. Sun and J. Heitman.
Should Y stay or should Y go: The evolution of non-recombining
  sex chromosomes.
BioEssays, 3420 (11):0 938–942, 2012.

[Sunnaker et al.(2013)Sunnaker, Busetto, Numminen, Corander, Foll, and
  Dessimoz]sunnaker
M. Sunnaker, A.G. Busetto, E. Numminen, J. Corander, M. Foll, and C. Dessimoz.
Approximate Bayesian Computation.
PLoS Comput. Biol., 9(1):0 e1002803, 2013.

[Tavaré et al.(1997)Tavaré, Balding, Griffiths, and
  Donnelly]tbgd1997
S. Tavaré, D.J. Balding, R.C. Griffiths, and P. Donnelly.
Inferring coalescence times from DNA sequence data.
Genetics, 145(2):0 505–518, 1997.



