Expectation-Maximization Algorithm for Determining Natural Selection of Y-Linked Genes Through Two-Sex Branching Processes

1. Introduction

Recent research has shown the importance of certain genes linked to the Y chromosome in populations of humans (Hughes et al., 2005) and other animals (Charlesworth et al., 2005). This chromosome has the particularity of being male-specific (the SRY gene is responsible for maleness) and haploid, and of having a region that escapes recombination (the nonrecombining region, NRY, which is 95% of the chromosome in humans—see, for example, Graves, 2006). The unique properties of the Y chromosome have major consequences for its population genetics: The NRY region passes down from father to son largely unchanged, preserving the paternal genetic legacy, and is therefore very useful for studying how populations have evolved. A history of paternal lineages can be reproduced by examining the differences (such as DNA polymorphisms) among modern Y chromosomes. There have been many studies in this sense in the context of populations of humans (e.g., The Y Chromosome Consortium studies or Rosa et al., 2007) and other species (Hellborg et al., 2005). In human populations, the surname can also be regarded as a Y-linked characteristic, and there have been studies aimed at determining its relationship with Y-chromosome lineages (King et al., 2006).

Another singular question associated with the Y chromosome is that of the microdeletions of this chromosome's long arm (Yq). The Yq deletion is associated with males who have fertility problems (Krausz et al., 2003), but many cases have been reported in which the natural transmission of this genetic defect from fathers to sons has occurred (Kuhnert et al., 2004). Obviously, determining the evolution of the number of males with this genetic defect in a human population is an important medical problem (Fitch et al., 2005), but it has also been investigated in other species (Toure et al., 2004). Moreover, there is evidence that the Y chromosome plays a role in skeletal growth, germ-cell tumorigenesis, and graft rejection, and that its genes might also influence gender-specific differences in disease susceptibility.

Appropriate mathematical models are needed to understand the evolution of Y chromosome lineages (for instance, to help solve the problem of Y-chromosomal Adam—the theoretical male who is the most recent common patrilineal ancestor of all living humans; estimations of the date of this common ancestor is an important problem), Yq deletions, or other Y-linked genes.

Many models used in population genetics are based on the Wright-Fisher model, although branching processes naturally come to mind in this context and represent a clear alternative approach. These processes are stochastic models, which arise in the description of population dynamics, being of particular use in describing the extinction/growth of populations (Haccou et al., 2005). Branching models have been applied to many biological problems in such fields as epidemiology, genetics, and cell kinetics. Examples include the evolution of infectious diseases (Garske and Rhodes, 2008), population genetics (Iwasa et al., 2005), and stem cells (Yakovlev and Yanev 2006). Further examples are reviewed in the recent monographs of Kimmel and Axelrod (2002), Pakes (2003), and González et al. (2010). A comparison between Wright-Fisher and branching models can be found in the recent article of Cyran and Kimmel (2010).

The simplest branching models are the Galton-Watson and the Markov branching processes. They have been used to model Y-chromosome lineages and their female analogues—mitochondrial DNA lineages (Neves and Moreira, 2006; Cyran and Kimmel, 2010). But more accurate models are needed in which all the phases of sexual reproduction can be considered, including the interaction between females and males in producing offspring. Recently, two models (González et al., 2006; González et al., 2009b) have been presented to describe the evolution of the number of carriers of the two alleles of a Y-linked gene (so that there are two types of male, each carrying one of these alleles) in a two-sex monogamic population. In the first, it was considered that the characters controlled by such a gene can influence the mating process of the species, with females having a preference for males carrying one of the alleles of the gene (Pidancier et al., 2006). It was shown (González et al., 2008) that this preference can sometimes be definitive in determining the survival of the different genotypes in the population. This model was denominated a Y-linked bisexual branching process (Y-linked BBP) with preference. And in the second, González et al. (2009b), it was considered that females choose their mates without caring about their genotype that is, each female makes a blind choice of the genotype of her mate. This model was called Y-linked BBP with blind choice.

The focus of the present article is the first model, a Y-linked BBP with preference, to pattern the evolution of the number of carriers of each allele of a Y-linked gene or of Y chromosome lineages in a two-sex monogamic population, assuming that this gene influences the mating process. In González et al. (2006) and González et al. (2008), it was shown that the behavior of genes that fit the pattern of a Y-linked BBP with preference is strongly related to the reproduction laws of each genotype, that is, those that model natural selection. In practice, these offspring distributions are usually unknown and need to be estimated to guarantee the applicability of these models. In this article, we deal with the problem of estimating the offspring distribution of each genotype of a Y-linked gene (as well as some related parameters, such as their mean values and future population sizes). We consider a frequentist and nonparametric framework.

First, we obtain the maximum likelihood estimators (MLEs) of the parameters when the complete family tree is observed up to some fixed generation. The limiting behavior (consistency and asymptotic normality) of these estimators is also studied. Since it is usually impossible in practice to observe the entire family tree, we then consider the problem of estimating the main parameters of the model, using only the sample given by the numbers of females and of the two different types of male in each generation, which are more easily observed. We approach this problem as an incomplete data estimation problem. This leads us to apply an expectation-maximization (EM) algorithm (McLachlan and Krishnan, 2008) in order to obtain the MLEs. (For a review on the use of this kind of algorithm in genetics see Laird, 2010).

Besides the Introduction, this article consists of four additional sections. In Section 2, we provide the definition of the Y-linked BBP with preference. In Section 3, we obtain the MLEs, assuming the complete and incomplete sampling schemes indicated above, and present the development of the EM algorithm. The accuracy of this algorithm is illustrated in Section 4 by means of a simulated example (figures are included as Supplementary Material, available online at www.liebertpub.com/cmb). Some concluding remarks are provided in Section 5. Finally, the proofs of some theoretical results related to the asymptotic properties of the MLEs based on the complete family tree sample are shown in the Supplementary Material.

2. The probability model

The probability model we deal with is the Y-linked BBP with preference that was introduced in González et al. (2006). This model describes the evolution of the number of carriers of a Y-linked gene generation-by-generation. It is assumed that the gene has a pair of alleles, denoted by R and r, which are expressed in the male phenotype (r can model the absence of R). We are thus assuming a population formed by females and males, where two types of male can be observed depending on the allele they carry. Males with R allele are denoted R males, while males with r allele are denoted r males. Hence, two types of (male–female) couple are formed—those consisting of one female and one R male (resp. r male) are denoted R (resp. r) couples.

Assuming nonoverlapping generations, and having fixed the number of couples of each type at the initial (n = 0) generation, the population size is determined in each generation according to two phases: reproduction and mating. In the reproduction phase, each couple is assumed to randomly produce offspring independently of the other couples. The probability distributions of these variables are the same for all the couples with a given genotype. Moreover, following the inheritance rules, R couples can generate females and R males, while r couples can generate females and r males (no mutation is assumed). More formally, we consider two independent sequences \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\{ ( FR_{ni}, MR_{ni} ) : i = 1, 2, \ldots; n = 0, 1, \ldots \} \ { \rm and } \ \{ ( Fr_{nj}, Mr_{nj} ) : j = 1, 2, \ldots; n = 0, 1, \ldots \} \end{align*}\end{document}

of independent, identically distributed, non-negative and integer-valued bivariate random vectors, where (FR_ni, MR_ni) (resp. (Fr_nj, Mr_nj)) represents the number of females and males generated by the i-th R couple (resp. j-th r couple) in generation n.

In general, (FR_ni, MR_ni) and (Fr_nj, Mr_nj) may have different distributions, modeling the natural selection between genotypes (i.e., their possibly different reproductive abilities). In particular, the total number of offspring generated by an R couple (resp. r couple) is specified by a probability distribution \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p^R = \{ p_k^R \} _{k \in S^R}$$\end{document} (resp. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p^r = \{ p_l^r \} _{l \in S^r}$$\end{document} ), where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p_k^R = P ( FR_{ni} + MR_{ni} = k ), k \in S^R \ ( { \rm resp}. \ p_l^r = P ( Fr_{nj} + Mr_{nj} = l ), l \in S^r )$$\end{document} , with S^R (resp. S^r) being the support of the distribution that is considered finite. This probability distribution is called the reproduction law of the R genotype (resp. r genotype). Moreover, we denote by m_R (resp. m_r) the average number of offspring per R couple (resp. r couple).

In order to model the sex designation, we consider that each offspring will be female with probability α, 0 < α < 1, or male with probability 1 − α (i.e., a binomial reproduction scheme). These sex designations are made independently among the offspring of any couple, and it is assumed that the genotype has no influence on sex determination, so that α is the same for both genotypes. Then, given an R couple (resp. r couple) that has produced k (resp. l) offspring, the number of females among these—i.e., FR_ni (resp. Fr_nj)—follows a binomial distribution of size k (resp. l) and probability α. Hence, the average number of females and males per R couple (resp. r couple) is, respectively, αm_R and (1 − α)m_R (resp. αm_r and (1 − α)m_r).

Therefore, for a generation n with total numbers of R and r couples ZR_n and Zr_n, respectively, one obtains the total number of females in generation n + 1, as well as the number of males stemming from R couples (resp. r couples) in generation n + 1, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}F_{n + 1} = \sum_{i = 1}^{ZR_n}FR_{ni} + \sum_{j = 1}^{Zr_n}Fr_{nj}, \quad MR_{n + 1} = \sum_{i = 1}^{ZR_n}MR_{ni} \quad \quad \bigg( { \rm resp.} \ Mr_{n + 1} = \sum_{j = 1}^{Zr_n}Mr_{nj}\bigg). \tag{1}\end{align*}\end{document}

Once the total numbers of females and males of each type in generation (n + 1) are known (i.e., F_n+1,MR_n+1, and Mr_n+1), one deals with the mating phase. To determine the total number of couples of each type, we assume perfect fidelity mating (i.e., each individual mates with only one individual of the other sex provided that some of them are still available), and also that females choose their partner with a preference for R males. Hence, R males are chosen first, so the total number of R couples is determined by the minimum of the number of females and the number of R males: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}ZR_{n + 1} = \min \{ F_{n + 1}, MR_{n + 1} \}. \tag{2}\end{align*}\end{document}

Therefore, the number of females that do not mate with R males is max{0,F_n+1 -MR_n+1}. These females (if any) mate with r males, and the assumption of perfect fidelity implies that the number of r couples is \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Zr_{n + 1} = \min \{ \max \{ 0, F_{n + 1} - MR_{n + 1} \}, Mr_{n + 1} \}. \tag{3}\end{align*}\end{document}

Notice that the number of couples of each type in the (n + 1)-th generation is given deterministically once the total numbers of females and of males of each type in this generation are known.

From the definition of the model, the number of couples of each genotype in the next generation depends only on the present number of mating units and not on the number of ancestors that belonged to past generations. Furthermore, since each reproduction law remains the same over the generations, the transitions from one generation to another are homogeneous. The process {(ZR_n, Zr_n)}_n≥0 is therefore a homogeneous two-type Markov chain.

Some basic properties of this model are established in González et al. (2006). Among them, particularly worthy of note is that each genotype presents the dual behavior typical of branching processes: either it becomes extinct or the number of couples of this genotype eventually reaches arbitrarily large values. The latter event is known as the explosion or indefinite growth of this particular genotype. Consequently, the whole population also presents this duality. Although this property might seem unrealistic, it merely expresses what would be the ideal long-term evolution of a population when its development is not constrained by any external bound.

In González et al. (2009a), conditions for the survival/fixation of one genotype and the extinction/survival of the whole population are reviewed. These conditions depend on α and on the reproduction laws p^R and p^r through their mean values m_R and m_r, respectively. These values also determine the asymptotic behavior of the genotypes (as was proved in González et al., 2008). Since in practice these parameters are usually unknown, in order for these models to be applicable it is necessary to develop the estimation theory for the above parameters, including the reproduction laws. Then, knowing these estimators, predictions about the number of individuals and couples in future generations can also be established.

3. Maximum likelihood estimators with complete and incomplete data

In this section, we shall study the MLEs of the parameters α, p^R, and p^r. We shall also derive from them the MLEs for the reproduction means m_R and m_r. First we consider that the entire family tree up to some generation N is observed. This is the set of random vectors \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\{ ( FR_{ni}, MR_{ni} ), ( Fr_{nj}, Mr_{nj} ), i = 1, \ldots, ZR_n; j = 1, \ldots, Zr_n; n = 0, \ldots, N - 1 \}.\end{align*}\end{document}

From these random vectors, assuming that (ZR₀, Zr₀) is known and using Equations (1)–(3), one can obtain the sets \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \cal FM}_N = \{ Z R_0, Zr_{0}, F_n, MR_n, Mr_n, n = 1, \ldots, N \}$$\end{document} , containing the initial number of couples of each type and the total number of females and the total number of males of each type until generation N; and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \cal Z}_{N} = \{ ZR_n ( k ), k \in S^R, Zr_n ( l ), l \in S^r, n = 0, \ldots, N - 1 \},$$\end{document} where, with I _A denoting the indicator function of the set \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$A, ZR_n ( k ) = \sum\nolimits_{i = 1}^{ZR_n}I_{ \{ FR_{ni} + MR_{ni} = k \} }$$\end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$Zr_n ( l ) = \sum\nolimits_{j = 1}^{Zr_n}I_{ \{ Fr_{nj} + Mr_{nj} = l \} }$$\end{document} represent the total number of couples of each type that have generated, respectively, k and l individuals in the generation n.

Therefore, taking into account the binomial scheme and that mating units reproduce independently, it is not hard to obtain that the complete likelihood function based on the sample \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$( { \cal Z}_N, { \cal FM}_N )$$\end{document} is given by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}L ( p^R, p^r, \alpha \mid { \cal Z}_{N}, { \cal FM}_{N} ) \propto \prod_{n = 0}^{N - 1} \alpha^{F_{n + 1}} ( 1 - \alpha ) ^{ ( MR_{n + 1} + Mr_{n + 1} ) } \prod_{k \in S^R} ( p_k^R ) ^{ZR_n ( k ) } \prod_{l \in S^r} ( p_l^r ) ^{Zr_n ( l ) }. \tag{4}\end{align*}\end{document}

From this expression, and adapting some classical procedures of estimation theory in branching processes (see Supplementary Material, Theorem 1), one can obtain that the MLEs for α, p^R, and p^r based on the sample \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$( { \cal Z}_N, { \cal FM}_N )$$\end{document} are given by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\hat { \alpha } = \frac { \sum_ { n = 0 } ^ { N - 1 } F_ { n + 1 } } { \sum_ { n = 0 } ^ { N - 1 } ( F_ { n + 1 } + MR_ { n + 1 } + Mr_ { n + 1 } ) }, \tag { 5 } \end{align*}\end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\hat { p } _k^R = \frac { \sum_ { n = 0 } ^ { N - 1 } ZR_n ( k ) } { \sum_ { n = 0 } ^ { N - 1 } ZR_n }, k \in S^R, \ \ { \rm and } \ \ \hat { p } _l^r = \frac { \sum_ { n = 0 } ^ { N - 1 } Zr_n ( l ) } { \sum_ { n = 0 } ^ { N - 1 } Zr_n }, l \in S^r.\end{align*}\end{document}

The estimator for α is intuitively very reasonable, since it is obtained by means of the proportion of females among all observed individuals. The estimator for \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p_k^R$$\end{document} with \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$k \in S^R$$\end{document} (resp. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p_l^r$$\end{document} with \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$l \in S^R$$\end{document} ) is obtained as the total number of R couples (resp. r couples) that have generated k (resp. l) offspring as a fraction of the total number of R couples (resp. r couples).

From the estimators of p^R and p^r, one deduces that the MLEs for m_R and m_r based on the sample \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$( { \cal Z}_N, { \cal FM}_N )$$\end{document} are \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\hat { m } _R = \frac { \sum_ { n = 1 } ^ { N } ( FR_n + MR_n ) } { \sum_ { n = 0 } ^ { N - 1 } ZR_n } \quad { \rm and } \quad \hat { m } _r = \frac { \sum_ { n = 1 } ^ { N } ( Fr_n + Mr_n ) } { \sum_ { n = 0 } ^ { N - 1 } Zr_n },\end{align*}\end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$FR_n = \sum\nolimits_{i = 1}^{ZR_{n - 1}}FR_{n - 1i}$$\end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$Fr_n = \sum\nolimits_{j = 1}^{Zr_{n - 1}}Fr_{n - 1j}$$\end{document} , for all \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$n = 1, \ldots, N$$\end{document} , are the total numbers of females generated by each type of couple. Notice that for \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$n = 1, \ldots, N$$\end{document} , FR_n and Fr_n are derived from \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$( { \cal Z}_N, { \cal FM}_N )$$\end{document} , since \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$MR_n + FR_n = \sum\nolimits_{k \in S^R}kZR_{n - 1} ( k )$$\end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$Mr_n + Fr_n = \sum\nolimits_{k \in S^r}lZr_{n - 1} ( l )$$\end{document} .

All of these estimators verify some properties related to their asymptotic behavior. Specifically, on the nonextinction set, each estimator is strongly consistent, and, suitably normalized, converges in distribution to a standard normal distribution (see Supplementary Material).

Notice that the above estimators depend on the sample \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \cal Z}_N$$\end{document} which, in most real situations, is impossible to observe. Usually, only the total number of individuals of each type can be observed (recall that the Y-linked genes present different phenotypes). Thus, there arises an interesting estimation problem from assuming that only the sample \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \cal FM}_N$$\end{document} is observed. Since, from the definition of the model [Eq. (2) and (3)], one obtains ZR_n and Zr_n deterministically, knowing the total number of females and the total numbers of males of each type, one can insert the variables ZR_n and Zr_n into the sample \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \cal FM}_N$$\end{document} . Hence, writing \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$ZFM_n = \{ ZR_n, Zr_n, F_{n + 1}, MR_{n + 1}, Mr_{n + 1} \}, n = 0, \ldots, N - 1$$\end{document} , one is considering that the sample observed is \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \cal FM}_N = \{ ZFM_0, \ldots, ZFM_{N - 1} \} $$\end{document} .

Assuming that \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \cal Z}_N$$\end{document} is unknown and only the total number of individuals and of couples are observed, one is faced with an incomplete data estimation problem. In such a case, it seems appropriate to use an expectation-maximization (EM) algorithm (McLachlan and Krishnan, 2008), extensively used to deal with maximum likelihood calculations when there are missing or incomplete data. This algorithm is an iterative method that consists of two steps. In the first step (E step), the expectation of the complete log-likelihood is calculated using the distribution of the unobserved data. The second step (M step) consists of finding the values of the parameters that maximize the expectation that had been calculated in the E step. The E and M steps are repeated until convergence is attained. In our case, starting with initial values (p^R(0), p^r(0), α⁽⁰⁾), we shall obtain a sequence {(p^R(i), p^r(i), α⁽ⁱ⁾)}_i≥0 that is updated in each iteration of the method, as will be described in the following section.

4. Simulation study

The method presented in the previous section will now be applied using the R statistical computing language and environment (R Development Core Team, 2012) to estimate the parameters of a Y-linked BBP with preference using simulated data. To this end, we consider a process with the following parameters: the probability to be female is α = 0.4 and the reproduction laws of each type of couple are \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p^R = ( p_0^R, p_1^R, p_2^R ) = ( 0.0225, 0.2550, 0.7225 )$$\end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p^r = ( p_0^r, p_1^r, p_2^r, p_3^r ) = ( 0.0025, 0.0462, 0.3004, 0.6509 )$$\end{document} .

Note that we have chosen the sex-ratio to be less than a half since, in most populations, the sex-ratio is different from 0.5, and the analysis of Y-linked gene evolution turns out to be more interesting when α < 0.5 (González et al., 2006 and González et al., 2008). Also, the average number of individuals generated by each type of couple are m_R = 1.7 and m_r = 2.6, respectively, reflecting the possible difference between the reproductive capacity of mating units of each type that exists in nature.

For this model, we simulated 20 generations starting with (ZR₀, Zr₀) = (3,10). Table 1 lists the sample fm₂₀ formed by the total numbers of females and of males of each type obtained in these generations. The relatively small amount of sample information that this represents would make it difficult to determine at first sight anything about the future behavior of the Y-linked character on the basis of these observations.

Let us now apply the EM algorithm using the above sample, fm₂₀. To start the algorithm, we need the initial values (p^R(0), p^r(0), α⁽⁰⁾). Assuming the lack of information, we choose the values p^R(0) and p^r(0) according to uniform distributions of sizes 3 and 4, respectively. Thus, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p^{R ( 0 ) }_k = 1 / 3$$\end{document} , with k = 0, 1, 2 and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$p^{r ( 0 ) }_l = 1 / 4$$\end{document} , with l = 0, 1, 2, 3. The best option to initialize α⁽⁰⁾ is the MLE of α based on the entire family tree, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{ \alpha}$$\end{document} (see Eq. (5)—recall that \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{ \alpha}$$\end{document} only depends on the values recorded in fm₂₀). Therefore, as was indicated in the previous section, the sequence \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\{ \alpha^{ ( i ) } \} _{i \geq 0}$$\end{document} is constant, and of value \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{ \alpha}_{EM} = \hat{ \alpha}$$\end{document} — in the example, equal to 0.416.

We ran the EM algorithm and observed the sequence \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\{ ( p^{R ( i ) }, p^{r ( i ) }, \alpha^{ ( i ) } ) \} _{i \geq 0}, { \rm with} \ p^{R ( i ) } = \{ p^{R ( i ) }_k \} _{k = 0, 1, 2}, p^{r ( i ) } = \{ p^{r ( i ) }_l \} _{l = 0, 1, 2, 3}$$\end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\alpha^{ ( i ) } = \hat{ \alpha}_{EM}, i \geq 0$$\end{document} , to converge from iteration 500 onwards (with the difference between consecutive elements of the sequence being less than 10⁻⁵). The values obtained in the last iteration were taken to be the expectation-maximization ML estimates. A discrete sensitivity analysis applied to study the influence of the initial values (p^R(0), p^r(0), α⁽⁰⁾) on the convergence of the method showed that the procedure is stable with respect to the initial values, with there being no changes in the limit.

From \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\{ ( p^{R ( i ) }, p^{r ( i ) } ) \} _{i \geq 0}$$\end{document} , it is straightforward to obtain the sequence \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\{ ( m^{ ( i ) }_R, m^{ ( i ) }_r ) \} _{i \geq 0}$$\end{document} with the means of the distributions p^R(i) and p^r(i), respectively, in each iteration of the method. This last sequence converges to the expectation-maximization MLEs for m_R and m_r, denoted by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{m}^{EM}_R$$\end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{m}^{EM}_r$$\end{document} , respectively. From the values of the sequence \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\{ ( m^{ ( i ) }_R, m^{ ( i ) }_r ) \} _{i = 1, \ldots, 500}$$\end{document} , one obtains that they are quite stable from iteration 200 onwards, with the resulting expectation-maximization ML estimates of m_R and m_r based on fm₂₀ being \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{m}^{EM}_R = 1.724$$\end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{m}^{EM}_r = 2.605$$\end{document} , respectively.

In order to analyze the consistency of the expectation-maximization MLEs, we next applied the EM algorithm by varying the number of generations observed, that is, we applied the algorithm 20 times, taking the sample to be fm_N, with \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$N = 1, \ldots, 20$$\end{document} . Each of these times, we performed 500 iterations of the method and saved the estimates given in the last iteration, taking them to be the expectation-maximization ML estimates of the corresponding parameters. At the end of the process, we thus had a sequence \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$( \hat{ \alpha}_{EM, N}, \hat{m}^{EM}_{R, N}, \hat{m}^{EM}_{r, N} )$$\end{document} with \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$N = 1, \ldots, 20$$\end{document} . The three components of this sequence are plotted in Supplementary Figure S1. The more generations one has, the closer the estimate approaches the true value of the parameter (dashed line). Actually, under weak general conditions, the EM method leads to consistent estimates (McLachlan and Krishnan, 2008), as in the case of the usual MLEs based on complete data samples (see Supplementary Material).

To approximate the sampling distributions of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{p}^{R}_{EM, 20}, \hat{p}^{r}_{EM, 20}$$\end{document} , and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{ \alpha}_{EM, 20}$$\end{document} , we applied a bootstrap procedure, making use of the EM estimates obtained on the sample fm₂₀ (i.e., the values of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{p}^{R}_{EM, 20}, \hat{p}^{r}_{EM, 20}$$\end{document} , and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{ \alpha}_{EM, 20}$$\end{document} . These values were used as parameters to perform a Monte-Carlo simulation of 2,000 processes until generation 20. For each of these bootstrap samples, we applied the EM method thus obtaining bootstrap approximations to the sampling distributions of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{p}^{R}_{EM, 20}, \hat{p}^{r}_{EM, 20}$$\end{document} , and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{ \alpha}_{EM, 20}$$\end{document} . Obviously, from them it is straightforward to obtain a bootstrap sample of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{m}^{EM}_{R, 20}$$\end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{m}^{EM}_{r, 20}$$\end{document} . Supplementary Figure S2 illustrates the bootstrap approximation to these sampling distributions. One observes that the variability associated with the distribution of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{m}^{EM}_R$$\end{document} was greater than that of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{m}^{EM}_r$$\end{document} . This may have been because there were fewer R males recorded in each generation than r males.

An interesting question applied is to predict, on the basis of the observed data, whether or not the process will survive over time. Theoretically, it is known that the condition α < 0.5 and 1 < (1 − α)m_R < αm_r ensures that there exists a positive probability for both genotypes to grow without limit over time (González et al., 2008). From the bootstrap approximation to the sampling distributions of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{m}^{EM}_R, \hat{m}^{EM}_r$$\end{document} , and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{\alpha}_{EM}$$\end{document} , we calculated the proportion of samples in generation 20 that satisfied \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\hat{ \alpha}_{EM} < 0.5$$\end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$1 < ( 1 - \hat{ \alpha}_{EM} ) \hat{m}_R^{EM} < \hat{ \alpha}_{EM} \hat{m}_r^{EM}$$\end{document} , finding the approximate value 0.886. The high value of this calculated proportion is indicative that the theoretical condition might be satisfied. Indeed, the true values of the parameters do satisfy this condition, and therefore there exists a positive probability that both genotypes grow over the course of the generations.

Finally Supplementary Figure S3 illustrates the predictive distribution of the total numbers of females and of each type of male in the 21st generation. The predicted behavior is in keeping with the fact that there exists a positive probability that both genotypes survive over time.

5. Concluding remarks

In order to study the natural selection of Y-linked genes, the estimation of the main parameters of the Y-linked BBP with preference has been considered in a general nonparametric context. The model assumes males can be distinguished by certain genetic characteristics linked to the Y chromosome, characteristics which they either do or do not possess. The females choose their mates preferentially according to whether or not this characteristic is present. Firstly, we assumed that the entire family tree can be observed up to some generation and obtained the corresponding MLEs, studying their asymptotic properties—consistency and limiting normality. The procedure applied represented a methodological adaptation to the Y-linked models of some classical estimation theory procedures used in branching processes. Secondly, we considered the problem of estimating the main parameters of the model using only the sample information that is usually more plausibly observable in practice—that given simply by the number of females and of the two different types of male in each generation. We approached this problem as an incomplete data estimation problem, applying the expectation-maximization method that has proven very effective in solving it. How well this estimation procedure works was illustrated by means of a simulated example, in which we also showed the consistency of the estimates, obtained bootstrap approximations to their sampling distributions, and inferred the behavior of the process for future generations. This second procedure represents the principal objective of the present communication, allowing the use of these Y-linked models in applied problems under realistic assumptions.

We also showed that, when the only observable data are the total number of mating units of each genotype, the expectation-maximization method cannot be relied on to operate appropriately in estimating the probability of an individual being female, the reason being that the incomplete data likelihood function may not be unimodal. We concluded that it is necessary to observe, as a minimum, the numbers of females and of both male genotypes in each generation to guarantee the validity of the method.

A line for future research is the question of inferences for the two-sex branching model introduced in González et al. (2009b), in which it is considered that Y-linked genes are not expressed in the phenotype of males, so that females mate following a blind choice. In this framework, the total number of mating units of each type is not determined one-to-one from the total number of females and males of each type, and a random component underlies the mating process. Computationally, therefore, sampling the branching tree latent vector, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\cal Z}_N$$\end{document} , is more difficult and needs to be studied in some specific way. This complexity will probably lead to estimators whose sampling distributions will have large variances.