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.

— Kenneth Boulding

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:

  1. Deterministic projections of a specific time length
  2. Ordered or cyclical projections of a specific time length
  3. Replicated simulations assuming temporal stochasticity
  4. Density dependent deterministic projections
  5. Ordered or cyclical density dependent simulations
  6. 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.000

Now 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:

  1. projection - A list with two levels of elements. The top-level elements correspond to the patches or population-patch combinations in the input lefkoMat object, 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 (option nreps, 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.
  2. stage_dist - A list showing the stage distribution projected per time, with the same order and structure as projection. Only estimated if growthonly = FALSE.
  3. 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.
  4. 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).
  5. labels - A data frame giving the population, patch, and year of each matrix in the lefkoMat object used as input, in order.
  6. ahstages - The stageframe used in analysis, though modified and edited per rules set for function-based matrix estimators.
  7. hstages - A data frame matrix showing the pairing of ahistorical stages used to create historical stage pairs (only if historical).
  8. agestages - A data frame showing the order of age-stage pairs (only if age-by-stage).
  9. control - An integer vector indicating the number of replicates and the number of time steps.
  10. 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)
Ahistorical (a) vs. historical (b) deterministic projection

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)
Historical projection allowing integer population size only

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   100

In 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] NA

We 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)
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

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] NA

Let’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)
Ordered projections can be cyclical. (a) A 20 year ordered projection. (b) A 10,000 year projection cycling through 20 ordered years

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] NA

Now 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)
Cyclical ahistorical (a) and historical (b) projections

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")
Patch-level projections for patch A (top left), patch B (top right), patch C (bottom left), and the patch-weighted population (bottom right)

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     0

10.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))
Stochastic projection with exploding population size

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] NA

We 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)
Stochastic ahistorical (a) vs. historical (b) projection

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+00

Note 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] NA

Let’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-02

This 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")
Stochastic patch-level projections

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)
Stochastic projections using equal (a) or different (b) matrix weights

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] NA

Let’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)
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.

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)
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.

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] NA

Let’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)
Ordered,density dependent projections can be cyclical. (a) A 20 year ordered projection. (b) A 10,000 year projection cycling through 20 ordered years

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] NA

Let’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)
Cyclical density dependent projections using the Ricker function

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] NA

Once 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)
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.

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)
Stochastic density dependent projections with replicates assuming the two-parameter Ricker function. (a) alpha = 0.05, beta = 0; and (b) alpha = 1, beta = 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")
Predicted annual precipitation for Cypripedium candidum

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:

  1. 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.
  2. 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.
  3. 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 = FALSE and repvalue = TRUE.
  4. 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).
  5. ahstages - The stageframe used in analysis, though modified and edited per rules set for function-based MPM estimators.
  6. hstages - A data frame matrix showing the pairing of ahistorical stages used to create historical stage pairs (only if historical).
  7. agestages - A data frame showing the order of age-stage pairs (only if age-by-stage).
  8. labels - A data frame showing the order of population and patch identities used.
  9. control - An integer vector indicating the number of replicates and the number of time steps.
  10. density - The data frame input under the density option. This data frame shows density dependence per matrix element, if such information was provided.
  11. density_vr - The data frame input under the density_vr option. 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] NA

Let’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] NA

Now 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)
Deterministic function-based projections for (a) 2010-2030 and (b) 2070-2090

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)
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

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] NA

Here 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] NA

Let’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)
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).

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] NA

No 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)
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.

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] NA

Let’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] NA

Let’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)
Stochastic function-based projections testing the impacts of natural fecundity values (a) vs. doubled fecundity (b)

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] NA

Our 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 = "")
Stochastic function-based projections testing doubled fecundity with replicates

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] NA

Let’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)
Stochastic function-based projections assuming (a) equal weights for all years, and (b) user-defined weights for different years

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)
Density independent (a) vs. dependent (b) function-based projections

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)
Density independent (a) vs. dependent (b) function-based projections with density dependent vital rates

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] NA

Here 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] NA

As 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] NA

Now 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)
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.

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)
Ordered function-based projection with vital rate density dependence.

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] NA

Here 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] NA

Let’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 = "")
Stochastic density dependent function-based projections, assuming two different sets of values for the Ricket function

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

  1. 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() and f_projection3().
  2. 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.
  3. 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.
  4. 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.
  5. Users may set density dependence in matrix elements or in vital rates, using the density_input() and density_vr() functions, respectively.

References

Caswell, H. (2001). Matrix population models: Construction, analysis, and interpretation. Second edition. Sinauer Associates, Inc., Sunderland, Massachusetts, USA.
Crone, E.E., Ellis, M.M., Morris, W.F., Stanley, A., Bell, T., Bierzychudek, P., et al. (2013). Ability of matrix models to explain the past and predict the future of plant populations. Conservation Biology, 27, 968–978.
Jensen, A.L. (1995). Simple density-dependent matrix model for population projection. Ecological Modelling, 77, 43–48.
Leslie, P.H. (1959). The properties of a certain lag type of population growth and the influence of an external random factor on a number of such populations. Physiological Zoology, 32, 151–159.
Shefferson, R.P., Mizuta, R. & Hutchings, M.J. (2017). Predicting evolution in response to climate change: The example of sprouting probability in three dormancy-prone orchid species. Royal Society Open Science, 4, 160647.