matpopmod.model¶
This module implements the class MPM
used to represent matrix
population models. It provides a flexible way to create matrix population
models and a high-level interface to study their mathematical properties and
biological descriptors (growth rate, stable distribution, reproductive
values, etc). The descriptors are computed when they are first needed, and
then stored as immutable attributes.
>>> m = MPM(A = "0.1 2; 0.4 0.8")
>>> m.primitive # is the projection matrix primitive?
True
>>> m.lmbd # asymptotic growth rate
1.4104686356149274
>>> m.w # stable distribution
array([0.60414407, 0.39585593])
>>> m.lmbd = 0 # trying to change a descriptor throws an error.
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
AttributeError: can't set attribute
In order to be computed, many descriptors require knowing which entries of the
projection matrix correspond to survival probabilities and which correspond to
fertilities, i.e. its decomposition into the survival matrix S and the
fertility matrix F. Trying to access such descriptors without having
specified the matrices S and F when creating the model will raise
NotAvailable
:
>>> m1 = MPM(A = "0.1 2; 0.4 0.8")
>>> m1.lmbd
1.4104686356149274
>>> m1.R0
Traceback (most recent call last):
...
matpopmod.utils.NotAvailable: 'A = S + F decomposition required.'
>>> m2 = MPM(S = "0.1 0; 0.4 0.8", F = "0 2; 0 0")
>>> m2.lmbd
1.4104686356149274
>>> m2.R0
4.4444444444444455
In what follows, we detail how to create MPM
objects and we
give the complete list of descriptors that can be computed with them.
- class matpopmod.model.MPM(A=None, S=None, F=None, fertilities=None, delimiter=None, metadata=None)¶
A matrix population model is defined by its projection matrix A. Some descriptors (such the growth rate, the stable distribution or the reproductive values) can be computed directly from this projection matrix. However, other ones (such as the net reproductive rate or the generation time) can be computed only if the A = S + F decomposition of the projection matrix into its survival and fertility components is known. Thus, there are two ways to create a matrix population model:
Using the projection matrix A only:
MPM(A = "0.1 2; 0.4 0.8")
If the model is created that way, trying to access descriptor that require knowing S and F will raise
NotAvailable
.By specifying any two of the three matrices A, S and F, for instance:
MPM(S = "0.1 0; 0.4 0.8", F = "0 2; 0 0")
If the three matrices are given, A will be silently ignored (that is, only S and F will be used, without checking whether A = S + F).
For convenience, instead of specifying two matrices it is also possible to provide A and the list of entries that correspond to reproduction. This is done with the fertility argument:
MPM(A = "0.1 2; 0.4 0.8", fertilities = [(0, 1)])
However this is less flexible, since this does not allow for overlapping entries between S and F. Note that the indexing of the entries starts at 0.
In addition to this, there are several ways to specify any of the matrices in the arguments:
Using a string, as we saw in the examples above. Rows have to be separated by a semicolon and columns by a space and/or a comma.
Using Python iterables, such as a list of lists:
MPM(A = [[0.1, 2], [0.4, 0.8]])
Valid arguments are the same as for NumPy’s
array()
function. In particular, NumPy arrays can be used.Reading the matrix from a text file. In that case, the name of the file containing the matrix must be passed as a
Path
from Python’spathlib
module.file_A = pathlib.Path("~/myfolder/myfile.txt") MPM(A = file_A)
Valid files contain one row of the matrix per line, with columns separated by whitespace (the character separating columns can be specified using the delimiter argument – use
","
to load CSV files). Lines starting with#
,%
or//
are ignored.
Finally, each model has a
metadata
attribute. It is simply a dictionary that can be used to store relevant information about the model.Complete list of available descriptors
Unless a specific reference is given, the formulas used are classic and can be found in textbooks such as [Casw00] or [KeCa05].
- Descriptors of the projection matrix
- Descriptors of the survival and fertility matrices
- Functions of the initial population
I. Descriptors of the projection matrix
I.1. Structural properties
I.1.a. The matrix A
- A¶
The projection matrix \(\mathbf{A} = (a_{ij})\), where \(a_{ij}\) is the expected per-capita contribution of individuals of class j to the abundance of class i at the next time-step.
- dim¶
The number of classes of the model.
I.1.b. Irreducibility and related quantities
- irreducible¶
Whether the projection matrix A is irreducible, meaning for any two classes i and j, there exists a directed path from i to j in the life-cycle graph. Equivalently, a non-negative matrix A is irreducible if and only if for all i, j there exists k such that the (i, j)-th entry of \(\mathbf{A}^k\) is positive.
See
is_irreducible()
from the modulemathtools
for a function that can be used on any NumPy array and for details about the implementation.
- irreducible_components¶
The decomposition of the projection matrix A into irreducible components, as a couple (submatrices, classes) where submatrix[k] is the submatrix of A that corresponds to the k-th irreducible component and classes[k][i] is the class of A that corresponds to the i-th class of that component. Example:
>>> m = MPM(A = "1 0 0 0; 2 0 3 0; 0 4 0 0; 0 0 5 0") >>> m.A array([[1, 0, 0, 0], [2, 0, 3, 0], [0, 4, 0, 0], [0, 0, 5, 0]]) >>> m.irreducible_components ([array([[1]]), array([[0, 3], [4, 0]])], [[0], [1, 2]])
Note that classes are indexed from 0 to n - 1, where n = dim(A), and that classes that do not belong to any directed cycle (as is the case above with class
3
) are not part of any irreducible component.See
strongly_connected_components()
from the modulemathtools
for a function that can be used on any NumPy array and for details about the implementation. See alsonormal_form()
for an alternative representation of the decomposition that retains the transitions between the classes.
- normal_form¶
A normal form of the projection matrix A is obtained by reordering its classes so as to get a projection matrix of the form
\[\begin{split}\begin{pmatrix} \mathbf{A}_{11} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{A}_{21} & \mathbf{A}_{22} & \ddots & \vdots \\ \vdots & \vdots & \ddots & \mathbf{0} \\ \mathbf{A}_{2m} & \mathbf{A}_{2m} & \cdots & \mathbf{A}_{mm} \end{pmatrix} , \end{split}\]where each of the diagonal blocks \(\mathbf{A}_{11}, \ldots, \mathbf{A}_{mm}\) is either an irreducible matrix or a null matrix – see e.g. Chapter XIII, §4 of [Gant59].
The new matrix is returned with a NumPy array permutation representing the reordering of the classes: permutation[i] is the index of the row/column of A that corresponds to the i-th row/column in the normal form. Indexing starts at 0. Example:
>>> m = MPM(A = "1 0 2 0; 3 0 0 4; 5 0 0 0; 0 6 7 0") >>> m.A array([[1, 0, 2, 0], [3, 0, 0, 4], [5, 0, 0, 0], [0, 6, 7, 0]]) >>> m.normal_form (array([[1, 2, 0, 0], [5, 0, 0, 0], [3, 0, 0, 4], [0, 7, 6, 0]]), array([0, 2, 1, 3]))
Note that the normal form is not unique, as we can reorder the classes inside each irreducible component to get a different one (if some of the blocks under the diagonal are null, we may also be able to permutate the irreducible blocks). This implementation tries to preserve the original order of the classes. Also note that several authors (e.g, [Varg00] and [HoJo13]) define the normal form to be block upper-triangular, rather than block lower-triangular. However, this is at odds with the usage in matrix population models, where projection matrices are typically lower-triangular (e.g, Leslie matrices).
See
irreducible_components
for a function that extracts the submatrices corresponding to the irreducible components.
- quasi_irreducible¶
Whether the projection matrix A is quasi-irreducible, i.e. whether it has only one irreducible component whose dominant eigenvalue is equal to λ.
Unlike irreducibility, quasi-irreducibility is not a standard notion. However, it is useful to treat models with post-reproductive classes or metapopulation models with uni-directional migration: indeed, the vast majority of matrix population are quasi-irreducible, and many of the results for irreducible matrices are easily adapted to quasi-irreducible ones. See
quasi_primitive
for more information.
I.1.c. Aperiodicity and related quantities
- periods¶
The vector of periods of the projection matrix A. The period of the class i is the greatest common divisor of the lengths of the cycles goings through the corresponding vertex in the life-cycle graph (if there is no such cycle, then it is
None
). See alsoindex_of_imprimitivity
for the greatest common divisor of the periods the classes.See
periods()
from themathtools
module for a function that can be used on any NumPy array and for details about the implementation.
- index_of_imprimitivity¶
The index of imprimitivity of the projection matrix A, also known as its period.
When A is irreducible, the coefficient of imprimitivity is both the greatest common divisor of the lengths of the cycles in the life-cycle graph and the number of eigenvalues whose modulus is equal to the asymptotic growth rate λ. It corresponds to the period of the asymptotic oscillations of the population structure.
When A is reducible, the coefficient of imprimitivity is defined as the least common multiple of the periods of its dominant components (see the documentation of the function
index_of_imprimitivity()
from the modulemathtools
for details). It corresponds to the longest possible period of the asymptotic oscillations: if we start from any¹ population vector with individuals in every classes, there will be asymptotic oscillations whose period is the coefficient of imprimitivity; but some initial population vectors can produce asymptotic oscillations with a period that is a divisor of the index of imprimitivity.When the index of imprimitivity is equal to 1, the population structure can never exhibit asymptotic oscillations and A is said to be
aperiodic
.
- aperiodic¶
Whether the projection matrix A is aperiodic, i.e. whether the index of imprimitivity is equal to 1. Aperiodicity is a necessary condition for the convergence of the population structure: otherwise, it can exhibit asymptotic oscillations (see
index_of_imprimitivity
for details).
I.1.d. Primitivity and quasi-primitivity
- primitive¶
Whether the projection matrix A is primitive, meaning that there exists an integer k such that \(\mathbf{A}^k\) is positive. Equivalently, a non-negative matrix is primitive if and only if it is irreducible and aperiodic.
Primitivity is a central notion of the Perron-Frobenius theory. It ensures that, as \(t \to \infty\),
\[\frac{\mathbf{A}^t}{\lambda^t} \to \mathbf{w v} , \]with λ the dominant eigenvalue of A, v and w the corresponding left and right eigenvectors (scaled so that vw = 1), and wv the matrix \((w_i v_j)\). As a result, for any¹ initial population vector \(\mathbf{n}(0)\), \(\mathbf{n}(t) / \lambda^t \to c \mathbf{w}\), where \(c = \mathbf{vn}(0)\). Note that these conclusions also hold if A is
quasi_primitive
.See
is_primitive()
from the modulemathtools
for a function that can be used on any NumPy array and for details about the implementation.
- quasi_primitive¶
Whether the projection matrix A is quasi-primitive, i.e. quasi-irreducible and aperiodic. Equivalently, a non-negative matrix is quasi-primitive if and only if it has a single eigenvalue of maximum modulus and this eigenvalue has algebraic multiplicity 1.
Quasi-primitivity is not a standard notion, but it is a useful relaxation of primitivity when working with matrix population models. Indeed, like primitivity it ensures that \(\mathbf{A}^t / \lambda^t \to \mathbf{w v}\) as \(t \to \infty\). The only difference with the primitive case is that the dominant eigenvectors v and w are allowed to have entries that are equal to zero (corresponding to “sink” classes for v and to “source” ones for w).
Further generalizations are possible – see e.g. [Roth81] – but they are more technical and typically not useful for matrix population models, where the structure of matrices is subject to biological constraints. See Section 9.3 of [Hogb06] for an overview of Perron-Frobenius theory for reducible matrices.
I.2. Eigen-elements and related quantities
I.2.a. Asymptotic dynamics
- lmbd¶
The asymptotic growth rate λ. For large times t, the population size will grow (or decrease, if \(\lambda < 1\)) like \(\lambda^t\) times a constant that depends on the initial composition of the population.
Mathematically, λ is the spectral radius of the projection matrix A, i.e. the maximum of the modulus of its eigenvalues. Note that the spectral radius of a non-negative matrix is always an eigenvalue of that matrix (see e.g. Section 8.3 of [Meye00]). However, there can be several eigenvalues of maximum modulus; in fact, λ is the only eigenvalue of maximum modulus if and only if A is aperiodic – see Fact 3.(a) in Section 9.3 of [Hogb06]. When λ is the only eigenvalue of maximum modulus, we refer to it as the dominant eigenvalue.
- w¶
The stable class distribution vector w. If the projection matrix A is primitive, then for any¹ initial population structure \(\mathbf{n}(0)\) the population vector n will satisfy
\[\frac{\mathbf{n}(t)\;}{\|\mathbf{n}(t)\|_1} \to \mathbf{w} \quad\text{as}\quad t \to \infty , \]where \(\|\mathbf{n}(t)\|_1 = \sum_i |n_i(t)|\) denotes the total population size. This convergence will also hold if A is quasi-primitive and we start with at least one individual in the dominant component (or one of its sources); see the documentation of
quasi_primitive
for details.Mathematically, w is the right eigenvector associated to the asymptotic growth rate λ. If A is not quasi-primitive, then w is not well-defined (the asymptotic population structure may either depend on the initial one, due to reducibility; or not converge, due to periodicity). In that case, a vector of
nan
will be returned.
- v¶
The reproductive values vector v, scaled so that \(\mathbf{vw} = \sum_i v_i w_i = 1\). The term “reproductive value” is justified by the fact that, for quasi-primitive projection matrices, as \(t \to \infty\) the population size is asymptotically equivalent to \(c\lambda^t\), where
\[c \;=\; \mathbf{v} \mathbf{n}(0) \;=\; \sum_i v_i n_i(0)\,.\]Thus, \(v_i\) is the relative contribution of an individual of class i to the composition of the population in the distant future.
Mathematically, v is the left eigenvector associated to the asymptotic growth rate λ. If A is not quasi-primitive, then v is not well-defined. In that case, a vector of
nan
will be returned.Note that even though NumPy does not make a distinction between row vectors and column vectors (both correspond to 1-D arrays), in our mathematical expressions we consider v to be a row vector, This explains why we write vw for the scalar \(\sum_i v_i w_i\) and wv for the matrix \((v_i w_j)\).
I.2.b. Second order dynamics
- damping_ratio¶
The damping ratio of the projection matrix A, defined as the growth rate λ divided by the maximum of the modulus of the subdominant eigenvalues, that is,
\[\frac{\lambda}{\tau} \quad\text{where}\quad \tau = \max\big\{|\mu| : \mu \text{ is an eigenvalue with } |\mu| < \lambda\big\},\]where \(\max \varnothing = 0\) – i.e. when there is no eigenvalue \(\mu\) such that \(0 < |\mu| < \lambda\), the damping ratio is infinite and
inf
is returned.The damping ratio corresponds to the rate of convergence to the steady regime. Indeed, for primitive projection matrices, for any \(r > \tau\),
\[\mathbf{A}^t \;=\; \lambda^t \mathbf{wv} \;+\; O(r^t)\]see e.g. Theorem 1.2 in [Sene06] for a slightly more precise statement.
Note that the name “damping ratio” is now well-established in the literature on matrix population models (see e.g. [Casw00]), but that it is slightly unfortunate since \(\lambda / \tau\) does not correspond to the notion of damping ratio in classical physics, which also depends on the period of the second-order oscillations (see
second_order_period
). The damping ratio from matrix population models is actually a geometric rate of convergence; in all likeliness the confusion comes from the interpretation of the word ratio (which, in physics, refers to the ratio of the actual damping to the critical one).
- second_order_period¶
The period of the second-order oscillations,
\[\mathcal{P} \;=\; \left|\frac{2 \pi}{\arg \lambda_2}\right|, \]where \(\lambda_2\) and its complex conjugate \(\lambda_2^*\) are the only eigenvalues with largest modulus strictly less than λ (note that these two eigenvalues give the same \(\mathcal{P}\) so it does not matter which one we choose in the definition). When \(\lambda_2 = \lambda_2^* \geq 0\), then there are no oscillations associated to the eigenvalue of second-largest modulus and
nan
is returned. Note that there could still be oscillations associated to eigenvalues of smaller modulus.In the unlikely case where there are more than two eigenvalues of second-largest modulus, then the corresponding oscillations are periodic if and only if the corresponding periods \(\mathcal{P}_1, \ldots, \mathcal{P}_q\) are commensurable, i.e. if there exists positive integers \(k_1, \ldots, k_q\) such that \(k_1 \mathcal{P}_1 = \cdots = k_q \mathcal{P}_q\). The period is then the smallest possible such common value of \(k_i \mathcal{P}_i\). Of course, it does not make sense to ask if real numbers are commensurable when they are represented by floating point numbers. Thus, in that case we test if there are small integers such that the relation above holds up to numerical precision; otherwise a
CannotCompute
error is raised.The values that we compute may differ from those calculated by other software. Whenever we identified discrepancies, they were due to the fact that other software chose their arguments in \([0, 2\pi]\) or failed to take into account the possibility of more than two eigenvalues of second largest modulus, which are errors of implementation.
I.2.c. Other eigen-elements
- eigenvalues¶
The eigenvalues of the projection matrix A, in decreasing order of their modulus and repeated according to their algebraic multiplicities.
See
right_eigenvectors
andleft_eigenvectors
for the corresponding eigenvectors, andeigen_elements()
from the modulemathtools
for details about the implementation.
- right_eigenvectors¶
Right eigenvectors of the projection matrix A, that is, each \(\mathbf{w}^{(k)}\) = right_eigenvectors[k] is a solution of
\[\mathbf{A}\mathbf{w}^{(k)} = \lambda^{(k)}\,\mathbf{w}^{(k)},\]where \(\lambda^{(k)}\) = eigenvalues[k]. Note that there are always some small errors (typically of the order of the machine epsilon, although there are no strict guarantees on that) in the obtention of those solutions:
>>> m = MPM(A = "0.1 2; 0.4 0.8") >>> lmbd0, w0 = m.eigenvalues[0], m.right_eigenvectors[0] >>> m.A @ w0 array([0.85212626, 0.55834237]) >>> lmbd0 * w0 array([0.85212626, 0.55834237]) >>> m.A @ w0 - lmbd0 * w0 array([-1.11022302e-16, 0.00000000e+00])
Each vector is scaled so that
\[\|\mathbf{w}^{(k)}\|_1 \;:=\; \sum_i |w_i^{(k)}| \;=\; 1, \]and is chosen to be non-negative when possible.
Importantly, note that unless A is diagonalizable, some of its eigenvectors will be linearly dependent (that includes the possibility of repeated or null eigenvectors). Thus, although
right_eigenvectors
will always return exactlydim
eigenvectors (because the eigenvalues are repeated according to their algebraic multiplicities ineigenvalues
), this should not be interpreted as meaning that A “hasdim
eigenvectors” (an abusive formulation that implicitly refers to linearly independent eigenvectors). If needed, functions to compute the geometric multiplicities and to test linear independence are provided in the modulemathtools
.
- left_eigenvectors¶
Left eigenvectors of the projection matrix A, that is, each \(\mathbf{v}^{(k)}\) = left_eigenvectors[k] is a solution of
\[\mathbf{v}^{(k)} \mathbf{A} = \lambda^{(k)}\,\mathbf{v}^{(k)},\]where \(\lambda^{(k)}\) = eigenvalues[k].
Note that some authors define left eigenvectors as column vectors, and will sometimes introduce them as “the right eigenvectors of the conjugate transpose matrix \(\mathbf{A}^*\)”. In the context of matrices with real entries, all definitions are equivalent up to transposition, since \(\mathbf{A}^* = \mathbf{A}^\top\) and \((\mathbf{A}^\top \mathbf{x})^\top = \mathbf{x A}\).
Although in our mathematical expressions we consider left eigenvectors as row vectors, NumPy does not make a distinction between row vectors and column vectors and all vectors in this library are represented using NumPy 1-D arrays.
The left eigenvectors returned here are all scaled so that
\[\|\mathbf{v}^{(k)}\|_1 \;:=\; \sum_i |v_i^{(k)}| \;=\; 1, \]however note that a different scaling is used for the vector of reproductive values
v
.Finally, the comments on the accuracy of the determination of the eigenvectors made in the documentation of
right_eigenvectors
also apply here.
I.2.d. Derivatives (sensitivities and elasticities)
- sensitivities¶
The matrix of sensitivities of the growth rate λ to the entries of the projection matrix A, that is,
\[s_\lambda(a_{ij}) \;=\; \frac{\partial \lambda}{\partial a_{ij}} \,.\]If A is primitive, then
\[s_\lambda(a_{ij}) \;=\; v_i w_j, \]where the vector of reproductive values v and the stable class distribution w are such that vw = 1. In fact, this formula is valid whenever A is quasi-irreducible; however for this v and w have to be defined as the only non-negative eigenvectors associated to λ, which cannot necessarily be interpreted as the reproductive values and stable distribution.
If A is not quasi-irreducible, then λ is not a differentiable function of its entries. In that case, a warning is issued and a matrix of
nan
is returned.The standard reference for the formula \(s_\lambda(a_{ij}) = v_i w_j\) is [Casw78], which was the first time it appeared in the biological literature. However, this paper does not contain a mathematical proof. The first proof seems to have been published in the earlier paper [Vahr76]. As of June 2021, this paper was cited only 6 times.
See Chapter 3 of [KiNe19] for a more general yet very accessible treatment.
- elasticities¶
The matrix of elasticities of the growth rate λ to the entries of the projection matrix A, that is,
\[e_\lambda(a_{ij}) \;=\; \frac{a_{ij}}{\lambda} \frac{\partial \lambda}{\partial a_{ij}} \;=\; \frac{\partial \log \lambda}{\partial \log a_{ij}} \,.\]Contrary to the sensitivities, which measure the absolute change in λ in response to an additive perturbation of the entries of A, the elasticities measure the relative change in λ in response to a multiplicative perturbation of the entries of A.
The elasticities also have a simple interpretation as the asymptotic frequencies of the transitions of the life cycle that we encounter as we follow the lineage of a focal individual back in time; see [BALM17] for more on this.
If A is primitive, then
\[e_\lambda(a_{ij}) \;=\; v_i\, a_{ij}\, w_j \,/\, \lambda,\]where the vector of reproductive values v and the stable class distribution w are such that vw = 1. See the documentation of
sensitivities
for details about the validity of this formula and for references.
I.3. The genealogical Markov chain
I.3.a Transition matrix and stationary distribution
- P¶
The transition of matrix of the genealogical Markov chain, \(\mathbf{P} = (p_{ij})\) where
\[p_{ij} = \frac{a_{ij} w_j}{\lambda w_i},\]with λ the asymptotic growth rate and w the stable distribution.
The genealogical Markov chain describes the sequence of classes that we encounter as we go “up” the genealogy of the population by following the lineage of a focal individual backwards in time: \(p_{ij}\) is the probability that an individual in class i at time t was in class j at time t-1 (or had its mother in class j, if the individual was not born yet at time t-1).
The genealogical Markov chain was first introduced in [Deme74] to define population entropy, but has recently found its use in a variety of problems. See [BALM17] for a complete introduction.
See also
Ps
andPf
for the survival and fertility components of P, andpi
for its stationary distribution.
- pi¶
The vector of class reproductive values \(\pi_i = v_i w_i\). This is also the stationary distribution of the Markov chains corresponding to
P
, see [BALM17].
I.3.b Entropy rate and other statistics
- entropy_rate¶
The entropy rate of the genealogical Markov chain
P
, that is\[H \;=\; -\sum_{i, j} \pi_{i} p_{ij} \log_2 p_{ij} .\]This measures the information density of the trajectories of the genealogical chain, and can therefore be seen as a measure of the complexity of the life cycle (although, when comparing life cycles with different number of classes, remember that H takes values in \([0, \log_2 n]\), where n is the number of classes). Mathematically, if we let \(X_1, X_2, \ldots\) be a trajectory of the genealogical chain and \(H(X_1, \ldots, X_t)\) denote the joint entropy of \(X_1, \ldots, X_t\), then
\[H \;=\; \lim_{t \to \infty} \frac{1}{t} H(X_1, \ldots, X_t)\]The first use of the entropy rate of P in matrix population models is due to Demetrius in [Deme74], where it is refered to as the population entropy. The same terminology is used in [Tulj82], but Demetrius later changed it²; nowadays, \(H\) is often referred to as the evolutionary entropy, the term population entropy being reserved to the quantity \(S = HT_a\), where \(T_a\) is the generation time.
- hitting_times¶
(coming in version 0.2.0)
- kemeny_constant¶
(coming in version 0.2.0)
II. Descriptors of the survival and fertility matrices
II.1. The A = S + F decomposition
II.1.a Matrices
- split¶
Whether the A = S + F decomposition of the projection matrix into survival probabilities vs fertilities is available. Note that many descriptors of matrix population models require this decomposition in order to be computed.
- S¶
The survival matrix \(\mathbf{S} = (s_{ij})\), where \(s_{ij}\) is the probability that an individual of class j survives and goes to class i at the next time-step. See also
F
andsplit
.
- F¶
The fertility matrix \(\mathbf{F} = (f_{ij})\), where \(f_{ij}\) is the expected number of offspring of class i of an individual in class j. See also
S
andsplit
.
- mixed_transitions¶
Whether the model has “mixed transitions”, i.e. transitions in the life-cycle that can correspond to both survival and reproduction. In other words, whether there are overlapping entries between the matrices S and F.
- fundamental_matrix¶
The fundamental matrix associated to S, that is,
\[\mathbf{N} = (\mathbf{I} - \mathbf{S})^{-1} \,.\]The (i, j)-th entry of N is the expected number of times that an individual from class j will visit the class i during its remaining lifetime.
See e.g. Chapter III of [KeSn76] for more on the mathematical properties of N.
- Ps¶
Survival component of the genealogical matrix
P
: \(p_s(i, j)\) is the probability that an individual in class i at time t was alive in class j at time t-1.
- Pf¶
Fertility component of the genealogical matrix
P
: \(p_f(i, j)\) is the probability that an individual in class i has just been born to a mother from class j.
- fundamental_matrix_Ps¶
The fundamental matrix associated to \(\mathbf{P}_s\), that is
\[(\mathbf{I} - \mathbf{P}_s)^{-1} \,.\]The (i, j)-th entry of this matrix is the expected number of times that an individual from class i has been in class j in the past.
II.1.b Leslie and Usher models
- leslie¶
Whether the model is a Leslie model, i.e. an age-classified model with a projection matrix is of the form
\[\begin{split}\mathbf{A} \;= \; \begin{pmatrix} \; f_1 & f_2 & \cdots & f_{n-1} & f_n \\[1ex] \; s_1 & 0 & \cdots & 0 & 0 \\ \; 0 & s_2 & \ddots & & \vdots \\ \; \vdots & & \ddots & 0 & 0 \\ \; 0 & 0 & \cdots & s_{n-1} & s_{n} \\ \end{pmatrix}\end{split}\]where \(s_{i}\) is the probability of surviving from class i to class i + 1 and \(f_i\) is the fertility of class i. The last class regroups all individuals older than a certain age.
Note that here the classes have to be labeled by increasing age, otherwise the model will not be considered a Leslie model, even if the life-cycle graph is isomorphic to a Leslie model. You can therefore rely on the fact that any model such that
leslie == True
is in the exact form above. To test whether a model is a “relabeled” Leslie model, seerelabeled_leslie
.Also note that for a model to be considered a Leslie model, the survival and fertility matrices have to have been provided when creating the model. Otherwise there simply is no way to know whether this is really a Leslie model and this could lead to incorrect results.
>>> m0 = MPM(A = "0.5 1; 0.3 0.8") >>> m0.leslie Traceback (most recent call last): ... matpopmod.utils.NotAvailable: A = S + F decomposition required. >>> m1 = MPM(A = "0.5 1; 0.3 0.8", F = "0.5 1; 0 0") >>> m1.leslie True >>> m1.R0 2.0 >>> m2 = MPM(A = "0.5 1; 0.3 0.8", F = "0 1; 0 0") >>> m2.leslie False >>> m2.R0 # Different from the corresponding Leslie model m1. 3.0000000000000004
- usher¶
Whether the model is an Usher model (also sometimes called “Lefkovitch model” or “standard size-based model”), i.e. has a projection matrix of the form
\[\begin{split}\mathbf{A} \;= \; \begin{pmatrix} \; f_1 + r_1 & f_2 & \cdots & f_{n-1} & f_n \\[1ex] \; g_1 & r_2 & \cdots & 0 & 0 \\ \; 0 & g_2 & \ddots & & \vdots \\ \; \vdots & & \ddots & r_{n-1} & 0 \\ \; 0 & 0 & \cdots & g_{n-1} & r_{n} \\ \end{pmatrix}\end{split}\]where \(f_i\) is the fertility of class i, \(g_i\) the probability of “growing” from class i to class i + 1, and \(r_i\) the probability of surviving and remaining in class i. Of course the term “growing” is an image, since we have no way to know whether the classes are actually based on size.
Note that, mathematically, Usher models are a generalization of Leslie models. Thus, every model such that
usher == True
is also such thatleslie == True
.See also
relabeled_usher
for whether the model can be obtained by reordering the classes of an Usher model.
- relabeled_leslie¶
(coming in version 0.2.0)
- relabeled_usher¶
(coming in version 0.2.0)
II.1.c Newborn and reproductive classes
- newborn_classes¶
The vector whose i-th entry is equal to:
\(0\) if no individual in class i is a newborn, i.e. if \(f_{ij} = 0\) for all j.
\(1\) if every individual in class i is a newborn, i.e. if \(s_{ij} = 0\) for all j (and there exists j such that \(f_{ij} > 0\)).
nan
otherwise. In that case, class i may or may not be considered newborn class, since some but not all of its individuals are newborns.
If you want to use a different convention, you can use NumPy’s
nan_to_num()
function:>>> m = MPM(S = "0.1 0; 0.4 0.8", F = "0 2; 0 0") >>> m.newborn_classes array([nan, 0.]) >>> numpy.nan_to_num(m.newborn_classes, nan = 0) array([0., 0.]) >>> numpy.nan_to_num(m.newborn_classes, nan = 1) array([1., 0.])
(this requires NumPy ≥1.17)
See also
proportion_newborns
for the proportion of newborns in each stage, andnu
andclass_of_birth
for related quantities.
- unique_newborn_class¶
Whether the model has a single, well-defined newborn class, i.e. whether the two following conditions are met: (1) all individuals are born in the same class and (2) all individuals in that class are newborns.
- reproductive_classes¶
The vector whose i-th entry is equal to \(1\) if i is a reproductive class and \(0\) otherwise.
- postreproductive_classes¶
The vector whose i-th entry is equal to \(1\) if i is a post-reproductive class and \(0\) otherwise.
Note that if i is post-reproductive then its reproductive value \(v_i\) is equal to zero, but that the converse is not true: in a meta-population model, the classes of a sink patch can have \(v_i = 0\) without being post-reproductive.
- proba_maturity¶
(coming in version 0.2.0)
II.1.d Steady state distributions
- proportion_newborns¶
The vector whose i-th entry is the proportion of the individuals of class i that have just been born, assuming the population is at its stable class distribution. That is,
\[\sum_{j} \frac{f_{ij} w_j}{\lambda w_i} \, .\]If the proportion of newborns in the class is not defined (e.g, because it depends on the composition of the population and there is no stable class distribution), the corresponding entry is set to
nan
.
- nu¶
The vector \(\boldsymbol\nu = (\nu_i)\) giving the fraction of offspring produced in a given year that are born in class i, assuming the population is at its stable class distribution. That is,
\[\boldsymbol\nu \;=\; \frac{\mathbf{Fw}}{\|\mathbf{Fw}\|_1} \,, \]where \(\|\cdot\|_1\) denotes the \(L_1\)-norm (sum of the absolute values of the components). If there is no stable class distribution but all individuals are born in class k, then the vector \(\boldsymbol\nu = (\nu_i)\), where \(\nu_i = 1\) if and only if i = k, is returned.
- class_of_birth¶
The matrix \(\mathbf{B} = (b_{ij})\) whose (i, j)-th entry gives the probability that an individual currently in class j was born in class i, that is,
\[\begin{split}b_{ij} \;&=\; \mathbb{P}(\text{born in } i \mid \text{in } j \text{ at }t = 0)\\ \;&=\; \sum_{t \geq 0} \mathbb{P}\big(\text{born at } -t \,\big|\, \text{in } i \text{ at }-t\big) \, \mathbb{P}\big(\text{in } i \text{ at } -t \,\big|\, \text{in } j \text{ at }t = 0\big) \\ \;&=\; \sum_{t \geq 0} \sum_k p_f(i, k) \, p_s^{(t)}(j, i) \,, \end{split}\]with \(\mathbf{P}_{\!f} = (p_f(i, j))\) and \(\mathbf{P}_{\!s}^t = (p_s^{(t)}(i, j))\) the genealogical matrices. In matrix form,
\[\mathbf{B} \;=\; (\mathbf{I} - \mathbf{P}_{\!s}^\top)^{-1} \otimes\mathbf{P}_{\!f} \, \mathbf{E}\,\]where E is a square matrix of ones and \(\otimes\) denotes the entry-wise product.
- class_of_death¶
(coming in version 0.2.0)
II.1.d Fertility and survival excesses
- fertility_excess¶
The fertility excess (non-standard terminology), defined as the constant c > 0 such that \(\rho(\mathbf{S} + c^{-1} \mathbf{F}) = 1\), where \(\rho\) denotes the spectral radius (i.e. the asymptotic growth rate λ, in the case of a primitive matrix). In other words, the fertility excess is the factor by which all fertilities must be rescaled to get a stationary projection matrix.
- survival_excess¶
The survival excess (non-standard terminology), defined as the constant c > 0 such that \(\rho(c^{-1}\mathbf{S} + \mathbf{F}) = 1\), where \(\rho\) denotes the spectral radius (i.e. the asymptotic growth rate λ, in the case of a primitive matrix). In other words, the survival excess, if it exists, is the factor by which all fertilities must be rescaled to get a stationary projection matrix.
When λ > 1, if the spectral radius of F is greater than 1, then c is not defined and
NaN
will be returned. When λ < 1, the matrix \(c^{-1}\mathbf{S}\) can have columns that sum to more than 1. This means that the original model is such that it is impossible to get a sationary projection matrix by increasing survival by a common factor, even though c is mathematically well-defined.
- theta¶
(coming in version 0.2.0)
II.2. Biological descriptors
II.2.a Reproductive output
- G¶
The next generation matrix
\[\mathbf{G} = \mathbf{F} (\mathbf{I} - \mathbf{S})^{-1} \,,\]whose (i, j)-th entry is the expected number of offspring of class i that an individual of class j will produce in its remaining lifetime. See
R0
for the dominant eigenvalue of this matrix, a commonly used measure of the net reproductive rate.
- w_G¶
Cohort stable distribution (dominant right-eigenvector of
G
, normalized so that its entries sum to 1).
- v_G¶
Cohort reproductive values (dominant left-eigenvector of
G
, normalized so that \(v_c \cdot w_c = 1\)).
- total_reproductive_output¶
The matrix \(\mathbf{R} = (r_{ij})\) such that \(r_{ij}\) is the expected number of offspring of class i that an individual currently in class j is expected to produce over the course of its life. This includes offspring produced before the present:
\[\begin{split} \mathbf{R} \;&=\; \underbrace{\sum_{t \geq 1} \mathbf{F} (\mathbf{P}_{\!s}^\top)^t}_{\text{before}} \;+\; \underbrace{\sum_{t \geq 0} \mathbf{F} \mathbf{S}^t}_{\text{after}} \\ \;&=\; \mathbf{F} (\mathbf{I} - \mathbf{P}_{\!s}^\top)^{-1} \mathbf{P}_{\!s}^\top \;+\; \mathbf{G}\,.\end{split}\]where \(\mathbf{P}_{\!s}\) is the survival component of the genealogical matrix
P
.
- proba_repro¶
The vector whose i-th entry is the probability that an individual currently in class i will reproduce during its remaining lifetime. Use
nu @ proba_repro
to get the proportion of newborns that reproduce before dying in the stable population.See
life_expectancy_repro
for details on how this is calculated.
II.2.b Net reproductive rate
- R0¶
The measure of the net reproductive rate (also called net reproduction number) introduced by Cushing and Zhou in [CuZh94], namely the dominant eigenvalue of the next generation matrix \(\mathbf{G} = \mathbf{F}(\mathbf{I} - \mathbf{S})^{-1}\).
Although this is a relevant measure of the net reproductive rate, one should note that if there are several newborn classes then this does not necessarily correspond to the expected number of offspring that newborn will produce over the course of its life, in the usual sense. This is because there are two, non-equivalent ways to choose a newborn:
According to the dominant eigenvector of G, as done here.
Uniformly among the individuals born in a given year, as done in
cohort_R0
.
These two measures can differ significantly:
>>> mpm.examples.thalia_democratica.R0 4.001730883164397 >>> mpm.examples.thalia_democratica.cohort_R0 5.222416673531082
- cohort_R0¶
The net reproductive rate, as the average number of offspring that a newborn is expected to produce over its lifetime in the stable population. This is given by
\[R_0^* \;=\; \|\mathbf{G} \boldsymbol \nu \|_1\]with G the next-generation matrix and \(\boldsymbol\nu = (\nu_i)\), where \(\nu_i\) is the fraction of newborns that are born in class i when the population is at its stable distribution (see
nu
).
II.2.c Generation time
- T_a¶
The generation time, as the mean age of mothers in the stable population (that is, when there is a birth in the stable population we record the age of the mother).
As suggested by [CoEl92], all births are not counted equally: instead, they are weighted by the reproductive value of the newborns. This yields a mathematically more tractable measure. In particular, [BiLe15] showed that \(T_a\) can be interpreted as the expected time between two reproductive events along the ancestral lineage of an individual, and is given by
\[T_a = \frac{\lambda}{\mathbf{vFw}} \, .\]See [Elln18] for a detailed discussion.
- T_R0¶
The \(R_0\) generation time, i.e. the time it takes for the population to grow by a factor of its net reproductive rate \(R_0\):
\[T_{R_0} \;=\; \frac{\log R_0}{\log\lambda} \, .\]When \(\lambda \to 1\), \(T_{R_0} \to T_a\), the average age of mothers in the stable population; see [Elln18]. Thus we can extend \(T_{R_0}\) by continuity at \(\lambda = 1\).
For coherence with the literature and other libraries, the dominant eigenvalue of the next generation matrix G is used to compute the net reproductive rate \(R_0\). However, the interpretation of this quantity is not as straightforward as is usually assumed, see
R0
for details.
- T_G¶
A variant of the classic “cohort generation time”
mu1
. Here births are weighted by the cohort reproductive valuesv_G
instead of equally, as suggested by [StTC14]. This yields the simple formula\[T_G = \mathbf{v}_G (\mathbf{I} - \mathbf{S})^{-1} \mathbf{w}_G,\]as shown by [Elln18]. However the interpretation of this measure is subject to the same difficulties as that of
mu1
, and this should not be used as a measure of the typical age at which an individual is expected to reproduce.
- mu1¶
The mean age of the mothers when considering all offspring produced by a cohort. This is sometimes referred to as the cohort generation time and is given by the formula
\[\mu_1 \;=\; \frac{\mathbf{e F} (\mathbf{I} - \mathbf{S})^{-2} \mathbf{F w}}{\mathbf{e F} (\mathbf{I} - \mathbf{S})^{-1} \mathbf{F w}},\]where e is a vector of ones. See [CoEl92] and [Elln18].
Despite the name “mean age of mothers in a cohort”, this is not a good measure of the typical age at which an individual is expected to reproduce. In particular, \(\mu_1\) can be greater than the life expectancy conditional on reproduction:
>>> astrocaryum = mpm.examples.astrocaryum_mexicanum >>> astrocaryum.mu1 275.1594139265308 >>> astrocaryum.life_expectancy_repro 232.18934046629715
See [Bien19] for a detailed explanation of the problems associated with \(\mu_1\), and
mean_age_repro()
for an alternative measure of the typical age at reproduction.
- mean_age_repro(n=1000, reproduction='poisson', ini=None, target_err=0.01, return_err=False)¶
Note: unlike most descriptors, this function will return a different value every time it is called. This is because the numerical estimation is based on individual-based simulations, not an analytical formula. The last computed value is stored in
self._LAST_COMPUTED_MEAN_AGE_REPRO
.Be aware that the computing time can be very long, especially when few newborns get to reproduce.
Numerical estimate of the mean age at which a typical mother produces offspring, using individual-based simulations. Formally, this is the expected age at offspring production, conditional on reproduction; see [Bien19].
n
The number of replicates.
reproduction
The distribution of the number of offspring. Use
"bernoulli"
for monotocous species and"poisson"
(the default) otherwise.ini
The birth class of the focal individual, as a vector of weights. The default is to use
nu
, the stable distribution of newborns.target_err
if
None
, this parameter will be ignored and n replicates will be used. Otherwise, the number of replicates will start at n and double that number untilerr < target_err * mean
(see below).return_err
Whether to return the number of replicates used and \(\mathsf{err}=\mathrm{t}_{n-1} \hat{\sigma} / \sqrt{n}\), where \(\mathrm{t}_{n-1}\) is the 0.975 quantile of the Student distribution with n-1 degrees of freedom. Thus,
err
is the half-width of the 95% confidence interval for the mean of normally distributed variables. If return_err is set, then a tuple(mean, reps, err)
is returned; otherwise onlymean
is returned.
II.2.d Life expectancy
Important remark: an unavoidable problem when using matrix population models to compute age-specific quantities is that the exactitude of the formulas depend on the specific assumptions made when building the model.
For instance, in pre-breeding models, “newborns” are \(1 - \varepsilon\) years old, with \(\varepsilon\) small, whereas in post-breeding models they are \(\varepsilon\) years old. Similarly, if an individual is alive at time t and not alive at time t + 1, then we know that it died between times t and t + 1; but we have no way of knowing when.
As a result, when interpreting quantities such as the life expectancy or the mean age in a class, one should remember that “age” is counted in number of observations of the individuals, rather than absolute time. One may also want to add/subtract a correction term, e.g, based on knowledge on when deaths are most likely to occur in the projection interval.
- survivorship(t)¶
The survivorship function \(\ell\), which gives the probability that a newborn reaches age t in the stable population:
\[\ell(t) \;=\; \mathbf{e\,S}^t \boldsymbol\nu \,, \]where S is the survival matrix; e is a vector of ones; and \(\boldsymbol \nu\) is the distribution of newborns in the stable population (see
nu
).
- class_survivorship(t)¶
The class-structured survivorship, which gives the probability that an individual survives at least t time steps, as a function of the current class of that individual. That is,
\[\ell_i(t) \;=\; \mathbb{P}(T_i \geq t)\,, \]where \(T_i\) is the remaining lifespan of an individual from class i. For instance, \(\ell_i(1)\) is the probability that the individual survives to the next time step. We have:
\[\ell(t) \;=\; \mathbf{e\,S}^t\,, \]where S is the survival matrix and e is a vector of ones.
- life_expectancy¶
The expected number of projection intervals that an individual born in the stable population lives. The convention used here is to return 1 for annuals. For instance,
>>> m = MPM(S = "0", F = "1.2") >>> m.life_expectancy 1.0
Note that whether this corresponds to the expected age at death depends on the type of census (pre-breeding, post-breeding, birth-flow…), as well as on when deaths occur in the projection interval. It is also possible to get different results for models that describe the same population, if the models are build differently:
>>> eg = matpopmod.examples >>> eg.passerine_prebreeding.life_expectancy 1.7 >>> eg.passerine_postbreeding.life_expectancy 1.3399999999999999
This is normal: individuals who die shortly after their birth are not taken into account in the pre-breeding model, which biases the life expectancy upwards.
With our convention, the life expectancy is given by
\[\sum_{t \geq 0} \ell(t) \;=\; \mathbf{e}(\mathbf{I} - \mathbf{S})^{-1} \boldsymbol\nu \,, \]where \(\ell\) is the
survivorship
, e is a vector of ones and \(\boldsymbol \nu\) is the distribution of newborns in the stable population (seenu
).
- remaining_life_expectancy¶
The remaining life expectancy of an individual, as a function of its current class. See
life_expectancy
for a discussion of what this quantity represents. It is given by\[\sum_{t \geq 0} \ell_i(t) \;=\; \mathbf{e}(\mathbf{I} - \mathbf{S})^{-1}\,, \]where \(\ell_i\) is the survivorship of class i and e is a vector of ones.
- conditional_life_expectancy¶
(coming in version 0.2.0)
- life_expectancy_repro¶
The life expectancy conditional on reproducing, assuming that the fertilities are Poisson-distributed. See
life_expectancy
for a discussion of what we mean by “life expectancy”.To compute this quantity, we build a Markov chain tracking the movement of individuals in the life-cycle, with two absorbing states: corresponding to dead individuals that reproduced and dead individuals that did not reproduce. This requires knowing the probability of reproducing in each class of the model, which explains that we need to assume that the fertilities are Poisson-distributed.
Once we have this extended Markov chain, we condition it on absorption in the “died after reproducing” state (see e.g. Section 11.1.2 of [KeCa05], or Chapter III of [KeSn76] for a more detailed treatment) and compute the expected time to absorption when starting from the distribution of newborns in the stable population (see
nu
).
- remaining_life_expectancy_repro¶
The remaining life expectancy of an individual, as a function of its current class and conditional on this individual reproducing in the future. See
life_expectancy_repro
for details on how this is calculated.
- variance_lifespan¶
(coming in version 0.2.0)
- variance_remaining_lifespan¶
(coming in version 0.2.0)
II.2.e Age structure
- mean_age_class¶
The vector whose i-th component is the mean age in class i, when the population is at its stable structure. It is given by
\[\mathbf{y} \;=\; \big((\mathbf{I} - \lambda^{-1} \mathbf{S})^{-2} \boldsymbol \nu \big) \oslash \big((\mathbf{I} - \lambda^{-1} \mathbf{S})^{-1} \boldsymbol \nu \big)\,, \]where \(\boldsymbol\nu\) is the distribution of newborns in the stable population and \(\oslash\) denotes the entry-wise division – see e.g. Equation (23) in [CoEl92].
For post-breeding models, one may want to substract 1 from this quantity; here individuals in newborn classes are 1 year old.
- mean_age_population¶
The mean age in the stable population, given by \(\mathbf{y} \cdot \mathbf{w}\), where y is the vector whose i-th entry is the mean age in class i (see
mean_age_class
).For post-breeding models, one may want to substract 1 from this quantity.
- variance_age_class¶
(coming in version 0.2.0)
- variance_age_population¶
(coming in version 0.2.0)
- age_specific_fertility(t)¶
The matrix \(\mathbf{\Phi}(t)\), whose (i, j)-th entry gives the expected number of class-i offspring produced at age t by an individual in class j at age 0. It is given by
\[\mathbf{\Phi(t)} \;=\; \mathbf{F} \tilde{\mathbf{S}}(t) , \]where the (i, j)-th entry \(\tilde{s}_{ij}(t)\) of \(\tilde{\mathbf{S}}(t) = (\tilde{s}_{ij}(t))\) is the probability that an individual in class j at time 0 is in class i at time t, conditional on having survived to that time:
\[\tilde{s}_{ij}(t) \;=\; \frac{\mathbf{S}^t(i, j)}{\sum_{k} \mathbf{S}^t(k, j)}.\]See e.g. Section 11.3.2 of [KeCa05].
- mean_age_maturity¶
(coming in version 0.2.0)
- mean_age_first_repro¶
(coming in version 0.2.0)
II.2.f Entropies
See
entropy_rate
above for the entropy rate H of the genealogical Markov chain.- population_entropy¶
The population entropy, defined by Demetrius as²
\[S = H \cdot T_a\]where H is the
entropy_rate
of the model and \(T_a\) is the generation time, as given by the mean age of mothers in the stable population.For Leslie models, S has several interpretations; it corresponds to both:
The Shannon entropy of the difference of age between mothers and daughters in the stable population.
The Shannon entropy of a random trajectory in the genealogy of the population, which we refer to as the genealogical entropy.
Moreover, for those models S has the following simple expression:
\[S \;=\; -\sum_i \phi_i \lambda^{-i} \log_2 \phi_i \lambda^{-i},\]where \(\phi_i\) is the age-specific net fertility and \(\lambda\) is the asymptotic growth rate – see [Deme74], but note that in this article Demetrius refers to H, not S, as the entropy of the population.
For general projection matrices, there is no clear interpretation of S as an entropy in the usual sense.
- birth_entropy¶
(coming in version 0.2.0)
- genealogical_entropy¶
(coming in version 0.2.0)
- lifetable_entropy¶
The life-table entropy, also known as Keyfitz’s entropy to distinguish it from the population entropy introduced by Demetrius. Introduced by Keyfitz in [Keyf77], is defined as
\[- \frac{1}{L} \sum_{t \geq 0} \ell(t) \log_2 \ell(t)\,, \]where \(\ell\) is the
survivorship
, that is, the probability of reaching age t and \(L = \sum_{t \geq 0} \ell(t)\) is the life expectancy.Note that despite its accepted name, the life-table entropy is not an entropy in the usual sense. This makes its interpretation challenging – all the more so given that in the literature it is often conflated with Demetrius’ population entropy (see e.g. Section 4.3.1 of [KeCa05]).
III. Functions of the initial population
III.1 Distance to the stable population
- keyfitz_delta(n0)¶
Keyfitz’s \(\Delta\), which corresponds to the total variation distance between the stable class distribution w and the initial class distribution \(\mathbf{x} = \mathbf{n}(0) / \|\mathbf{n}(0)\|_1\). It is given by
\[\Delta \;=\; \frac{1}{2} \|\mathbf{x} - \mathbf{w}\|_1 \;=\; \frac{1}{2} \sum_{i} |x_i - w_i| \,.\]The factor 1/2 ensures that \(\Delta \in [0, 1]\).
- cohen_D1(n0)¶
Cohen’s \(D_1\) index of cumulative distance, introduced in [Cohe79]:
\[D_1 = \Bigg\|\sum_{t \geq 0} \big(\mathbf{n}(t) / \lambda^t - \mathbf{vn}(0) \mathbf{w}\big)\Bigg\|_1 \,.\]This corresponds to accumulating the difference between \(\mathbf{n}(t) / \lambda^t\) and its limit along the trajectory of the population, then taking sum of the absolute values of the entries of the resulting vector.
What makes \(D_1\) noteworth is that when the projection matrix A is (quasi-)primitive we have
\[D_1 = \big( (\mathbf{I} + \mathbf{wv} - \lambda^{-1} \mathbf{A}) - \mathbf{wv} \big) \, \mathbf{n}(0), \]which makes it straightforward to compute \(D_1\) in practice. However, one of the major drawbacks of \(D_1\) is that when \(n_i(t) / \lambda^t\) fluctuates around \(\mathbf{vn}(0)\, w_i\), the terms of the time series tend to cancel out. For instance, consider the following situation:
>>> dp = mpm.examples.dipsacus_sylvestris >>> dp.w array([0.63771991, 0.26395418, 0.01215485, 0.06931315, 0.01224112, 0.00461679]) >>> n0 = [0.2635, 0.6420, 0.0157, 0.0483, 0.0013, 0.0292] >>> dp.cohen_D1(n0) 0.037129817800633416
Here, \(D_1\) is small even though the initial population structure is very different from the stable one, and the fluctuations of \(\mathbf{n}(t) / \lambda^t\) around \(\mathbf{vn}(0) \mathbf{w}\) are initially much larger than \(D_1\):
>>> traj = dp.trajectory(n0, t_max=5) >>> numpy.sum(traj.rescaled_Y - (dp.v @ n0) * dp.w, axis=1) array([-1.43193488, 2.16076214, -0.48105167, -0.9132861 , 1.02463032, -0.31244266])
This makes the relevance of \(D_1\) questionable.
- population_momentum(n0, auto=True)¶
Population momentum is the name given by Keyfitz in [Keyf71] to the following phenomenon: when a population undergoes a demographic transition, it can keep growing even after its asymptotic growth rate has become equal to 1. The reason for this is that the population is not at its stable structure, and that the increase / decrease that we see corresponds to the transient regime.
The standard way to quantify this is to use the following quantity, often referred to as the momentum of the population:
\[M = \lim_{t \to +\infty} \frac{\|\mathbf{n}(t)\|_1}{\|\mathbf{n}(0)\|_1}, \]where t = 0 is the time at which the asymptotic growth rate of the population becomes equal to 1. For (quasi-)primitive, stationary projection matrices, this limit always exists and is given by
\[M = \frac{\mathbf{vn}(0)}{\|\mathbf{n}(0)\|_1}, \]where v is the vector of reproductive values of the projection matrix, rescaled so that vw = 1.
Computing the population momentum requires knowing two things: (1) the population vector at the time t = 0 of the transitions and (2) the projection matrix after the transition; specifying the projection matrix before the transition is not sufficient. Therefore, this function makes the following assumptions:
If λ = 1, we assume that the projection matrix of the model describes the dynamics of the population after the demographic transition, and we return \(\mathbf{vn}(0) / \|\mathbf{n}(0)\|_1\).
If λ ≠ 1, we do not know what projection matrix should be used after the transition, so
NotAvailable
is raised. However, if auto isTrue
, then instead of raising an error we assume that the projection matrix after the transition is\[\mathbf{A}' = \mathbf{S} + c^{-1} \mathbf{F}\]where c is the
fertility_excess
, i.e. the constant such that \(\lambda(\mathbf{S} + c^{-1} \mathbf{F}) = 1\). This corresponds to Keyfitz’s original idea of assuming that all fertilities get reduced by the same factor during the demographic transition.
III.2 Trajectories
- trajectory(n0, t_max, step=1)¶
The trajectory of the population vector associated to the initial condition n0, up to time t = t_max. This is the solution of
\[\mathbf{n}(t+1) \;=\; \mathbf{A} \mathbf{n}(t)\,, \quad t = 0, \ldots, t_{\max}.\]The population vector will only be computed for multiples of the optional argument step. For instance, use step = 3 to get the population vector at times \(t= 0, 3, 6, \ldots\)
Returns a
Trajectory
object. Trajectories have a set of methods that facilitate their study and can be plotted using the functiontrajectory()
from the moduleplot
.
- stochastic_trajectories(n0, t_max, reps, reproduction='poisson', Z=None)¶
Returns a list of reps independent stochastic trajectories corresponding to the matrix population model. Each trajectory starts at time t = 0 from the integer-valued population vector n0, and ends at time t = t_max.
In ecological terms, those stochastic trajectories account for demographic stochasticity – that is, the randomness stemming from the variation in the realizations of the demographic rates of the individuals (as opposed to randomness in the environment affecting the demographic rates population-wide).
Formally, the trajectories are realizations of a multitype Galton–Watson process whose types correspond to the classes of the model. At each time-step, each class-j individual x:
produces a random number \(F_{ij}(x)\) of class-i individuals, where \(\mathbb{E}(F_{ij}(x)) = f_{ij}\), the (i, j)-th entry of the fertility matrix F;
survives and goes to class i with probability \(s_{ij}\), the (i, j)-th entry of the survival matrix S.
What happens to individuals is independent of what happens to other individuals or of anything that happened in the past. This crucial assumption is the cornerstone of the link between multitype Galton–Watson processes and matrix population models. It ensures that \(\mathbf{n}(t)\), the expected value of the (integer-valued) vector tracking the number of individuals in each class of the stochastic process, satisfies the linear equation
\[\mathbf{n}(t+1) = \mathbf{A} \mathbf{n}(t)\]with A = S + F is the projection matrix of the matrix population model.
Each Galton–Watson process corresponds to a single matrix population model, but there are an infinity of Galton–Watson processes for a given matrix population model. To fully specify a Galton–Watson process, one needs to provide the joint law of the family \((S_{ij}, F_{ij})_i\), for every j. Two standard sets of assumptions can be selected via the argument reproduction:
"poisson"
— for all j, \((F_{ij})_i\) are independent Poisson variables with respective parameters \(f_{ij}\), and are independent of the family \((S_{ij})_i\). If the survival and fertility matrices are not available, instead of the model described above we simply take the number of class-i individuals that a class-j individual leaves at the next time-step to follow a \(\mathrm{Poisson}(a_{ij})\) distribution (instead of \(\mathrm{Bernoulli}(s_{ij}) + \mathrm{Poisson}(f_{ij})\) in the more detailed model).This is the assumption that we recommend using in the absence of additional information.
"bernoulli"
— for all j, \((F_{ij})_i\) are incompatible Bernoulli variables that are independent of \((S_{ij})_i\). In other words, each class-j individual produces an offspring with probability \(\sum_k f_{kj} \leq 1\), and that offspring is produced in class i with probability \(f_{ij} / \sum_k f_{kj}\). If the survival and fertility matrices are not available or if \(\sum_k f_{kj} > 1\), an error will be raised.This assumption is what we recommend using for monotocous species.
Alternatively, one can use a custom Galton–Watson process by specifying the argument Z. In that case, Z should be a random Python function such that the dynamics of the population vector \(N(t)\) is defined by:
\[N(t+1) = Z\big(\mathbf{F}, \mathbf{S}, N(t)\big)\]if S and F are known, and \(N(t+1) = Z(\mathbf{A}, \mathrm{None}, N(t))\) otherwise.
Because this function focuses on efficiency, it is not possible to specify \((S_{ij}, F_{ij})_i\) only; the user has to use them to work out how to sample Z. See the module
ibm
to generate stochastic trajectories directly from \((S_{ij}, F_{ij})_i\).
Footnotes
When we say that the population satisfies \(\mathbf{n}(t) \sim \mathbf{vn}(0) \lambda^t \mathbf{w}\) for “any” initial population vector, in fact this is true if and only if \(\mathbf{vn}(0) \neq 0\). For instance, if \(\mathbf{n}(0)\) were colinear with another right-eigenvector of the projection matrix – say \(\mathbf{w}'\) with associated eigenvalue \(\lambda'\) – then we would have \(\mathbf{n}(t) \sim \mathbf{n}(0) (\lambda')^t\).
However, the set of vectors x such that vx = 0 is a vector space of dimension n - 1 in the n-dimensional space of population vectors. Thus, when choosing an initial population vector according to any “reasonable” probability distribution, the probability of getting a vector such that \(\mathbf{vn}(0) = 0\) is equal to 0, just as there is a probability zero of falling exactly on a given line when throwing a point in the plane. As a result, the case \(\mathbf{vn}(0) = 0\) can be considered a degenerate case with no biological relevance.
The terminology regarding the quantities H and S has not been very consistent; in fact Demetrius himself has used both population entropy and evolutionary entropy to refer to either of the two. The convention we use here corresponds to a distinction explicitly made by Demetrius in a correspondence with FB, where he referred to H as the evolutionary entropy and to S as the population entropy. However, we favor the term entropy rate over evolutionary entropy because it is more explicit.
While the definition and interpretation of S are very clear and natural in the context of Leslie models, this is not the case for general matrix population models. Based on the same correspondence with Demetrius, we adopt the point of view that S should be defined so as to ensure that \(S = H / T_a\). This, however, means that we have to renounce to the probabilistic interpretations that S has in the context of Leslie models.