Chapter 10 Population Projection III: Projection Simulations
Anyone who believes in indefinite growth of anything physical on a physically finite planet is either a madman or an economist.
So far we have covered a number of standard approaches to matrix projection analysis. Basic deterministic analyses and projection assuming temporal environmental stochasticity dominate the MPM literature. However, these analysis types have strong assumptions that are often not met in real life. Indeed, lack of concordance between these assumptions and reality may be behind the general failure of matrix-based population viability analyses to predict successfully population dynamics beyond just a time step or two into the future (Crone et al. 2013).
Package lefko3 provides some general projection options for users interested in more complex analyses. For example, users can project function-based matrices using a particular set of individual covariate values, or conduct density-dependent analyses using a variety of density dependence functions, or even run cyclical projections or replicate stochastic simulations. Further, package lefko3 provides the added capability to enforce substochasticity in all projections, thus preventing analysis from involving impossible rates and probabilities.
Chapter 9 introduced an equation showing how matrices are projected under the assumption of temporal environmental stochasticity:
\[\begin{equation} \tag{10.1} \mathbf{n_t} = \mathbf{A_t A_{t-1} \cdots A_1 n_0} \end{equation}\]
This is really a general projection equation, where each \(\mathbf{A_i}\) simply represents the matrix to be used at time \(i\). We can project forward any number of times, perhaps only a small number to see how a particular management regime might impact population dynamics in the immediate term, or perhaps a very large number to see what long-term patterns we find, such as cycles, plateaus, or chaos. These matrices can be chosen randomly or set to a specific order.
Package lefko3 provides two functions for general projection: projection3() and f_projection3(). The former is used to project matrices that have already been developed and exist within a lefkoMat object. The latter is used to take vital rate models and build new function-based matrices for each time step. Both functions have been developed to run as binaries in R, making them run as fast as possible.
Many kinds of projections are possible with these functions, including deterministic, ordered / cyclical, and stochastic, both with and without density dependence. Let’s now to turn to density dependence, in particular, since we have not discussed the topic before.
10.1 Density dependence
Density dependent matrix projections involve the modification of matrix elements by some function of population size. Many different functions can be applied to matrix elements in order to incorporate density dependence into matrix projection. The oldest typically applied a function such as the logistic function to all matrix elements, generally multiplying each matrix element by some function of density at either the current time or some time in the past (e.g., Leslie 1959). More recent approaches typically separate density dependent functions by vital rate, and alter only certain specific matrix elements (e.g., Jensen 1995).
To provide the greatest flexibility, lefko3 currently implements density dependence using the density_input() function, which allows users to list all matrix elements that should be subject to density dependence, and to stipulate the characteristics of the density dependence. Shorthand abbreviations can be used, such as all for all stages, rep for reproductive stages, nrep for mature but non-reproductive stages, etc. A default time delay of one time step is applied, and this can be set to any positive integer. Currently, lefko3 includes four density dependence functions that can be chosen, with projections capable of including multiple functions and time delays. The first is the two-parameter Ricker function, given as
\[\begin{equation} \tag{10.2} \phi_{t+1} = \phi_t \times \alpha e^{-\beta n_t} \end{equation}\]
where \(\alpha\) and \(\beta\) are the density dependence parameters, and \(\beta\) in particular gives the strength of density dependence. This equation is among the most interesting in the density dependence literature, because it is capable of producing single or multi-period oscillations, as well as chaos.
The second density dependence function in lefko3 is the two-parameter Beverton-Holt function, given as
\[\begin{equation} \tag{10.3} \phi_{t+1} = \phi_t \times \frac{\alpha}{1 + \beta n_t} \end{equation}\]
where \(\alpha\) and \(\beta\) are, once again, the density dependence parameters, and \(\beta\) in particular gives the strength of density dependence. This function generally asymptotes at an equilibrium density, and so is not capable of producing some of the complex patterns that the Ricker function is capable of.
The third function implemented is the Usher function, given as
\[\begin{equation} \tag{10.4} \phi_{t+1} = \phi_t \times \frac{1}{1 + e^{\alpha n_t + \beta}} \end{equation}\]
where both \(\alpha\) and \(\beta\) give the strength of density dependence, though the former via an interaction with density and the latter via an addition to the exponential effect.
Finally, lefko3 also implements a form of the logistic function, given as
\[\begin{equation} \tag{10.5} \phi_{t+1} = \phi_t \times (1 - \frac{n_t}{K}) \end{equation}\]
which classically takes only one parameter, the carrying capacity \(K\). For this function, the user may also stipulate whether \(K\) is a hard limit, or whether the time lag in density can result in overshooting \(K\).
10.1.1 Substochasticity
Functions projection3() and f_projection3() allow the user to stipulate whether the projection should force matrices to remain substochastic. Substochasticity refers to the condition that estimated probabilities stay within the bounds of probabilities, meaning that they should not equal values below 0.0 or above 1.0. Values outside this range may occur when using some density dependence functions, for example when survival transitions are subject to the Ricker function.
Two forms of substochastic enforcement are included in these functions: mild and hard. In the mild form, matrix elements used as survival-transition probabilities are prevented from moving outside of the interval [0, 1]. In the hard form, column sums in survival-transition matrices (equivalent to the survival probabilities of stages, stage-pairs, or age-stages) are kept within the interval [0, 1] by adjusting matrix elements by proportionate amounts. This latter form is due to the fact that the column sums of each survival-transition matrix, \(\mathbf{U}\), must equal the overall time-step survival of each stage, age-stage combination, stage-pair, or whatever aspect of the life history model that each column represents. Both forms of substochastic enforcement also prevent fecundity from becoming negative. The default setting does not enforce substochasticity, but we encourage users to try enforcing substochasticity when experimenting with density dependence, as logically impossible values may occur, and others have warned that the consequences of incorporating logically impossible values on analyses of population dynamics are largely unknown (Caswell 2001).
Let’s consider what happens mathematically. Both forms of substochastic enforcement start by identifying which elements are negative, and changing them to 0. This is done first within the \(\mathbf{U}\) matrix, and then in the \(\mathbf{F}\) matrix. Thus, if \(U _{ij}\) is the survival transition element in the \(i\)th row and the \(j\)th column, then
\[\begin{equation} \tag{10.6} U _{ij} = \left\{\begin{array}{cc} U _{ij} & U _{ij} > 0 \\ 0 & U _{ij} < 0 \end{array}\right. \end{equation}\]
for all \(U _{ij}\). The equation is the same for all \(F _{ij}\).
Once this operation is complete, mild substochastic enforcement changes all survival-transition elements above 1.0 to 1.0, as in
\[\begin{equation} \tag{10.7} U _{ij} = \left\{\begin{array}{cc} U _{ij} & U _{ij} < 1 \\ 1 & U _{ij} > 1 \end{array}\right. \end{equation}\]
for all \(U _{ij}\). Under hard substochastic enforcement, we instead identify survival-transition column sums greater than 1.0 and alter survival-transitions, such that
\[\begin{equation} \tag{10.8} U _{ij} = \left\{\begin{array}{cc} U _{ij} & \sum _i ^m U _{ij} < 1 \\ \frac{U _{ij}}{\sum _i ^m U _{ij}} & \sum _i ^m U _{ij} > 1 \end{array}\right. \end{equation}\]
for all \(U _{ij}\). Here, \(m\) refers to the number of rows (and hence columns, since the matrices must be square).
10.2 Projecting existing MPMs
Let’s explore general projection simulations with projection3() using our raw MPMs developed with the Cypripedium candidum dataset. Here we see the life history model that we will use, which was originally introduced in chapter 3.
Figure 10.1: Life history model of Cypripedium candidum for use in raw MPMs
As it has been several chapters since we have seen that code, we show the full code to build our raw Cypripedium MPMs below.
rm(list=ls(all=TRUE))
library(lefko3)
data(cypdata)
sizevector <- c(0, 0, 0, 0, 0, 0, 1, 3, 6, 11, 19.5)
stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "XSm", "Sm", "Md", "Lg",
"XLg")
repvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
obsvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
matvector <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)
immvector <- c(0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0)
propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
indataset <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)
binvec <- c(0, 0, 0, 0, 0, 0.5, 0.5, 1.5, 1.5, 3.5, 5)
comments <- c("Dormant seed", "1st yr protocorm", "2nd yr protocorm",
"3rd yr protocorm", "Seedling", "Dormant adult",
"Extra small adult (1 shoot)", "Small adult (2-4 shoots)",
"Medium adult (5-7 shoots)", "Large adult (8-14 shoots)",
"Extra large adult (>14 shoots)")
cypframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector,
repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
propstatus = propvector, immstatus = immvector, indataset = indataset,
binhalfwidth = binvec, comments = comments)
cypraw_v1 <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004,
patchidcol = "patch", individcol = "plantid", blocksize = 4,
sizeacol = "Inf2.04", sizebcol = "Inf.04", sizeccol = "Veg.04",
repstracol = "Inf.04", repstrbcol = "Inf2.04", fecacol = "Pod.04",
stageassign = cypframe_raw, stagesize = "sizeadded", NAas0 = TRUE,
NRasRep = TRUE)
seeds_per_fruit <- 5000
sl_mult <- 0.7
cypsupp2_raw <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "SL", "D",
"XSm", "Sm", "SD", "P1"),
stage2= c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "SL", "rep", "rep"),
eststage3 = c(NA, NA, NA, NA, NA, NA, "D", "XSm", "Sm", NA, NA),
eststage2 = c(NA, NA, NA, NA, NA, NA, "XSm", "XSm", "XSm", NA, NA),
givenrate = c(0.08, 0.10, 0.10, 0.10, 0.05, 0.05, NA, NA, NA, NA, NA),
multiplier = c(NA, NA, NA, NA, NA, NA, sl_mult, sl_mult, sl_mult,
0.5 * seeds_per_fruit, 0.5 * seeds_per_fruit),
type =c(1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3),
stageframe = cypframe_raw, historical = FALSE)
cypsupp3_raw <- supplemental(stage3 = c("SD", "SD", "P1", "P1", "P2", "P2",
"P3", "SL", "SL", "SL", "D", "XSm", "Sm", "D", "XSm", "Sm", "mat", "mat",
"mat", "SD", "P1"),
stage2 = c("SD", "SD", "SD", "SD", "P1", "P1", "P2", "P3", "SL", "SL", "SL",
"SL", "SL", "SL", "SL", "SL", "D", "XSm", "Sm", "rep", "rep"),
stage1 = c("SD", "rep", "SD", "rep", "SD", "rep", "P1", "P2", "P3", "SL",
"SL", "SL", "SL", "P3", "P3", "P3", "SL", "SL", "SL", "mat", "mat"),
eststage3 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "D", "XSm", "Sm", "D",
"XSm", "Sm", "mat", "mat", "mat", NA, NA),
eststage2 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "XSm", "XSm", "XSm",
"XSm", "XSm", "XSm", "D", "XSm", "Sm", NA, NA),
eststage1 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "XSm", "XSm", "XSm",
"XSm", "XSm", "XSm", "XSm", "XSm", "XSm", NA, NA),
givenrate = c(0.08, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.05, 0.05, 0.05, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA),
multiplier = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, sl_mult, sl_mult,
sl_mult, sl_mult, sl_mult, sl_mult, 1, 1, 1, 0.5 * seeds_per_fruit,
0.5 * seeds_per_fruit),
type =c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3),
stageframe = cypframe_raw, historical = TRUE)
cypmatrix2r <- rlefko2(data = cypraw_v1, stageframe = cypframe_raw,
year = "all", stages = c("stage3", "stage2"),
size = c("size3added", "size2added"), supplement = cypsupp2_raw,
yearcol = "year2", patchcol = "patchid", indivcol = "individ")
cypmatrix2rp <- rlefko2(data = cypraw_v1, stageframe = cypframe_raw,
year = "all", patch = "all", stages = c("stage3", "stage2"),
size = c("size3added", "size2added"), supplement = cypsupp2_raw,
yearcol = "year2", patchcol = "patchid", indivcol = "individ")
cypmatrix3r <- rlefko3(data = cypraw_v1, stageframe = cypframe_raw,
year = "all", stages = c("stage3", "stage2", "stage1"),
size = c("size3added", "size2added", "size1added"), supplement = cypsupp3_raw,
yearcol = "year2", patchcol = "patchid", indivcol = "individ")
cypmatrix3rp <- rlefko3(data = cypraw_v1, stageframe = cypframe_raw,
year = "all", patch = "all", stages = c("stage3", "stage2", "stage1"),
size = c("size3added", "size2added", "size1added"), supplement = cypsupp3_raw,
yearcol = "year2", patchcol = "patchid", indivcol = "individ")
We see that we have four MPMs. The first two, cypmatrix2r and cypmatrix2rp, are ahistorical MPMs at the population and patch level, respectively. The next two, cypmatrix3r and cypmatrix3rp, are historical MPMs at the population and patch level, respectively. Let’s see their summaries to jumpstart our memories about them.
summary(cypmatrix2r)
>
> This ahistorical lefkoMat object contains 5 matrices.
>
> Each matrix is square with 11 rows and columns, and a total of 121 elements.
> A total of 120 survival transitions were estimated, with 24 per matrix.
> A total of 40 fecundity transitions were estimated, with 8 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 5 time steps.
>
> The dataset contains a total of 74 unique individuals and 320 unique transitions.
>
> Survival probability sum check (each matrix represented by column in order):
> [,1] [,2] [,3] [,4] [,5]
> Min. 0.000 0.050 0.050 0.000 0.050
> 1st Qu. 0.100 0.140 0.140 0.100 0.140
> Median 0.689 0.870 0.864 0.610 0.882
> Mean 0.552 0.629 0.629 0.528 0.627
> 3rd Qu. 1.000 1.000 1.000 0.960 1.000
> Max. 1.000 1.000 1.000 1.000 1.000
summary(cypmatrix2rp)
>
> This ahistorical lefkoMat object contains 15 matrices.
>
> Each matrix is square with 11 rows and columns, and a total of 121 elements.
> A total of 266 survival transitions were estimated, with 17.733 per matrix.
> A total of 70 fecundity transitions were estimated, with 4.667 per matrix.
> This lefkoMat object covers 1 population, 3 patches, and 5 time steps.
>
> The dataset contains a total of 74 unique individuals and 320 unique transitions.
>
> Survival probability sum check (each matrix represented by column in order):
> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
> Min. 0.000 0.000 0.000 0.000 0.000 0.000 0.050 0.050 0.000 0.050 0.000 0.000
> 1st Qu. 0.075 0.025 0.075 0.025 0.075 0.075 0.140 0.140 0.100 0.140 0.100 0.100
> Median 0.180 0.100 0.180 0.100 0.180 0.180 0.909 0.778 0.686 0.857 0.750 0.575
> Mean 0.457 0.361 0.471 0.328 0.417 0.464 0.631 0.611 0.530 0.631 0.562 0.523
> 3rd Qu. 0.955 0.769 1.000 0.592 0.781 1.000 1.000 1.000 0.955 1.000 1.000 1.000
> Max. 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
> [,13] [,14] [,15]
> Min. 0.000 0.000 0.000
> 1st Qu. 0.075 0.075 0.100
> Median 0.180 0.180 0.750
> Mean 0.432 0.450 0.562
> 3rd Qu. 0.875 1.000 1.000
> Max. 1.000 1.000 1.000
summary(cypmatrix3r)
>
> This historical lefkoMat object contains 4 matrices.
>
> Each matrix is square with 121 rows and columns, and a total of 14641 elements.
> A total of 242 survival transitions were estimated, with 60.5 per matrix.
> A total of 52 fecundity transitions were estimated, with 13 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 4 time steps.
>
> The dataset contains a total of 74 unique individuals and 320 unique transitions.
>
> Survival probability sum check (each matrix represented by column in order):
> [,1] [,2] [,3] [,4]
> Min. 0.000 0.000 0.000 0.000
> 1st Qu. 0.000 0.000 0.000 0.000
> Median 0.000 0.000 0.000 0.000
> Mean 0.173 0.179 0.166 0.198
> 3rd Qu. 0.100 0.100 0.100 0.100
> Max. 1.000 1.000 1.000 1.000
summary(cypmatrix3rp)
>
> This historical lefkoMat object contains 12 matrices.
>
> Each matrix is square with 121 rows and columns, and a total of 14641 elements.
> A total of 516 survival transitions were estimated, with 43 per matrix.
> A total of 68 fecundity transitions were estimated, with 5.667 per matrix.
> This lefkoMat object covers 1 population, 3 patches, and 4 time steps.
>
> The dataset contains a total of 74 unique individuals and 320 unique transitions.
>
> Survival probability sum check (each matrix represented by column in order):
> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
> Min. 0.000 0.0000 0.0000 0.000 0.000 0.000 0.00 0.000 0.000 0.0000 0.000
> 1st Qu. 0.000 0.0000 0.0000 0.000 0.000 0.000 0.00 0.000 0.000 0.0000 0.000
> Median 0.000 0.0000 0.0000 0.000 0.000 0.000 0.00 0.000 0.000 0.0000 0.000
> Mean 0.107 0.0945 0.0851 0.101 0.158 0.158 0.14 0.169 0.119 0.0851 0.119
> 3rd Qu. 0.000 0.0000 0.0000 0.000 0.100 0.100 0.05 0.100 0.000 0.0000 0.000
> Max. 1.000 1.0000 1.0000 1.000 1.000 1.000 1.00 1.000 1.000 1.0000 1.000
> [,12]
> Min. 0.000
> 1st Qu. 0.000
> Median 0.000
> Mean 0.144
> 3rd Qu. 0.000
> Max. 1.000
Users of lefko3 might wish to conduct several styles of projections, including the following:
- Deterministic projections of a specific time length
- Ordered or cyclical projections of a specific time length
- Replicated simulations assuming temporal stochasticity
- Density dependent deterministic projections
- Ordered or cyclical density dependent simulations
- Replicated simulations assuming temporal stochasticity and density dependence
Projection styles 1, 2, 4, and 5 assume no randomness and so we assume that users are not interested in producing replicates for those styles. However, both projection3() and f_projection3() allow replication for all styles. Let’s begin with the first style of projection.
10.2.1 Deterministic projections of a specific time length
One of the more common situations we might find ourselves in is wishing to conduct a deterministic projection using an arithmetic mean matrix. Perhaps we wish to project the mean forward by ten time steps simply to see the short-term predicted trend, or perhaps we wish to project many thousands of time steps in order to see the asymptotic properties of the matrix. Let’s start off by building arithmetic mean matrices for the ahistorical and historical population-level MPMs.
cypmean2 <- lmean(cypmatrix2r)
cypmean3 <- lmean(cypmatrix3r)
summary(cypmean2)
>
> This ahistorical lefkoMat object contains 1 matrix.
>
> Each matrix is square with 11 rows and columns, and a total of 121 elements.
> A total of 33 survival transitions were estimated, with 33 per matrix.
> A total of 10 fecundity transitions were estimated, with 10 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 0 time steps.
>
> The dataset contains a total of 74 unique individuals and 320 unique transitions.
>
> Survival probability sum check (each matrix represented by column in order):
> [,1]
> Min. 0.050
> 1st Qu. 0.140
> Median 0.800
> Mean 0.593
> 3rd Qu. 0.921
> Max. 1.000
summary(cypmean3)
>
> This historical lefkoMat object contains 1 matrix.
>
> Each matrix is square with 121 rows and columns, and a total of 14641 elements.
> A total of 99 survival transitions were estimated, with 99 per matrix.
> A total of 28 fecundity transitions were estimated, with 28 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 0 time steps.
>
> The dataset contains a total of 74 unique individuals and 320 unique transitions.
>
> Survival probability sum check (each matrix represented by column in order):
> [,1]
> Min. 0.000
> 1st Qu. 0.000
> Median 0.000
> Mean 0.179
> 3rd Qu. 0.200
> Max. 1.000Now let’s project these two forward. Since they have only single matrices, we do not need to indicate any specific year or patch to project.
cypproj2 <- projection3(cypmean2, times = 10)
> Warning: This function takes annual matrices as input. This lefkoMat object
> appears to be a set of mean matrices, and may lack annual matrices. Will project
> only the mean.
cypproj2
> $projection
> $projection[[1]]
> $projection[[1]][[1]]
> [,1] [,2] [,3] [,4] [,5] [,6]
> [1,] 1 1.025278e+04 8976.9428602 7579.2162919 6559.5257069 5791.2902872
> [2,] 1 1.025280e+04 9181.9984786 7758.7551491 6711.1100327 5922.4808013
> [3,] 1 1.000000e-01 1025.2800919 918.1998479 775.8755149 671.1110033
> [4,] 1 1.000000e-01 0.0100000 102.5280092 91.8199848 77.5875515
> [5,] 1 1.000000e-01 0.0100000 0.0010000 5.1264505 4.8473218
> [6,] 1 2.326441e-01 0.2158192 0.2037291 0.1945073 0.4237303
> [7,] 1 1.496794e+00 1.3301888 1.2290474 1.1663551 2.9546675
> [8,] 1 1.804320e+00 1.8827350 1.8582091 1.7999057 2.7390456
> [9,] 1 8.390066e-01 0.8995651 0.9136912 0.8912658 0.8505446
> [10,] 1 1.075116e+00 0.9570211 0.8349439 0.7303300 0.6449596
> [11,] 1 5.952381e-01 0.4000111 0.2911504 0.2250937 0.1821021
> [,7] [,8] [,9] [,10] [,11]
> [1,] 5777.0586419 6330.3074391 7144.4860609 8023.4195510 8864.9356819
> [2,] 5892.8844477 6445.8486120 7271.0922097 8166.3092722 9025.4040729
> [3,] 592.2480801 589.2884448 644.5848612 727.1092210 816.6309272
> [4,] 67.1111003 59.2248080 58.9288445 64.4584861 72.7109221
> [5,] 4.1217437 3.5616422 3.1393225 3.1034083 3.3780947
> [6,] 0.5946490 0.6972683 0.7576925 0.7932213 0.8270709
> [7,] 4.0352491 4.6021723 4.9080956 5.0742834 5.2639516
> [8,] 3.7343588 4.5180195 5.0926985 5.5058964 5.8572084
> [9,] 0.9816220 1.2047260 1.4431232 1.6600921 1.8434275
> [10,] 0.6345045 0.6788768 0.7530879 0.8390462 0.9253580
> [11,] 0.1524758 0.1366669 0.1329884 0.1382168 0.1490176
>
>
>
> $stage_dist
> $stage_dist[[1]]
> [1] 0
>
>
> $rep_value
> $rep_value[[1]]
> [1] 0
>
>
> $pop_size
> $pop_size[[1]]
> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
> [1,] 11 20511.92 19189.93 16364.03 14148.47 12475.11 12343.56 13440.07
> [,9] [,10] [,11]
> [1,] 15135.32 16998.41 18797.93
>
>
> $labels
> pop patch
> 1 1 1
>
> $ahstages
> stage_id stage original_size original_size_b original_size_c min_age max_age
> 1 1 SD 0.0 NA NA 0 NA
> 2 2 P1 0.0 NA NA 0 NA
> 3 3 P2 0.0 NA NA 0 NA
> 4 4 P3 0.0 NA NA 0 NA
> 5 5 SL 0.0 NA NA 0 NA
> 6 6 D 0.0 NA NA 0 NA
> 7 7 XSm 1.0 NA NA 0 NA
> 8 8 Sm 3.0 NA NA 0 NA
> 9 9 Md 6.0 NA NA 0 NA
> 10 10 Lg 11.0 NA NA 0 NA
> 11 11 XLg 19.5 NA NA 0 NA
> repstatus obsstatus propstatus immstatus matstatus entrystage indataset
> 1 0 0 1 0 0 1 0
> 2 0 0 0 1 0 1 0
> 3 0 0 0 1 0 0 0
> 4 0 0 0 1 0 0 0
> 5 0 0 0 1 0 0 0
> 6 0 0 0 0 1 0 1
> 7 1 1 0 0 1 0 1
> 8 1 1 0 0 1 0 1
> 9 1 1 0 0 1 0 1
> 10 1 1 0 0 1 0 1
> 11 1 1 0 0 1 0 1
> binhalfwidth_raw sizebin_min sizebin_max sizebin_center sizebin_width
> 1 0.0 0.0 0.0 0.0 0
> 2 0.0 0.0 0.0 0.0 0
> 3 0.0 0.0 0.0 0.0 0
> 4 0.0 0.0 0.0 0.0 0
> 5 0.0 0.0 0.0 0.0 0
> 6 0.5 -0.5 0.5 0.0 1
> 7 0.5 0.5 1.5 1.0 1
> 8 1.5 1.5 4.5 3.0 3
> 9 1.5 4.5 7.5 6.0 3
> 10 3.5 7.5 14.5 11.0 7
> 11 5.0 14.5 24.5 19.5 10
> binhalfwidthb_raw sizebinb_min sizebinb_max sizebinb_center sizebinb_width
> 1 NA NA NA NA NA
> 2 NA NA NA NA NA
> 3 NA NA NA NA NA
> 4 NA NA NA NA NA
> 5 NA NA NA NA NA
> 6 NA NA NA NA NA
> 7 NA NA NA NA NA
> 8 NA NA NA NA NA
> 9 NA NA NA NA NA
> 10 NA NA NA NA NA
> 11 NA NA NA NA NA
> binhalfwidthc_raw sizebinc_min sizebinc_max sizebinc_center sizebinc_width
> 1 NA NA NA NA NA
> 2 NA NA NA NA NA
> 3 NA NA NA NA NA
> 4 NA NA NA NA NA
> 5 NA NA NA NA NA
> 6 NA NA NA NA NA
> 7 NA NA NA NA NA
> 8 NA NA NA NA NA
> 9 NA NA NA NA NA
> 10 NA NA NA NA NA
> 11 NA NA NA NA NA
> group comments alive almostborn
> 1 0 Dormant seed 1 0
> 2 0 1st yr protocorm 1 0
> 3 0 2nd yr protocorm 1 0
> 4 0 3rd yr protocorm 1 0
> 5 0 Seedling 1 0
> 6 0 Dormant adult 1 0
> 7 0 Extra small adult (1 shoot) 1 0
> 8 0 Small adult (2-4 shoots) 1 0
> 9 0 Medium adult (5-7 shoots) 1 0
> 10 0 Large adult (8-14 shoots) 1 0
> 11 0 Extra large adult (>14 shoots) 1 0
>
> $hstages
> NA
> 1 NA
>
> $agestages
> NA
> 1 NA
>
> $control
> [1] 1 10
>
> attr(,"class")
> [1] "lefkoProj"
Here we can see the structure of the resulting object, which is of class lefkoProj. Here are the elements:
- projection - A list of with number of elements equal to the number of patches or populations. Each element is also a list with number of elements equal to the number of replicates (option
nreps, which defaults to 1). Each element of this lower-level list is a matrix with number of rows equal to the number of stages (if ahistorical; paired stages if historical; ages if Leslie; and age-stages if age-by-stage), and the number of columns equal to the projected time steps plus 1 (the first time is time 0). The elements of the matrix show the projected number of individuals in that stage in that time. - stage_dist - A list showing the stage distribution projected per time, with the same order and structure as projection. Only estimated if
growthonly = FALSE. - rep_value - A list showing the reproductive value of each stage projected per time, with the same order and structure as projection. Only estimated if
growthonly = FALSE. - pop_size - A list with number of elements equal to the number of projected populations or patches. Each element of this list is a matrix showing the total population size at each time (column) per replicate (row).
- ahstages - The stageframe used in analysis, though modified and edited per rules set for function-based matrix estimators.
- hstages - A data frame matrix showing the pairing of ahistorical stages used to create historical stage pairs (only if historical).
- agestages - A data frame showing the order of age-stage pairs (only if age-by-stage).
- labels - A data frame giving the population, patch, and year of each matrix in the
lefkoMatobject used as input, in order. - control - An integer vector indicating the number of replicates and the number of time steps.
- density - The data frame input under the density option. Only provided if input by the user.
Let’s also conduct a historical projection, and then compare summaries.
cypproj3 <- projection3(cypmean3, times = 10)
> Warning: This function takes annual matrices as input. This lefkoMat object
> appears to be a set of mean matrices, and may lack annual matrices. Will project
> only the mean.
summary(cypproj2)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 10 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 3 1
> 6 1
> 8 1
> 11 1
summary(cypproj3)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 10 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 3 1
> 6 1
> 8 1
> 11 1
The summary() calls give us a reasonable amount of information, showing the number of population-patches, projected time steps and number of replicates. We also see how many replicates still have population sizes above 1.0 at 0%, 25%, 50%, 75%, and 100% of the way through the projection. A number of options exist for summary() calls with projection3() and f_projection3(), and we encourage the user to explore upon reading the help pages for function summary.lefkoProj.
Let’s also plot the projections for comparison (figure 10.2).
par(mfrow=c(1,2))
plot(cypproj2)
title("a)", adj = 0)
plot(cypproj3, ylab = "")
title("b)", adj = 0)
Figure 10.2: Ahistorical (a) vs. historical (b) deterministic projection
These are fascinating results suggesting that accounting for history reduces the population growth rate in this example. Indeed, the ahistorical population projection shows growth, while the historical projection suggests a steeper initial decline followed by more gradual population growth.
A close look at the population sizes projected at each time shows that decimal values are allowed. Naturally, partial individuals cannot exist in a population. Although the default in lefko3 is to project fractional individuals, we may allow only integers, which results in the projected number of individuals to be rounded down to the nearest whole number. We can set this with the integeronly option, as below (figure 10.3).
cypproj3i <- projection3(cypmean3, times = 10, integeronly = TRUE)
> Warning: This function takes annual matrices as input. This lefkoMat object
> appears to be a set of mean matrices, and may lack annual matrices. Will project
> only the mean.
plot(cypproj3i)
Figure 10.3: Historical projection allowing integer population size only
We will not set integeronly = TRUE in the rest of this chapter, mostly because we are interested predominantly in the overall patterns of population dynamics, and forcing the projection of whole individuals makes extinction more and more likely with smaller population sizes. However, users interested in absolute population numbers should use this setting.
Before moving on to ordered / cyclical projections, let’s look at another setting of function projection3(): the start vector. Function projection3() assumes an initial starting vector of one individual in each stage, stage pair (if historical), or age-stage (if age-by-stage). We can see this by taking a look at the first row of the core ahistorical projection.
cypproj2$projection[[1]][[1]][,1]
> [1] 1 1 1 1 1 1 1 1 1 1 1
We see that we are starting with one individual of each stage. But what if we wished to use a different start vector? In that circumstance, we have two options. The first is to use the start_vec option to load a vector of starting numbers of each stage, stage pair, or age-stage. Here is an example where we start with 1000 seeds and 100 first-year protocorms.
cypproj2_1000 <- projection3(cypmean2, times = 10,
start_vec = c(1000, 100, 0, 0, 0, 0, 0, 0, 0, 0, 0))
> Warning: This function takes annual matrices as input. This lefkoMat object
> appears to be a set of mean matrices, and may lack annual matrices. Will project
> only the mean.
summary(cypproj2_1000)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 10 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 3 1
> 6 1
> 8 1
> 11 1
The start_vec approach is very useful, but it requires a vector as long as the number of rows in the projection matrix. This can make it very difficult to create the right start vector for a historical matrix, an age-by-stage matrix, or an IPM. For the latter options, we can use the start_input() function together with the start_frame option. We use start_input() to create a lefkoStart object, which is a data frame that contains everything we need to know for the start vector. Let’s see this approach in action, starting the population with 1000 3rd year protocorms and 100 dormant adults.
c2m_sv <- start_input(cypmean2, stage2 = c("P3", "D"), value = c(1000, 100))
c2m_sv
> stage2 stage_id_2 stage1 stage_id_1 age2 row_num value
> 1 P3 4 NA NA NA 4 1000
> 2 D 6 NA NA NA 6 100In the chunk above, we told R that we wished to start with 1000 3rdyear protocorms and 100 dormant adults, and zero individuals of all other stages. The result is this start frame, which is a data frame showing the stages with non-zero starting individuals.
Now that we have our start frame, let’s use it in a projection.
cypproj2_s1000 <- projection3(cypmean2, times = 10, start_frame = c2m_sv)
> Warning: This function takes annual matrices as input. This lefkoMat object
> appears to be a set of mean matrices, and may lack annual matrices. Will project
> only the mean.
summary(cypproj2_s1000)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 10 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 3 1
> 6 1
> 8 1
> 11 1We can compare these projections via plots (figure 10.4).
par(mfrow = c(1, 3))
plot(cypproj2)
title("a)", adj = 0)
plot(cypproj2_1000, ylab = "")
title("b)", adj = 0)
plot(cypproj2_s1000, ylab = "")
title("c)", adj = 0)
Figure 10.4: Using different start frames in projection. Projections started with (a) 1 of each stage, (b) 1000 dormant seeds and 100 seedlings set via start_vec, and (c) 1000 dormant seeds and 100 seedlings set via start_frame
The impact of a different start vector can be quite profound. The top left plot shows a population with a single individual of each stage, and so the population is able to grow right away. In contrast, the top right and bottom left plots show populations that initially contain no reproductive individuals. In the top left plot’s case, it takes a few years for the dormant seeds and 1st year protocorms to mature and the population to start growing, while in the bottom left plot, the population starts growing once some of the dormant individuals become reproductive, which only takes one or two years.
10.2.2 Ordered / cyclical projections of a specific time length
Situations may arise in which users will want to project an MPM with a specific order of matrices in mind. Typically these situations arise when a specific order of annual matrices is meant to correspond to assumptions of predicted conditions. For example, we might have several years worth of matrices, and may believe due to climate forecasts that the next ten years should be relatively rainy while the ten years after that should be drier. Or we may have a situation in which a hurricane has affected our population, and we wish to project the population forward avoiding matrices in the immediate aftermath of the hurricane, since those matrices likely reflect dramatic changes as the population adjusts to the mass mortality and other issues that might have destabilized population dynamics during the hurricane itself.
Let’s project the population forward, but this time let’s force the projection to cycle through the matrices for 2007 and 2008 for ten years, and then to cycle between 2004 and 2005 for another ten years (2004 and 2005 were relatively dry years, while 2007 and 2008 were relatively wet years). Our projection will be 20 years in total.
year_order <- c(2007, 2008, 2007, 2008, 2007, 2008, 2007, 2008, 2007, 2008,
2004, 2005, 2004, 2005, 2004, 2005, 2004, 2005, 2004, 2005)
cypproj2_ord1 <- projection3(cypmatrix2r, year = year_order, times = 20,
start_frame = c2m_sv)
summary(cypproj2_ord1)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 20 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 6 1
> 11 1
> 16 1
> 21 1
If we have an MPM with a set of annual matrices, we might wish to project the population forward cyclically through those matrices. If the number of matrices in the lefkoMat object is shorter than the number of time steps set, then the default behavior of lefko3’s projection functions is to cycle through the MPMs. For example, if we leave the times option out from the above projection, then lefko3 will assume the default of 10,000 time steps, and cycle through our 20 defined years over and over, as below.
cypproj2_ord2 <- projection3(cypmatrix2r, year = year_order,
start_frame = c2m_sv)
summary(cypproj2_ord2)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 10000 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 2501 1
> 5001 1
> 7501 1
> 10001 1Let’s see plots of these projections (figure 10.5).
par(mfrow=c(1,2))
plot(cypproj2_ord1)
title("a)", adj = 0)
plot(cypproj2_ord2, ylab = "")
title("b)", adj = 0)
Figure 10.5: Ordered projections can be cyclical. (a) A 20 year ordered projection. (b) A 10,000 year projection cycling through 20 ordered years
Of course, the scale of the y axis differs dramatically, because the time scale differs so dramatically. The population is generally increasing, and so reaches much greater numbers in 10,000 time steps than in 20.
If the year option is not used and stochastic = FALSE, then the default behavior of this function is to cycle through all matrices in a lefkoMat object in order. Here, we will cycle through all of the matrices in the ahistorical and historical population-level MPMs for 100 time steps each.
cypproj2_c100 <- projection3(cypmatrix2r, times = 100)
cypproj3_c100 <- projection3(cypmatrix3r, times = 100)
summary(cypproj2_c100)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 100 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
summary(cypproj3_c100)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 100 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1Now let’s plot the projections for comparison (figure 10.6).
par(mfrow=c(1,2))
plot(cypproj2_c100)
title("a)", adj = 0)
plot(cypproj3_c100, ylab = "")
title("b)", adj = 0)
Figure 10.6: Cyclical ahistorical (a) and historical (b) projections
Here we see that individual history appears to reduce the growth rate relative to ahistorical projection, although clearly there are years in which the population growth rate appears quite high in both cases. The cyclical nature of the projection is obvious in both plots.
We might also wish to see what happens among different patches. Let’s try a projection with patch structure, here using cypmatrix3rp. This MPM contains four matrices each for three patches.
cypproj3p_c100 <- projection3(cypmatrix3rp, times = 100)
summary(cypproj3p_c100)
>
> The input lefkoProj object covers 4 population-patches.
> It includes 100 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 A 1 B 1 C 1 0
> 1 1 1 1 1
> 26 1 1 1 1
> 51 1 1 1 1
> 76 1 1 1 1
> 101 0 1 1 1
Function projection3() has projected the matrix sets from each patch and the patch-weighted population in an ordered cycle for 100 time steps each. The summary shows us that patch A may go extinct some time between time step 76 and time step 100. The population-level projection is the rightmost, and shows no extinction.
Let’s now plot the results.
par(mfrow=c(2,2))
plot(cypproj3p_c100, patch = "all")
Figure 10.7: Patch-level projections for patch A (top left), patch B (top right), patch C (bottom left), and the patch-weighted population (bottom right)
The plots are made in the order of the labels element in the lefkoProj object, which we see below. The plots are propagated by row, meaning that patch A is in the top-left, patch B is in the top-right, patch C is in the bottom-left, and the patch-weighted population mean is in the bottom-right, as shown in the labels element below. These results differ from the population-level projection because the population-level results from a patch-stratified MPM involve averaging the patches equally, whereas a population-level projection weights individual transitions as equal regardless of source patch.
cypproj3p_c100$labels
> pop patch
> 1 1 A
> 2 1 B
> 3 1 C
> 4 1 010.2.3 Replicated simulations assuming temporal stochasticity
Let’s now turn our attention to stochastic simulations. Package lefko3 currently implements temporal environmental stochasticity through its two projection functions. In projection3(), this means that the matrices in a MPM can be shuffled. Let’s conduct a simple population-level stochastic simulation with 100 replicates and 10,000 time steps (note that 10,000 time steps is the default). We will then summarize the results.
set.seed(42)
cypstoch2 <- projection3(cypmatrix2r, nreps = 100, stochastic = TRUE)
summary(cypstoch2)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 10000 projected steps per replicate and 100 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 100
> 2501 100
> 5001 100
> 7501 100
> 10001 0
The summary shows us that we have one population-patch covered in this projection, and that it includes 10,000 time steps and a single replicate. The final portion of the summary shows the number of replicates yielding a population size above a threshold value at particular milepost times in the projection. The default threshold population value is 1.0, and the default milepost times are the times once 0%, 25%, 50%, 75%, and 100% of the projection has been run. We see that extinction may have occurred in every replicate some time after time step 7,501. However, this is unlikely, as we know from the cyclical analysis that the ahistorical MPM suggests a growing population in the long-term. Let’s take a look at what might be going on by examining the pop_size elements in the first ten replicates right after time step 7,500.
cypstoch2$pop_size[[1]][c(1:10), c(7501:7503)]
> [,1] [,2] [,3]
> [1,] 8.321329e+252 3.753757e+252 1.269240e+253
> [2,] 1.612221e+253 2.310869e+253 1.539012e+253
> [3,] 3.266403e+257 3.040670e+257 3.963822e+257
> [4,] 4.595511e+257 3.245229e+257 4.734050e+257
> [5,] 1.276522e+256 1.461003e+256 6.718161e+255
> [6,] 7.869332e+253 3.224921e+253 2.112759e+254
> [7,] 3.699333e+256 5.988414e+256 1.572965e+257
> [8,] 1.049118e+254 1.465183e+254 1.135257e+254
> [9,] 1.265153e+258 5.958920e+257 6.857992e+257
> [10,] 1.443963e+259 1.299965e+259 3.667100e+258
We see that the population has exploded in size, and so the zeros in the milepost section of the summary occur because the population sizes projected became too high for our computers to handle (to see the highest number you can handle, type .Machine$double.xmax at the prompt). To see this, let’s view a plot (note the ylim values - exceeding the maximum allowed number in our computers means that we need to designate appropriate limits to the y-axis or else the plot() function will produce an error).
plot(cypstoch2, ylim = c(0, 10000000000000000000000000000000000000000000000000))
Figure 10.8: Stochastic projection with exploding population size
How does this compare to a historical projection? Let’s limit the number of time steps to 1000 to keep this quick and then take a look at a summary.
cypstoch3 <- projection3(cypmatrix3r, nreps = 100, times = 1000,
stochastic = TRUE)
summary(cypstoch3)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 1000 projected steps per replicate and 100 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 100
> 251 100
> 501 100
> 751 100
> 1001 100We see neither extinction, nor computationally problematic population growth. Let’s plot this projection against the preceding projection, but standardizing the scales of the x and y axes to compare properly (figure 10.9).
par(mfrow=c(1, 2))
plot(cypstoch2, xlim = c(0, 1001), ylim = c(0, 35000000000))
title("a)", adj = 0)
plot(cypstoch3, ylim = c(0, 35000000000), ylab = "")
title("b)", adj = 0)
Figure 10.9: Stochastic ahistorical (a) vs. historical (b) projection
All of our replicates survived and increased in the historical case, but more moderately than in the ahistorical case.
Let’s try another projection, using the patch-level historical MPM. We will also set growthonly = FALSE to allow the stage structure and reproductive values to be tracked.
cypstoch2p <- projection3(cypmatrix2rp, stochastic = TRUE, nreps = 100,
times = 200, growthonly = FALSE)
summary(cypstoch2p)
>
> The input lefkoProj object covers 4 population-patches.
> It includes 200 projected steps per replicate and 100 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 A 1 B 1 C 1 0
> 1 100 100 100 100
> 51 96 100 100 100
> 101 82 100 100 100
> 151 79 100 100 100
> 201 68 100 100 100
Our summary shows us that we have three patches and the full (patch-weighted) population projected. There are 100 replicates of each. We can call each patch using the first set of double brackets for the projection, stage_dist, rep_value, and pop_size elements of the resulting lefkoProj object, and each replicate within each patch using the second set of brackets. Here, we will take a look at the overall stage distribution at the population level in the first replicate at the 100th time by focusing on the first element within the fourth list, and the 100th column of that element. If we wanted to see the equivalent stage distribution for the 100th replicate, we would adjust the second set of brackets accordingly.
writeLines("Replicate 1:")
> Replicate 1:
cypstoch2p$stage_dist[[4]][[1]][,100]
> [1] 4.669882e-01 4.774991e-01 5.292660e-02 1.939510e-03 2.011212e-04
> [6] 2.086144e-05 1.942072e-04 1.782106e-04 5.180013e-05 3.739271e-07
> [11] 0.000000e+00
writeLines("\nReplicate 100:")
>
> Replicate 100:
cypstoch2p$stage_dist[[4]][[100]][,100]
> [1] 4.778352e-01 4.847863e-01 3.508813e-02 1.678213e-03 5.298717e-05
> [6] 2.896109e-05 2.657260e-04 2.070924e-04 5.473344e-05 2.624975e-06
> [11] 0.000000e+00These are fractional individuals, suggesting extinction.
Let’s conduct one more projection, this time as a simulation study of the historical patch-level MPM with 100 replicates and 1000 time steps each. Note that this might take a little while to run.
cypstoch3p <- projection3(cypmatrix3rp, nreps = 100, times = 1000,
stochastic = TRUE, growthonly = FALSE)Let’s see a summary.
summary(cypstoch3p)
>
> The input lefkoProj object covers 4 population-patches.
> It includes 1000 projected steps per replicate and 100 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 A 1 B 1 C 1 0
> 1 100 100 100 100
> 251 0 1 7 0
> 501 0 0 0 0
> 751 0 0 0 0
> 1001 0 0 0 0
The new output shows us how different predictions can be between historical vs. ahistorical MPMs, as well as between population-level and patch-level MPMs. The ahistorical population-level MPMs showed explosively growing population size, while the historical version showed moderate, fluctuating growth, often still dropping quite low. When run with patch-level MPMs, extinction was predicted either way. To see how quickly extinction happens, let’s change the milepost settings in the summary.lefkoProj() function, as below.
summary(cypstoch3p, milepost= c(1, 10, 20, 25, 30, 50, 100, 150, 200, 250, 300))
>
> The input lefkoProj object covers 4 population-patches.
> It includes 1000 projected steps per replicate and 100 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 A 1 B 1 C 1 0
> 1 100 100 100 100
> 10 95 100 100 100
> 20 68 100 100 100
> 25 55 100 100 100
> 30 55 100 100 100
> 50 16 100 99 95
> 100 0 83 83 1
> 150 0 24 48 0
> 200 0 6 28 0
> 250 0 1 8 0
> 300 0 0 4 0Let’s zoom in on the historical population-level’s 100th replicate, and take a peek at the first 100 population sizes.
cypstoch3p$pop_size[[4]][100,c(1:100)]
> [1] 121.000000 46179.502222 29023.695555 22837.294912 7033.629303
> [6] 8990.143963 11966.585249 3319.758097 3808.467930 1781.964874
> [11] 10204.204434 7451.053601 5378.446872 4460.233720 1504.415848
> [16] 7100.823513 2037.357258 1445.202800 4973.835482 1714.938486
> [21] 2290.734278 2402.992611 882.279525 429.399209 635.994069
> [26] 3496.124971 831.284762 634.363016 678.903569 569.063124
> [31] 2335.717319 1749.165617 1257.248280 369.835533 894.208457
> [36] 320.150773 454.546602 795.335980 570.544055 183.458250
> [41] 162.247241 146.381632 184.258993 962.803243 316.705579
> [46] 919.381047 764.787634 217.100369 86.073363 340.554620
> [51] 447.166229 138.589790 61.565752 26.710117 12.080205
> [56] 102.601484 23.751709 13.044780 53.367519 45.137515
> [61] 52.692054 53.175603 13.297251 9.005166 10.275935
> [66] 5.789382 25.903322 33.408648 32.518226 29.207798
> [71] 8.125889 28.556519 7.356006 30.927875 10.193891
> [76] 3.595580 2.848797 4.257870 19.892382 19.694320
> [81] 5.359228 18.696500 5.240442 10.490754 3.513380
> [86] 13.683151 9.544740 2.760140 2.312355 2.830370
> [91] 1.283799 6.145632 5.408657 1.914644 1.364086
> [96] 1.244728 3.346333 4.676133 3.440483 1.125806
Each patch or population’s pop_size element is a matrix with rows equal to the number of replicates and columns equal to the number of times projected. In contrast, the projection, stage_dist, and rep_value elements are lists nested within lists, where the top level is the patch / population, and the bottom level is the replicate. Within the replicate, there is a matrix with rows equal to the number of stages or stage-pairs, and the columns equal to the number of times projected. Let’s take a look at the first five time steps in the 100th replicate of the population level, for each of these three elements.
cypstoch3p$projection[[4]][[100]][,c(1:5)]
> [,1] [,2] [,3] [,4] [,5]
> [1,] 1 0.5800000 2.308247e+03 1.289094e+03 1.035021e+03
> [2,] 1 0.6000000 2.308259e+03 1.335259e+03 1.060803e+03
> [3,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [4,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [5,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [6,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [7,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [8,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [9,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [10,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [11,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [12,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [13,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [14,] 1 0.6000000 2.308261e+03 1.335260e+03 1.065419e+03
> [15,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [16,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [17,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [18,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [19,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [20,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [21,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [22,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [23,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [24,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [25,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [26,] 1 0.1000000 6.000000e-02 2.308261e+02 1.335260e+02
> [27,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [28,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [29,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [30,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [31,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [32,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [33,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [34,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [35,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [36,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [37,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [38,] 1 0.0500000 5.000000e-03 3.000000e-03 1.154131e+01
> [39,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [40,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [41,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [42,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [43,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [44,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [45,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [46,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [47,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [48,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [49,] 1 0.1000000 7.500000e-03 6.250000e-04 1.812500e-04
> [50,] 1 0.1444444 1.083333e-02 9.027778e-04 0.000000e+00
> [51,] 1 1.0333333 7.750000e-02 6.458333e-03 1.127778e-03
> [52,] 1 0.1555556 1.166667e-02 9.722222e-04 1.409722e-03
> [53,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [54,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [55,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [56,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [57,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [58,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [59,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [60,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [61,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [62,] 1 0.0000000 0.000000e+00 0.000000e+00 3.316353e-01
> [63,] 1 0.6666667 1.169312e-01 1.157111e-01 1.187400e-01
> [64,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [65,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [66,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [67,] 1 119.0476190 2.650227e+02 3.723882e+02 0.000000e+00
> [68,] 1 119.0476190 2.650227e+02 3.723882e+02 0.000000e+00
> [69,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [70,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [71,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [72,] 1 0.2063492 3.363001e-01 3.307325e-01 0.000000e+00
> [73,] 1 2.2261905 3.128061e+00 3.153509e+00 2.012200e+00
> [74,] 1 0.3888889 5.226631e-01 5.311735e-01 2.109449e+00
> [75,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [76,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [77,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [78,] 1 2546.2962963 2.515432e+03 2.314603e+03 0.000000e+00
> [79,] 1 2546.2962963 2.515432e+03 2.314603e+03 0.000000e+00
> [80,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [81,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [82,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [83,] 1 0.6666667 1.814815e-01 1.781099e-01 0.000000e+00
> [84,] 1 0.9629630 1.050000e+00 1.015436e+00 3.547638e-02
> [85,] 1 2.1666667 2.232716e+00 2.046362e+00 1.847207e+00
> [86,] 1 0.5000000 2.777778e-01 1.481481e-01 4.629739e-01
> [87,] 1 0.0000000 0.000000e+00 0.000000e+00 1.219136e-01
> [88,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [89,] 1 7083.3333333 7.152778e+03 6.631944e+03 1.859568e+03
> [90,] 1 7083.3333333 7.152778e+03 6.631944e+03 1.859568e+03
> [91,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [92,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [93,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [94,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [95,] 1 0.1666667 8.333333e-02 4.629630e-02 0.000000e+00
> [96,] 1 1.0000000 6.666667e-01 3.657407e-01 0.000000e+00
> [97,] 1 1.1666667 1.138889e+00 1.041667e+00 7.932099e-01
> [98,] 1 0.0000000 0.000000e+00 0.000000e+00 3.472222e-01
> [99,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [100,] 1 5833.3333333 5.555556e+02 0.000000e+00 0.000000e+00
> [101,] 1 5833.3333333 5.555556e+02 0.000000e+00 0.000000e+00
> [102,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [103,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [104,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [105,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [106,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [107,] 1 0.3333333 1.111111e-01 0.000000e+00 0.000000e+00
> [108,] 1 0.6666667 1.111111e-01 0.000000e+00 0.000000e+00
> [109,] 1 0.3333333 0.000000e+00 0.000000e+00 0.000000e+00
> [110,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [111,] 1 7500.0000000 5.555556e+02 0.000000e+00 0.000000e+00
> [112,] 1 7500.0000000 5.555556e+02 0.000000e+00 0.000000e+00
> [113,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [114,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [115,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [116,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [117,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [118,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [119,] 1 0.0000000 0.000000e+00 0.000000e+00 0.000000e+00
> [120,] 1 0.3333333 1.111111e-01 0.000000e+00 0.000000e+00
> [121,] 1 0.3333333 0.000000e+00 0.000000e+00 0.000000e+00
cypstoch3p$stage_dist[[4]][[100]][,c(1:5)]
> [,1] [,2] [,3] [,4] [,5]
> [1,] 0.008264463 1.255968e-05 7.952976e-02 5.644689e-02 1.471532e-01
> [2,] 0.008264463 1.299278e-05 7.953016e-02 5.846836e-02 1.508187e-01
> [3,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [4,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [5,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [6,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [7,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [8,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [9,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [10,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [11,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [12,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [13,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [14,] 0.008264463 1.299278e-05 7.953023e-02 5.846841e-02 1.514751e-01
> [15,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [16,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [17,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [18,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [19,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [20,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [21,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [22,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [23,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [24,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [25,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [26,] 0.008264463 2.165463e-06 2.067276e-06 1.010742e-02 1.898394e-02
> [27,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [28,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [29,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [30,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [31,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [32,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [33,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [34,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [35,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [36,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [37,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [38,] 0.008264463 1.082731e-06 1.722730e-07 1.313641e-07 1.640875e-03
> [39,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [40,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [41,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [42,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [43,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [44,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [45,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [46,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [47,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [48,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [49,] 0.008264463 2.165463e-06 2.584095e-07 2.736751e-08 2.576906e-08
> [50,] 0.008264463 3.127891e-06 3.732582e-07 3.953085e-08 0.000000e+00
> [51,] 0.008264463 2.237645e-05 2.670232e-06 2.827976e-07 1.603408e-07
> [52,] 0.008264463 3.368498e-06 4.019704e-07 4.257169e-08 2.004260e-07
> [53,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [54,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [55,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [56,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [57,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [58,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [59,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [60,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [61,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [62,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 4.714995e-05
> [63,] 0.008264463 1.443642e-05 4.028819e-06 5.066762e-06 1.688175e-05
> [64,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [65,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [66,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [67,] 0.008264463 2.577932e-03 9.131252e-03 1.630614e-02 0.000000e+00
> [68,] 0.008264463 2.577932e-03 9.131252e-03 1.630614e-02 0.000000e+00
> [69,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [70,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [71,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [72,] 0.008264463 4.468416e-06 1.158709e-05 1.448212e-05 0.000000e+00
> [73,] 0.008264463 4.820733e-05 1.077761e-04 1.380859e-04 2.860828e-04
> [74,] 0.008264463 8.421245e-06 1.800815e-05 2.325904e-05 2.999090e-04
> [75,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [76,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [77,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [78,] 0.008264463 5.513910e-02 8.666822e-02 1.013519e-01 0.000000e+00
> [79,] 0.008264463 5.513910e-02 8.666822e-02 1.013519e-01 0.000000e+00
> [80,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [81,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [82,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [83,] 0.008264463 1.443642e-05 6.252873e-06 7.799082e-06 0.000000e+00
> [84,] 0.008264463 2.085261e-05 3.617734e-05 4.446395e-05 5.043823e-06
> [85,] 0.008264463 4.691836e-05 7.692735e-05 8.960614e-05 2.626250e-04
> [86,] 0.008264463 1.082731e-05 9.570724e-06 6.487115e-06 6.582290e-05
> [87,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 1.733296e-05
> [88,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [89,] 0.008264463 1.533870e-01 2.464461e-01 2.903997e-01 2.643824e-01
> [90,] 0.008264463 1.533870e-01 2.464461e-01 2.903997e-01 2.643824e-01
> [91,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [92,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [93,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [94,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [95,] 0.008264463 3.609105e-06 2.871217e-06 2.027223e-06 0.000000e+00
> [96,] 0.008264463 2.165463e-05 2.296974e-05 1.601506e-05 0.000000e+00
> [97,] 0.008264463 2.526373e-05 3.923997e-05 4.561252e-05 1.127739e-04
> [98,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 4.936601e-05
> [99,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [100,] 0.008264463 1.263187e-01 1.914145e-02 0.000000e+00 0.000000e+00
> [101,] 0.008264463 1.263187e-01 1.914145e-02 0.000000e+00 0.000000e+00
> [102,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [103,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [104,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [105,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [106,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [107,] 0.008264463 7.218210e-06 3.828290e-06 0.000000e+00 0.000000e+00
> [108,] 0.008264463 1.443642e-05 3.828290e-06 0.000000e+00 0.000000e+00
> [109,] 0.008264463 7.218210e-06 0.000000e+00 0.000000e+00 0.000000e+00
> [110,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [111,] 0.008264463 1.624097e-01 1.914145e-02 0.000000e+00 0.000000e+00
> [112,] 0.008264463 1.624097e-01 1.914145e-02 0.000000e+00 0.000000e+00
> [113,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [114,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [115,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [116,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [117,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [118,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [119,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [120,] 0.008264463 7.218210e-06 3.828290e-06 0.000000e+00 0.000000e+00
> [121,] 0.008264463 7.218210e-06 0.000000e+00 0.000000e+00 0.000000e+00
cypstoch3p$rep_value[[4]][[100]][,c(1:5)]
> [,1] [,2] [,3] [,4] [,5]
> [1,] 9.944209e-07 1.003321e-06 1.132723e-06 2.147589e-06 2.400209e-06
> [2,] 7.388011e-06 9.315816e-06 9.349685e-06 9.759336e-06 1.907842e-05
> [3,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [4,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [5,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [6,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [7,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [8,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [9,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [10,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [11,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [12,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [13,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [14,] 6.554559e-05 7.517480e-05 9.522544e-05 9.473648e-05 9.542433e-05
> [15,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [16,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [17,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [18,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [19,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [20,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [21,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [22,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [23,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [24,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [25,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [26,] 6.681412e-04 6.669422e-04 7.684302e-04 9.648799e-04 9.263094e-04
> [27,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [28,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [29,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [30,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [31,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [32,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [33,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [34,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [35,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [36,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [37,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [38,] 1.367068e-02 1.359700e-02 1.363485e-02 1.557237e-02 1.886870e-02
> [39,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [40,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [41,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [42,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [43,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [44,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [45,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [46,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [47,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [48,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [49,] 1.367068e-02 1.359700e-02 1.363485e-02 1.557237e-02 1.886870e-02
> [50,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [51,] 2.008882e-02 2.019836e-02 1.838060e-02 2.086784e-02 2.814545e-02
> [52,] 3.041762e-02 3.593027e-02 4.783299e-02 2.899659e-02 1.421092e-02
> [53,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [54,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [55,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [56,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [57,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [58,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [59,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [60,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [61,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [62,] 0.000000e+00 0.000000e+00 0.000000e+00 9.595070e-03 0.000000e+00
> [63,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 9.630127e-03
> [64,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [65,] 0.000000e+00 0.000000e+00 0.000000e+00 1.628174e-02 0.000000e+00
> [66,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [67,] 1.014142e-06 1.025484e-06 1.175112e-06 2.196685e-06 2.425730e-06
> [68,] 7.388011e-06 9.315816e-06 9.349685e-06 9.759336e-06 1.907842e-05
> [69,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [70,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [71,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [72,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [73,] 2.129873e-02 2.142412e-02 1.978531e-02 2.086784e-02 2.814545e-02
> [74,] 3.606384e-02 4.165045e-02 5.438832e-02 2.899659e-02 1.421092e-02
> [75,] 0.000000e+00 0.000000e+00 0.000000e+00 9.788790e-02 0.000000e+00
> [76,] 1.122580e-02 2.183984e-03 0.000000e+00 0.000000e+00 0.000000e+00
> [77,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [78,] 1.014142e-06 1.025484e-06 1.175112e-06 2.196685e-06 2.425730e-06
> [79,] 7.388011e-06 9.315816e-06 9.349685e-06 9.759336e-06 1.907842e-05
> [80,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [81,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [82,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [83,] 0.000000e+00 0.000000e+00 0.000000e+00 6.566015e-03 0.000000e+00
> [84,] 2.261354e-02 2.338477e-02 2.021562e-02 2.163602e-02 3.512915e-02
> [85,] 5.041651e-02 5.794661e-02 7.252931e-02 9.330865e-02 9.753824e-02
> [86,] 1.455068e-01 1.528321e-01 1.085502e-01 2.141032e-01 5.971789e-02
> [87,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [88,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [89,] 1.014142e-06 1.025484e-06 1.175112e-06 2.196685e-06 2.425730e-06
> [90,] 7.388011e-06 9.315816e-06 9.349685e-06 9.759336e-06 1.907842e-05
> [91,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [92,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [93,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [94,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [95,] 0.000000e+00 0.000000e+00 0.000000e+00 1.443972e-02 4.055609e-02
> [96,] 7.106409e-02 7.258836e-02 1.050841e-01 2.035843e-02 8.155862e-02
> [97,] 2.991879e-01 2.835810e-01 2.584208e-01 2.324307e-01 2.871370e-01
> [98,] 5.725843e-02 5.288144e-02 6.055035e-02 0.000000e+00 0.000000e+00
> [99,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [100,] 1.014142e-06 1.025484e-06 1.175112e-06 2.196685e-06 2.425730e-06
> [101,] 7.388011e-06 9.315816e-06 9.349685e-06 9.759336e-06 1.907842e-05
> [102,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [103,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [104,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [105,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [106,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [107,] 5.853599e-02 4.397815e-02 8.026688e-02 3.325177e-02 0.000000e+00
> [108,] 2.377941e-02 3.426758e-02 6.697348e-03 0.000000e+00 1.088939e-01
> [109,] 4.257138e-02 4.551928e-02 3.060485e-02 3.460920e-02 4.775964e-02
> [110,] 6.490698e-02 6.647482e-02 6.883086e-02 7.352662e-02 0.000000e+00
> [111,] 1.014142e-06 1.025484e-06 1.175112e-06 2.196685e-06 2.425730e-06
> [112,] 7.388011e-06 9.315816e-06 9.349685e-06 9.759336e-06 1.907842e-05
> [113,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [114,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [115,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [116,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [117,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [118,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [119,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [120,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 2.752000e-02
> [121,] 1.693866e-02 1.716052e-02 1.966596e-02 0.000000e+00 8.095853e-02This is a huge deal of data, and could be of great use to ecologists studying changes in population structure, natural selection, etc.
The output for our historical MPM suggests that the population will decline quite quickly. Let’s see how quickly by plotting the first 100 projected time steps of each of the 100 replicates, as below (figure 10.10).
par(mfrow = c(2,2))
plot(cypstoch3p, xlim = c(1, 100), patch = "all")
Figure 10.10: Stochastic patch-level projections
We see very quick declines in patch A and the population as a whole, though patches B and C appear capable of surviving some time longer in some replicates.
Users running stochastic projections may wish to weight the probability of sampling each matrix differently. This can be done using the tweights option, which takes as its input a vector giving the respective weight of each matrix. Users may fill this vector with probabilities, or may use integers or other numbers to provide a relative weighting. In the latter case, function projection3() will calculate the probability as the vector element divided by the vector sum. Let’s project our historical population-level MPM assuming that the last matrix is chosen 91% of the time while the other three matrices are chosen 3% of the time each.
cypstoch3_weights <- projection3(cypmatrix3r, nreps = 100, times = 200,
stochastic = TRUE, tweights = c(0.03, 0.03, 0.03, 0.91))Now let’s plot the results, comparing them to the original projection with equal weighting (figure 10.11).
par(mfrow = c(1,2))
plot(cypstoch3)
title("a)", adj = 0)
plot(cypstoch3_weights, ylab = "")
title("b)", adj = 0)
Figure 10.11: Stochastic projections using equal (a) or different (b) matrix weights
Choosing the first matrix 91% of the time may lead to greater population growth.
10.2.4 Density dependent deterministic projections
Let’s now move on to the simplest density dependent analysis: a single matrix projected forward and altered at each time step according to some function of population size, \(N\). This matrix will have particular elements that will be altered according to \(N\) at each time step. We might do this if we have an inkling that density has some role in population dynamics. For example, we might suspect that natural selection is at work in the population, and natural selection often works under the assumption of conspecific competition. Alternatively, we might believe that the stability of the population dynamics that we see might be due to a limit on the number of individuals that a community can support. In other cases we may see actual population dynamics that suggest a role for density, such as cycles of different periods.
Let’s start our analysis by projecting the mean ahistorical matrix forward using four different parameterizations of the Ricker function modifying germination. Our first step will be to create four density frames (lefkoDens objects), which will describe the relationships between the matrices and density.
c2d_1 <- density_input(cypmean2, stage3 = c("P1", "P1"), stage2= c("SD", "rep"),
style = 1, time_delay = 1, alpha = 0.05, beta = 0, type = c(2, 2))
c2d_2 <- density_input(cypmean2, stage3 = c("P1", "P1"), stage2= c("SD", "rep"),
style = 1, time_delay = 1, alpha = 1, beta = 0, type = c(2, 2))
c2d_3 <- density_input(cypmean2, stage3 = c("P1", "P1"), stage2= c("SD", "rep"),
style = 1, time_delay = 1, alpha = 0.05, beta = 0.0005, type = c(2, 2))
c2d_4 <- density_input(cypmean2, stage3 = c("P1", "P1"), stage2= c("SD", "rep"),
style = 1, time_delay = 1, alpha = 1, beta = 0.0005, type = c(2, 2))
c2d_1
> stage3 stage2 stage1 age2 style time_delay alpha beta type type_t12
> 1 P1 SD <NA> <NA> 1 1 0.05 0 2 NA
> 2 P1 XSm <NA> <NA> 1 1 0.05 0 2 NA
> 3 P1 Sm <NA> <NA> 1 1 0.05 0 2 NA
> 4 P1 Md <NA> <NA> 1 1 0.05 0 2 NA
> 5 P1 Lg <NA> <NA> 1 1 0.05 0 2 NA
> 6 P1 XLg <NA> <NA> 1 1 0.05 0 2 NA
c2d_2
> stage3 stage2 stage1 age2 style time_delay alpha beta type type_t12
> 1 P1 SD <NA> <NA> 1 1 1 0 2 NA
> 2 P1 XSm <NA> <NA> 1 1 1 0 2 NA
> 3 P1 Sm <NA> <NA> 1 1 1 0 2 NA
> 4 P1 Md <NA> <NA> 1 1 1 0 2 NA
> 5 P1 Lg <NA> <NA> 1 1 1 0 2 NA
> 6 P1 XLg <NA> <NA> 1 1 1 0 2 NA
c2d_3
> stage3 stage2 stage1 age2 style time_delay alpha beta type type_t12
> 1 P1 SD <NA> <NA> 1 1 0.05 5e-04 2 NA
> 2 P1 XSm <NA> <NA> 1 1 0.05 5e-04 2 NA
> 3 P1 Sm <NA> <NA> 1 1 0.05 5e-04 2 NA
> 4 P1 Md <NA> <NA> 1 1 0.05 5e-04 2 NA
> 5 P1 Lg <NA> <NA> 1 1 0.05 5e-04 2 NA
> 6 P1 XLg <NA> <NA> 1 1 0.05 5e-04 2 NA
c2d_4
> stage3 stage2 stage1 age2 style time_delay alpha beta type type_t12
> 1 P1 SD <NA> <NA> 1 1 1 5e-04 2 NA
> 2 P1 XSm <NA> <NA> 1 1 1 5e-04 2 NA
> 3 P1 Sm <NA> <NA> 1 1 1 5e-04 2 NA
> 4 P1 Md <NA> <NA> 1 1 1 5e-04 2 NA
> 5 P1 Lg <NA> <NA> 1 1 1 5e-04 2 NA
> 6 P1 XLg <NA> <NA> 1 1 1 5e-04 2 NA
Determining the values of \(\alpha\) and \(\beta\) to use for the alpha and beta options generally takes a great deal of extra work. We have done this work separately and present only four interesting combinations that produce different, interesting results in this case. Note also that the style, alpha, and beta terms can be vectors if necessary, giving different kinds of density dependence and different values for \(\alpha\) and \(\beta\) when needed.
Now let’s conduct some short simulations. We will only project forward 100 time steps, and we will enforce full substochasticity.
cypproj2_d1 <- projection3(cypmean2, times = 100, substoch = 2, density = c2d_1)
> Warning: This function takes annual matrices as input. This lefkoMat object
> appears to be a set of mean matrices, and may lack annual matrices. Will project
> only the mean.
cypproj2_d2 <- projection3(cypmean2, times = 100, substoch = 2, density = c2d_2)
> Warning: This function takes annual matrices as input. This lefkoMat object
> appears to be a set of mean matrices, and may lack annual matrices. Will project
> only the mean.
cypproj2_d3 <- projection3(cypmean2, times = 100, substoch = 2, density = c2d_3)
> Warning: This function takes annual matrices as input. This lefkoMat object
> appears to be a set of mean matrices, and may lack annual matrices. Will project
> only the mean.
cypproj2_d4 <- projection3(cypmean2, times = 100, substoch = 2, density = c2d_4)
> Warning: This function takes annual matrices as input. This lefkoMat object
> appears to be a set of mean matrices, and may lack annual matrices. Will project
> only the mean.
summary(cypproj2_d1)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 100 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
summary(cypproj2_d2)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 100 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
summary(cypproj2_d3)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 100 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
summary(cypproj2_d4)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 100 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1Let’s plot these projections together and take a look (figure 10.12).
par(mfrow = c(2, 2))
plot(cypproj2_d1)
title("a)", adj = 0)
plot(cypproj2_d2, ylab = "")
title("b)", adj = 0)
plot(cypproj2_d3)
title("c)", adj = 0)
plot(cypproj2_d4, ylab = "")
title("d)", adj = 0)
Figure 10.12: Density dependent projections using the two-parameter Ricker function. (a) alpha = 0.05, beta = 0; (b) alpha = 1, beta = 0; (c) alpha = 0.05, beta = 0.0005; and (d) alpha = 1, beta = 0.0005.
Each projection behaves a little differently. Notably, raising \(\alpha\) from 0 to 1 appears to keep the population from declining (b, d), and changing \(\beta\) from 0 to 0.0005 leads to a seemingly stabilized population size (d) rather than what appears to be extinction (c). The patterns are contingent on this particular MPM and to some extent on the start vector, meaning that other MPMs will yield different patterns with these \(\alpha\) and \(\beta\) values.
Let’s also try four parameterizations of the two-parameter Beverton-Holt function.
c2d_5 <- density_input(cypmean2, stage3 = c("P1", "P1"), stage2= c("SD", "rep"),
style = 2, time_delay = 1, alpha = 0.05, beta = 0, type = c(2, 2))
c2d_6 <- density_input(cypmean2, stage3 = c("P1", "P1"), stage2= c("SD", "rep"),
style = 2, time_delay = 1, alpha = 1, beta = 0, type = c(2, 2))
c2d_7 <- density_input(cypmean2, stage3 = c("P1", "P1"), stage2= c("SD", "rep"),
style = 2, time_delay = 1, alpha = 0.05, beta = 1, type = c(2, 2))
c2d_8 <- density_input(cypmean2, stage3 = c("P1", "P1"), stage2= c("SD", "rep"),
style = 2, time_delay = 1, alpha = 1, beta = 1, type = c(2, 2))
cypproj2_d5 <- projection3(cypmean2, times = 100, substoch = 2, density = c2d_5)
> Warning: This function takes annual matrices as input. This lefkoMat object
> appears to be a set of mean matrices, and may lack annual matrices. Will project
> only the mean.
cypproj2_d6 <- projection3(cypmean2, times = 100, substoch = 2, density = c2d_6)
> Warning: This function takes annual matrices as input. This lefkoMat object
> appears to be a set of mean matrices, and may lack annual matrices. Will project
> only the mean.
cypproj2_d7 <- projection3(cypmean2, times = 100, substoch = 2, density = c2d_7)
> Warning: This function takes annual matrices as input. This lefkoMat object
> appears to be a set of mean matrices, and may lack annual matrices. Will project
> only the mean.
cypproj2_d8 <- projection3(cypmean2, times = 100, substoch = 2, density = c2d_8)
> Warning: This function takes annual matrices as input. This lefkoMat object
> appears to be a set of mean matrices, and may lack annual matrices. Will project
> only the mean.Now let’s plot them (figure 10.13).
par(mfrow = c(2, 2))
plot(cypproj2_d5); title("a)", adj = 0)
plot(cypproj2_d6, ylab = ""); title("b)", adj = 0)
plot(cypproj2_d7); title("c)", adj = 0)
plot(cypproj2_d8, ylab = ""); title("d)", adj = 0)
Figure 10.13: Density dependent projections using the two-parameter Beverton-Holt function. (a) alpha = 0.05, beta = 0; (b) alpha = 1, beta = 0; (c) alpha = 0.05, beta = 1; and (d) alpha = 1, beta = 1.
The Beverton-Holt function is less likely to produce unusual shapes to population projections than the Ricker function. For example, it is generally not capable of producing chaos. This is both a strength and a weakness - its simplicity makes it attractive in many situations, but prevents it from matching reality in others. Users will need to explore the different density functions under different parameterizations to decide on an appropriate parameterization of density dependence for their research.
10.2.5 Ordered or cyclical density dependent simulations
We have already seen that we can set our projections to run through a specific sequence of years. This sequence can be combined with density dependence in any way that we wish. Here, we conduct a density dependent version of the previous projection in which we set the first ten years to repeatedly switch between two wet years (2007 and 2008), followed by ten years of switching between two drier years (2004 and 2005). To make things interesting in our 20 year projection, we utilize an Usher function with \(\alpha = 1\) and \(\beta = 0\), and our previous start frame with 1000 3rd year protocorms and 100 dormant adults.
year_order <- c(2007, 2008, 2007, 2008, 2007, 2008, 2007, 2008, 2007, 2008,
2004, 2005, 2004, 2005, 2004, 2005, 2004, 2005, 2004, 2005)
c2d_9 <- density_input(cypmean2, stage3 = c("P1", "P1"), stage2= c("SD", "rep"),
style = 3, time_delay = 1, alpha = 1, beta = 0, type = c(2, 2))
cypproj2_ord1d <- projection3(cypmatrix2r, year = year_order, times = 20,
start_frame = c2m_sv, density = c2d_9)
summary(cypproj2_ord1d)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 20 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 6 1
> 11 1
> 16 1
> 21 1
As before, if the sequence of years noted in the year option is shorter than the time steps designated in the times option, then lefko3 will cycle through the sequences of years until the projection is done. For example, if we leave the times option out from the above projection, then lefko3 will assume the default of 10,000 time steps, and cycle through our 20 defined years over and over, as below.
cypproj2_ord2d <- projection3(cypmatrix2r, year = year_order,
start_frame = c2m_sv, density = c2d_9)
summary(cypproj2_ord2d)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 10000 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 2501 0
> 5001 0
> 7501 0
> 10001 0Let’s see plots of these projections (figure 10.14).
par(mfrow=c(1,2))
plot(cypproj2_ord1d)
title("a)", adj = 0)
plot(cypproj2_ord2d, ylab = "")
title("b)", adj = 0)
Figure 10.14: Ordered,density dependent projections can be cyclical. (a) A 20 year ordered projection. (b) A 10,000 year projection cycling through 20 ordered years
It appears that this particular density dependence leads to extinction fairly early on.
Projecting a full MPM with multiple annual matrices defaults to cycling through the matrices in order, also when density dependence i set. Let’s try the first four Ricker projections again, but this time using the full ahistorical MPM rather than the arithmetic mean matrix to allow cycling of year sequences.
cypproj2r_d1 <- projection3(cypmatrix2r, times = 100, substoch = 2,
density = c2d_1)
cypproj2r_d2 <- projection3(cypmatrix2r, times = 100, substoch = 2,
density = c2d_2)
cypproj2r_d3 <- projection3(cypmatrix2r, times = 100, substoch = 2,
density = c2d_3)
cypproj2r_d4 <- projection3(cypmatrix2r, times = 100, substoch = 2,
density = c2d_4)
summary(cypproj2r_d1)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 100 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
summary(cypproj2r_d2)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 100 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
summary(cypproj2r_d3)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 100 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
summary(cypproj2r_d4)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 100 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1Let’s see what these projections look like (10.15).
par(mfrow=c(2, 2))
plot(cypproj2r_d1)
title("a)", adj = 0)
plot(cypproj2r_d2, ylab = "")
title("b)", adj = 0)
plot(cypproj2r_d3)
title("c)", adj = 0)
plot(cypproj2r_d4, ylab = "")
title("d)", adj = 0)
Figure 10.15: Cyclical density dependent projections using the Ricker function
We can see the obvious cycles all projections, but the population seems to be on the road to extinction in a) and c). Projection c) appears to show a growing population, while d) seems to be entering a stable oscillation.
10.2.6 Replicated simulations assuming temporal stochasticity and density dependence
Finally, we may wish to conduct projections assuming both temporal stochasticity and density dependence. Here, we will set up four projections, each assuming one of four density frames assuming different parameterizations of the Ricker function working on germination (we have defined these density frames previously).
set.seed(42)
cypproj2r_sd1 <- projection3(cypmatrix2r, stochastic = TRUE, times = 100,
substoch = 2, density = c2d_1)
cypproj2r_sd2 <- projection3(cypmatrix2r, stochastic = TRUE, times = 100,
substoch = 2, density = c2d_2)
cypproj2r_sd3 <- projection3(cypmatrix2r, stochastic = TRUE, times = 100,
substoch = 2, density = c2d_3)
cypproj2r_sd4 <- projection3(cypmatrix2r, stochastic = TRUE, times = 100,
substoch = 2, density = c2d_4)
summary(cypproj2r_sd1)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 100 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
summary(cypproj2r_sd2)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 100 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
summary(cypproj2r_sd3)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 100 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
summary(cypproj2r_sd4)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 100 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1Once again, let’s plot the projections for comparison.
par(mfrow=c(2, 2))
plot(cypproj2r_sd1)
title("a)", adj = 0)
plot(cypproj2r_sd2, ylab = "")
title("b)", adj = 0)
plot(cypproj2r_sd3)
title("c)", adj = 0)
plot(cypproj2r_sd4, ylab = "")
title("d)", adj = 0)
Figure 10.16: Stochastic density dependent projections assuming the two-parameter Ricker functon. (a) alpha = 0.05, beta = 0; (b) alpha = 1, beta = 0; (c) alpha = 0.05, beta = 0.0005; and (d) alpha = 1, beta = 0.0005.
We can clearly see the impact of temporal environmental stochasticity in each plot, and different density dependence relationships influences population trajectories strongly across the plots. Projections a) and c) appear to be moving to extinction, while projection b) is growing and projection d) appears to be fluctuating around a stable number. We can alter other settings, such as changing the weights of the matrices and using different density dependence relationships, as well as replicating our stochastic runs. Here, we re-run the first two projections with 100 replicates each and plot the results.
set.seed(42)
cypproj2r_sd1_100 <- projection3(cypmatrix2r, stochastic = TRUE, times = 100,
nreps = 100, substoch = 2, density = c2d_1)
cypproj2r_sd2_100 <- projection3(cypmatrix2r, stochastic = TRUE, times = 100,
nreps = 100, substoch = 2, density = c2d_2)
par(mfrow = c(1, 2))
plot(cypproj2r_sd1_100)
title("a)", adj = 0)
plot(cypproj2r_sd2_100, ylab = "")
title("b)", adj = 0)
Figure 10.17: Stochastic density dependent projections with replicates assuming the two-parameter Ricker function. (a) alpha = 0.05, beta = 0; and (b) alpha = 1, beta = 0.
The general trends are the same as before, but we see the degree of variability in these trends is large.
10.3 Projecting function-based MPMs
Function-based projections provide us with the flexibility to construct custom matrices at each time step given some set of conditions that we are interested in. These are quite useful because of the flexibility offered. The price paid is that it takes extra time to develop matrices at each time step in the projection, making standard MPM projection through function projection3() preferable except in cases where custom matrices are genuinely needed at each time step.
What sorts of situations might require custom matrices at each time step? One obvious situation is when matrices are to be built using set values of some individual or environmental covariate. For example, we may have developed a set of vital rate models that show a relationship with climatic variables that are likely to change going into the future. We might then wish to project the population forward under different climate change scenarios, with vectors of specific values of these climatic variables derived from climate prediction models (e.g. Shefferson et al. 2017). In this circumstance, using already made matrices in some sequence will not likely be of much help to us because matrices should change at each time step as a function of climate. Instead, we should create new matrices at each projected time, assuming the values of the climatic variables that are relevant.
Let’s explore this with the Cypripedium candidum dataset. We introduced the model and code for for our stageframe and supplements in Chapter 2, and showed how to build function-based MPMs with all of this material in Chapter 5. Let’s look again at the life history model that we will use (note that this is the same as figure 2.2).
Figure 10.18: Life history model of Cypripedium candidum for use in function-based MPMs
This is a fairly large life history model with two adult life stages for each number of sprouts - a reproductive stage that flowers and a non-reproductive stage that does not flower. Let’s load the raw dataset and set up the stageframe first.
data(cypdata)
stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "V1", "V2", "V3", "V4",
"V5", "V6", "V7", "V8", "V9", "V10", "V11", "V12", "V13", "V14", "V15", "V16",
"V17", "V18", "V19", "V20", "V21", "V22", "V23", "V24", "F1", "F2", "F3",
"F4", "F5", "F6", "F7", "F8", "F9", "F10", "F11", "F12", "F13", "F14", "F15",
"F16", "F17", "F18", "F19", "F20", "F21", "F22", "F23", "F24")
indataset <- c(0, 0, 0, 0, 0, rep(1, 49))
sizevector <- c(0, 0, 0, 0, 0, seq(from = 0, t = 24), seq(from = 1, to = 24))
repvector <- c(0, 0, 0, 0, 0, rep(0, 25), rep(1, 24))
obsvector <- c(0, 0, 0, 0, 0, 0, rep(1, 48))
matvector <- c(0, 0, 0, 0, 0, rep(1, 49))
immvector <- c(0, 1, 1, 1, 1, rep(0, 49))
propvector <- c(1, rep(0, 53))
comments <- c("Dormant seed", "Yr1 protocorm", "Yr2 protocorm", "Yr3 protocorm",
"Seedling", "Veg dorm", "Veg adult 1 stem", "Veg adult 2 stems",
"Veg adult 3 stems", "Veg adult 4 stems", "Veg adult 5 stems",
"Veg adult 6 stems", "Veg adult 7 stems", "Veg adult 8 stems",
"Veg adult 9 stems", "Veg adult 10 stems", "Veg adult 11 stems",
"Veg adult 12 stems", "Veg adult 13 stems", "Veg adult 14 stems",
"Veg adult 15 stems", "Veg adult 16 stems", "Veg adult 17 stems",
"Veg adult 18 stems", "Veg adult 19 stems", "Veg adult 20 stems",
"Veg adult 21 stems", "Veg adult 22 stems", "Veg adult 23 stems",
"Veg adult 24 stems", "Flo adult 1 stem", "Flo adult 2 stems",
"Flo adult 3 stems", "Flo adult 4 stems", "Flo adult 5 stems",
"Flo adult 6 stems", "Flo adult 7 stems", "Flo adult 8 stems",
"Flo adult 9 stems", "Flo adult 10 stems", "Flo adult 11 stems",
"Flo adult 12 stems", "Flo adult 13 stems", "Flo adult 14 stems",
"Flo adult 15 stems", "Flo adult 16 stems", "Flo adult 17 stems",
"Flo adult 18 stems", "Flo adult 19 stems", "Flo adult 20 stems",
"Flo adult 21 stems", "Flo adult 22 stems", "Flo adult 23 stems",
"Flo adult 24 stems")
cypframe_fb <- sf_create(sizes = sizevector, stagenames = stagevector,
repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
propstatus = propvector, immstatus = immvector, indataset = indataset,
comments = comments)Now let’s standardize the dataset into a historically formatted vertical (hfv) data frame. We will also need to introduce our individual covariate, total annual precipitation, into the data, as below, since the dataset does not currently include climatic data. These data were derived from the National Centers for Environmental Information and the U.S. National Oceanic and Atmospheric Administration. We will insert these climatic data as constants into the data frame holding the dataset.
cypdata_env <- cypdata
cypdata_env$prec.04 <- 92.2
cypdata_env$prec.05 <- 57.6
cypdata_env$prec.06 <- 96.0
cypdata_env$prec.07 <- 109.8
cypdata_env$prec.08 <- 111.9
cypdata_env$prec.09 <- 106.8
cypfb_env <- verticalize3(data = cypdata_env, noyears = 6, firstyear = 2004,
patchidcol = "patch", individcol = "plantid", blocksize = 4,
sizeacol = "Inf2.04", sizebcol = "Inf.04", sizeccol = "Veg.04",
repstracol = "Inf.04", repstrbcol = "Inf2.04", fecacol = "Pod.04",
indcovacol = c("prec.04", "prec.05", "prec.06", "prec.07", "prec.08",
"prec.09"),
stageassign = cypframe_fb, stagesize = "sizeadded", NAas0 = TRUE,
age_offset = 4)
summary(cypfb_env)
> rowid popid patchid individ year2
> Min. : 1.00 :320 A: 93 Length:320 Min. :2004
> 1st Qu.:21.00 B:154 Class :character 1st Qu.:2005
> Median :37.50 C: 73 Mode :character Median :2006
> Mean :38.45 Mean :2006
> 3rd Qu.:56.00 3rd Qu.:2007
> Max. :77.00 Max. :2008
> firstseen lastseen obsage obslifespan
> Min. :2004 Min. :2004 Min. :5.000 Min. :0.000
> 1st Qu.:2004 1st Qu.:2009 1st Qu.:6.000 1st Qu.:5.000
> Median :2004 Median :2009 Median :7.000 Median :5.000
> Mean :2004 Mean :2009 Mean :6.853 Mean :4.556
> 3rd Qu.:2004 3rd Qu.:2009 3rd Qu.:8.000 3rd Qu.:5.000
> Max. :2008 Max. :2009 Max. :9.000 Max. :5.000
> sizea1 sizeb1 sizec1 size1added
> Min. :0.000000 Min. : 0.0000 Min. : 0.0 Min. : 0.000
> 1st Qu.:0.000000 1st Qu.: 0.0000 1st Qu.: 0.0 1st Qu.: 0.000
> Median :0.000000 Median : 0.0000 Median : 1.0 Median : 2.000
> Mean :0.009375 Mean : 0.7469 Mean : 1.9 Mean : 2.656
> 3rd Qu.:0.000000 3rd Qu.: 1.0000 3rd Qu.: 3.0 3rd Qu.: 4.000
> Max. :1.000000 Max. :18.0000 Max. :13.0 Max. :21.000
> repstra1 repstrb1 feca1 indcova1
> Min. : 0.0000 Min. :0.000000 Min. :0.0000 Min. : 0.00
> 1st Qu.: 0.0000 1st Qu.:0.000000 1st Qu.:0.0000 1st Qu.: 57.60
> Median : 0.0000 Median :0.000000 Median :0.0000 Median : 92.20
> Mean : 0.7469 Mean :0.009375 Mean :0.2656 Mean : 70.64
> 3rd Qu.: 1.0000 3rd Qu.:0.000000 3rd Qu.:0.0000 3rd Qu.: 96.00
> Max. :18.0000 Max. :1.000000 Max. :7.0000 Max. :109.80
> juvgiven1 obsstatus1 repstatus1 fecstatus1
> Min. :0 Min. :0.0000 Min. :0.0000 Min. :0.0000
> 1st Qu.:0 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000
> Median :0 Median :1.0000 Median :0.0000 Median :0.0000
> Mean :0 Mean :0.7469 Mean :0.2875 Mean :0.1344
> 3rd Qu.:0 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:0.0000
> Max. :0 Max. :1.0000 Max. :1.0000 Max. :1.0000
> matstatus1 alive1 stage1 stage1index
> Min. :0.0000 Min. :0.0000 Length:320 Min. : 0.00
> 1st Qu.:1.0000 1st Qu.:1.0000 Class :character 1st Qu.: 6.00
> Median :1.0000 Median :1.0000 Mode :character Median : 8.00
> Mean :0.7688 Mean :0.7688 Mean :14.17
> 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:31.00
> Max. :1.0000 Max. :1.0000 Max. :51.00
> sizea2 sizeb2 sizec2 size2added
> Min. :0.000000 Min. : 0.0000 Min. : 0.000 Min. : 0.000
> 1st Qu.:0.000000 1st Qu.: 0.0000 1st Qu.: 1.000 1st Qu.: 1.000
> Median :0.000000 Median : 0.0000 Median : 2.000 Median : 2.000
> Mean :0.009375 Mean : 0.8969 Mean : 2.416 Mean : 3.322
> 3rd Qu.:0.000000 3rd Qu.: 1.0000 3rd Qu.: 3.000 3rd Qu.: 4.000
> Max. :1.000000 Max. :18.0000 Max. :13.000 Max. :24.000
> repstra2 repstrb2 feca2 indcova2
> Min. : 0.0000 Min. :0.000000 Min. :0.0000 Min. : 57.60
> 1st Qu.: 0.0000 1st Qu.:0.000000 1st Qu.:0.0000 1st Qu.: 92.20
> Median : 0.0000 Median :0.000000 Median :0.0000 Median : 96.00
> Mean : 0.8969 Mean :0.009375 Mean :0.2906 Mean : 92.77
> 3rd Qu.: 1.0000 3rd Qu.:0.000000 3rd Qu.:0.0000 3rd Qu.:109.80
> Max. :18.0000 Max. :1.000000 Max. :7.0000 Max. :111.90
> juvgiven2 obsstatus2 repstatus2 fecstatus2 matstatus2
> Min. :0 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :1
> 1st Qu.:0 1st Qu.:1.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:1
> Median :0 Median :1.0000 Median :0.0000 Median :0.0000 Median :1
> Mean :0 Mean :0.9531 Mean :0.3688 Mean :0.1562 Mean :1
> 3rd Qu.:0 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:1
> Max. :0 Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1
> alive2 stage2 stage2index sizea3
> Min. :1 Length:320 Min. : 6.00 Min. :0.000000
> 1st Qu.:1 Class :character 1st Qu.: 7.00 1st Qu.:0.000000
> Median :1 Mode :character Median :10.00 Median :0.000000
> Mean :1 Mean :18.17 Mean :0.009375
> 3rd Qu.:1 3rd Qu.:32.00 3rd Qu.:0.000000
> Max. :1 Max. :54.00 Max. :1.000000
> sizeb3 sizec3 size3added repstra3
> Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. : 0.000
> 1st Qu.: 0.000 1st Qu.: 1.000 1st Qu.: 1.000 1st Qu.: 0.000
> Median : 0.000 Median : 1.000 Median : 2.000 Median : 0.000
> Mean : 1.069 Mean : 2.209 Mean : 3.288 Mean : 1.069
> 3rd Qu.: 1.000 3rd Qu.: 3.000 3rd Qu.: 4.000 3rd Qu.: 1.000
> Max. :18.000 Max. :13.000 Max. :24.000 Max. :18.000
> repstrb3 feca3 indcova3 juvgiven3 obsstatus3
> Min. :0.000000 Min. :0.0000 Min. : 57.6 Min. :0 Min. :0.0
> 1st Qu.:0.000000 1st Qu.:0.0000 1st Qu.: 96.0 1st Qu.:0 1st Qu.:1.0
> Median :0.000000 Median :0.0000 Median :106.8 Median :0 Median :1.0
> Mean :0.009375 Mean :0.4562 Mean : 96.1 Mean :0 Mean :0.9
> 3rd Qu.:0.000000 3rd Qu.:0.0000 3rd Qu.:109.8 3rd Qu.:0 3rd Qu.:1.0
> Max. :1.000000 Max. :8.0000 Max. :111.9 Max. :0 Max. :1.0
> repstatus3 fecstatus3 matstatus3 alive3
> Min. :0.0 Min. :0.0000 Min. :1 Min. :0.0000
> 1st Qu.:0.0 1st Qu.:0.0000 1st Qu.:1 1st Qu.:1.0000
> Median :0.0 Median :0.0000 Median :1 Median :1.0000
> Mean :0.4 Mean :0.2219 Mean :1 Mean :0.9469
> 3rd Qu.:1.0 3rd Qu.:0.0000 3rd Qu.:1 3rd Qu.:1.0000
> Max. :1.0 Max. :1.0000 Max. :1 Max. :1.0000
> stage3 stage3index
> Length:320 Min. : 0.00
> Class :character 1st Qu.: 7.00
> Mode :character Median :10.00
> Mean :18.57
> 3rd Qu.:33.00
> Max. :54.00
Now we will load the ahistorical supplement table that we will use for the ahistorical projections. We will not work with the historical version in this chapter, just due to the increased amount of processing time required for the historical modelsearch() run.
seeds_per_fruit <- 5000
sl_mult <- 0.7
cypsupp2_fb <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "SL", "D",
"V1", "V2", "V3", "SD", "P1"),
stage2 = c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "SL", "SL", "rep",
"rep"),
eststage3 = c(NA, NA, NA, NA, NA, NA, "D", "V1", "V2", "V3", NA, NA),
eststage2 = c(NA, NA, NA, NA, NA, NA, "D", "D", "D", "D", NA, NA),
givenrate = c(0.08, 0.1, 0.1, 0.1, 0.05, 0.05, NA, NA, NA, NA, NA, NA),
multiplier = c(NA, NA, NA, NA, NA, NA, sl_mult, sl_mult, sl_mult, sl_mult,
0.5 * seeds_per_fruit, 0.5 * seeds_per_fruit),
type = c("S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "R", "R"),
stageframe = cypframe_fb, historical = FALSE)
cypsupp2_fb
> stage3 stage2 stage1 eststage3 eststage2 eststage1 givenrate multiplier
> 1 SD SD <NA> <NA> <NA> <NA> 0.08 1.0
> 2 P1 SD <NA> <NA> <NA> <NA> 0.10 1.0
> 3 P2 P1 <NA> <NA> <NA> <NA> 0.10 1.0
> 4 P3 P2 <NA> <NA> <NA> <NA> 0.10 1.0
> 5 SL P3 <NA> <NA> <NA> <NA> 0.05 1.0
> 6 SL SL <NA> <NA> <NA> <NA> 0.05 1.0
> 7 D SL <NA> D D <NA> NA 0.7
> 8 V1 SL <NA> V1 D <NA> NA 0.7
> 9 V2 SL <NA> V2 D <NA> NA 0.7
> 10 V3 SL <NA> V3 D <NA> NA 0.7
> 11 SD rep <NA> <NA> <NA> <NA> NA 2500.0
> 12 P1 rep <NA> <NA> <NA> <NA> NA 2500.0
> convtype convtype_t12
> 1 1 1
> 2 1 1
> 3 1 1
> 4 1 1
> 5 1 1
> 6 1 1
> 7 1 1
> 8 1 1
> 9 1 1
> 10 1 1
> 11 3 1
> 12 3 1
Let’s now determine the proper distributions for size and fecundity. As a reminder, size in the Cypripedium dataset is given as the numbers of sprouts, and our proxy for fecundity can be either the number of flowers or the number of fruits produced by an individual. So, size and fecundity in this dataset are all non-negative integers. Count variables generally break the assumptions of a Gaussian distribution because the mean and variance are strongly related, possible values are discrete, and distributions are often dramatically skewed. Additionally, because we have absorbed individuals with a size of zero into the dormant stage, there are no zeros left for the size model and so we will need a zero-truncated distribution. Fecundity, in contrast, is a count variable that includes zeros and so we should test whether we need a zero-inflated model. As for which specific distributions to use, we will use function sf_distrib() to make this determination.
sf_distrib(cypfb_env, sizea = c("size3added", "size2added"), obs3= "obsstatus3",
fec = c("feca3", "feca2"), repst = c("repstatus3", "repstatus2"),
zisizea = FALSE)
> Mean sizea is 3.653
>
> The variance in sizea is 13.41
>
> The probability of this dispersion level by chance assuming that
> the true mean sizea = variance in sizea,
> and an alternative hypothesis of overdispersion, is 3.721e-138
>
> Primary size is significantly overdispersed.
>
> Mean fec is 0.8036
>
> The variance in fec is 1.601
>
> The probability of this dispersion level by chance assuming that
> the true mean fec = variance in fec,
> and an alternative hypothesis of overdispersion, is 0.07088
>
> Dispersion level in fecundity matches expectation.
>
>
> Mean lambda in fec is 0.4477
> The actual number of 0s in fec is 65
> The expected number of 0s in fec under the null hypothesis is 50.15
> The probability of this deviation in 0s from expectation by chance is 1.742e-06
>
> Fecundity is significantly zero-inflated.The results show that size is significantly overdispersed. This means that we cannot use the Poisson distribution, and should instead use the negative binomial distribution. Fecundity is not significantly overdispersed, but it is significantly zero-inflated, so we will require a zero-inflated Poisson distribution for that vital rate.
Now we will run function modelsearch() and find the right parameterizations for our vital rates. Note that we need to include the option indcova = c("indcova3", "indcova2", "indcova1") in order to make sure that modelsearch() includes our individual covariate, which reflects total annual precipitation. If these covariates were categorical and random, then we would also need to stipulate random.indcova = TRUE (our individual covariate is actually quantitative and numerical, so we will use the default setting of FALSE, instead). To prevent this function from taking too long to run, we will drop the patch term, set suite = "main" to prevent the inclusion of interaction terms, and run only an ahistorical set.
cypmodels2_env1 <- modelsearch(cypfb_env, historical = FALSE, approach= "mixed",
vitalrates = c("surv", "obs", "size", "repst", "fec"), sizedist = "negbin",
size.trunc = TRUE, fecdist = "poisson", fec.zero = TRUE, suite = "main",
size = c("size3added", "size2added", "size1added"),
indcova = c("indcova3", "indcova2", "indcova1"), quiet = TRUE)Let’s see the model summaries.
summary(cypmodels2_env1)
> This LefkoMod object includes 5 linear models.
> Best-fit model criterion used: aicc&k
>
>
>
> Survival model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: alive3 ~ indcova2 + (1 | year2) + (1 | individ)
> Data: subdata
> AIC BIC logLik deviance df.resid
> 87.6075 102.6808 -39.8038 79.6075 316
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 552.9
> year2 (Intercept) 394.6
> Number of obs: 320, groups: individ, 74; year2, 5
> Fixed Effects:
> (Intercept) indcova2
> 476.174 -3.278
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
>
>
>
> Observation model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: obsstatus3 ~ size2added + (1 | year2) + (1 | individ)
> Data: subdata
> AIC BIC logLik deviance df.resid
> 118.2567 133.1117 -55.1284 110.2567 299
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 1.078e-05
> year2 (Intercept) 8.776e-01
> Number of obs: 303, groups: individ, 70; year2, 5
> Fixed Effects:
> (Intercept) size2added
> 2.4904 0.3134
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
>
>
>
> Size model:
> Formula:
> size3added ~ indcova2 + size2added + (1 | year2) + (1 | individ)
> Data: subdata
> AIC BIC logLik df.resid
> 1001.3778 1023.3555 -494.6889 282
> Random-effects (co)variances:
>
> Conditional model:
> Groups Name Std.Dev.
> year2 (Intercept) 0.06554
> individ (Intercept) 0.96060
>
> Number of obs: 288 / Conditional model: year2, 5; individ, 70
>
> Dispersion parameter for truncated_nbinom2 family (): 1.87e+12
>
> Fixed Effects:
>
> Conditional model:
> (Intercept) indcova2 size2added
> 0.070927 0.004918 0.023979
>
>
>
> Secondary size model:
> [1] 1
>
>
>
> Tertiary size model:
> [1] 1
>
>
>
> Reproductive status model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: repstatus3 ~ indcova2 + repstatus2 + size2added + (1 | year2) +
> (1 | individ)
> Data: subdata
> AIC BIC logLik deviance df.resid
> 330.0418 352.0196 -159.0209 318.0418 282
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 0.0002695
> year2 (Intercept) 0.2479120
> Number of obs: 288, groups: individ, 70; year2, 5
> Fixed Effects:
> (Intercept) indcova2 repstatus2 size2added
> -4.23766 0.02954 1.71091 0.17187
>
>
>
> Fecundity model:
> Formula: feca2 ~ size2added + (1 | year2) + (1 | individ)
> Zero inflation: ~size2added + (1 | year2) + (1 | individ)
> Data: subdata
> AIC BIC logLik df.resid
> 248.8609 271.0264 -116.4305 110
> Random-effects (co)variances:
>
> Conditional model:
> Groups Name Std.Dev.
> year2 (Intercept) 0.5760
> individ (Intercept) 0.1639
>
> Zero-inflation model:
> Groups Name Std.Dev.
> year2 (Intercept) 1.642e-06
> individ (Intercept) 3.089e-04
>
> Number of obs: 118 / Conditional model: year2, 5; individ, 51 / Zero-inflation model: year2, 5; individ, 51
>
> Fixed Effects:
>
> Conditional model:
> (Intercept) size2added
> -0.54014 0.06174
>
> Zero-inflation model:
> (Intercept) size2added
> 3.865 -1.574
>
>
> Juvenile survival model:
> [1] 1
>
>
>
> Juvenile observation model:
> [1] 1
>
>
>
> Juvenile size model:
> [1] 1
>
>
>
> Juvenile secondary size model:
> [1] 1
>
>
>
> Juvenile tertiary size model:
> [1] 1
>
>
>
> Juvenile reproduction model:
> [1] 1
>
>
>
> Juvenile maturity model:
> [1] 1
>
>
>
>
>
> Number of models in survival table: 8
>
> Number of models in observation table: 8
>
> Number of models in size table: 8
>
> Number of models in secondary size table: 1
>
> Number of models in tertiary size table: 1
>
> Number of models in reproduction status table: 8
>
> Number of models in fecundity table: 61
>
> Number of models in juvenile survival table: 1
>
> Number of models in juvenile observation table: 1
>
> Number of models in juvenile size table: 1
>
> Number of models in juvenile secondary size table: 1
>
> Number of models in juvenile tertiary size table: 1
>
> Number of models in juvenile reproduction table: 1
>
> Number of models in juvenile maturity table: 1
>
>
>
>
>
> General model parameter names (column 1), and
> specific names used in these models (column 2):
> parameter_names mainparams
> 1 time t year2
> 2 individual individ
> 3 patch patch
> 4 alive in time t+1 surv3
> 5 observed in time t+1 obs3
> 6 sizea in time t+1 size3
> 7 sizeb in time t+1 sizeb3
> 8 sizec in time t+1 sizec3
> 9 reproductive status in time t+1 repst3
> 10 fecundity in time t+1 fec3
> 11 fecundity in time t fec2
> 12 sizea in time t size2
> 13 sizea in time t-1 size1
> 14 sizeb in time t sizeb2
> 15 sizeb in time t-1 sizeb1
> 16 sizec in time t sizec2
> 17 sizec in time t-1 sizec1
> 18 reproductive status in time t repst2
> 19 reproductive status in time t-1 repst1
> 20 maturity status in time t+1 matst3
> 21 maturity status in time t matst2
> 22 age in time t age
> 23 density in time t density
> 24 individual covariate a in time t indcova2
> 25 individual covariate a in time t-1 indcova1
> 26 individual covariate b in time t indcovb2
> 27 individual covariate b in time t-1 indcovb1
> 28 individual covariate c in time t indcovc2
> 29 individual covariate c in time t-1 indcovc1
> 30 stage group in time t group2
> 31 stage group in time t-1 group1
>
>
>
>
>
> Quality control:
>
> Survival estimated with 74 individuals and 320 individual transitions.
> Survival accuracy is 1.
> Observation estimated with 70 individuals and 303 individual transitions.
> Observation accuracy is 0.95.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Primary size pseudo R-squared is 0.011.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Reproductive status accuracy is 0.719.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Fecundity pseudo R-squared is 0.32.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity probability not estimated.Total annual precipitation is a factor in the best-fit models of survival, primary size, and reproductive status, suggesting that it is quite important to the demography of the population. The accuracy of survival and observation status is really excellent, while the other models have poor to fair accuracy or pseudo-R2.
Now that we have developed our vital rate models, let’s also build vectors showing our expectation for shifts in the precipitation over the next 20 years from the year that data collection stopped. This 20 year period will be our main time window of interest for projection. Our climate values will be real predictions from Shefferson et al. (2017), which were extracted from a real climatological forecast model maintained by Japan’s National Meteorological Research Institute. We will use the values predicted for 2010 to 2030. Just out of interest, we will also create a vector of predicted values for 2070 to 2090.
pred_vals_2010to30 <- c(110.8152136, 98.87159542, 127.0623005, 110.0501797,
82.03005742, 103.6893125, 104.2549713, 97.86989854, 111.25054, 116.4481035,
125.9652504, 98.89823738, 105.7075192, 135.7032496, 149.6705465, 149.0110188,
115.3965537, 119.7017105, 99.5415225, 111.7898123, 90.51828146)
pred_vals_2070to90 <- c(101.7828506, 80.70954521, 106.760383, 142.3355673,
113.6519693, 122.3257958, 100.8745268, 138.814026, 111.1349475, 92.00440938,
120.2391638, 115.6253448, 102.6373673, 125.3819469, 98.25464504, 102.4367633,
118.6401779, 125.6021119, 101.9117747, 127.2149739, 109.244138)Let’s see what these values look like graphically (figure 10.19).
plot(x = c(1:21), y = pred_vals_2010to30, type = "l", lty = 1, lwd = 3,
ylab = "Annual precip (cm)", xlab = "Projected year", ylim = c(80, 170),
bty = "n")
lines(x = c(1:21), y = pred_vals_2070to90, lty = 2, lwd = 3)
legend("topright", c("2010-2030", "2070-2090"), lty = c(1, 2), lwd = 3,
bty = "n")
Figure 10.19: Predicted annual precipitation for Cypripedium candidum
Now let’s run our projections. Note that many of the options in f_projection3() work similarly to those of projection3(). For example, projection can be limited to whole individuals if integeronly = TRUE is set (we will not do so here in order to see the overall population patterns over time, and limit the influence of population size on predicted extinction risk). Users should experiment with the various options to see the impacts.
10.3.1 Deterministic projections of a specific time length
The first style of projection that we might try is a standard deterministic projection using our individual covariate values. We have two choices - we may either choose a single year to run (our vital rates all include random year terms) with our projected precipitation values, or we may try a specific sequence of years. There are different merits to both approaches. For instructional purposes, let’s choose a single year in the middle of the monitoring - 2007. We also need to format the individual covariate into a three-column data frame with as many rows as we wish to project, and set zero as the default value for unused covariates. Here, we develop our first projection using the 2010 to 2030 predictions. Note that f_projection3() does not calculate mean matrices, unlike function projection3(), so we cannot use that approach.
ind_frame_1 <- cbind.data.frame(inda = pred_vals_2010to30, indb = 0, indc = 0)
trial2f_1 <- f_projection3(cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, times = 21,
year = 2007, ind_terms = ind_frame_1)
> Warning: Option patch not set, so will set to first patch/population.
The resulting object is of class lefkoProj, just as in the case of function projection3() above. The object includes the following elements:
- projection - A list with a single element representing the one population or patch projected. This element is a list with number of elements equal to the number of replicates (option
nreps, which defaults to 1). Each element in this lower-level list is a matrix with number of rows equal to the number of stages (if ahistorical; paired stages if historical; ages if Leslie; and age-stages if age-by-stage), and number of columns equal to the projected time steps plus 1 (the start time is time 0). The elements of the matrix show the projected number of individuals in that stage in that time. - stage_dist - A list showing the stage distribution projected per time, with the same order and number of elements as projection. Only estimated if
growthonly = FALSE. - rep_value - A list showing the reproductive value of each stage projected per time, with the same order and number of elements as projection. Only estimated if
growthonly = FALSEandrepvalue = TRUE. - pop_size - A list with number of elements equal to the number of projected populations or patches (in
f_projection3(), this can only be one population or patch, so the length of the list will always be one for this function). This one element is a matrix showing the total population size at each time (column) per replicate (row). - ahstages - The stageframe used in analysis, though modified and edited per rules set for function-based matrix estimators.
- hstages - A data frame matrix showing the pairing of ahistorical stages used to create historical stage pairs (only if historical).
- agestages - A data frame showing the order of age-stage pairs (only if age-by-stage).
- labels - Currently, this shows only the population and patch identities used.
- control - An integer vector indicating the number of replicates and the number of time steps.
- density - The data frame input under the density option. Only provided if input by the user.
We can view a summary of this object, as follows.
summary(trial2f_1)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 21 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 6 1
> 11 1
> 16 1
> 22 1Let’s now also build a projection covering climatic values projected for 2070 to 2090.
ind_frame_2 <- cbind.data.frame(inda = pred_vals_2070to90, indb = 0, indc = 0)
trial2f_2 <- f_projection3(cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, times = 21,
year = 2007, ind_terms = ind_frame_2)
> Warning: Option patch not set, so will set to first patch/population.
summary(trial2f_2)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 21 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 6 1
> 11 0
> 16 0
> 22 0Now let’s plot our projections side-by-side (figure 10.20).
par(mfrow=c(1,2))
plot(trial2f_1, xlab = "Time from 2010")
title("a)", adj = 0)
plot(trial2f_2, xlab = "Time from 2070", ylab = "")
title("b)", adj = 0)
Figure 10.20: Deterministic function-based projections for (a) 2010-2030 and (b) 2070-2090
Clearly our different climate projections have dramatically different impacts on population dynamics. The population seems to grow in response to climate predicted for 2010 to 2030, but quickly plunges between 2070 and 2090.
As before, we can reset the starting numbers of individuals in different stages. Function f_projection3() assumes one individual per stage at the start by default, similarly to projection3(). As in projection3(), we can set the starting vector with the start_vec and start_frame options. Here we show an example of the latter, which is more powerful than the former. Here, we will start the population with 1000 dormant seeds and 100 1st year protocorms (figure 10.21).
c2m_sv <- start_input(cypmean2, stage2 = c("SD", "P1"), value = c(1000, 100))
trial2f_3 <- f_projection3(cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, times = 21,
year = 2007, ind_terms = ind_frame_1, start_frame = c2m_sv)
> Warning: Option patch not set, so will set to first patch/population.
par(mfrow = c(1, 2))
plot(trial2f_1)
title("a)", adj = 0)
plot(trial2f_3, ylab = "")
title("b)", adj = 0)
Figure 10.21: Deterministic function-based projections assuming different start frames. (a) Start with one individual of each stage. (b) Start with 1000 dormant seeds and 100 1st year protocorms
In the new projection, we start the population with dormant seeds and 1st year protocorms only, which means that there is no reproduction in the population for the first few years. This causes an immediate drop and plateau until some seeds and protocorms start reproducing. Once reproductive individuals enter the population, the slope of population size begins increasing.
10.3.2 Ordered / cyclical projections of a specific time length
The vital rate models developed with modelsearch() generally include year2 terms denoting time at time t. This is a categorical variable, and so each year has its own value. If we do not specify a specific year to project, then f_projection3() will cycle through the years in the dataset in order for the number of time steps set, using the appropriate year2 term in each vital rate model used to estimate each matrix at each time step. Alternatively, if the year option is set to a vector of years, then these particular years will cycle in order for the number of time steps set. Note, however, that this function only sets the values for the categorical year2 term in the vital rate models - any related terms, such as individual covariates, are unaffected by this.
Let’s take an example using three years in particular: 2006, 2007, and 2008. By choosing these years and loading this choice as a vector in the year option, we will have lefko3 cycle through them in order. However, because our vital rate models also include precipitation, we need to incorporate values for precipitation, and our choices will make a big difference. As in the year case, f_projection3() will cycle though whatever values we provide until the projection is finished. These cycles are independent of the year option cycles, and so care needs to be taken to load only values genuinely needed. Here, we will run three projections using our year vector of three years. Each one will assume a different precipitation input.
ind_real <- data.frame(inda = c(92.2, 57.6, 96.0, 109.8, 111.9, 106.8),
indb = 0, indc = 0)
short_real <- as.data.frame(ind_real[1,])
mean_real <- data.frame(inda = mean(ind_real$inda), indb = 0, indc = 0)
trial2f_4a <- f_projection3(cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, times = 101,
year = c(2006, 2007, 2008), ind_terms = ind_real)
> Warning: Option patch not set, so will set to first patch/population.
> Warning: Fewer values of individual covariates have been supplied than required,
> so input values will be cycled.
trial2f_4b <- f_projection3(cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, times = 101,
year = c(2006, 2007, 2008), ind_terms = short_real)
> Warning: Option patch not set, so will set to first patch/population.
> Warning: Fewer values of individual covariates have been supplied than required,
> so input values will be cycled.
trial2f_4c <- f_projection3(cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, times = 101,
year = c(2006, 2007, 2008), ind_terms = mean_real)
> Warning: Option patch not set, so will set to first patch/population.
> Warning: Fewer values of individual covariates have been supplied than required,
> so input values will be cycled.
summary(trial2f_4a)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 102 1
summary(trial2f_4b)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 102 1
summary(trial2f_4c)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 102 1Here is another example using the full set of years.
trial2f_5a <- f_projection3(cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, times = 101,
ind_terms = ind_real)
> Warning: Option year not set, so will cycle through existing years.
> Warning: Option patch not set, so will set to first patch/population.
> Warning: Fewer values of individual covariates have been supplied than required,
> so input values will be cycled.
trial2f_5b <- f_projection3(cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, times = 101,
ind_terms = short_real)
> Warning: Option year not set, so will cycle through existing years.
> Warning: Option patch not set, so will set to first patch/population.
> Warning: Fewer values of individual covariates have been supplied than required,
> so input values will be cycled.
trial2f_5c <- f_projection3(cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, times = 101,
ind_terms = mean_real)
> Warning: Option year not set, so will cycle through existing years.
> Warning: Option patch not set, so will set to first patch/population.
> Warning: Fewer values of individual covariates have been supplied than required,
> so input values will be cycled.
summary(trial2f_5a)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 102 1
summary(trial2f_5b)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 0
> 51 0
> 76 0
> 102 0
summary(trial2f_5c)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 0
> 51 0
> 76 0
> 102 0Let’s plot the results together (figure 10.22).
par(mfrow=c(3,2))
plot(trial2f_4a)
title("a)", adj = 0)
plot(trial2f_4b, ylab = "")
title("b)", adj = 0)
plot(trial2f_4c)
title("c)", adj = 0)
plot(trial2f_5a, ylab = "")
title("d)", adj = 0)
plot(trial2f_5b)
title("e)", adj = 0)
plot(trial2f_5c, ylab = "")
title("f)", adj = 0)
Figure 10.22: Cyclical function-based projections using 2007-2009 only (a, b, c), or all years (d, e, f,), assuming all 6 real precip values (a, d), precip for 2004 only (b, e), and mean precip over 2004-2009 (d, f).
The choice of which years and individual covariate values to use has profound implications.
The approach above can be used to program any particular order of years. The year option should simply be paired with a vector of whatever length showing the order of year terms desired. Here, we use a sequence of years initially switching between years 2007 and 2008, then settling into 2008 for nine years, and ending on year 2005. Let’s assume constant mean precipitation.
year_order <- c(2007, 2008, 2007, 2008, 2007, 2008, 2007, 2008, 2007, 2008,
2008, 2008, 2008, 2008, 2008, 2008, 2008, 2008, 2008, 2005)
trial2f_6 <- f_projection3(cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, year = year_order,
times = 20, ind_terms = mean_real)
> Warning: Option patch not set, so will set to first patch/population.
> Warning: Fewer values of individual covariates have been supplied than required,
> so input values will be cycled.
As before, if we set a larger number for times than we will cycle through the sequence.
trial2f_7 <- f_projection3(cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, year = year_order,
times = 100, ind_terms = mean_real)
> Warning: Option patch not set, so will set to first patch/population.
> Warning: Fewer values of individual covariates have been supplied than required,
> so input values will be cycled.Let’s take a look at plots of these projections (figure 10.23).
par(mfrow=c(1, 2))
plot(trial2f_6)
title("a)", adj= 0)
plot(trial2f_7, ylab = "")
title("b)", adj= 0)
Figure 10.23: Ordered progression of year terms used in projection. (a) Specific sequence of 20 years projected for 20 years. (b) Specific sequence of 20 years projected cyclically for 100 years.
The population appears to grow. When cycled, the year values cause dips at roughly 20 year intervals, with larger upswings as time goes on.
10.3.3 Replicated simulations assuming temporal stochasticity
Function f_projection3() can also handle temporally stochastic projections and replicated simulations. In this case, randomness is included through random draws of year2 terms in the vital rate models, meaning that all other terms are unaffected and can be run in whatever order the user deems fit. Here we conduct a simple stochastic projection for 101 time steps in a single replicate, using the mean precipitation over the study period.
set.seed(42)
trial2f_8 <- f_projection3(cypfb_env, format = 3, stochastic = TRUE,
stageframe = cypframe_fb, supplement = cypsupp2_fb,
modelsuite = cypmodels2_env1, times = 101, ind_terms = mean_real)
> Warning: Option patch not set, so will set to first patch/population.
> Warning: Fewer values of individual covariates have been supplied than required,
> so input values will be cycled.
summary(trial2f_8)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 0
> 102 0Let’s say that we wished to compare this particular scenario with one in which we double fecundity through artificial pollination. For this purpose, we will define a new supplement table in which we double our two routes of fecundity. Then we will run a projection with this new supplement table.
cypsupp2_fb_alt <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "SL",
"D", "V1", "V2", "V3", "SD", "P1"),
stage2 = c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "SL", "SL", "rep",
"rep"),
eststage3 = c(NA, NA, NA, NA, NA, NA, "D", "V1", "V2", "V3", NA, NA),
eststage2 = c(NA, NA, NA, NA, NA, NA, "D", "D", "D", "D", NA, NA),
givenrate = c(0.08, 0.1, 0.1, 0.1, 0.05, 0.05, NA, NA, NA, NA, NA, NA),
multiplier = c(NA, NA, NA, NA, NA, NA, sl_mult, sl_mult, sl_mult, sl_mult,
0.5 * seeds_per_fruit * 2, 0.5 * seeds_per_fruit * 2),
type = c("S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "R", "R"),
stageframe = cypframe_fb, historical = FALSE)
set.seed(42)
trial2f_9 <- f_projection3(cypfb_env, format = 3, stochastic = TRUE,
stageframe = cypframe_fb, supplement = cypsupp2_fb_alt,
modelsuite = cypmodels2_env1, times = 101, ind_terms = mean_real)
> Warning: Option patch not set, so will set to first patch/population.
> Warning: Fewer values of individual covariates have been supplied than required,
> so input values will be cycled.
summary(trial2f_9)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 102 0Let’s look at plots to compare the results (figure 10.24).
par(mfrow=c(1,2))
plot(trial2f_8)
title("a)", adj = 0)
plot(trial2f_9, ylab = "")
title("b)", adj = 0)
Figure 10.24: Stochastic function-based projections testing the impacts of natural fecundity values (a) vs. doubled fecundity (b)
It certainly appears that doubling fecundity has a positive impact, with higher spikes in numbers and more visible cycles. However, these are single stochastic replicates, and we have no idea what the level of uncertainty is. So, let’s repeat the projections, but using 100 replicates each.
set.seed(42)
trial2f_10 <- f_projection3(cypfb_env, format = 3, stochastic = TRUE,
stageframe = cypframe_fb, supplement = cypsupp2_fb,
modelsuite = cypmodels2_env1, nreps = 100, times = 101,
ind_terms = mean_real)
set.seed(42)
trial2f_11 <- f_projection3(cypfb_env, format = 3, stochastic = TRUE,
stageframe = cypframe_fb, supplement = cypsupp2_fb_alt,
modelsuite = cypmodels2_env1, nreps = 100, times = 101,
ind_terms = mean_real)Let’s see some summaries.
summary(trial2f_10)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 100 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 100
> 26 80
> 51 41
> 76 16
> 102 3
summary(trial2f_11)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 100 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 100
> 26 87
> 51 68
> 76 47
> 102 29Let’s compare the results (figure 10.25).
par(mfrow=c(1,2))
plot(trial2f_10)
plot(trial2f_11, ylab = "")
Figure 10.25: Stochastic function-based projections testing doubled fecundity with replicates
We see here that actually, the original simulation with natural fecundity leads to extinction in seemingly all cases, while some simulations still hold extant populations at 100 time steps when fecundity is doubled.
Let’s try one more stochastic projection. This time, we will weight the matrices differently, with the 2007 matrix being chosen 88% of the time while the other matrices are chosen 3% of the time each.
set.seed(42)
trial2f_12 <- f_projection3(cypfb_env, format = 3, stochastic = TRUE,
stageframe = cypframe_fb, supplement = cypsupp2_fb,
modelsuite = cypmodels2_env1, nreps = 100, times = 101,
ind_terms = mean_real, tweights = c(0.03, 0.03, 0.03, 0.88, 0.03))And the summary….
summary(trial2f_12)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 100 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 100
> 26 100
> 51 100
> 76 100
> 102 100Let’s compare the new projection with the equivalent projection with equal year weights (figure 10.26).
par(mfrow = c(1, 2))
plot(trial2f_10); title("a)", adj = 0)
plot(trial2f_12, ylab = ""); title("b)", adj = 0)
Figure 10.26: Stochastic function-based projections assuming (a) equal weights for all years, and (b) user-defined weights for different years
Clearly, the overweighting on year 2007 had a positive effect on the population, with seemingly many more replicates maintaining themselves for the entire projection.
10.3.4 Density dependent deterministic projections
Density dependent projections are run similarly to function projection3(). Users may use the two-parameter Ricker function, two-parameter Beverton-Holt function, the Usher function, and the logistic function, and need to specify which transitions to apply these to. Likewise, function f_projection3() can also enforce substochasticity, with both a mild setting forcing survival-transition elements to be in the range [0, 1], and a strict setting forcing survival-transition column sums to this range. Let’s see a simple example.
cyp_ex <- flefko2(year = 2004, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, data = cypfb_env)
c2d <- density_input(cyp_ex, stage3 = c("P1", "P1"),
stage2 = c("SD", "rep"), style = 1, time_delay = 1, alpha = 1, beta = 0.0005,
type = c(2, 2))
trial2f_13 <- f_projection3(cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, times = 101,
year = 2007, ind_terms = mean_real, density = c2d)
> Warning: Option patch not set, so will set to first patch/population.
> Warning: Fewer values of individual covariates have been supplied than required,
> so input values will be cycled.
summary(trial2f_13)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 102 1
We can compare the results to the previous projection assuming the year2 value for 2007, but without density dependence (figure 10.27). Note the extreme impact that density dependence has on the population trajectory, yielding a stable population size somewhere around 25,000.
par(mfrow=c(1,2))
plot(trial2f_2)
title("a)", adj = 0)
plot(trial2f_13, ylab = "")
title("b)", adj = 0)
Figure 10.27: Density independent (a) vs. dependent (b) function-based projections
10.3.5 Ordered / cyclical density dependent simulations
We can use f_projections() to run density dependent projections cycling through a specific set of year2 terms. Here is an example cycling through all terms, which is the default, while using the mean precipitation value for all.
trial2f_14 <- f_projection3(cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, times = 101,
ind_terms = mean_real, density = c2d)
> Warning: Option year not set, so will cycle through existing years.
> Warning: Option patch not set, so will set to first patch/population.
> Warning: Fewer values of individual covariates have been supplied than required,
> so input values will be cycled.
summary(trial2f_14)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 0
> 51 0
> 76 0
> 102 0Here is a further set cycling through just three particular years: 2005, 2006, and 2008.
trial2f_15 <- f_projection3(cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, times = 101,
year = c(2005, 2006, 2008), ind_terms = mean_real, density = c2d)
> Warning: Option patch not set, so will set to first patch/population.
> Warning: Fewer values of individual covariates have been supplied than required,
> so input values will be cycled.
summary(trial2f_15)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 26 0
> 51 0
> 76 0
> 102 0As before, we can also run a specific ordered sequence of years. Here we use the previously defined 20 year vector, initially switching for ten years between 2007 and 2008 before heading on to nine years of 2008 and ending on 2005.
trial2f_16 <- f_projection3(cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, times = 20,
year = year_order, ind_terms = mean_real, density = c2d)
> Warning: Option patch not set, so will set to first patch/population.
> Warning: Fewer values of individual covariates have been supplied than required,
> so input values will be cycled.
summary(trial2f_16)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 20 projected steps per replicate and 1 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 1
> 6 1
> 11 1
> 16 1
> 21 1Now let’s compare the results (figure 10.28).
par(mfrow = c(2, 2))
plot(trial2f_14)
title("a)", adj = 0)
plot(trial2f_15, ylab = "")
title("b)", adj = 0)
plot(trial2f_16)
title("c)", adj = 0)
Figure 10.28: Cyclical and ordered density dependent function-based projections. (a) Cycle of all 6 year values for 101 years; (b) cycle of year values for 2005, 2006, and 2008 for 101 years; and (c) ordered progression for 20 years.
Extinction happens quickly in our two cyclical projections, though we seem to get there after a much lower high population size when cycling through year terms for all years in order. In contrast, our ordered projection seems stable over the 20 projected years.
10.3.6 Replicated simulations assuming temporal stochasticity and density dependence
Finally, we might conduct a density dependent stochastic projection. In this case, we will probably wish to run a simulation with replicates. Here is one such projection using the previous density frame.
set.seed(42)
trial2f_17 <- f_projection3(cypfb_env, format = 3, stochastic = TRUE,
substoch = 2, stageframe = cypframe_fb, supplement = cypsupp2_fb,
modelsuite = cypmodels2_env1, times = 101, nreps = 100,
ind_terms = mean_real, density = c2d)The summary….
summary(trial2f_17)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 100 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 100
> 26 48
> 51 11
> 76 3
> 102 0Here is one further such projection, using a different density dependence relationship.
c2d1 <- density_input(cyp_ex, stage3 = c("P1", "P1"),
stage2 = c("SD", "rep"), style = 1, time_delay = 1, alpha = 1, beta = 0.05,
type = c(2, 2))
trial2f_18 <- f_projection3(cypfb_env, format = 3, stochastic = TRUE,
substoch = 2, stageframe = cypframe_fb, supplement = cypsupp2_fb,
modelsuite = cypmodels2_env1, times = 101, nreps = 100,
ind_terms = mean_real, density = c2d1)summary(trial2f_18)
>
> The input lefkoProj object covers 1 population-patches.
> It includes 101 projected steps per replicate and 100 replicates.
> The number of replicates with population size above the threshold size of 1 is as in
> the following matrix, with pop-patches given by column and milepost times given by row:
> 1 1
> 1 100
> 26 14
> 51 0
> 76 0
> 102 0Let’s take a peek at what the results look like (figure 10.29).
par(mfrow = c(1, 2))
plot(trial2f_17)
plot(trial2f_18, ylab = "")
Figure 10.29: Stochastic density dependent function-based projections, assuming two different sets of values for the Ricket function
The Ricker function has some interesting effects here. We encourage the user to experiment further with different settings.
10.4 Points to remember
- Although deterministic and stochastic analyses are very useful, sometimes more general tools are needed to assess population dynamics. Examples include density dependent analyses, projections under altered climatic regimes, and replicated simulations. For these situations, we use the general projection functions
projection3()andf_projection3(). - Function
projection3()is the quickest simulation tool for most situations, because it uses pre-defined matrices and so does not create matrices at every step.
- The greatest generality is available with function
f_projection3(), because it creates new matrices at every time step based on input vital rate models and user-defined conditions.
- Both projection functions include two substochasticity settings that can prevent impossible values, such as negative fecundity or survival probability above 1.0, from being incorporated into matrix projection.