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 analyses 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 population dynamics successfully beyond just a few time steps 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, shown below.
\[\begin{equation} \mathbf{n_t} = \mathbf{A_t A_{t-1} \cdots A_1 n_0} \tag{10.1} \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 this topic before.
10.1 Density dependence
Density dependent matrix projections involve the modification of matrix elements or vital rates by some function of population size. Many different functions can be applied to incorporate density dependence into matrix projection. The first projections 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 on matrix elements 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 (the function density_vr() can be used to incorporate density dependence characteristics on vital rates, if a function-based MPM is used). Density dependence is assigned to specific elements or sets of elements by noting the transitions to be targeted, using the same format as in function supplemental(). 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 below.
\[\begin{equation} \phi_{t+1} = \phi_t \times \alpha e^{-\beta n_t} \tag{10.2} \end{equation}\]
Here, \(\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 plateaus, single or multi-period oscillations, and even chaos.
The second density dependence function in lefko3 is the two-parameter Beverton-Holt function, shown below.
\[\begin{equation} \phi_{t+1} = \phi_t \times \frac{\alpha}{1 + \beta n_t} \tag{10.3} \end{equation}\]
Here, \(\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 below.
\[\begin{equation} \phi_{t+1} = \phi_t \times \frac{1}{1 + e^{\alpha n_t + \beta}} \tag{10.4} \end{equation}\]
Here, both \(\alpha\) and \(\beta\) give the strength of density dependence, 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 below.
\[\begin{equation} \phi_{t+1} = \phi_t \times (1 - \frac{n_t}{K}) \tag{10.5} \end{equation}\]
The logistic function classically takes only one parameter as input, 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\). These can both be set using alpha to set \(K\) and beta to set whether \(K\) is a hard limit in the density_input() function.
Users working with density dependence forms that involve user-provided exponents should be aware of the limits of their computers to handle particularly large and particularly small numbers. Generally speaking, the largest number capable of being handled by most computers is likely to be 1.7976931^{308}. This roughly corresponds to \(e^{709.7}\). The minimum limit of 2.2250739^{-308} corresponds to \(e^{-708.3}\) If these limits are crossed, then projections will likely still run but yield absurd results, such as switching to negative numbers of individuals in a single time step.
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]. However, simply adjusting matrix elements to the interval [0, 1] may still lead to total survival probabilities greater than 1.0, potentially leading to biased analyses. This happens because the column sums in 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. That is why we have also developed the hard form of substochasticity, in which survival-transition matrix column sums are kept within the interval [0, 1] by adjusting matrix elements by proportionate amounts. 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 the consequences of incorporating logically impossible values in analyses of population dynamics are largely unknown (Caswell 2001).
Let’s consider substochastic forcing 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 we have the following.
\[\begin{equation} U _{ij} = \left\{\begin{array}{cc} U _{ij} & U _{ij} > 0 \\ 0 & U _{ij} < 0 \end{array}\right. \tag{10.6} \end{equation}\]
The above applies for all \(U _{ij}\). The equation is the same for all \(F _{ij}\).
Once this operation is complete, mild substochastic forcing changes all survival-transition elements above 1.0 to 1.0, as in the following for all \(U _{ij}\).
\[\begin{equation} U _{ij} = \left\{\begin{array}{cc} U _{ij} & U _{ij} < 1 \\ 1 & U _{ij} > 1 \end{array}\right. \tag{10.7} \end{equation}\]
Under hard substochastic forcing, we instead identify survival-transition column sums greater than 1.0 and alter survival-transitions, as below for all \(U _{ij}\) and with \(m\) referring to the number of rows (and hence columns, since the matrices must be square).
\[\begin{equation} 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. \tag{10.8} \end{equation}\]
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. The blue box shows the stages that were actually monitored within the field study.
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 54 fecundity transitions were estimated, with 13.5 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 70 fecundity transitions were estimated, with 5.833 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 most common projection styles involves projecting an arithmetic mean matrix deterministically. 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 mean matrices forward. Since these are only single matrices, we do not need to indicate any specific year or patch to project. We will just instruct R to project the ahistorical mean ten time steps forward. Then, we will take a look at the resulting object.
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 with two levels of elements. The top-level elements correspond to the patches or population-patch combinations in the input
lefkoMatobject, with as many elements as there are population-patch combinations. Each of these top-level elements is also a list with number of elements equal to the number of replicates (optionnreps, which defaults to 1). Each element of this lower-level list is a matrix with rows corresponding to the MPM stages (if ahistorical; paired stages if historical; ages if Leslie; and age-stages if age-by-stage), and the number of columns corresponding to the projected times (total number of columns equal to the number of time steps plus 1, with the first time called time 0). The elements of the matrix show the projected number of individuals in that stage / age-stage / historical stage-pair 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 patches or population-patches. Each element of this list is a matrix showing the total population size at each time (column) per replicate (row).
- labels - A data frame giving the population, patch, and year of each matrix in the
lefkoMatobject used as input, in order. - 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).
- 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:
> $milepost_sums
> 1 1
> 1 1
> 3 1
> 6 1
> 8 1
> 11 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 3 1
> 6 1
> 8 1
> 11 1
>
> $extinction_times
> [1] NA
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 non-negative integer. 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, in the impacts of individual fate, or in the impacts of demographic stochasticity 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 in which we start with 1000 seeds and 100 first-year protocorms, and no individuals of any other stage.
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:
> $milepost_sums
> 1 1
> 1 1
> 3 1
> 6 1
> 8 1
> 11 1
>
> $extinction_times
> [1] NA
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 situations, 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 code above, we told R that we wished to start with 1000 3rdyear protocorms and 100 dormant adults. By not stating any further information, we implied zero individuals for 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:
> $milepost_sums
> 1 1
> 1 1
> 3 1
> 6 1
> 8 1
> 11 1
>
> $extinction_times
> [1] NAWe 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_input
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 leftmost plot, it takes a few years for the dormant seeds and 1st year protocorms to mature and the population to start growing. In the rightmost 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 we wish 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 the assumptions of predicted conditions. For example, we might have several years worth of matrices with some under high rain conditions and other under drought. We may believe due to climate forecasts that the next ten years should be relatively rainy while the ten years after that should be drier, and may wish to project our matrices forward in an order that matches these expectation. 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:
> $milepost_sums
> 1 1
> 1 1
> 6 1
> 11 1
> 16 1
> 21 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 2501 1
> 5001 1
> 7501 1
> 10001 1
>
> $extinction_times
> [1] NALet’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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
>
> $extinction_times
> [1] NANow 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 through comparison of the scales of the y axes 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:
> $milepost_sums
> 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
>
> $extinction_times
> [1] NA
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 summarized in the the rightmost column of the milepost_sums, and shows no extinction. Let’s plot these 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 above are produced 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:
> $milepost_sums
> 1 1
> 1 100
> 2501 100
> 5001 100
> 7501 100
> 10001 0
>
> $extinction_times
> [1] NA
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 extinction 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 your computer 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:
> $milepost_sums
> 1 1
> 1 100
> 251 100
> 501 100
> 751 100
> 1001 100
>
> $extinction_times
> [1] NAWe 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:
> $milepost_sums
> 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
>
> $extinction_times
> [1] NA
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+00Note that although these are different distributions, they are nonetheless extremely similar. This reflects the strong stochastic ergodic theorem, which predicts that that the stage structure should converge to a stationary distribution in a stochastic environment (Caswell 2001).
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:
> $milepost_sums
> 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
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 A 1 B 1 C 1 0
> 1 100 100 100 100
> 10 95 100 100 100
> 20 70 100 100 100
> 25 55 100 100 100
> 30 55 100 100 100
> 50 16 100 99 98
> 100 0 83 83 1
> 150 0 24 48 0
> 200 0 6 28 0
> 250 0 1 8 0
> 300 0 0 4 0
>
> $extinction_times
> [1] NALet’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 7643.197205
> [6] 9081.579148 11978.167039 3320.958945 3808.661359 2153.897008
> [11] 10310.177333 7497.822071 5409.348722 4487.853611 1680.404924
> [16] 7203.782417 2318.704951 1498.582135 5040.397381 1737.632999
> [21] 2319.786969 2434.495892 1008.608147 454.399196 754.616703
> [26] 3583.433577 850.056523 680.232056 697.419729 581.051819
> [31] 2399.858548 1797.896250 1292.585451 381.666618 922.366477
> [36] 423.705934 483.055480 822.205618 588.699106 207.467436
> [41] 170.084749 152.944574 219.276019 1011.397505 361.317706
> [46] 967.865891 801.997131 237.946495 92.049747 358.774818
> [51] 472.099252 154.537674 66.314238 28.414543 12.795816
> [56] 108.581301 25.138376 13.802560 56.575191 47.809178
> [61] 55.878980 56.392295 14.108342 10.042440 10.976212
> [66] 6.152709 27.485961 35.446991 34.501840 31.014467
> [71] 9.064451 30.411878 8.595342 32.981367 10.847486
> [76] 3.822600 5.375812 4.877180 21.269629 21.015718
> [81] 5.966223 19.995727 5.676440 11.413514 4.896152
> [86] 15.046584 10.421207 3.010737 2.865675 3.133895
> [91] 1.422278 6.759237 5.958552 2.370643 1.542012
> [96] 1.378010 3.707283 5.180411 3.811777 1.373851
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 population-patch, and the lower level is the replicate. Within the replicate, there is a matrix with rows corresponding to the stages or stage-pairs, and the columns corresponding 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 3.047840e+02
> [79,] 1 2546.2962963 2.515432e+03 2.314603e+03 3.047840e+02
> [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.354173e-01
> [2,] 0.008264463 1.299278e-05 7.953016e-02 5.846836e-02 1.387905e-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.393945e-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.746992e-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.510010e-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.371390e-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.475531e-07
> [52,] 0.008264463 3.368498e-06 4.019704e-07 4.257169e-08 1.844414e-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.338960e-05
> [63,] 0.008264463 1.443642e-05 4.028819e-06 5.066762e-06 1.553538e-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.632668e-04
> [74,] 0.008264463 8.421245e-06 1.800815e-05 2.325904e-05 2.759903e-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 3.987650e-02
> [79,] 0.008264463 5.513910e-02 8.666822e-02 1.013519e-01 3.987650e-02
> [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 4.641563e-06
> [85,] 0.008264463 4.691836e-05 7.692735e-05 8.960614e-05 2.416799e-04
> [86,] 0.008264463 1.082731e-05 9.570724e-06 6.487115e-06 6.057333e-05
> [87,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 1.595060e-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.432971e-01
> [90,] 0.008264463 1.533870e-01 2.464461e-01 2.903997e-01 2.432971e-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.037799e-04
> [98,] 0.008264463 0.000000e+00 0.000000e+00 0.000000e+00 4.542892e-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.883159e-07 9.885868e-07 1.106955e-06 2.096488e-06 2.389431e-06
> [2,] 7.341332e-06 9.215603e-06 9.155864e-06 9.523130e-06 1.897771e-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.515367e-05 7.432931e-05 9.360616e-05 9.264027e-05 9.488778e-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.644859e-04 6.596660e-04 7.549893e-04 9.471198e-04 9.230610e-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.359207e-02 1.345554e-02 1.340092e-02 1.527817e-02 1.887407e-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.359207e-02 1.345554e-02 1.340092e-02 1.527817e-02 1.887407e-02
> [50,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [51,] 1.996688e-02 1.997815e-02 1.806562e-02 2.047292e-02 2.813370e-02
> [52,] 3.024501e-02 3.557370e-02 4.709978e-02 2.851265e-02 1.421137e-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.411851e-03 0.000000e+00
> [63,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 9.604091e-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.589258e-02 0.000000e+00
> [66,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [67,] 1.007844e-06 1.010383e-06 1.148395e-06 2.144449e-06 2.414708e-06
> [68,] 7.341332e-06 9.215603e-06 9.155864e-06 9.523130e-06 1.897771e-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.116926e-02 2.118584e-02 1.943840e-02 2.047292e-02 2.813370e-02
> [74,] 3.585611e-02 4.120960e-02 5.350609e-02 2.851265e-02 1.421137e-02
> [75,] 0.000000e+00 0.000000e+00 0.000000e+00 9.630101e-02 0.000000e+00
> [76,] 1.117750e-02 4.104774e-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.007844e-06 1.010383e-06 1.148395e-06 2.144449e-06 2.414708e-06
> [79,] 7.341332e-06 9.215603e-06 9.155864e-06 9.523130e-06 1.897771e-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.425907e-03 0.000000e+00
> [84,] 2.247720e-02 2.313248e-02 1.987204e-02 2.122501e-02 3.513143e-02
> [85,] 5.012610e-02 5.733413e-02 7.136429e-02 9.179451e-02 9.775391e-02
> [86,] 1.460718e-01 1.521127e-01 1.155540e-01 2.104936e-01 6.001083e-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.007844e-06 1.010383e-06 1.148395e-06 2.144449e-06 2.414708e-06
> [90,] 7.341332e-06 9.215603e-06 9.155864e-06 9.523130e-06 1.897771e-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.416613e-02 4.052185e-02
> [96,] 7.080763e-02 7.328743e-02 1.034554e-01 3.796758e-02 8.167958e-02
> [97,] 2.993882e-01 2.826221e-01 2.568144e-01 2.283851e-01 2.878609e-01
> [98,] 5.753146e-02 5.211191e-02 5.917720e-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.007844e-06 1.010383e-06 1.148395e-06 2.144449e-06 2.414708e-06
> [101,] 7.341332e-06 9.215603e-06 9.155864e-06 9.523130e-06 1.897771e-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.849602e-02 4.637516e-02 7.895477e-02 3.270260e-02 0.000000e+00
> [108,] 2.412809e-02 3.395090e-02 1.250809e-02 0.000000e+00 1.083988e-01
> [109,] 4.327867e-02 4.692303e-02 2.999251e-02 3.378413e-02 4.750575e-02
> [110,] 6.448293e-02 6.548406e-02 6.726622e-02 7.181279e-02 0.000000e+00
> [111,] 1.007844e-06 1.010383e-06 1.148395e-06 2.144449e-06 2.414708e-06
> [112,] 7.341332e-06 9.215603e-06 9.155864e-06 9.523130e-06 1.897771e-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.735823e-02
> [121,] 1.683328e-02 1.690771e-02 1.921892e-02 0.000000e+00 8.059011e-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, as we reach higher numbers more quickly under that assumption (note the difference in scale of the x axes).
10.2.4 Density dependent deterministic projections
Let’s now move on to the simplest density dependent analysis: a single matrix projected forward with some elements altered at each time step according to some function of population size, \(N\). We might conduct such an analysis 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 trends 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 alpha beta time_delay type type_t12
> 1 P1 SD NA NA 1 0.05 0 1 2 1
> 2 P1 XSm NA NA 1 0.05 0 1 2 1
> 3 P1 Sm NA NA 1 0.05 0 1 2 1
> 4 P1 Md NA NA 1 0.05 0 1 2 1
> 5 P1 Lg NA NA 1 0.05 0 1 2 1
> 6 P1 XLg NA NA 1 0.05 0 1 2 1
c2d_2
> stage3 stage2 stage1 age2 style alpha beta time_delay type type_t12
> 1 P1 SD NA NA 1 1 0 1 2 1
> 2 P1 XSm NA NA 1 1 0 1 2 1
> 3 P1 Sm NA NA 1 1 0 1 2 1
> 4 P1 Md NA NA 1 1 0 1 2 1
> 5 P1 Lg NA NA 1 1 0 1 2 1
> 6 P1 XLg NA NA 1 1 0 1 2 1
c2d_3
> stage3 stage2 stage1 age2 style alpha beta time_delay type type_t12
> 1 P1 SD NA NA 1 0.05 5e-04 1 2 1
> 2 P1 XSm NA NA 1 0.05 5e-04 1 2 1
> 3 P1 Sm NA NA 1 0.05 5e-04 1 2 1
> 4 P1 Md NA NA 1 0.05 5e-04 1 2 1
> 5 P1 Lg NA NA 1 0.05 5e-04 1 2 1
> 6 P1 XLg NA NA 1 0.05 5e-04 1 2 1
c2d_4
> stage3 stage2 stage1 age2 style alpha beta time_delay type type_t12
> 1 P1 SD NA NA 1 1 5e-04 1 2 1
> 2 P1 XSm NA NA 1 1 5e-04 1 2 1
> 3 P1 Sm NA NA 1 1 5e-04 1 2 1
> 4 P1 Md NA NA 1 1 5e-04 1 2 1
> 5 P1 Lg NA NA 1 1 5e-04 1 2 1
> 6 P1 XLg NA NA 1 1 5e-04 1 2 1
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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
>
> $extinction_times
> [1] NALet’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 to extinction (b, d). 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:
> $milepost_sums
> 1 1
> 1 1
> 6 1
> 11 1
> 16 1
> 21 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 2501 0
> 5001 0
> 7501 0
> 10001 0
>
> $extinction_times
> [1] NALet’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 form of density dependence leads to extinction fairly early on, though only after some strong fluctuations.
Projecting a full MPM with multiple annual matrices defaults to cycling through the matrices in order, even when density dependence is 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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
>
> $extinction_times
> [1] NALet’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 b) 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 with 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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
>
> $extinction_times
> [1] NAOnce 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 influence 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 some 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 projection3() allows users to project matrices forward in time, but is limited to the case in which we have already developed our MPM and wish to project using only those matrices. What if we wished to create a new matrix at each time step, perhaps because of some changing conditions that should alter basic vital rates? Function-based projection through function f_projection3() provides us with the flexibility to construct custom matrices at each time step given some set of conditions that we are interested in. This approach is 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 an individual or environmental covariate. For example, we may have developed a set of vital rate models that have relationships 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.3).
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. 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_hfv(cypfb_env)
>
> This hfv dataset contains 320 rows, 57 variables, 1 population,
> 3 patches, 74 individuals, and 5 time steps.
> 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 # Mean number of seeds per fruit
sl_mult <- 0.7 # Product of assumed germination probability and survival to 1st year
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.
hfv_qc(data = cypfb_env, vitalrates = c("surv", "obs", "size", "repst", "fec"),
size = c("size3added", "size2added", "size1added"),
indcova = c("indcova3", "indcova2", "indcova1"))
> Survival:
>
> Data subset has 58 variables and 320 transitions.
>
> Variable alive3 has 0 missing values.
> Variable alive3 is a binomial variable.
>
>
> Observation status:
>
> Data subset has 58 variables and 303 transitions.
>
> Variable obsstatus3 has 0 missing values.
> Variable obsstatus3 is a binomial variable.
>
>
> Primary size:
>
> Data subset has 58 variables and 288 transitions.
>
> Variable size3added has 0 missing values.
> Variable size3added appears to be an integer variable.
>
> Variable size3added is fully positive, lacking even 0s.
>
> Overdispersion test:
> Mean size3added is 3.653
> The variance in size3added is 13.41
> The probability of this dispersion level by chance assuming that
> the true mean size3added = variance in size3added,
> and an alternative hypothesis of overdispersion, is 3.721e-138
> Variable size3added is significantly overdispersed.
>
> Zero-inflation and truncation tests:
> Mean lambda in size3added is 0.02592
> The actual number of 0s in size3added is 0
> The expected number of 0s in size3added under the null hypothesis is 7.465
> The probability of this deviation in 0s from expectation by chance is 0.9964
> Variable size3added is not significantly zero-inflated.
>
> Variable size3added does not include 0s, suggesting that a zero-truncated distribution may be warranted.
>
>
> Reproductive status:
>
> Data subset has 58 variables and 288 transitions.
>
> Variable repstatus3 has 0 missing values.
> Variable repstatus3 is a binomial variable.
>
>
> Fecundity:
>
> Data subset has 58 variables and 118 transitions.
>
> Variable feca2 has 0 missing values.
> Variable feca2 appears to be an integer variable.
>
> Variable feca2 is fully non-negative.
>
> Overdispersion test:
> Mean feca2 is 0.7881
> The variance in feca2 is 1.536
> The probability of this dispersion level by chance assuming that
> the true mean feca2 = variance in feca2,
> and an alternative hypothesis of overdispersion, is 0.1193
> Dispersion level in feca2 matches expectation.
>
> Zero-inflation and truncation tests:
> Mean lambda in feca2 is 0.4547
> The actual number of 0s in feca2 is 68
> The expected number of 0s in feca2 under the null hypothesis is 53.65
> The probability of this deviation in 0s from expectation by chance is 5.904e-06
> Variable feca2 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. All other response variables appear to match the expectations of the binomial distribution, and so are fine.
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 model estimated with 74 individuals and 320 individual transitions.
> Survival model accuracy is 1.
> Observation status model estimated with 70 individuals and 303 individual transitions.
> Observation status model accuracy is 0.95.
> Primary size model estimated with 70 individuals and 288 individual transitions.
> Primary size model R-squared is 0.829.
> Secondary size model not estimated.
> Tertiary size model not estimated.
> Reproductive status model estimated with 70 individuals and 288 individual transitions.
> Reproductive status model accuracy is 0.719.
> Fecundity model estimated with 51 individuals and 118 individual transitions.
> Fecundity model R-squared is 0.549.
> Juvenile survival model not estimated.
> Juvenile observation status model not estimated.
> Juvenile primary size model not estimated.
> Juvenile secondary size model not estimated.
> Juvenile tertiary size model not estimated.
> Juvenile reproductive status model not estimated.
> Juvenile maturity status model 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 reasonable accuracy or 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 will 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(data = 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 itself 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 in element 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 in element projection. Only estimated if
growthonly = FALSEandrepvalue = TRUE. - pop_size - A list with number of elements equal to the number of projected population-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 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 MPM 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 showing the order of 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
densityoption. This data frame shows density dependence per matrix element, if such information was provided. - density_vr - The data frame input under the
density_vroption. This data frame shows density dependence per vital rate model, if such information was provided.
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:
> $milepost_sums
> 1 1
> 1 1
> 6 1
> 11 1
> 16 1
> 22 1
>
> $extinction_times
> [1] NALet’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(data = 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:
> $milepost_sums
> 1 1
> 1 1
> 6 1
> 11 0
> 16 0
> 22 0
>
> $extinction_times
> [1] NANow 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 the more powerful approach. We will start the population with 1000 dormant seeds and 100 1st year protocorms, using the 2010 to 2030 climate projection (figure 10.21).Note that prior to using this option, we need to create a new MPM in the same format as our projection, as the start_input() function will need to look at the structure of the matrices and row / column indices to determine the proper vector.
trial_mpm <- flefko2(data = cypfb_env, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, year = 2007)
c2m_sv <- start_input(trial_mpm, stage2 = c("SD", "P1"), value = c(1000, 100))
trial2f_3 <- f_projection3(data = 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 in population size until some seeds and protocorms start reproducing. Once reproductive individuals enter the population, 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 t. These are terms within 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: the first is a vector of six precipitation values, the second is just the first of those values, and the third is the mean of those six values.
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(data = 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(data = 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(data = 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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 102 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 102 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 102 1
>
> $extinction_times
> [1] NAHere is another example using the full set of years, though assuming the same three precipitation vectors as the last runs.
trial2f_5a <- f_projection3(data = 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(data = 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(data = 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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 102 1
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 26 0
> 51 0
> 76 0
> 102 0
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 1
> 26 0
> 51 0
> 76 0
> 102 0
>
> $extinction_times
> [1] NALet’s plot the results of all six runs 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. Limiting the runs to using years 2006, 2007, and 2008 led to projections that grew through the projection, though with a small single-period cycle, no what matter which precipitation projection was used. Using all years made the choice of precipitation projection crucial, as only the first projection assuming a six-year precipitation cycle led to population growth while all others led to extinction.
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(data = 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.
summary(trial2f_6)
>
> 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:
> $milepost_sums
> 1 1
> 1 1
> 6 1
> 11 1
> 16 1
> 21 1
>
> $extinction_times
> [1] NA
As before, if we set a larger number for times than we will cycle through the sequence.
trial2f_7 <- f_projection3(data = 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.
summary(trial2f_7)
>
> 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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 101 1
>
> $extinction_times
> [1] NANo extinction happened in either case. Let’s take a look at plots of these projections for further details (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 large upswings as time goes on.
10.3.3 Replicated simulations assuming temporal stochasticity
Function f_projection3() can also run 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(data = 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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 0
> 102 0
>
> $extinction_times
> [1] NALet’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(data = 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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 102 0
>
> $extinction_times
> [1] NALet’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. Graphically, it appears that the population also goes extinct at a later time when fecundity is doubled. 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(data = 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(data = 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:
> $milepost_sums
> 1 1
> 1 100
> 26 80
> 51 41
> 76 16
> 102 3
>
> $extinction_times
> [1] NA
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:
> $milepost_sums
> 1 1
> 1 100
> 26 87
> 51 68
> 76 47
> 102 29
>
> $extinction_times
> [1] NAOur summaries now include information on the numbers of replicates still alive at each milepost time. They did before as well, but we only ran single replicates before. Let’s compare these results graphically (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(data = 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:
> $milepost_sums
> 1 1
> 1 100
> 26 100
> 51 100
> 76 100
> 102 100
>
> $extinction_times
> [1] NALet’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 can run in two ways using function f_projection3(). The first involves setting density dependence relationships in matrix elements, just as is done with 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. Remember that, as in start_input(), we need to provide density_input() with a MPM in with the same characteristics as the projection that we will run, so that the function can infer the proper stage distribution and produce the proper output vector.
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(data = 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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 102 1
>
> $extinction_times
> [1] NA
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
Function f_projection3() can also run density dependent simulations assuming that vital rates are density dependent. This situation allows us to modify the vital rates used to estimate the matrix elements. As before, it is possible to use the two-parameter Ricker function, the two-parameter Beverton-Holt function, the Usher function, or the logistic function. Let’s try a simulation in which we make survival and fecundity density dependent, using the Ricker function for both. We will use the density_vr() function for this purpose.
vr1 <- density_vr(density_yn = c(T, F, F, F, F, F, T, F, F, F, F, F, F, F),
style = c(1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0),
alpha = c(1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0),
beta = c(-0.0000000005, 0, 0, 0, 0, 0, 0.0000000003, 0, 0, 0, 0, 0, 0, 0))
trial2f_14 <- f_projection3(data = cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, times = 101,
year = 2007, ind_terms = mean_real, density_vr = vr1, substoch = 2)
> 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:
> $milepost_sums
> 1 1
> 1 1
> 26 1
> 51 1
> 76 1
> 102 1
>
> $extinction_times
> [1] NA
In our density_vr() call, we need to set the density dependence relationships for 14 different vital rates that are used by the function-based matrix estimators in lefko3. The 14 models, in order, are: survival, observation, primary size, secondary size, tertiary size, reproduction status, and fecundity, followed by juvenile survival, juvenile observation, juvenile primary size, juvenile secondary size, juvenile tertiary size, juvenile reproduction status, and juvenile maturity status. The density_yn option stipulates with TRUE and FALSE statements which of these will be modified, and the style, alpha, beta, and time_delay option give the settings to use (last option not shown). Let’s plot the results.
par(mfrow=c(1,2))
plot(trial2f_2)
title("a)", adj = 0)
plot(trial2f_14, ylab = "")
title("b)", adj = 0)
Figure 10.28: Density independent (a) vs. dependent (b) function-based projections with density dependent vital rates
We see a dramatic difference in projected population dynamics, resulting from two density dependent vital rates.
10.3.5 Ordered / cyclical density dependent simulations
We can also use f_projection() 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_15 <- f_projection3(data = 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_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:
> $milepost_sums
> 1 1
> 1 1
> 26 0
> 51 0
> 76 0
> 102 0
>
> $extinction_times
> [1] NAHere is a further set cycling through just three particular years: 2005, 2006, and 2008.
trial2f_16 <- f_projection3(data = 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_16)
>
> 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:
> $milepost_sums
> 1 1
> 1 1
> 26 0
> 51 0
> 76 0
> 102 0
>
> $extinction_times
> [1] NAAs 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_17 <- f_projection3(data = 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_17)
>
> 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:
> $milepost_sums
> 1 1
> 1 1
> 6 1
> 11 1
> 16 1
> 21 1
>
> $extinction_times
> [1] NANow let’s compare the results (figure 10.29).
par(mfrow = c(2, 2))
plot(trial2f_15)
title("a)", adj = 0)
plot(trial2f_16, ylab = "")
title("b)", adj = 0)
plot(trial2f_17)
title("c)", adj = 0)
Figure 10.29: 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.
As before, we can use density dependence in vital rates even with a cyclical or ordered simulation. In that circumstance, the density option simply needs to be replaced with an appropriate density_vr input.
trial2f_18 <- f_projection3(data = cypfb_env, format = 3, stageframe = cypframe_fb,
supplement = cypsupp2_fb, modelsuite = cypmodels2_env1, times = 20,
year = year_order, ind_terms = mean_real, density_vr = vr1)
> 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.
plot(trial2f_18)
Figure 10.30: Ordered function-based projection with vital rate density dependence.
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_19 <- f_projection3(data = 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_19)
>
> 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:
> $milepost_sums
> 1 1
> 1 100
> 26 48
> 51 11
> 76 3
> 102 0
>
> $extinction_times
> [1] NAHere 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_20 <- f_projection3(data = 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_20)
>
> 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:
> $milepost_sums
> 1 1
> 1 100
> 26 14
> 51 0
> 76 0
> 102 0
>
> $extinction_times
> [1] NALet’s take a peek at what the results look like (figure 10.30).
par(mfrow = c(1, 2))
plot(trial2f_19)
plot(trial2f_20, ylab = "")
Figure 10.31: 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.
- Users may set density dependence in matrix elements or in vital rates, using the
density_input()anddensity_vr()functions, respectively.