Asymptotic Distribution of Density-Dependent Stage-Grouped Population Dynamics Models

Abstract

Matrix models are widely used in biology to predict the temporal evolution of stage-structured populations. One issue related to matrix models that is often disregarded is the sampling variability. As the sample used to estimate the vital rates of the models are of finite size, a sampling error is attached to parameter estimation, which has in turn repercussions on all the predictions of the model. In this study, we address the question of building confidence bounds around the predictions of matrix models due to sampling variability. We focus on a density-dependent Usher model, the maximum likelihood estimator of parameters, and the predicted stationary stage vector. The asymptotic distribution of the stationary stage vector is specified, assuming that the parameters of the model remain in a set of the parameter space where the model admits one unique equilibrium point. Tests for density-dependence are also incidentally provided. The model is applied to a tropical rain forest in French Guiana.

Appendices

Appendix A: Differential of the Application $\theta\mapsto{\mathbf{N}}_{\infty}$

Assuming that it is an application, h maps the vector of parameters θ onto the asymptotic stage distribution ${\mathbf N}_{\infty}.$ We hereafter focus on the simplified model (3). Thus, $h:{\mathbf R}^{10}\rightarrow{\mathbf R}^m,$ and the differential d_θ h is a m × 10 matrix whose element on the ith row and jth column is $\partial N_i/\partial\theta_j,$ where N _i is the ith element of ${\mathbf N}_{\infty}$ and θ _j is the jth element of θ. The function h is implicitly defined as Φ[h(θ),θ] = 0 where $\Phi:{\mathbf R}^{m}\times{\mathbf R}^{10}\rightarrow{\mathbf R}^{m}$ is given by (4). Let Φ _i be the ith element of Φ (1 ≤ i ≤ m), let d_θΦ be the partial differential of Φ represented by the m × 10 matrix $(\partial\Phi_i/\partial\theta_j)_{1\leq i\leq m,1\leq j\leq 10},$ and let $\hbox{d}_{\mathbf N}\Phi$ be the partial differential of Φ represented by the m × m matrix $(\partial\Phi_i/\partial N_j)_{1\leq i,j\leq m}.$ Then, by the implicit function theorem,

$$ \hbox{d}_{\theta}h=-\left(\hbox{d}_{\mathbf N}\Phi)^{-1}\right|_{(h(\theta),\theta)} \cdot\left.\hbox{d}_{\theta}\Phi\right|_{(h(\theta),\theta)} $$

We are now left with computing the partial derivatives of Φ.

For the simplified model (3) and $\forall i,j,$

$$ \frac{\partial p_i}{\partial N_j}=-\xi a_j\,p_i({\mathbf{N}}), \quad\frac{\partial q_i}{\partial N_j}=-\kappa a_j\,q_i({\mathbf{N}}), \quad\frac{\partial f}{\partial N_j}=-\beta a_j\,f({\mathbf{N}}) $$

where a _j is the jth element of a. Using the developed form of the matrix relationship,

$$ \Phi_i({\mathbf{N}},\theta)=\left\{ \begin{array}{lr} N_1-p_1({\mathbf{N}};\theta)\,N_1-f({\mathbf{N}};\theta)\,\sum_{j=1}^mN_j & (i=1)\\ N_i-p_i({\mathbf{N}};\theta)\,N_i-q_{i-1}({\mathbf{N}};\theta)\,N_{i-1} & (2\leq i\leq m) \end{array} \right. $$

The partial derivatives of Φ with respect to ${\mathbf N}$ are thus:

• for the first element of Φ,

  $$ \frac{\partial\Phi_1}{\partial N_j}=\left\{ \begin{array}{lr} 1-p_1+\xi a_1p_1N_1+\beta a_1f\sum_{k=1}^mN_k-f & (j=1)\\ \xi a_jp_1N_1+\beta a_jf\sum_{k=1}^mN_k-f &(2\leq j\leq m) \end{array} \right. $$

• for 2 ≤ i ≤ m,

  $$ \frac{\partial\Phi_i}{\partial N_j}=\left\{ \begin{array}{lr} 1-p_i+\xi a_ip_iN_i+\kappa a_iq_{i-1}N_{i-1} & (j=i)\\ -q_{i-1}+\xi a_{i-1}p_iN_i+\kappa a_{i-1}q_{i-1}N_{i-1} & (j=i-1)\\ \xi a_jp_iN_i+\kappa a_jq_{i-1}N_{i-1} & (j\neq i,i-1) \end{array} \right. $$

Moreover, given that the parameters are in the order θ = (μ₀, μ₁, μ₂, ξ, ν₀, ν₁, ν₂, κ, α, β), the partial derivatives of Φ with respect to θ are:

• for 1 ≤ j ≤ 4,

  $$ \frac{\partial\Phi_i}{\partial\theta_j}= -\frac{\partial p_i}{\partial\theta_j}N_i\quad(1\leq i\leq m) $$

• for 5 ≤ j ≤ 8,

  $$ \frac{\partial\Phi_i}{\partial\theta_j}=\left\{ \begin{array}{lr} 0 & (i=1)\\ -(\partial q_{i-1}/\partial\theta_j)N_{i-1} & (2\leq i\leq m) \end{array} \right. $$

• for j = 9, 10,

  $$ \frac{\partial\Phi_i}{\partial\theta_j}=\left\{ \begin{array}{lr} -(\partial f/\partial\theta_j)\sum_{k=1}^mN_k & (i=1)\\ 0 & (2\leq i\leq m) \end{array} \right. $$

where the partial derivatives of the vital rates with respect to the parameters θ are:

$$ \begin{aligned} \frac{\partial p_i}{\partial\theta_j} &= i^{j-1}\exp(-\xi{\mathbf{N}}.{\mathbf{a}})\hbox{ for }1\leq j\leq 3,\quad\frac{\partial p_i}{\partial\theta_4}=-{\mathbf{N}}.{\mathbf{a}}\,p_i\\ \frac{\partial p_i}{\partial\theta_j} &= 0\hbox{ for }5\leq j\leq 10\\ \frac{\partial q_i}{\partial\theta_j} &= i^{j-5}\exp(-\kappa{\mathbf{N}}.{\mathbf{a}})\hbox{ for }5\leq j\leq 7,\quad\frac{\partial q_i}{\partial\theta_8}=-{\mathbf{N}}.{\mathbf{a}}\,q_i\\ \frac{\partial q_i}{\partial\theta_j} &= 0\hbox{ for }1\leq j\leq 4\hbox{ or }j=9,10\\ \frac{\partial f}{\partial\theta_j} &= 0\hbox{ for }1\leq j\leq 8,\quad\frac{\partial f}{\partial\theta_9}=\frac{f}{\alpha},\quad \frac{\partial f}{\partial\theta_{10}}=-{\mathbf{N}}.{\mathbf{a}}\,f \end{aligned} $$

Finally, to compute the differential of h at θ it is necessary to compute ${\mathbf N}_{\infty}=h(\theta),$ that is to say to solve $\Phi({\mathbf N},\theta)=0$ for N. To this end, we used the numerical method proposed by Caswell (2001, Sect. 16.2.1 p. 518).

Appendix B: Fisher Information Matrix of θ

Computations are done for the simplified model (3). The 10 parameters are taken in the order θ = (μ₀, μ₁, μ₂, ξ, ν₀, ν₁, ν₂, κ, α, β). The Fisher information matrix I(θ) for a sample of size one is a 10 × 10 matrix with element:

$$ [I(\theta)]_{ij}=\hbox{E}\left[\left(\frac{\partial \ln \ell_{\theta}(X)} {\partial\theta_i}\right) \left(\frac{\partial \ln \ell_{\theta}(X)}{\partial\theta_j}\right)\right] $$

where ℓ_θ is given by (7) and the expectation is over the distribution (5) of the observation X. As the estimators of α and β are independent from the other parameters, I(θ) is actually of the form: where $\theta^{\ast}$ is the restriction of θ to the eight first parameters. The logarithm of ℓ_θ is:

$$ \begin{aligned} \ln \ell_{\theta}(X) &= \sum_{i=1}^m\sum_{k=1}^K \ln d_{ik}{\mathbf{1}}_{X=(i,\bullet,k)}+ \sum_{k=1}^K[ \ln (1-f_k^{\ast}){\mathbf{1}}_{X=(\bullet,\bullet,k)} + \ln f_k^{\ast}{\mathbf{1}}_{X=(0,1,k)}] \\ &+ \sum_{i=1}^m\sum_{k=1}^K[ \ln p_{ik}{\mathbf{1}}_{X=(i,i,k)}+ \ln q_{ik}{\mathbf{1}}_{X=(i,i+1,k)}+ \ln (1-p_{ik}-q_{ik}) {\mathbf{1}}_{X=(i,\dag,k)}] \\ &+ \sum_{k=1}^K \ln \rho_k{\mathbf{1}}_{X=(.,.,k)} \end{aligned} $$

where the small dot denotes any stage (including 0), the big dot denotes any stage except 0 and, by convention, q _mk  = 0. The partial derivatives of ln ℓ_θ with respect to θ then are:

$$ \frac{\partial \ln \ell_{\theta}(X)}{\partial\theta_j}=\left\{ \begin{array}{lr} \sum \limits_{i=1}^m\sum\limits_{k=1}^K\frac{\partial p_{ik}}{\partial\theta_j}\left[\frac{{\mathbf{1}}_{X=(i,i,k)}}{p_{ik}} -\frac{{\mathbf{1}}_{X=(i,\dag,k)}}{1-p_{ik}-q_{ik}}\right] & (1\leq j\leq 4)\\ \sum\limits_{i=1}^{m-1}\sum\limits_{k=1}^K\frac{\partial q_{ik}}{\partial\theta_j}\left[\frac{{\mathbf{1}}_{X=(i,i+1,k)}}{q_{ik}} -\frac{{\mathbf{1}}_{X=(i,\dag,k)}}{1-p_{ik}-q_{ik}}\right] & (5\leq j\leq 8)\\ \sum\limits_{k=1}^K\frac{\partial f_k^{\ast}}{\partial\theta_j}\left[ \frac{{\mathbf{1}}_{X=(0,1,k)}}{f_k^{\ast}}-\frac{{\mathbf{1}}_{X=( \bullet,\bullet,k)}}{1-f_k^{\ast}}\right] & (9\leq j\leq 10) \end{array} \right. $$

where, using (6),

$$ \frac{\partial f_k^{\ast}}{\partial\theta_j}=\frac{1}{(1+f_k)^2} \frac{\partial f_k}{\partial\theta_j} $$

and the partial derivatives of the vital rates with respect to the parameters θ are given in Appendix A. The non-null elements of the Fisher information matrix then are:

• for 1 ≤ j,l ≤ 4,

  $$ [I(\theta)]_{jl}=\sum_{i=1}^m\sum_{k=1}^K\frac{\partial p_{ik}}{\partial\theta_j}\frac{\partial p_{ik}}{\partial\theta_l} \rho_k(1-f_k^{\ast})d_{ik}\left(\frac{1}{p_{ik}}+\frac{1}{1-p_{ik}-q_{ik}} \right) $$

• for 5 ≤ j,l ≤ 8,

  $$ [I(\theta)]_{jl}=\sum_{i=1}^{m-1}\sum_{k=1}^K\frac{\partial q_{ik}}{\partial\theta_j}\frac{\partial q_{ik}}{\partial\theta_l} \rho_k(1-f_k^{\ast})d_{ik}\left(\frac{1}{q_{ik}}+\frac{1}{1-p_{ik}-q_{ik}} \right) $$

• for 1 ≤ j ≤ 4 and 5 ≤ l ≤ 8,

  $$ [I(\theta)]_{jl}=[I(\theta)]_{lj}=\sum_{i=1}^m\sum_{k=1}^K\frac{\partial p_{ik}}{\partial\theta_j}\frac{\partial q_{ik}}{\partial\theta_l} \rho_k(1-f_k^{\ast})d_{ik}\frac{1}{1-p_{ik}-q_{ik}} $$

• for 9 ≤ j,l ≤ 10,

  $$ \begin{aligned} {[I(\theta)]_{9,9}} &= \frac{1}{\alpha^2}\sum_{k=1}^K\rho_k\frac{f_k}{(1+f_k)^2}, \quad [I(\theta)]_{10,10}=\sum_{k=1}^K\rho_k({\mathbf{N}}_k.{\mathbf{a}})^2 \frac{f_k}{(1+f_k)^2},\\ {[I(\theta)]_{9,10}} &= [I(\theta)]_{10,9}=-\frac{1}{\alpha}\sum_{k=1}^K\rho_k({\mathbf{N}}_k.{\mathbf{a}}) \frac{f_k}{(1+f_k)^2} \end{aligned} $$