Chapter 5 Matrix Models II: Function-based MPMs

“There is a simple rule here, a rule of legislation, a rule of business, a rule of life: beyond a certain point, complexity is fraud.”

— P.J. O’Rourke

Matrix projection models, whether historical or ahistorical, can be categorized as either raw or function-based. The oldest and most common in the literature is the raw MPM, and chapter 4 discusses that approach in detail. In this chapter, we will cover the function-based approach.

Function-based MPMs are created differently than raw MPMs. Although they can look the same to the untrained eye, and although the matrix elements themselves still represent survival-transition probabilities and fecundity rates, matrix elements are calculated differently. Rather than estimating matrix elements as proportions of individuals moving between stages at consecutive times or as products of average offspring production from the original dataset, matrix elements in function-based MPMs are estimated using functions that represent the vital rates or demographic processes underlying the dynamics of the population. These functions are linked together via kernels that treat them typically as conditional rates and probabilities.

How are these vital rate functions estimated? Let’s consider a simple function-based size-classified model. In this model, the survival-transition probability is estimated as the product of the probability that an individual in stage \(j\) at occasion t survives to occasion t+1, and the probability that that same individual becomes stage \(k\) in time t+1 assuming survival to that time. The latter probability is conditional on the first, and so can be treated as independent and simply multiplied by it. Thus, we have the following product, where \(\sigma _{j}\) is the survival probability of an individual in stage \(j\) at time t to time t+1, and \(\psi _{kj}\) is the probability of transition from stage \(j\) in time t to stage \(k\) in time t+1 assuming survival to time t+1:

\[\begin{equation} \tag{5.1} a _{kj} = \sigma _{j} \psi _{kj} \end{equation}\]

Continuing with this simple case, fecundity elements assuming a pre-breeding model are given as

\[\begin{equation} \tag{5.2} a _{kj} = f_{kj} \sigma _{k.} \end{equation}\]

and assuming a post-breeding model are given as

\[\begin{equation} \tag{5.3} a _{kj} = \sigma _{.j} f_{kj} \end{equation}\]

where \(f_{kj}\) is the production of stage \(k\) offspring by a reproductive individual in stage \(j\) in time t, \(\sigma _{k.}\) is the survival probability of newborns to age 1 at time t+1 when they are first censused (assuming a pre-breeding model), and \(\sigma _{.j}\) is the survival of the parent from the occasion preceding reproduction to the time of reproduction (assuming a post-breeding model). Note that reproductive stages are defined differently in the post-breeding model than in the pre-breeding model to account for the need to incorporate parental survival in the former but not in the latter (see Chapter 2 for further details).

Complications in the life history model generally increase the terms to be multiplied in each product to produce the corresponding element. For example, if we are dealing with an herbaceous plant that undergoes vegetative dormancy, then we might change the survival-transition kernel to two different forms, with the equation used dependent on whether the stage we are transitioning to in time t+1 is dormancy or not. A dormant plant has no aboveground size, and so the equation used for dormancy would be the following instead of equation 5.1, where \(\rho _{.}\) is the probability of sprouting in time t+1 assuming survival from time t and stage \(D\) is the dormant stage:

\[\begin{equation} \tag{5.4} a _{Dj} = \sigma _{j} (1 - \rho _{.}) \end{equation}\]

The probability of transitioning to a non-dormant stage \(k\) is given as

\[\begin{equation} \tag{5.5} a _{kj} = \sigma _{j} \rho _{.} \psi _{kj} \end{equation}\]

\(\psi _{kj}\) is the transition probability from stage \(j\) in time t to stage \(k\) in time t+1 assuming both survival to time t+1 and sprouting in time t+1. This makes it a little different than \(\psi _{kj}\) in equation 5.1, where we assumed a life history that did not involve the possibility of unobservable life stages.

Size in function-based MPMs is handled using size classes, which often have size bins set to standardized widths. Such widths may be set to 1.0, as might be typical of count variables, but they can also be set narrower or wider. The function-based approach is very flexible, and allows stages and their associated size classes to be cut quite finely by using different vital rate functions on a case-by-case basis. For example, we may wish to create two sets of stages that have the same size classes, but differ in reproductive status. In the herbaceous perennial plant case that we have been working with, such a matrix would be composed of survival-transition probabilities determined by three functions. The first would be given by equation 5.4, for the dormant plant case. The second would be the following equation, corresponding to the transition to non-reproductive stage \(k\), where \(\gamma _{k}\) is the probability of becoming reproductive in time t+1 of an individual in size class \(k\) in that time:

\[\begin{equation} \tag{5.6} a _{kj} = \sigma _{j} \rho _{.} \psi _{kj} (1 - \gamma _{k}) \end{equation}\]

Finally, we have the following equation giving the survival-transition probability to reproductive stage \(k\):

\[\begin{equation} \tag{5.7} a _{kj} = \sigma _{j} \rho _{.} \psi _{kj} \gamma _{k} \end{equation}\]

The functions used in these equations can be developed in several ways. In lefko3, we utilize the linear modeling approach. In this approach, we use hfv-formatted datasets to analyze the main vital rates via either generalized linear models or generalized linear mixed models. The best-fit models are then used as inputs into general functions created for the job, with a great deal of flexibility provided for the assumptions inherent in the system being analyzed.

5.1 Developing vital rate models

Package lefko3 employs two methods to estimate the vital rates used to create function-based MPMs. The first, main method is automated vital rate modeling. Alternatively, users may create their own vital rate models separately and incorporate them into MPM modeling. In lefko3, the function modelsearch() is the main tool provided to estimate vital rate functions via the former method. It also allows users to explore whether history influences component vital rates within a rigorous statistical framework. The results can be used not only to decide whether a historical MPM is justified, but also to develop function-based hMPMs and ahMPMs, including age-by-stage MPMs and even integral projection models (IPMs). Package lefko3 can estimate linear models to estimate up to 14 different vital rates:

  1. Survival probability - This is the probability of surviving from occasion t to occasion t+1, given that the individual is in stage \(j\) in occasion t (and, if historical, in stage \(l\) in occasion t-1). In lefko3, this parameter may be modeled as a function of up to three size metrics, reproductive status, patch, year, age, and individual identity, and a number of individual or environmental covariates in occasions t and t-1. This parameter is required in all function-based matrices.
  2. Observation probability - This is the probability of observation in occasion t+1 of an individual in stage \(k\) given survival from occasion t to occasion t+1. This parameter is only used when at least one stage is technically not observable. For example, some plants are capable of vegetative dormancy, in which case they are alive but do not necessarily sprout in all years. In these cases, the probability of sprouting may be estimated as the observation probability. Note that this probability does not refer to observer effort, and so should only be used to differentiate completely unobservable stages where the observation status refers to an important biological phenomenon, such as when individuals may be alive but have a size of zero. In lefko3, this parameter may be modeled as a function of up to three size metrics, reproductive status, patch, year, age, and individual identity, and a number of individual or environmental covariates in occasions t and t-1.
  3. Primary size transition probability - This is the probability of becoming size \(k\) in occasion t+1 assuming survival from occasion t to occasion t+1 and observation in that time. If multiple size metrics are used, then this refers only to the first of these, which we may refer to as the primary size variable. In lefko3, this parameter may be modeled as a function of up to three size metrics, reproductive status, patch, year, age, and individual identity, and a number of individual or environmental covariates in occasions t and t-1. This parameter is required in all function-based size-classified matrices.
  4. Secondary size transition probability - This is the probability of becoming size \(k\) in occasion t+1 assuming survival from occasion t to occasion t+1 and observation in that time, within a second size metric used for classification in addition to the primary metric. In lefko3, this parameter may be modeled as a function of up to three size metrics, reproductive status, patch, year, age, and individual identity, and a number of individual or environmental covariates in occasions t and t-1.
  5. Tertiary size transition probability - This is the probability of becoming size \(k\) in occasion t+1 assuming survival from occasion t to occasion t+1 and observation in that time, within a third size metric used for classification in addition to the primary and secondary metrics. In lefko3, this parameter may be modeled as a function of up to three size metrics, reproductive status, patch, year, age, and individual identity, and a number of individual or environmental covariates in occasions t and t-1.
  6. Reproduction probability - This is the probability of becoming reproductive in occasion t+1 given survival from occasion t to occasion t+1, and observation in that time. Note that this should be used only if the researcher wishes to separate breeding from non-breeding mature stages. If all adult stages are potentially reproductive and no separation of reproducing from non-reproducing adults is required by the life history model, then this parameter should not be estimated. In lefko3, this parameter may be modeled as a function of up to three size metrics, reproductive status, patch, year, age, and individual identity, and a number of individual or environmental covariates in occasions t and t-1.
  7. Fecundity rate - Under the default setting, this is the rate of successful production of offspring in occasion t by individuals alive, observable, and reproductive in that time, and, if assuming a pre-breeding model and sufficient information is provided in the dataset, the survival of those offspring into occasion t+1 in whatever juvenile class is possible. Thus, the fecundity rate of seed-producing plants might be split into seedlings, which are plants that germinated within a year of seed production, and dormant seeds. Alternatively, it may be given only as produced fruits or seeds, with the survival and germination of seeds provided elsewhere in the MPM development process, such as within a supplement table. An additional setting allows fecundity rate to be estimated using data provided for occasion t+1 instead of occasion t. In lefko3, this parameter may be modeled as a function of up to three size metrics, reproductive status, patch, year, age, and individual identity, and a number of individual or environmental covariates in occasions t and t-1.
  8. Juvenile survival probability - This is the probability of surviving from juvenile stage \(j\) in occasion t to occasion t+1. It is used only when the user wishes to model vital rates for a single size-unclassified juvenile period separately from adults. In lefko3, this parameter may be modeled as a function of up to three size metrics, patch, year, and individual identity, and a number of individual or environmental covariates in occasions t and t-1.
  9. Juvenile observation probability - This is the probability of observation in occasion t+1 of an individual in juvenile stage \(j\) in occasion t given survival from occasion t to occasion t+1. It is used only when the user wishes to model vital rates for a single size-unclassified juvenile period separately from adults. In lefko3, this parameter may be modeled as a function of up to three size metrics, patch, year, and individual identity, and a number of individual or environmental covariates in occasions t and t-1.
  10. Juvenile primary size transition probability - This is the probability of becoming size-classified stage \(k\) in occasion t+1 assuming survival from juvenile stage \(j\) in occasion t to occasion t+1 and observation in that time. It is in terms of a single size metric. It is used only when the user wishes to model vital rates for a single size-unclassified juvenile period separately from adults. In lefko3, this parameter may be modeled as a function of up to three size metrics, patch, year, and a number of individual or environmental covariates in occasions t and t-1, and individual identity.
  11. Juvenile secondary size transition probability - This is the probability of becoming size-classified stage \(k\) in occasion t+1 assuming survival from juvenile stage \(j\) in occasion t to occasion t+1 and observation in that time, in a secondary size metric in addition to the primary size metric. It is used only when the user wishes to model vital rates for a single size-unclassified juvenile period separately from adults. In lefko3, this parameter may be modeled as a function of up to three size metrics, patch, year, and individual identity, and a number of individual or environmental covariates in occasions t and t-1.
  12. Juvenile tertiary size transition probability - This is the probability of becoming size-classified stage \(k\) in occasion t+1 assuming survival from juvenile stage \(j\) in occasion t to occasion t+1 and observation in that time, in a tertiary size metric in addition to the primary and secondary size metrics. It is used only when the user wishes to model vital rates for a single size-unclassified juvenile period separately from adults. In lefko3, this parameter may be modeled as a function of up to three size metrics, patch, year, and individual identity, and a number of individual or environmental covariates in occasions t and t-1.
  13. Juvenile reproduction probability - This is the probability of reproducing in mature stage \(k\) in occasion t+1 given survival from juvenile stage \(j\) in occasion t to occasion t+1, and observation in that time. It is used only when the user wishes to model vital rates for a single size-unclassified juvenile period separately from adults. In lefko3, this parameter may be modeled as a function of up to three size metrics, patch, year, and individual identity, and a number of individual or environmental covariates in occasions t and t-1.
  14. Juvenile maturity probability - This is the probability of becoming mature in occasion t+1 given survival from juvenile stage \(j\) in occasion t to occasion t+1. It is used only when the user wishes to model vital rates for a single size-unclassified juvenile period separately from adults. In lefko3, this parameter may be modeled as a function of up to three size metrics, patch, year, and individual identity, and a number of individual or environmental covariates in occasions t and t-1. Note that this parameter denotes transition to maturity.

Of these fourteen vital rates, most users will estimate at least parameters (1) survival probability, (3) size transition probability, and (7) fecundity. These three are the default set for function modelsearch(). Parameters (2) observation probability and (6) reproduction probability may be used when some stages are included that are completely unobservable (and so do not have any size), or that are mature but non-reproductive, respectively. Parameters (4) secondary size transition and (5) tertiary size transition should only be used when size classification involves more than one size variable. Parameters (8) through (14) should only be added if the dataset contains juvenile individuals transitioning to maturity, and these juveniles live essentially as a single juvenile stage for some amount of time before transitioning to maturity, or before transitioning to a stage that is size-classified in the same manner as adult stages are. If juveniles can be classified by size similarly to adults (or at least on the same scale), then only vital rates (1) through (7) should be used. If multiple juvenile stages exist on a different size classification system than adults, then we advise the user to model only rates (1) through (7) and to use stage groups to stratify vital rate models properly.

5.1.1 Options in modeling vital rates with modelsearch()

Function modelsearch() handles the entire modeling process using an exhaustive model building process combined with information theoretical model selection using model AICc and the number of estimated parameters in each model. It first develops global models, and then develops all nested models via exhaustive model building. Finally, it determines and selects the best-fit models. Users may provide this function with information about the following:

  1. Individual history - Are the matrices to be built historical or ahistorical? If the former, then the state of the individual in occasion t-1 will be included in modeling.
  2. Modeling approach - Should the models be estimated as generalized linear models (GLMs) or generalized linear mixed models (GLMMs, the default)? While most function-based matrix models are estimated as the former, the latter approach can account for repeated observations of the same individual by including individual identity as a random factor. Mixed models also allow time and patch to be treated as random variables, which allows for broader and more theoretically sound inference.
  3. Suite of factors - Should both size and reproductive status be tested as causal factors? Should two-way interactions be included? Should only the y-intercept be estimated?
  4. Suite of vital rates - Which adult demographic parameters should be estimated? The defaults are (1) survival, (3) primary size, and (7) fecundity. Should (2) observation status, (4) secondary size, (5) tertiary size, or (6) reproductive status also be modeled?
  5. Juvenile vital rate estimation - Should juvenile parameters (8) through (14) also be modeled?
  6. Best-fit criterion - If a model with fewer parameters exists with \(\Delta AICc \leq 2.0\), then should this model be used as the best-fit model (the default), or should the model with the lowest AICc always be chosen?
  7. Size distribution - Should size be modeled as a continuous variable under a Gaussian or a gamma distribution, or as a count variable under either a Poisson or negative binomial distribution? If as a count variable, then should the distribution be zero-inflated to account for excess zeroes, zero-truncated to account for a lack of zeroes, or left unaltered?
  8. Fecundity distribution - Should fecundity be modeled as a continuous variable under a Gaussian or a gamma distribution, or as a count variable under either a Poisson or negative binomial distribution? If a count variable, then should the distribution be zero-inflated to account for excess zeroes, zero-truncated to account for a lack of zeroes, or left unaltered?
  9. Timing of fecundity - Function modelsearch() assumes that linear models of fecundity use a metric counted or measured in occasion t as the response. This applies well with most herbaceous plant datasets, where flowers or seeds produced in one year might be the fecundity response measured. However, users not wishing to follow this default behavior can use the fectime option to stipulate a fecundity metric measured in occasion t+1.
  10. Age - Is an age-classified (Leslie) or age-by-stage MPM the main goal?
  11. Individual covariates and environmental state variables - Should the effects of any individual or environmental covariates on vital rates also be tested? If so, should they be treated as fixed, quantitative terms or as random, categorical terms?
  12. Stage groups - If stage groups are noted in the stageframe, then should they also be used as fixed, independent categorical predictors in modeling?
  13. Censoring - Should data points marked as questionable be used in or excluded from modeling?
  14. Variable names - The names of all relevant variables in the dataset need to be specified. Note that the default behavior assumes variable names produced via the historicalize3() or verticalize3() functions, which produce standardized historically-formatted vertical (hfv format) datasets.
  15. Global model only - Should a full exhaustive model building and selection exercise be performed (the default), or should only the global model be built and parameterized?
  16. Accuracy - Should the accuracy of binomial models and pseudo-R2 of size and fecundity models be evaluated (the default)?

Once all decisions have been made and associated input terms have been provided, modelsearch() goes to work. The result is a lefkoMod object, which is a list in which the first 14 elements are the best-fit models developed for each vital rate. These are followed by an equivalent number of elements showing the full model tables developed and tested, followed by a table of parameter names, an element detailing the best-fit criterion used, and ending on quality control data showing the number of individuals and the number of unique transitions used in the estimation of each model, and the accuracy or pseudo-R2 of each model. Depending on user choices, linear modeling is handled through the lm() and glm() functions in the stats package, the lmer() and glmer() functions in lme4, the glmmTMB() function in glmmTMB, the glm.nb() function in package MASS, the vglm() function in package VGAM, or the zeroinfl() function in package pscl. Exhaustive model building proceeds through the dredge() function in package MuMIn (Bartoń 2014). Model selection is handled through assessment of AICc and the number of parameters estimated per model by default (see 6. Best-fit criterion above).

If modelsearch() is set for historical analysis (historical = TRUE, the default), then the decision of whether to develop a historical MPM can be made on the basis of whether any best-fit vital rate model includes size or reproductive status in occasion t-1. If at least one vital rate does, then a historical MPM is justified. Otherwise, it is not. Regardless, the output can be used to create a function-based MPM in the next step.

5.2 Assessing vital rates in Cypripedium candidum

Let’s now work on the Cypripedium candidum case study. We introduced the model and code for for our stageframe and supplements in Chapter 2. However, as it has likely been a while since we have viewed that chapter, let’s take a look again at the life history model that we will use (note that this is the same as figure 2.2).

Figure 5.1: 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 first. Then we will set up our stageframe and standardize the data into a historically formatted vertical (hfv) data frame.

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)

cypfb_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_fb, stagesize = "sizeadded", NAas0 = TRUE,
  age_offset = 4)

Now we will load the supplement tables that we will use for the ahistorical and historical MPMs. Note that these are slightly different from the previous chapter (chapter 4), because of the different life history model that we are using.

seeds_per_fruit <- 5000
sl_mult <- 0.7

cypsupp2_fb <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "SL", "D",
    "V1", "V2", "V3", "SD", "P1"), 
  stage2 = c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "SL", "SL", "rep",
    "rep"), 
  eststage3 = c(NA, NA, NA, NA, NA, NA, "D", "V1", "V2", "V3", NA, NA), 
  eststage2 = c(NA, NA, NA, NA, NA, NA, "D", "D", "D", "D", NA, NA), 
  givenrate = c(0.08, 0.1, 0.1, 0.1, 0.05, 0.05, NA, NA, NA, NA, NA, NA), 
  multiplier = c(NA, NA, NA, NA, NA, NA, sl_mult, sl_mult, sl_mult, sl_mult,
    0.5 * seeds_per_fruit, 0.5 * seeds_per_fruit),
  type = c("S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "R", "R"),
  stageframe = cypframe_fb, historical = FALSE)

cypsupp3_fb <- supplemental(stage3 = c("SD", "SD", "P1", "P1", "P2", "P3", "SL",
    "SL", "SL", "D", "V1", "V2", "V3", "D", "V1", "V2", "V3", "mat", "mat",
    "mat", "mat", "SD", "P1"), 
  stage2 = c("SD", "SD", "SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "SL", 
    "SL", "SL", "SL", "SL", "SL", "SL", "D", "V1", "V2", "V3", "rep", "rep"), 
  stage1 = c("SD", "rep", "SD", "rep", "SD", "P1", "P2", "P3", "SL", "P3", "P3",
    "P3", "P3", "SL", "SL", "SL", "SL", "SL", "SL", "SL", "SL", "mat", "mat"), 
  eststage3 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, "D", "V1", "V2", "V3", "D",
    "V1", "V2", "V3", "mat", "mat", "mat", "mat", NA, NA), 
  eststage2 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, "D", "D", "D", "D", "D", 
    "D", "D", "D", "D", "V1", "V2", "V3", NA, NA), 
  eststage1 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, "D", "D", "D", "D", "D", 
    "D", "D", "D", "V1", "V1", "V1", "V1", NA, NA), 
  givenrate = c(0.08, 0.08, 0.1, 0.1, 0.1, 0.1, 0.1, 0.05, 0.05, NA, NA, NA, NA, 
    NA, NA, NA, NA, NA, NA, NA, NA, NA, NA),
  multiplier = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, sl_mult, sl_mult,
    sl_mult, sl_mult, sl_mult, sl_mult, sl_mult, sl_mult, 1, 1, 1, 1,
    0.5 * seeds_per_fruit, 0.5 * seeds_per_fruit),
  type = c("S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S",
    "S", "S", "S", "S", "S", "S", "S", "R", "R"), 
  type_t12 = c("S", "F", "S", "F", "S", "S", "S", "S", "S", "S", "S", "S", "S",
    "S", "S", "S", "S", "S", "S", "S", "S", "S", "S"),
  stageframe = cypframe_fb, historical = TRUE)

Before we can develop our vital rate models, we need to determine the proper distributions for size and fecundity. As a reminder, size in the Cypripedium dataset is given as a series of numbers of sprouts, and our proxy for size 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, or count variables. 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. Package lefko3 includes a function that can help in determining which distributions to use: sf_distrib(). Fecundity, in contrast, is a count variable that includes zeros and so we should test whether we need a zero-inflated model.

sf_distrib(cypfb_v1, sizea = c("size3added", "size2added"), obs3 = "obsstatus3",
  fec = c("feca3", "feca2"), repst = c("repstatus3", "repstatus2"),
  zisizea = FALSE)
> Mean sizea is 3.653
> 
> The variance in sizea is 13.41
> 
> The probability of this dispersion level by chance assuming that
> the true mean sizea = variance in sizea,
> and an alternative hypothesis of overdispersion, is 3.721e-138
> 
> Primary size is significantly overdispersed.
> 
> Mean fec is 0.8036
> 
> The variance in fec is 1.601
> 
> The probability of this dispersion level by chance assuming that
> the true mean fec = variance in fec,
> and an alternative hypothesis of overdispersion, is 0.07088
> 
> Dispersion level in fecundity matches expectation.
> 
> 
> Mean lambda in fec is 0.4477
> The actual number of 0s in fec is 65
> The expected number of 0s in fec under the null hypothesis is 50.15
> The probability of this deviation in 0s from expectation by chance is 1.742e-06
> 
> Fecundity is significantly zero-inflated.

The results show that size is significantly overdispersed. This means that we cannot use the Poisson distribution, and should instead use the negative binomial distribution. Fecundity is not significantly overdispersed, so we can use the Poisson distribution. However, fecundity is also significantly zero-inflated, so this will require a zero-inflated Poisson distribution.

5.3 Setting up function modelsearch()

Next, we will create a full suite of vital rate models for the Cypripedium candidum dataset using function modelsearch(). Before proceeding, we need to decide on the linear model building strategy, the correct vital rates to model, the proper statistical distributions for estimated vital rates, the proper parameterization for each vital rate, and the strategy for determination of the best-fit models.

First, we must determine the model building strategy. In most cases, the best procedure will be through mixed linear models in which monitoring occasion and individual identity are random terms. We will set monitoring occasion as random because we wish to make inferences for the population as a whole and do not wish to restrict ourselves to inference only for the years monitored (i.e. our distribution of monitoring occasions sampled is itself a sample of the population in time). We will set individual identity as random because many or most of the individuals that we have sampled to produce our dataset yield multiple observation data points across time. Thus, we will set approach = "mixed". To make sure that time and individual identity are treated as random, we will set the proper variable names for indiv and year, corresponding to individual identity (individ by default), and to occasion t (year2 by default). The year.as.random option is set to random by default. Setting year.as.random to FALSE would make time a fixed categorical variable in all cases. The patch.as.random option functions similarly with respect to patches / subpopulations.

The mixed modeling approach is usually preferable to the GLM approach because of its elimination of pseudoreplication. However, a mixed modeling strategy results in lower statistical power and a greater time used to estimate models (more correctly, it yields truer statistical power while the GLM approach inflates Type I error). Users of package lefko3 wishing to use a standard generalized linear modeling strategy should set approach = "glm". In this case, individual identity is not used, time is a fixed categorical factor (as is patch, if used), and all observed transitions are treated as independent.

Next, we must determine which vital rates to model. The default settings for modelsearch estimate 1) survival probability, 3) primary size distribution, and 7) fecundity, which are the minimum three vital rates required for a full MPM. Observation probability (option obs in vitalrates) should only be included when a life history stage or size exists that cannot be observed. For example, in the case of a plant with vegetative dormancy, the observation probability can be thought of as the sprouting probability, which is a biologically meaningful vital rate (Shefferson et al. 2001). Further, reproduction status (option repst in vitalrates) should only be modeled if size classification needs to be stratified by the ability to reproduce, as when some individuals with no fecundity occur within stages that also include individuals producing offspring. Since Cypripedium candidum is capable of long bouts of vegetative dormancy, since we wish to stratify the population into reproductive and non-reproductive stages of the same size classes, and since we have no data derived from juvenile individuals, we will set vitalrates = c("surv", "obs", "size", "repst", "fec") and we will not set juvestimate and juvsize.

Third, we need to set the proper statistical distribution for each parameter. Survival probability, observation probability, and reproductive status are all modeled as binomial variables, and this cannot be changed. In the case of this population of Cypripedium candidum, size was measured as the number of stems and so is a count variable. Likewise, fecundity is estimated as the number of fruits produced per plant, and so is also a count variable. We have already performed tests for overdispersion and zero-inflation, and we are also aware that size in observed stages cannot be zero, requiring zero-truncation in that parameter. So we will set size to the zero-truncated negative binomial distribution, and fecundity to the zero-inflated Poisson distribution.

Fourth, we need the proper model parameterizations for each vital rate, using the suite option. The default, suite = "main", under the mixed model setting (approach = "mixed") starts with the estimation of global models that include size and reproductive status in occasions t and t-1 as fixed factors, with individual identity and time in occasion t (year t) set as random categorical terms. Setting suite = "full" will yield global models that also include all two-way interactions among fixed terms ("full" is the only setting with interaction terms). If the population is not stratified by reproductive status, then suite = "size" will eliminate reproductive status terms and use all others in the global model. If size is not important, then suite = "rep" will eliminate size but keep reproductive status and all other terms. Finally, suite = "cons" will result in a global model in which neither reproductive status nor size are considered. Other terms can be specified, including individual covariates and age.

Finally, we should determine the proper strategy for the determination of the best-fit model. Model building proceeds through the dredge() function in package MuMIn (Bartoń 2014), and each model has an associated AICc value. The default setting in lefko3 (bestfit = "AICc&k") will compare all models within 2.0 AICc units of the model with \(\Delta AICc = 0\), and choose the one with the lowest degrees of freedom. This approach is generally better than the alternative, which simply uses the model with \(\Delta AICc = 0\) (bestfit = "AICc"), as all models within 2.0 AICc units of that model are equally parsimonious and so fewer degrees of freedom result from fewer parameters estimated (Burnham & Anderson 2002).

In the model building exercise below, we will use the suite = "main" option to run all main effects only. Normally we would set to suite = "full", but running all effects including their two-way interactions will likely tie up our computers for a little too long for this tutorial. While this runs, feel free to make yourself a coffee or a cappuccino, or perhaps something stronger (no judgement here…).

cypmodels3p <- modelsearch(cypfb_v1, historical = TRUE, approach = "mixed", 
  vitalrates = c("surv", "obs", "size", "repst", "fec"), patch = "patchid",
  sizedist = "negbin", size.trunc = TRUE, fecdist = "poisson", fec.zero = TRUE,
  suite = "main", size = c("size3added", "size2added", "size1added"),
  quiet = TRUE)

Once done, we will summarize the output with the summary() function.

summary(cypmodels3p)
> 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 ~ size2added + (1 | year2) + (1 | patchid) + (1 | individ)
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
> 130.1321 148.9737 -60.0660 120.1321      315 
> Random effects:
>  Groups  Name        Std.Dev. 
>  individ (Intercept) 1.199e+00
>  year2   (Intercept) 5.161e-05
>  patchid (Intercept) 1.113e-05
> Number of obs: 320, groups:  individ, 74; year2, 5; patchid, 3
> Fixed Effects:
> (Intercept)   size2added  
>      2.0356       0.6343  
> 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 | patchid) + (1 |  
>     individ)
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
> 120.2567 138.8254 -55.1284 110.2567      298 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 0.0000  
>  year2   (Intercept) 0.8776  
>  patchid (Intercept) 0.0000  
> Number of obs: 303, groups:  individ, 70; year2, 5; patchid, 3
> 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 ~ size2added + (1 | year2) + (1 | patchid) + (1 |  
>     individ)
> Data: subdata
>       AIC       BIC    logLik  df.resid 
> 1009.3168 1031.2946 -498.6584       282 
> Random-effects (co)variances:
> 
> Conditional model:
>  Groups  Name        Std.Dev.
>  year2   (Intercept) 0.1215  
>  patchid (Intercept) 0.2052  
>  individ (Intercept) 0.9497  
> 
> Number of obs: 288 / Conditional model: year2, 5; patchid, 3; individ, 70
> 
> Dispersion parameter for truncated_nbinom2 family (): 1.23e+09 
> 
> Fixed Effects:
> 
> Conditional model:
> (Intercept)   size2added  
>     0.54038      0.02191  
> 
> 
> 
> 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 ~ repstatus2 + size2added + (1 | year2) + (1 | patchid) +  
>     (1 | individ)
>    Data: subdata
>       AIC       BIC    logLik  deviance  df.resid 
>  333.4037  355.3815 -160.7019  321.4037       282 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 0.1776  
>  year2   (Intercept) 0.6636  
>  patchid (Intercept) 0.3501  
> Number of obs: 288, groups:  individ, 70; year2, 5; patchid, 3
> Fixed Effects:
> (Intercept)   repstatus2   size2added  
>     -1.3836       1.5543       0.1788  
> 
> 
> 
> Fecundity model:
> Formula:          
> feca2 ~ size1added + size2added + (1 | year2) + (1 | patchid) +  
>     (1 | individ)
> Zero inflation:         
> ~size2added + (1 | year2) + (1 | patchid) + (1 | individ)
> Data: subdata
>       AIC       BIC    logLik  df.resid 
>  252.3835  282.8610 -115.1917       107 
> Random-effects (co)variances:
> 
> Conditional model:
>  Groups  Name        Std.Dev. 
>  year2   (Intercept) 6.043e-01
>  patchid (Intercept) 2.130e-01
>  individ (Intercept) 5.213e-10
> 
> Zero-inflation model:
>  Groups  Name        Std.Dev. 
>  year2   (Intercept) 4.136e-12
>  patchid (Intercept) 1.460e-13
>  individ (Intercept) 2.951e-04
> 
> Number of obs: 118 / Conditional model: year2, 5; patchid, 3; individ, 51 / Zero-inflation model: year2, 5; patchid, 3; individ, 51
> 
> Fixed Effects:
> 
> Conditional model:
> (Intercept)   size1added   size2added  
>    -0.51758     -0.03543      0.08305  
> 
> Zero-inflation model:
> (Intercept)   size2added  
>       3.857       -1.575  
> 
> 
> 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: 16
> 
> Number of models in observation table: 16
> 
> Number of models in size table: 16
> 
> Number of models in secondary size table: 1
> 
> Number of models in tertiary size table: 1
> 
> Number of models in reproduction status table: 16
> 
> Number of models in fecundity table: 240
> 
> Number of models in juvenile survival table: 1
> 
> Number of models in juvenile observation table: 1
> 
> Number of models in juvenile size table: 1
> 
> Number of models in juvenile secondary size table: 1
> 
> Number of models in juvenile tertiary size table: 1
> 
> Number of models in juvenile reproduction table: 1
> 
> Number of models in juvenile maturity table: 1
> 
> 
> 
> 
> 
> General model parameter names (column 1), and 
> specific names used in these models (column 2): 
>                       parameter_names mainparams
> 1                              time t      year2
> 2                          individual    individ
> 3                               patch      patch
> 4                   alive in time t+1      surv3
> 5                observed in time t+1       obs3
> 6                   sizea in time t+1      size3
> 7                   sizeb in time t+1     sizeb3
> 8                   sizec in time t+1     sizec3
> 9     reproductive status in time t+1     repst3
> 10              fecundity in time t+1       fec3
> 11                fecundity in time t       fec2
> 12                    sizea in time t      size2
> 13                  sizea in time t-1      size1
> 14                    sizeb in time t     sizeb2
> 15                  sizeb in time t-1     sizeb1
> 16                    sizec in time t     sizec2
> 17                  sizec in time t-1     sizec1
> 18      reproductive status in time t     repst2
> 19    reproductive status in time t-1     repst1
> 20        maturity status in time t+1     matst3
> 21          maturity status in time t     matst2
> 22                      age in time t        age
> 23                  density in time t    density
> 24   individual covariate a in time t   indcova2
> 25 individual covariate a in time t-1   indcova1
> 26   individual covariate b in time t   indcovb2
> 27 individual covariate b in time t-1   indcovb1
> 28   individual covariate c in time t   indcovc2
> 29 individual covariate c in time t-1   indcovc1
> 30              stage group in time t     group2
> 31            stage group in time t-1     group1
> 
> 
> 
> 
> 
> Quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Survival accuracy is 0.947.
> Observation estimated with 70 individuals and 303 individual transitions.
> Observation accuracy is 0.95.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Primary size pseudo R-squared is NA.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Reproductive status accuracy is 0.74.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Fecundity pseudo R-squared is 0.344.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity probability not estimated.

The above code may take some time to run, depending on how powerful a computer it is run on. It is sufficiently complicated that R might crash if it is run while the user does other things in R or R Studio, so please leave your R or R Studio session alone while it runs. As the function runs, it will likely produce some text related to the modeling. Some of the text is meant to provide guideposts to where the modeling stands. Since the vital rates run in order, the user can get a handle on which vital rate best-fit models are finished and which are still waiting to be developed. Other text relates to notes, warnings, and other messages from the modeling functions used, which are found in other CRAN-based packages. We have set quiet = TRUE to prevent this text from getting very long.

The final object is a class lefkoMod object. This is an S3 object, meaning that it is a list, and it is composed of 31 elements. The first fourteen elements are the best-fit models, and are in the formats of the functions used to develop them. The next fourteen elements are the model tables, showing all of the models tested for each vital rate in order of \(\Delta AICc\) (\(\Delta AICc\) is given as the difference between the AICc values of the current model and the model with the lowest AICc). The final three elements are quality control data used in summary() calls.

The summary of this lefkoMod object is quite informative. We start off by seeing the best-fit model chosen for each vital rate. In the case of vital rates that are not estimated, they are set to constants, either 1 or 0. Below that, we see how many models were built and compared for each vital rate, as well as data on the accuracy of binomial models and the pseudo-R2 values of size and fecundity models. The next section is something of a crib sheet, showing us the names of different parameters used in the modeling. Finally, we see a section on quality control. This last section includes information on the numbers of individuals as well as individual transitions used to parameterize each model, and also shows their accuracy or pseudo-R2 values. We see very high accuracy in survival and observation status, but poorer accuracy in reproduction status.

Looking over the models, we can also see from the summary that historical size is included in the best-fit model of fecundity. This suggests that a historical MPM is a more parsimonious choice of MPM than an ahistorical MPM. We also see that the survival and observation status models are very accurate, but that the reproductive status model is less so. Regardless, in order to compare MPMs for educational purposes, we will also create an ahistorical model set. We will set quiet = TRUE to see less text output as the function works, but you may remove that option to view the guidepost text. Note that a vital rate model set that includes historical terms should never be used to make an ahistorical MPM.

cypmodels2p <- modelsearch(cypfb_v1, historical = FALSE, approach = "mixed", 
  vitalrates = c("surv", "obs", "size", "repst", "fec"), patch = "patchid",
  sizedist = "negbin", size.trunc = TRUE, fecdist = "poisson", fec.zero = TRUE,
  suite = "main", size = c("size3added", "size2added", "size1added"),
  quiet = TRUE)

Let’s see a summary of this lefkoMod object.

summary(cypmodels2p)
> 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 ~ size2added + (1 | year2) + (1 | patchid) + (1 | individ)
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
> 130.1321 148.9737 -60.0660 120.1321      315 
> Random effects:
>  Groups  Name        Std.Dev. 
>  individ (Intercept) 1.199e+00
>  year2   (Intercept) 5.161e-05
>  patchid (Intercept) 1.113e-05
> Number of obs: 320, groups:  individ, 74; year2, 5; patchid, 3
> Fixed Effects:
> (Intercept)   size2added  
>      2.0356       0.6343  
> 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 | patchid) + (1 |  
>     individ)
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
> 120.2567 138.8254 -55.1284 110.2567      298 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 0.0000  
>  year2   (Intercept) 0.8776  
>  patchid (Intercept) 0.0000  
> Number of obs: 303, groups:  individ, 70; year2, 5; patchid, 3
> 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 ~ size2added + (1 | year2) + (1 | patchid) + (1 |  
>     individ)
> Data: subdata
>       AIC       BIC    logLik  df.resid 
> 1009.3168 1031.2946 -498.6584       282 
> Random-effects (co)variances:
> 
> Conditional model:
>  Groups  Name        Std.Dev.
>  year2   (Intercept) 0.1215  
>  patchid (Intercept) 0.2052  
>  individ (Intercept) 0.9497  
> 
> Number of obs: 288 / Conditional model: year2, 5; patchid, 3; individ, 70
> 
> Dispersion parameter for truncated_nbinom2 family (): 1.23e+09 
> 
> Fixed Effects:
> 
> Conditional model:
> (Intercept)   size2added  
>     0.54038      0.02191  
> 
> 
> 
> 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 ~ repstatus2 + size2added + (1 | year2) + (1 | patchid) +  
>     (1 | individ)
>    Data: subdata
>       AIC       BIC    logLik  deviance  df.resid 
>  333.4037  355.3815 -160.7019  321.4037       282 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 0.1776  
>  year2   (Intercept) 0.6636  
>  patchid (Intercept) 0.3501  
> Number of obs: 288, groups:  individ, 70; year2, 5; patchid, 3
> Fixed Effects:
> (Intercept)   repstatus2   size2added  
>     -1.3836       1.5543       0.1788  
> 
> 
> 
> Fecundity model:
> Formula:          feca2 ~ (1 | year2) + (1 | patchid) + (1 | individ)
> Zero inflation:         
> ~size2added + (1 | year2) + (1 | patchid) + (1 | individ)
> Data: subdata
>       AIC       BIC    logLik  df.resid 
>  256.6460  281.5822 -119.3230       109 
> Random-effects (co)variances:
> 
> Conditional model:
>  Groups  Name        Std.Dev.
>  year2   (Intercept) 0.5908  
>  patchid (Intercept) 0.1970  
>  individ (Intercept) 0.3863  
> 
> Zero-inflation model:
>  Groups  Name        Std.Dev. 
>  year2   (Intercept) 2.458e-08
>  patchid (Intercept) 3.189e-07
>  individ (Intercept) 2.445e-04
> 
> Number of obs: 118 / Conditional model: year2, 5; patchid, 3; individ, 51 / Zero-inflation model: year2, 5; patchid, 3; individ, 51
> 
> Fixed Effects:
> 
> Conditional model:
> (Intercept)  
>     -0.2146  
> 
> Zero-inflation model:
> (Intercept)   size2added  
>       4.081       -1.597  
> 
> 
> 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: 4
> 
> Number of models in observation table: 4
> 
> Number of models in size table: 4
> 
> Number of models in secondary size table: 1
> 
> Number of models in tertiary size table: 1
> 
> Number of models in reproduction status table: 4
> 
> Number of models in fecundity table: 13
> 
> Number of models in juvenile survival table: 1
> 
> Number of models in juvenile observation table: 1
> 
> Number of models in juvenile size table: 1
> 
> Number of models in juvenile secondary size table: 1
> 
> Number of models in juvenile tertiary size table: 1
> 
> Number of models in juvenile reproduction table: 1
> 
> Number of models in juvenile maturity table: 1
> 
> 
> 
> 
> 
> General model parameter names (column 1), and 
> specific names used in these models (column 2): 
>                       parameter_names mainparams
> 1                              time t      year2
> 2                          individual    individ
> 3                               patch      patch
> 4                   alive in time t+1      surv3
> 5                observed in time t+1       obs3
> 6                   sizea in time t+1      size3
> 7                   sizeb in time t+1     sizeb3
> 8                   sizec in time t+1     sizec3
> 9     reproductive status in time t+1     repst3
> 10              fecundity in time t+1       fec3
> 11                fecundity in time t       fec2
> 12                    sizea in time t      size2
> 13                  sizea in time t-1      size1
> 14                    sizeb in time t     sizeb2
> 15                  sizeb in time t-1     sizeb1
> 16                    sizec in time t     sizec2
> 17                  sizec in time t-1     sizec1
> 18      reproductive status in time t     repst2
> 19    reproductive status in time t-1     repst1
> 20        maturity status in time t+1     matst3
> 21          maturity status in time t     matst2
> 22                      age in time t        age
> 23                  density in time t    density
> 24   individual covariate a in time t   indcova2
> 25 individual covariate a in time t-1   indcova1
> 26   individual covariate b in time t   indcovb2
> 27 individual covariate b in time t-1   indcovb1
> 28   individual covariate c in time t   indcovc2
> 29 individual covariate c in time t-1   indcovc1
> 30              stage group in time t     group2
> 31            stage group in time t-1     group1
> 
> 
> 
> 
> 
> Quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Survival accuracy is 0.947.
> Observation estimated with 70 individuals and 303 individual transitions.
> Observation accuracy is 0.95.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Primary size pseudo R-squared is NA.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Reproductive status accuracy is 0.74.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Fecundity pseudo R-squared is 0.362.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity probability not estimated.

The ahistorical set of models is rather similar to the historical case, since historical terms mostly dropped out of the historical analysis. However, the fecundity model is obviously different, since object cypmodels3p includes history only in fecundity, and historical terms were not tested for cypmodels2p. Further, since historical terms were not included in the global models in cypmodels2p, the number of models tested per vital rate is much lower.

Function modelsearch() has many options that can lead to very complicated global models and searches for vital rate models. For example, up to three size metrics may be used for classification, and these may assume different distributions. Also, up to three individual/environmental covariates coded in the hfv dataset may also be included in models, and these covariates may be quantitative (used as fixed factors), or random and categorical (options indcova, indcovb, and indcovc, with logical settings random.indcova, random.indcovb, and random.indcovc). If stage groups are used and were coded into the stageframe, then stage group can be included in models as a categorical fixed factor (logical setting test.group). Model tables may be eliminated from the lefkoMod output using show.model.tables = FALSE, and model selection can be eliminated with the global model alone being estimated with the option global.only = TRUE.

5.4 Size and fecundity distributions

Before we create MPMs, let’s discuss response distributions a bit further. The probabilities of survival, observation, and reproductive status are automatically set to the binomial distribution in function modelsearch(), and this cannot be altered. However, the probability of size transition and the fecundity rate can be set to the Gaussian, gamma, Poisson, or negative binomial distributions, with zero-inflated and zero-truncated versions of the Poisson and negative binomial also available. If size or fecundity rate is a continuous variable (i.e., not an integer or count variable), then it should be set to the Gaussian distribution if reasonably symmetric. If there is a strong right skew, then a gamma distribution may be warranted. This can be explored by using the plot() function in conjunction with density(). For example, here we look at the primary size distribution used in the Cypripedium candidum case (figure 5.2). Please note that there is no established test to determine whether a Gaussian distribution is better than the gamma, so this decision should be made at the user’s discretion.

plot(density(cypfb_v1$size2added))
Density distribution of size in Cypripedium candidum

Figure 5.2: Density distribution of size in Cypripedium candidum

If size or fecundity is a count variable, then it should be set to the Poisson distribution if the mean equals the variance, and this can be tested with the sf_distrib() function. The negative binomial distribution is provided in cases where the assumption that the mean equals the variance is clearly broken. We do not encourage the use of the negative binomial except in such cases, as the extra parameters estimated for the negative binomial distribution reduce the power of the modeling exercises conducted.

The Poisson and the negative binomial distributions both predict specific numbers of zeros in the response variable. If excess zeros occur within the dataset even after including the observation status and reproductive status as vital rates to absorb zeros, then a zero-inflated Poisson or negative binomial distribution may be used. These modeling approaches work by parameterizing a binomial model, typically with a logit link, to predict zero responses. The Poisson or negative binomial is then used to predict non-zero responses. In effect, a zero-inflated model is actually two models in which zeros are assumed to be predicted under potentially different processes than the remaining counts. Users should be aware that, because an extra model is built to cover zeros, zero-inflated models are much more complex and can include many more parameters than their non-inflated counterparts. The principle of parsimony suggests that they should only be used when there are significantly more zeros than expected. This may be tested with sf_distrib().

Cases may arise in which zeros do not occur in either size or fecundity. For these situations, we provide zero-truncated distributions. This may occur in size if all cases of size = 0 are absorbed by observation status, leaving only positive integers for the size of observed individuals. For example, if an unobservable stage such as vegetative dormancy occurs and absorbs all cases of size = 0, then a zero-truncated Poisson or negative binomial distribution will be more appropriate than the equivalent distribution without zero-truncation. It can also occur if all cases of fecundity = 0 are absorbed by reproductive status. Such distributions only involve the estimation of single, conditional models, and so are simpler than zero-inflated models.

We have assumed that users will only use a single size variable in classification. However, cases may arise in which stages are classified under two or three size variables, and these two or three size variables are also used in vital rate modeling. These size variables can be treated as independent for the purposes of modeling, meaning that they can also have completely different distributions. Function sf_distrib() can assess up to three size variables plus one fecundity variable in a single run.

Let’s now build some function-based MPMs.

5.5 Using lefkoMod objects to create function-based MPMs

MPM creation can be accomplished with eight different functions:

  1. rlefko2() - Creates raw ahistorical MPMs given a dataset, a stageframe, and either a supplement table or a reproductive matrix and an overwrite table.
  2. rlefko3() - Creates raw historical MPMs given a dataset, a stageframe, and either a supplement table or a reproductive matrix and an overwrite table.
  3. arlefko2() - Creates raw ahistorical age x stage MPMs given a dataset, a stageframe, and either a supplement table or a reproductive matrix and an overwrite table.
  4. rleslie() - Creates raw age-based (Leslie) MPMs given a dataset.
  5. flefko2() - Creates function-based ahistorical MPMs given a dataset, a set of models, a stageframe, and either a supplement table or a reproductive matrix and an overwrite table.
  6. flefko3() - Creates function-based historical MPMs given a dataset, a set of models, a stageframe, and either a supplement table or a reproductive matrix and an overwrite table.
  7. aflefko2() - Creates function-based ahistorical age x stage MPMs given a dataset, a set of models, a stageframe, and either a supplement table or a reproductive matrix and an overwrite table.
  8. fleslie() - Creates function-based age-based (Leslie) MPMs given a dataset and a set of models.

These functions incorporate binary kernels developed to handle the estimation of matrix elements quickly and efficiently. A single run of flefko3(), for example, should be able to yield all annual matrices for all patches for the Cypripedium candidum dataset provided with lefko3 in under a minute on most machines (14s or so on the author’s 2019 MacBook Pro with 2.3 GHz 8-Core Intel Core i9). Parallel computing should not be necessary even with the slowest of current machines, provided that the machine is current enough to handle at least R 3.6.3.

The output for each of these functions is a lefkoMat object, which is an S3 object and is described in detail in Chapter 4. These lists include elements A, U, and F, which are themselves lists of complete projection matrices, survival-transition matrices, and fecundity matrices, respectively. For example, code such as matobject$A[[1]] would access the first complete projection matrix in a lefkoMat object named matobject. They are followed by elements referred to as ahstages, hstages, and agestages, which provide data frames describing the stages (technically showing the edited stageframe), the historical stage-pairs, and the age-stage pairs shown in the order in which they occur within the matrices, respectively. The labels element is a data frame giving a description of each matrix in the order in which it occurs within the A, U, and F elements, including the population, patch, and monitoring time designations. The final elements are quality control outputs that vary in content depending on whether the output matrices are raw or function-based, but nonetheless always include at least the numbers of individuals and individual-transitions used for estimation.

Let’s create function-based MPMs with the lefkoMod objects that we created in this chapter. We will start with the ahistorical function-based MPM. We will also look at a summary of the resulting lefkoMat object.

cypmatrix2fp <- flefko2(stageframe = cypframe_fb, supplement = cypsupp2_fb, 
  modelsuite = cypmodels2p, data = cypfb_v1)

summary(cypmatrix2fp)
> 
> This ahistorical lefkoMat object contains 15 matrices.
> 
> Each matrix is square with 54 rows and columns, and a total of 2916 elements.
> A total of 36165 survival transitions were estimated, with 2411 per matrix.
> A total of 720 fecundity transitions were estimated, with 48 per matrix.
> This lefkoMat object covers 1 population, 3 patches, and 5 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Observation estimated with 70 individuals and 303 individual transitions.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
> 
> 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.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050
> 1st Qu. 0.991 0.991 0.991 0.991 0.991 0.991 0.991 0.991 0.991 0.991 0.991 0.991
> Median  0.999 1.000 1.000 1.000 0.999 0.998 1.000 0.999 0.999 0.999 0.997 0.999
> Mean    0.916 0.918 0.917 0.917 0.914 0.915 0.917 0.917 0.916 0.914 0.913 0.916
> 3rd Qu. 1.000 1.000 1.000 1.000 1.000 0.999 1.000 1.000 1.000 0.999 0.999 1.000
> Max.    1.000 1.000 1.000 1.000 1.000 0.999 1.000 1.000 1.000 1.000 0.999 1.000
>         [,13] [,14] [,15]
> Min.    0.050 0.050 0.050
> 1st Qu. 0.991 0.991 0.991
> Median  0.999 0.998 0.997
> Mean    0.916 0.915 0.912
> 3rd Qu. 0.999 0.999 0.999
> Max.    1.000 0.999 0.999

Here we have fifteen matrices. Since there are six years of data, there are a total of five pairs of years with which to estimate transitions. There are also three patches, yielding fifteen total matrices. We also have \(54^{2} = 2916\) elements per matrix, of which 2,459 elements are estimated (the rest are zeros). Let’s look at the first function-based matrix.

writeLines("\nFirst matrix in function-based ahMPM:")
> 
> First matrix in function-based ahMPM:
print(cypmatrix2fp$A[[1]], digits = 3)
>       [,1] [,2] [,3] [,4]   [,5]     [,6]     [,7]     [,8]     [,9]    [,10]
>  [1,] 0.08  0.0  0.0 0.00 0.0000 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [2,] 0.10  0.0  0.0 0.00 0.0000 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [3,] 0.00  0.1  0.0 0.00 0.0000 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [4,] 0.00  0.0  0.1 0.00 0.0000 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [5,] 0.00  0.0  0.0 0.05 0.0500 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [6,] 0.00  0.0  0.0 0.00 0.0329 4.70e-02 3.69e-02 2.81e-02 2.10e-02 1.56e-02
>  [7,] 0.00  0.0  0.0 0.00 0.1660 2.37e-01 2.40e-01 2.34e-01 2.23e-01 2.09e-01
>  [8,] 0.00  0.0  0.0 0.00 0.1044 1.49e-01 1.52e-01 1.50e-01 1.44e-01 1.36e-01
>  [9,] 0.00  0.0  0.0 0.00 0.0656 9.37e-02 9.63e-02 9.55e-02 9.25e-02 8.80e-02
> [10,] 0.00  0.0  0.0 0.00 0.0000 5.89e-02 6.11e-02 6.10e-02 5.96e-02 5.71e-02
> [11,] 0.00  0.0  0.0 0.00 0.0000 3.71e-02 3.87e-02 3.90e-02 3.84e-02 3.71e-02
> [12,] 0.00  0.0  0.0 0.00 0.0000 2.33e-02 2.45e-02 2.49e-02 2.47e-02 2.41e-02
> [13,] 0.00  0.0  0.0 0.00 0.0000 1.47e-02 1.56e-02 1.59e-02 1.59e-02 1.56e-02
> [14,] 0.00  0.0  0.0 0.00 0.0000 9.22e-03 9.86e-03 1.02e-02 1.02e-02 1.01e-02
> [15,] 0.00  0.0  0.0 0.00 0.0000 5.80e-03 6.25e-03 6.50e-03 6.60e-03 6.58e-03
> [16,] 0.00  0.0  0.0 0.00 0.0000 3.64e-03 3.96e-03 4.16e-03 4.25e-03 4.27e-03
> [17,] 0.00  0.0  0.0 0.00 0.0000 2.29e-03 2.51e-03 2.66e-03 2.74e-03 2.77e-03
> [18,] 0.00  0.0  0.0 0.00 0.0000 1.44e-03 1.59e-03 1.70e-03 1.76e-03 1.80e-03
>          [,11]    [,12]    [,13]    [,14]    [,15]    [,16]    [,17]    [,18]
>  [1,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [2,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [3,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [4,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [5,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [6,] 1.15e-02 8.47e-03 6.21e-03 4.55e-03 3.33e-03 2.44e-03 1.78e-03 1.30e-03
>  [7,] 1.93e-01 1.77e-01 1.60e-01 1.43e-01 1.27e-01 1.12e-01 9.81e-02 8.52e-02
>  [8,] 1.26e-01 1.16e-01 1.06e-01 9.58e-02 8.57e-02 7.61e-02 6.70e-02 5.86e-02
>  [9,] 8.26e-02 7.67e-02 7.04e-02 6.41e-02 5.78e-02 5.16e-02 4.58e-02 4.03e-02
> [10,] 5.40e-02 5.05e-02 4.67e-02 4.28e-02 3.89e-02 3.50e-02 3.13e-02 2.77e-02
> [11,] 3.53e-02 3.33e-02 3.10e-02 2.86e-02 2.62e-02 2.38e-02 2.14e-02 1.91e-02
> [12,] 2.31e-02 2.19e-02 2.06e-02 1.92e-02 1.76e-02 1.61e-02 1.46e-02 1.31e-02
> [13,] 1.51e-02 1.45e-02 1.37e-02 1.28e-02 1.19e-02 1.09e-02 9.97e-03 9.03e-03
> [14,] 9.89e-03 9.53e-03 9.08e-03 8.57e-03 8.01e-03 7.42e-03 6.81e-03 6.21e-03
> [15,] 6.47e-03 6.28e-03 6.03e-03 5.73e-03 5.39e-03 5.03e-03 4.65e-03 4.27e-03
> [16,] 4.23e-03 4.14e-03 4.00e-03 3.83e-03 3.63e-03 3.41e-03 3.18e-03 2.94e-03
> [17,] 2.77e-03 2.73e-03 2.66e-03 2.56e-03 2.45e-03 2.32e-03 2.17e-03 2.02e-03
> [18,] 1.81e-03 1.80e-03 1.76e-03 1.71e-03 1.65e-03 1.57e-03 1.48e-03 1.39e-03
>          [,19]    [,20]    [,21]    [,22]    [,23]    [,24]    [,25]    [,26]
>  [1,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [2,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [3,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [4,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [5,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [6,] 9.54e-04 6.97e-04 5.10e-04 3.73e-04 2.72e-04 1.99e-04 1.46e-04 1.06e-04
>  [7,] 7.35e-02 6.30e-02 5.38e-02 4.56e-02 3.86e-02 3.25e-02 2.72e-02 2.28e-02
>  [8,] 5.09e-02 4.39e-02 3.77e-02 3.22e-02 2.74e-02 2.32e-02 1.96e-02 1.65e-02
>  [9,] 3.53e-02 3.06e-02 2.65e-02 2.28e-02 1.95e-02 1.66e-02 1.41e-02 1.19e-02
> [10,] 2.44e-02 2.14e-02 1.86e-02 1.61e-02 1.38e-02 1.19e-02 1.02e-02 8.65e-03
> [11,] 1.69e-02 1.49e-02 1.30e-02 1.14e-02 9.85e-03 8.50e-03 7.31e-03 6.26e-03
> [12,] 1.17e-02 1.04e-02 9.15e-03 8.02e-03 7.00e-03 6.08e-03 5.26e-03 4.54e-03
> [13,] 8.11e-03 7.24e-03 6.42e-03 5.67e-03 4.98e-03 4.35e-03 3.79e-03 3.29e-03
> [14,] 5.62e-03 5.05e-03 4.51e-03 4.00e-03 3.54e-03 3.11e-03 2.73e-03 2.38e-03
> [15,] 3.89e-03 3.52e-03 3.16e-03 2.83e-03 2.51e-03 2.23e-03 1.96e-03 1.72e-03
> [16,] 2.69e-03 2.45e-03 2.22e-03 2.00e-03 1.79e-03 1.59e-03 1.41e-03 1.25e-03
> [17,] 1.87e-03 1.71e-03 1.56e-03 1.41e-03 1.27e-03 1.14e-03 1.02e-03 9.04e-04
> [18,] 1.29e-03 1.19e-03 1.09e-03 9.96e-04 9.03e-04 8.15e-04 7.32e-04 6.54e-04
>          [,27]    [,28]    [,29]    [,30]    [,31]    [,32]    [,33]    [,34]
>  [1,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 2.31e+03 1.77e+03 8.24e+02 2.26e+02
>  [2,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 2.31e+03 1.77e+03 8.24e+02 2.26e+02
>  [3,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [4,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [5,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [6,] 7.78e-05 5.69e-05 4.16e-05 3.04e-05 3.69e-02 2.81e-02 2.10e-02 1.56e-02
>  [7,] 1.90e-02 1.58e-02 1.31e-02 1.09e-02 1.19e-01 1.09e-01 9.72e-02 8.53e-02
>  [8,] 1.38e-02 1.16e-02 9.69e-03 8.09e-03 7.55e-02 6.96e-02 6.26e-02 5.53e-02
>  [9,] 1.01e-02 8.50e-03 7.14e-03 6.00e-03 4.79e-02 4.45e-02 4.03e-02 3.59e-02
> [10,] 7.35e-03 6.23e-03 5.27e-03 4.44e-03 3.04e-02 2.84e-02 2.60e-02 2.33e-02
> [11,] 5.35e-03 4.56e-03 3.88e-03 3.30e-03 1.92e-02 1.82e-02 1.67e-02 1.51e-02
> [12,] 3.90e-03 3.34e-03 2.86e-03 2.44e-03 1.22e-02 1.16e-02 1.08e-02 9.82e-03
> [13,] 2.84e-03 2.45e-03 2.11e-03 1.81e-03 7.73e-03 7.41e-03 6.93e-03 6.37e-03
> [14,] 2.07e-03 1.80e-03 1.55e-03 1.34e-03 4.90e-03 4.74e-03 4.46e-03 4.14e-03
> [15,] 1.51e-03 1.32e-03 1.15e-03 9.95e-04 3.11e-03 3.03e-03 2.88e-03 2.68e-03
> [16,] 1.10e-03 9.65e-04 8.45e-04 7.38e-04 1.97e-03 1.93e-03 1.85e-03 1.74e-03
> [17,] 8.01e-04 7.07e-04 6.23e-04 5.47e-04 1.25e-03 1.24e-03 1.19e-03 1.13e-03
> [18,] 5.83e-04 5.18e-04 4.59e-04 4.06e-04 7.91e-04 7.90e-04 7.68e-04 7.34e-04
>          [,35]    [,36]    [,37]    [,38]    [,39]    [,40]    [,41]    [,42]
>  [1,] 4.94e+01 1.02e+01 2.06e+00 4.18e-01 8.47e-02 1.71e-02 3.47e-03 7.03e-04
>  [2,] 4.94e+01 1.02e+01 2.06e+00 4.18e-01 8.47e-02 1.71e-02 3.47e-03 7.03e-04
>  [3,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [4,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [5,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [6,] 1.15e-02 8.47e-03 6.21e-03 4.55e-03 3.33e-03 2.44e-03 1.78e-03 1.30e-03
>  [7,] 7.40e-02 6.36e-02 5.43e-02 4.61e-02 3.89e-02 3.27e-02 2.75e-02 2.30e-02
>  [8,] 4.84e-02 4.19e-02 3.60e-02 3.08e-02 2.62e-02 2.22e-02 1.88e-02 1.58e-02
>  [9,] 3.16e-02 2.76e-02 2.39e-02 2.06e-02 1.77e-02 1.51e-02 1.28e-02 1.09e-02
> [10,] 2.07e-02 1.82e-02 1.59e-02 1.38e-02 1.19e-02 1.02e-02 8.75e-03 7.47e-03
> [11,] 1.35e-02 1.20e-02 1.05e-02 9.21e-03 8.01e-03 6.93e-03 5.98e-03 5.14e-03
> [12,] 8.85e-03 7.90e-03 7.00e-03 6.16e-03 5.40e-03 4.70e-03 4.08e-03 3.54e-03
> [13,] 5.79e-03 5.21e-03 4.65e-03 4.12e-03 3.63e-03 3.19e-03 2.79e-03 2.43e-03
> [14,] 3.79e-03 3.43e-03 3.08e-03 2.76e-03 2.45e-03 2.16e-03 1.91e-03 1.67e-03
> [15,] 2.48e-03 2.26e-03 2.05e-03 1.84e-03 1.65e-03 1.47e-03 1.30e-03 1.15e-03
> [16,] 1.62e-03 1.49e-03 1.36e-03 1.23e-03 1.11e-03 9.96e-04 8.90e-04 7.91e-04
> [17,] 1.06e-03 9.81e-04 9.02e-04 8.24e-04 7.48e-04 6.76e-04 6.08e-04 5.44e-04
> [18,] 6.92e-04 6.47e-04 5.99e-04 5.51e-04 5.04e-04 4.58e-04 4.15e-04 3.74e-04
>          [,43]    [,44]    [,45]    [,46]    [,47]    [,48]    [,49]    [,50]
>  [1,] 1.42e-04 2.88e-05 5.84e-06 1.18e-06 2.39e-07 4.85e-08 9.81e-09 1.99e-09
>  [2,] 1.42e-04 2.88e-05 5.84e-06 1.18e-06 2.39e-07 4.85e-08 9.81e-09 1.99e-09
>  [3,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [4,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [5,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [6,] 9.54e-04 6.97e-04 5.10e-04 3.73e-04 2.72e-04 1.99e-04 1.46e-04 1.06e-04
>  [7,] 1.92e-02 1.59e-02 1.32e-02 1.10e-02 9.11e-03 7.54e-03 6.23e-03 5.15e-03
>  [8,] 1.33e-02 1.11e-02 9.30e-03 7.76e-03 6.47e-03 5.39e-03 4.49e-03 3.73e-03
>  [9,] 9.19e-03 7.75e-03 6.52e-03 5.48e-03 4.60e-03 3.86e-03 3.23e-03 2.70e-03
> [10,] 6.36e-03 5.40e-03 4.58e-03 3.87e-03 3.27e-03 2.76e-03 2.32e-03 1.96e-03
> [11,] 4.41e-03 3.77e-03 3.21e-03 2.74e-03 2.33e-03 1.97e-03 1.67e-03 1.42e-03
> [12,] 3.05e-03 2.63e-03 2.25e-03 1.93e-03 1.65e-03 1.41e-03 1.20e-03 1.03e-03
> [13,] 2.11e-03 1.83e-03 1.58e-03 1.36e-03 1.18e-03 1.01e-03 8.67e-04 7.43e-04
> [14,] 1.46e-03 1.28e-03 1.11e-03 9.64e-04 8.35e-04 7.22e-04 6.24e-04 5.38e-04
> [15,] 1.01e-03 8.90e-04 7.79e-04 6.81e-04 5.94e-04 5.17e-04 4.49e-04 3.90e-04
> [16,] 7.02e-04 6.20e-04 5.47e-04 4.81e-04 4.22e-04 3.70e-04 3.23e-04 2.82e-04
> [17,] 4.86e-04 4.32e-04 3.84e-04 3.40e-04 3.00e-04 2.64e-04 2.33e-04 2.04e-04
> [18,] 3.37e-04 3.01e-04 2.69e-04 2.40e-04 2.13e-04 1.89e-04 1.67e-04 1.48e-04
>          [,51]    [,52]    [,53]    [,54]
>  [1,] 4.02e-10 8.15e-11 1.65e-11 3.34e-12
>  [2,] 4.02e-10 8.15e-11 1.65e-11 3.34e-12
>  [3,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [4,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [5,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00
>  [6,] 7.78e-05 5.69e-05 4.16e-05 3.04e-05
>  [7,] 4.25e-03 3.51e-03 2.89e-03 2.39e-03
>  [8,] 3.10e-03 2.57e-03 2.13e-03 1.77e-03
>  [9,] 2.26e-03 1.88e-03 1.57e-03 1.31e-03
> [10,] 1.64e-03 1.38e-03 1.16e-03 9.72e-04
> [11,] 1.20e-03 1.01e-03 8.54e-04 7.20e-04
> [12,] 8.73e-04 7.42e-04 6.30e-04 5.34e-04
> [13,] 6.36e-04 5.43e-04 4.64e-04 3.96e-04
> [14,] 4.63e-04 3.98e-04 3.42e-04 2.94e-04
> [15,] 3.37e-04 2.92e-04 2.52e-04 2.18e-04
> [16,] 2.46e-04 2.14e-04 1.86e-04 1.61e-04
> [17,] 1.79e-04 1.57e-04 1.37e-04 1.20e-04
> [18,] 1.30e-04 1.15e-04 1.01e-04 8.87e-05
>  [ reached getOption("max.print") -- omitted 36 rows ]

Notice that this matrix is bigger and more full of non-zero values than the associated raw matrix. This is due to the fact that function-based matrices use kernels to estimate every estimable element, while raw matrices are populated only by proportions of individuals making actual transitions. Thus, if no individuals make a particular transition, then a raw matrix will assign that survival-transition probability a zero, while a function-based matrix will provide a non-zero estimate based on the chosen vital rate linear models.

Let’s now estimate the function-based historical MPM.

cypmatrix3fp <- flefko3(stageframe = cypframe_fb, supplement = cypsupp3_fb, 
  modelsuite = cypmodels3p, data = cypfb_v1)

summary(cypmatrix3fp)
> 
> This historical lefkoMat object contains 15 matrices.
> 
> Each matrix is square with 2916 rows and columns, and a total of 8503056 elements.
> A total of 1768620 survival transitions were estimated, with 117908 per matrix.
> A total of 35280 fecundity transitions were estimated, with 2352 per matrix.
> This lefkoMat object covers 1 population, 3 patches, and 5 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Observation estimated with 70 individuals and 303 individual transitions.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
> 
> 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.000 0.000 0.000 0.000 0.000 0.000
> 1st Qu. 0.965 0.965 0.965 0.965 0.965 0.964 0.965 0.965 0.965 0.964 0.964 0.965
> Median  0.999 1.000 1.000 1.000 0.999 0.998 1.000 0.999 0.999 0.998 0.997 0.999
> Mean    0.820 0.820 0.820 0.820 0.820 0.819 0.820 0.820 0.820 0.819 0.819 0.820
> 3rd Qu. 1.000 1.000 1.000 1.000 1.000 0.999 1.000 1.000 0.999 0.999 0.998 1.000
> Max.    1.000 1.000 1.000 1.000 1.000 0.999 1.000 1.000 1.000 1.000 0.999 1.000
>         [,13] [,14] [,15]
> Min.    0.000 0.000 0.000
> 1st Qu. 0.965 0.964 0.964
> Median  0.999 0.998 0.997
> Mean    0.819 0.819 0.819
> 3rd Qu. 0.999 0.999 0.999
> Max.    1.000 0.999 0.999

Note the similarities and differences with the ahistorical MPMs. First, we see 15 matrices produced again. However, in the raw hMPM, we only saw 12 matrices, because of the fact that we need three monitoring occasions of data to parameterize each raw transition (six monitoring occasions total means four sets of three monitoring occasions, or four total time steps, in the historical case). Since this is a function-based MPM, we can use our functions to estimate transitions in the first year even without a full set of three years of data. The result is five time steps that we can estimate transitions for.

Second, these matrices are utterly huge. They have \(54^2 = 2916\) rows and columns, yielding \(2916^2 = 8,503,056\) total elements per matrix. And we have 120,260 elements estimated per matrix. This is vastly more than in the ahistorical case. However, it is also just a small fraction of the total number of elements in the matrix. In fact, in an unreduced hMPM estimated in Ehrlén format, the total number of elements that can be estimated is equal to \(\frac{r \times c}{m}\), where \(r\) and \(c\) are the numbers of rows and columns, respectively, and \(m\) is the number of stages in the stageframe. For example, if there are 10 stages in the stageframe, then only 1000 of the total 10,000 elements in the hMPM are estimable (10% of the elements).

We can also make deVries format hMPMs. Here, we just add a single option regarding format. The resulting summary shows bigger matrices reflecting the addition of the prior stage for newborns.

cypmatrix3fp_deV <- flefko3(stageframe = cypframe_fb, supplement = cypsupp3_fb, 
  modelsuite = cypmodels3p, data = cypfb_v1, format = "deVries")

summary(cypmatrix3fp_deV)
> 
> This historical lefkoMat object contains 15 matrices.
> 
> Each matrix is square with 2970 rows and columns, and a total of 8820900 elements.
> A total of 1767900 survival transitions were estimated, with 117860 per matrix.
> A total of 35280 fecundity transitions were estimated, with 2352 per matrix.
> This lefkoMat object covers 1 population, 3 patches, and 5 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Observation estimated with 70 individuals and 303 individual transitions.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
> 
> 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.000 0.000 0.000 0.000 0.000 0.000
> 1st Qu. 0.965 0.965 0.965 0.965 0.965 0.964 0.965 0.965 0.965 0.964 0.964 0.965
> Median  0.999 1.000 1.000 1.000 0.999 0.998 0.999 0.999 0.999 0.998 0.997 0.999
> Mean    0.803 0.804 0.804 0.804 0.803 0.803 0.804 0.803 0.803 0.803 0.802 0.803
> 3rd Qu. 1.000 1.000 1.000 1.000 1.000 0.999 1.000 1.000 0.999 0.999 0.998 1.000
> Max.    1.000 1.000 1.000 1.000 1.000 0.999 1.000 1.000 1.000 1.000 0.999 1.000
>         [,13] [,14] [,15]
> Min.    0.000 0.000 0.000
> 1st Qu. 0.965 0.964 0.964
> Median  0.998 0.998 0.997
> Mean    0.803 0.803 0.802
> 3rd Qu. 0.999 0.999 0.999
> Max.    1.000 0.999 0.999

Now that we have created our MPMs, we might wish to create element-wise arithmetic mean matrices to aid inference and further analysis. For example, we might be interested in developing patch-level means and an overall population mean, but one in which the element means treat each patch and each year as equal in proportional effect. For this purpose, we can use the lmean() function.

cyp2fp_mean <- lmean(cypmatrix2fp)
cyp3fp_mean <- lmean(cypmatrix3fp)

summary(cyp2fp_mean)
> 
> This ahistorical lefkoMat object contains 4 matrices.
> 
> Each matrix is square with 54 rows and columns, and a total of 2916 elements.
> A total of 9644 survival transitions were estimated, with 2411 per matrix.
> A total of 192 fecundity transitions were estimated, with 48 per matrix.
> This lefkoMat object covers 1 population, 4 patches, and 0 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Observation estimated with 70 individuals and 303 individual transitions.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
> 
> Survival probability sum check (each matrix represented by column in order):
>          [,1]  [,2]  [,3]  [,4]
> Min.    0.050 0.050 0.050 0.050
> 1st Qu. 0.991 0.991 0.991 0.991
> Median  1.000 0.999 0.998 0.999
> Mean    0.916 0.916 0.914 0.916
> 3rd Qu. 1.000 0.999 0.999 0.999
> Max.    1.000 1.000 0.999 1.000
summary(cyp3fp_mean)
> 
> This historical lefkoMat object contains 4 matrices.
> 
> Each matrix is square with 2916 rows and columns, and a total of 8503056 elements.
> A total of 471632 survival transitions were estimated, with 117908 per matrix.
> A total of 9408 fecundity transitions were estimated, with 2352 per matrix.
> This lefkoMat object covers 1 population, 4 patches, and 0 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Observation estimated with 70 individuals and 303 individual transitions.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
> 
> 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.965 0.965 0.964 0.965
> Median  1.000 0.999 0.998 0.999
> Mean    0.820 0.820 0.819 0.820
> 3rd Qu. 1.000 0.999 0.999 0.999
> Max.    1.000 1.000 0.999 1.000

The function-based mean matrices have the same number of estimated transitions as their constituent matrices, because all elements that are estimable are actually estimated using the vital rate models supplied. We have four mean matrices in each case - the first three are patch-level means for each of three patches, and the final matrix is the population mean matrix. The labels element in each lefkoMat object shows this.

In the preceding vital rate model sets, we created best-fit models that necessarily included a patch term. This will result necessarily in the creation of matrices at the patch level. To create matrices for the population level only, we remove the patch option from the preceding function calls, as below.

cypmodels3 <- modelsearch(cypfb_v1, historical = TRUE, 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"),
  quiet = TRUE)

cypmodels2 <- modelsearch(cypfb_v1, 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"),
  quiet = TRUE)

Let’s see some summaries.

summary(cypmodels3)
> 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 ~ size2added + (1 | year2) + (1 | individ)
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
> 128.1324 143.2057 -60.0662 120.1324      316 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 1.198361
>  year2   (Intercept) 0.008826
> Number of obs: 320, groups:  individ, 74; year2, 5
> Fixed Effects:
> (Intercept)   size2added  
>      2.0352       0.6344  
> 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 ~ (1 | year2) + (1 | individ)
> Data: subdata
>       AIC       BIC    logLik  df.resid 
> 1008.2748 1022.9266 -500.1374       284 
> Random-effects (co)variances:
> 
> Conditional model:
>  Groups  Name        Std.Dev.
>  year2   (Intercept) 0.1109  
>  individ (Intercept) 1.0562  
> 
> Number of obs: 288 / Conditional model: year2, 5; individ, 70
> 
> Dispersion parameter for truncated_nbinom2 family (): 9.81e+08 
> 
> Fixed Effects:
> 
> Conditional model:
> (Intercept)  
>       0.576  
> 
> 
> 
> 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 ~ repstatus2 + size2added + (1 | year2) + (1 | individ)
>    Data: subdata
>       AIC       BIC    logLik  deviance  df.resid 
>  333.6176  351.9324 -161.8088  323.6176       283 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 0.1829  
>  year2   (Intercept) 0.6250  
> Number of obs: 288, groups:  individ, 70; year2, 5
> Fixed Effects:
> (Intercept)   repstatus2   size2added  
>     -1.4630       1.6457       0.1715  
> 
> 
> 
> 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: 16
> 
> Number of models in observation table: 16
> 
> Number of models in size table: 16
> 
> Number of models in secondary size table: 1
> 
> Number of models in tertiary size table: 1
> 
> Number of models in reproduction status table: 16
> 
> Number of models in fecundity table: 239
> 
> Number of models in juvenile survival table: 1
> 
> Number of models in juvenile observation table: 1
> 
> Number of models in juvenile size table: 1
> 
> Number of models in juvenile secondary size table: 1
> 
> Number of models in juvenile tertiary size table: 1
> 
> Number of models in juvenile reproduction table: 1
> 
> Number of models in juvenile maturity table: 1
> 
> 
> 
> 
> 
> General model parameter names (column 1), and 
> specific names used in these models (column 2): 
>                       parameter_names mainparams
> 1                              time t      year2
> 2                          individual    individ
> 3                               patch      patch
> 4                   alive in time t+1      surv3
> 5                observed in time t+1       obs3
> 6                   sizea in time t+1      size3
> 7                   sizeb in time t+1     sizeb3
> 8                   sizec in time t+1     sizec3
> 9     reproductive status in time t+1     repst3
> 10              fecundity in time t+1       fec3
> 11                fecundity in time t       fec2
> 12                    sizea in time t      size2
> 13                  sizea in time t-1      size1
> 14                    sizeb in time t     sizeb2
> 15                  sizeb in time t-1     sizeb1
> 16                    sizec in time t     sizec2
> 17                  sizec in time t-1     sizec1
> 18      reproductive status in time t     repst2
> 19    reproductive status in time t-1     repst1
> 20        maturity status in time t+1     matst3
> 21          maturity status in time t     matst2
> 22                      age in time t        age
> 23                  density in time t    density
> 24   individual covariate a in time t   indcova2
> 25 individual covariate a in time t-1   indcova1
> 26   individual covariate b in time t   indcovb2
> 27 individual covariate b in time t-1   indcovb1
> 28   individual covariate c in time t   indcovc2
> 29 individual covariate c in time t-1   indcovc1
> 30              stage group in time t     group2
> 31            stage group in time t-1     group1
> 
> 
> 
> 
> 
> Quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Survival accuracy is 0.947.
> Observation estimated with 70 individuals and 303 individual transitions.
> Observation accuracy is 0.95.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Primary size pseudo R-squared is 0.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Reproductive status accuracy is 0.715.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Fecundity pseudo R-squared is 0.32.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity probability not estimated.
summary(cypmodels2)
> 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 ~ size2added + (1 | year2) + (1 | individ)
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
> 128.1324 143.2057 -60.0662 120.1324      316 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 1.198361
>  year2   (Intercept) 0.008826
> Number of obs: 320, groups:  individ, 74; year2, 5
> Fixed Effects:
> (Intercept)   size2added  
>      2.0352       0.6344  
> 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 ~ (1 | year2) + (1 | individ)
> Data: subdata
>       AIC       BIC    logLik  df.resid 
> 1008.2748 1022.9266 -500.1374       284 
> Random-effects (co)variances:
> 
> Conditional model:
>  Groups  Name        Std.Dev.
>  year2   (Intercept) 0.1109  
>  individ (Intercept) 1.0562  
> 
> Number of obs: 288 / Conditional model: year2, 5; individ, 70
> 
> Dispersion parameter for truncated_nbinom2 family (): 9.81e+08 
> 
> Fixed Effects:
> 
> Conditional model:
> (Intercept)  
>       0.576  
> 
> 
> 
> 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 ~ repstatus2 + size2added + (1 | year2) + (1 | individ)
>    Data: subdata
>       AIC       BIC    logLik  deviance  df.resid 
>  333.6176  351.9324 -161.8088  323.6176       283 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 0.1829  
>  year2   (Intercept) 0.6250  
> Number of obs: 288, groups:  individ, 70; year2, 5
> Fixed Effects:
> (Intercept)   repstatus2   size2added  
>     -1.4630       1.6457       0.1715  
> 
> 
> 
> 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: 4
> 
> Number of models in observation table: 4
> 
> Number of models in size table: 4
> 
> Number of models in secondary size table: 1
> 
> Number of models in tertiary size table: 1
> 
> Number of models in reproduction status table: 4
> 
> Number of models in fecundity table: 15
> 
> Number of models in juvenile survival table: 1
> 
> Number of models in juvenile observation table: 1
> 
> Number of models in juvenile size table: 1
> 
> Number of models in juvenile secondary size table: 1
> 
> Number of models in juvenile tertiary size table: 1
> 
> Number of models in juvenile reproduction table: 1
> 
> Number of models in juvenile maturity table: 1
> 
> 
> 
> 
> 
> General model parameter names (column 1), and 
> specific names used in these models (column 2): 
>                       parameter_names mainparams
> 1                              time t      year2
> 2                          individual    individ
> 3                               patch      patch
> 4                   alive in time t+1      surv3
> 5                observed in time t+1       obs3
> 6                   sizea in time t+1      size3
> 7                   sizeb in time t+1     sizeb3
> 8                   sizec in time t+1     sizec3
> 9     reproductive status in time t+1     repst3
> 10              fecundity in time t+1       fec3
> 11                fecundity in time t       fec2
> 12                    sizea in time t      size2
> 13                  sizea in time t-1      size1
> 14                    sizeb in time t     sizeb2
> 15                  sizeb in time t-1     sizeb1
> 16                    sizec in time t     sizec2
> 17                  sizec in time t-1     sizec1
> 18      reproductive status in time t     repst2
> 19    reproductive status in time t-1     repst1
> 20        maturity status in time t+1     matst3
> 21          maturity status in time t     matst2
> 22                      age in time t        age
> 23                  density in time t    density
> 24   individual covariate a in time t   indcova2
> 25 individual covariate a in time t-1   indcova1
> 26   individual covariate b in time t   indcovb2
> 27 individual covariate b in time t-1   indcovb1
> 28   individual covariate c in time t   indcovc2
> 29 individual covariate c in time t-1   indcovc1
> 30              stage group in time t     group2
> 31            stage group in time t-1     group1
> 
> 
> 
> 
> 
> Quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Survival accuracy is 0.947.
> Observation estimated with 70 individuals and 303 individual transitions.
> Observation accuracy is 0.95.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Primary size pseudo R-squared is 0.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Reproductive status accuracy is 0.715.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Fecundity pseudo R-squared is 0.32.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity probability not estimated.

A close look at the model output shows that the best-fit models have changed in structure. First, patch is no longer included - it had been included previously as a random factor in all models. The lack of patch as a random factor has also led to many other changes in the best-fit models. Let’s use these models to create new MPMs.

cypmatrix2f <- flefko2(stageframe = cypframe_fb, supplement = cypsupp2_fb, 
  modelsuite = cypmodels2, data = cypfb_v1)
cypmatrix3f <- flefko3(stageframe = cypframe_fb, supplement = cypsupp3_fb, 
  modelsuite = cypmodels3, data = cypfb_v1)

summary(cypmatrix3f)
> 
> This historical lefkoMat object contains 5 matrices.
> 
> Each matrix is square with 2916 rows and columns, and a total of 8503056 elements.
> A total of 589540 survival transitions were estimated, with 117908 per matrix.
> A total of 11760 fecundity transitions were estimated, with 2352 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 5 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Observation estimated with 70 individuals and 303 individual transitions.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
> 
> Survival probability sum check (each matrix represented by column in order):
>          [,1]  [,2]  [,3]  [,4]  [,5]
> Min.    0.000 0.000 0.000 0.000 0.000
> 1st Qu. 0.965 0.965 0.965 0.965 0.965
> Median  1.000 1.000 1.000 1.000 1.000
> Mean    0.820 0.820 0.820 0.820 0.820
> 3rd Qu. 1.000 1.000 1.000 1.000 1.000
> Max.    1.000 1.000 1.000 1.000 1.000
summary(cypmatrix2f)
> 
> This ahistorical lefkoMat object contains 5 matrices.
> 
> Each matrix is square with 54 rows and columns, and a total of 2916 elements.
> A total of 12055 survival transitions were estimated, with 2411 per matrix.
> A total of 240 fecundity transitions were estimated, with 48 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 5 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Observation estimated with 70 individuals and 303 individual transitions.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
> 
> Survival probability sum check (each matrix represented by column in order):
>          [,1]  [,2]  [,3]  [,4]  [,5]
> Min.    0.050 0.050 0.050 0.050 0.050
> 1st Qu. 0.991 0.991 0.991 0.991 0.991
> Median  1.000 1.000 1.000 1.000 1.000
> Mean    0.916 0.917 0.917 0.917 0.914
> 3rd Qu. 1.000 1.000 1.000 1.000 1.000
> Max.    1.000 1.000 1.000 1.000 1.000

We see that the main difference here is that there are only five matrices in each MPM, since the matrices are now annual matrices representing the whole population rather than each patch.

5.6 Using stage groups for complex MPMs

Let’s imagine a more complex life history. For example, imagine that our study organism had a life history in which the juvenile morphology inevitably led to stages that were non-reproductive but otherwise just like the reproductive adult stages, and that once the organism became reproductive it stayed so until death. The resulting life history model might look like this (figure 5.3).

Figure 5.3: Imaginary life history model with permanent transitions from non-reproductive to reproductive adult stages

Intriguingly, we do not need to alter our code much to account for this subtle but nonetheless major change. Indeed, we really need to alter just two steps. First, we will add stage groups to our stage frame. We will assign the dormant seed stage, three protocorm stages, and seedling stage to a juvenile stage group; the dormant and non-reproductive adult stages to an adult stage group; and the reproductive adult stages to a second adult stage group. This only requires adding one further vector coding these groups. Please bear in mind that stage groups must be numbered sequentially, and the first stage group must be group 0.

groupvec <- c(rep(0, 5), rep(1, 25), rep(2, 24))

cypframe_fb_group <- sf_create(sizes = sizevector, stagenames = stagevector, 
  repstatus = repvector, obsstatus = obsvector, matstatus = matvector, 
  propstatus = propvector, immstatus = immvector, indataset = indataset,
  group = groupvec, comments = comments)

cypframe_fb_group
>    stage size size_b size_c min_age max_age repstatus obsstatus propstatus
> 1     SD    0     NA     NA      NA      NA         0         0          1
> 2     P1    0     NA     NA      NA      NA         0         0          0
> 3     P2    0     NA     NA      NA      NA         0         0          0
> 4     P3    0     NA     NA      NA      NA         0         0          0
> 5     SL    0     NA     NA      NA      NA         0         0          0
> 6      D    0     NA     NA      NA      NA         0         0          0
> 7     V1    1     NA     NA      NA      NA         0         1          0
> 8     V2    2     NA     NA      NA      NA         0         1          0
> 9     V3    3     NA     NA      NA      NA         0         1          0
> 10    V4    4     NA     NA      NA      NA         0         1          0
> 11    V5    5     NA     NA      NA      NA         0         1          0
> 12    V6    6     NA     NA      NA      NA         0         1          0
> 13    V7    7     NA     NA      NA      NA         0         1          0
> 14    V8    8     NA     NA      NA      NA         0         1          0
> 15    V9    9     NA     NA      NA      NA         0         1          0
> 16   V10   10     NA     NA      NA      NA         0         1          0
> 17   V11   11     NA     NA      NA      NA         0         1          0
> 18   V12   12     NA     NA      NA      NA         0         1          0
> 19   V13   13     NA     NA      NA      NA         0         1          0
> 20   V14   14     NA     NA      NA      NA         0         1          0
> 21   V15   15     NA     NA      NA      NA         0         1          0
> 22   V16   16     NA     NA      NA      NA         0         1          0
> 23   V17   17     NA     NA      NA      NA         0         1          0
> 24   V18   18     NA     NA      NA      NA         0         1          0
> 25   V19   19     NA     NA      NA      NA         0         1          0
> 26   V20   20     NA     NA      NA      NA         0         1          0
> 27   V21   21     NA     NA      NA      NA         0         1          0
> 28   V22   22     NA     NA      NA      NA         0         1          0
> 29   V23   23     NA     NA      NA      NA         0         1          0
> 30   V24   24     NA     NA      NA      NA         0         1          0
> 31    F1    1     NA     NA      NA      NA         1         1          0
> 32    F2    2     NA     NA      NA      NA         1         1          0
> 33    F3    3     NA     NA      NA      NA         1         1          0
> 34    F4    4     NA     NA      NA      NA         1         1          0
>    immstatus matstatus indataset binhalfwidth_raw sizebin_min sizebin_max
> 1          0         0         0              0.5        -0.5         0.5
> 2          1         0         0              0.5        -0.5         0.5
> 3          1         0         0              0.5        -0.5         0.5
> 4          1         0         0              0.5        -0.5         0.5
> 5          1         0         0              0.5        -0.5         0.5
> 6          0         1         1              0.5        -0.5         0.5
> 7          0         1         1              0.5         0.5         1.5
> 8          0         1         1              0.5         1.5         2.5
> 9          0         1         1              0.5         2.5         3.5
> 10         0         1         1              0.5         3.5         4.5
> 11         0         1         1              0.5         4.5         5.5
> 12         0         1         1              0.5         5.5         6.5
> 13         0         1         1              0.5         6.5         7.5
> 14         0         1         1              0.5         7.5         8.5
> 15         0         1         1              0.5         8.5         9.5
> 16         0         1         1              0.5         9.5        10.5
> 17         0         1         1              0.5        10.5        11.5
> 18         0         1         1              0.5        11.5        12.5
> 19         0         1         1              0.5        12.5        13.5
> 20         0         1         1              0.5        13.5        14.5
> 21         0         1         1              0.5        14.5        15.5
> 22         0         1         1              0.5        15.5        16.5
> 23         0         1         1              0.5        16.5        17.5
> 24         0         1         1              0.5        17.5        18.5
> 25         0         1         1              0.5        18.5        19.5
> 26         0         1         1              0.5        19.5        20.5
> 27         0         1         1              0.5        20.5        21.5
> 28         0         1         1              0.5        21.5        22.5
> 29         0         1         1              0.5        22.5        23.5
> 30         0         1         1              0.5        23.5        24.5
> 31         0         1         1              0.5         0.5         1.5
> 32         0         1         1              0.5         1.5         2.5
> 33         0         1         1              0.5         2.5         3.5
> 34         0         1         1              0.5         3.5         4.5
>    sizebin_center sizebin_width binhalfwidthb_raw sizebinb_min sizebinb_max
> 1               0             1                NA           NA           NA
> 2               0             1                NA           NA           NA
> 3               0             1                NA           NA           NA
> 4               0             1                NA           NA           NA
> 5               0             1                NA           NA           NA
> 6               0             1                NA           NA           NA
> 7               1             1                NA           NA           NA
> 8               2             1                NA           NA           NA
> 9               3             1                NA           NA           NA
> 10              4             1                NA           NA           NA
> 11              5             1                NA           NA           NA
> 12              6             1                NA           NA           NA
> 13              7             1                NA           NA           NA
> 14              8             1                NA           NA           NA
> 15              9             1                NA           NA           NA
> 16             10             1                NA           NA           NA
> 17             11             1                NA           NA           NA
> 18             12             1                NA           NA           NA
> 19             13             1                NA           NA           NA
> 20             14             1                NA           NA           NA
> 21             15             1                NA           NA           NA
> 22             16             1                NA           NA           NA
> 23             17             1                NA           NA           NA
> 24             18             1                NA           NA           NA
> 25             19             1                NA           NA           NA
> 26             20             1                NA           NA           NA
> 27             21             1                NA           NA           NA
> 28             22             1                NA           NA           NA
> 29             23             1                NA           NA           NA
> 30             24             1                NA           NA           NA
> 31              1             1                NA           NA           NA
> 32              2             1                NA           NA           NA
> 33              3             1                NA           NA           NA
> 34              4             1                NA           NA           NA
>    sizebinb_center sizebinb_width binhalfwidthc_raw sizebinc_min sizebinc_max
> 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
> 12              NA             NA                NA           NA           NA
> 13              NA             NA                NA           NA           NA
> 14              NA             NA                NA           NA           NA
> 15              NA             NA                NA           NA           NA
> 16              NA             NA                NA           NA           NA
> 17              NA             NA                NA           NA           NA
> 18              NA             NA                NA           NA           NA
> 19              NA             NA                NA           NA           NA
> 20              NA             NA                NA           NA           NA
> 21              NA             NA                NA           NA           NA
> 22              NA             NA                NA           NA           NA
> 23              NA             NA                NA           NA           NA
> 24              NA             NA                NA           NA           NA
> 25              NA             NA                NA           NA           NA
> 26              NA             NA                NA           NA           NA
> 27              NA             NA                NA           NA           NA
> 28              NA             NA                NA           NA           NA
> 29              NA             NA                NA           NA           NA
> 30              NA             NA                NA           NA           NA
> 31              NA             NA                NA           NA           NA
> 32              NA             NA                NA           NA           NA
> 33              NA             NA                NA           NA           NA
> 34              NA             NA                NA           NA           NA
>    sizebinc_center sizebinc_width group           comments
> 1               NA             NA     0       Dormant seed
> 2               NA             NA     0      Yr1 protocorm
> 3               NA             NA     0      Yr2 protocorm
> 4               NA             NA     0      Yr3 protocorm
> 5               NA             NA     0           Seedling
> 6               NA             NA     1           Veg dorm
> 7               NA             NA     1   Veg adult 1 stem
> 8               NA             NA     1  Veg adult 2 stems
> 9               NA             NA     1  Veg adult 3 stems
> 10              NA             NA     1  Veg adult 4 stems
> 11              NA             NA     1  Veg adult 5 stems
> 12              NA             NA     1  Veg adult 6 stems
> 13              NA             NA     1  Veg adult 7 stems
> 14              NA             NA     1  Veg adult 8 stems
> 15              NA             NA     1  Veg adult 9 stems
> 16              NA             NA     1 Veg adult 10 stems
> 17              NA             NA     1 Veg adult 11 stems
> 18              NA             NA     1 Veg adult 12 stems
> 19              NA             NA     1 Veg adult 13 stems
> 20              NA             NA     1 Veg adult 14 stems
> 21              NA             NA     1 Veg adult 15 stems
> 22              NA             NA     1 Veg adult 16 stems
> 23              NA             NA     1 Veg adult 17 stems
> 24              NA             NA     1 Veg adult 18 stems
> 25              NA             NA     1 Veg adult 19 stems
> 26              NA             NA     1 Veg adult 20 stems
> 27              NA             NA     1 Veg adult 21 stems
> 28              NA             NA     1 Veg adult 22 stems
> 29              NA             NA     1 Veg adult 23 stems
> 30              NA             NA     1 Veg adult 24 stems
> 31              NA             NA     2   Flo adult 1 stem
> 32              NA             NA     2  Flo adult 2 stems
> 33              NA             NA     2  Flo adult 3 stems
> 34              NA             NA     2  Flo adult 4 stems
>  [ reached 'max' / getOption("max.print") -- omitted 20 rows ]

We now need to add an extra line to the supplement table to set transitions from group 2 to group 1 to a value of 0.0. We will do this with our ahistorical supplement table and create only ahistorical MPMs, but this may also be done with the historical supplement table to create historical MPMs.

cypsupp2_fb_group <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "SL",
    "D", "V1", "V2", "V3", "SD", "P1", "group1"), 
  stage2 = c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "SL", "SL", "rep",
    "rep", "group2"), 
  eststage3 = c(NA, NA, NA, NA, NA, NA, "D", "V1", "V2", "V3", NA, NA, NA), 
  eststage2 = c(NA, NA, NA, NA, NA, NA, "D", "D", "D", "D", NA, NA, NA), 
  givenrate = c(0.08, 0.1, 0.1, 0.1, 0.05, 0.05, NA, NA, NA, NA, NA, NA, 0.0), 
  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, NA),
  type = c("S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "R", "R", "S"),
  stageframe = cypframe_fb_group, historical = FALSE)

cypsupp2_fb_group
>    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
> 13 group1 group2   <NA>      <NA>      <NA>      <NA>      0.00        1.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
> 13        1            1

Finally, let’s create our new ahistorical MPM, and take a look at the first matrix.

cypmatrix2fp_group <- flefko2(stageframe = cypframe_fb_group, data = cypfb_v1,
  supplement = cypsupp2_fb_group, modelsuite = cypmodels2p)
cypmatrix2fp_group$A[[1]]
>       [,1] [,2] [,3] [,4]       [,5]         [,6]         [,7]         [,8]
>  [1,] 0.08  0.0  0.0 0.00 0.00000000 0.000000e+00 0.000000e+00 0.000000e+00
>  [2,] 0.10  0.0  0.0 0.00 0.00000000 0.000000e+00 0.000000e+00 0.000000e+00
>  [3,] 0.00  0.1  0.0 0.00 0.00000000 0.000000e+00 0.000000e+00 0.000000e+00
>  [4,] 0.00  0.0  0.1 0.00 0.00000000 0.000000e+00 0.000000e+00 0.000000e+00
>  [5,] 0.00  0.0  0.0 0.05 0.05000000 0.000000e+00 0.000000e+00 0.000000e+00
>  [6,] 0.00  0.0  0.0 0.00 0.03291476 4.702108e-02 3.686949e-02 2.809354e-02
>  [7,] 0.00  0.0  0.0 0.00 0.16596277 2.370897e-01 2.397022e-01 2.339950e-01
>  [8,] 0.00  0.0  0.0 0.00 0.10435961 1.490852e-01 1.519502e-01 1.495185e-01
>  [9,] 0.00  0.0  0.0 0.00 0.06562272 9.374674e-02 9.632315e-02 9.553964e-02
> [10,] 0.00  0.0  0.0 0.00 0.00000000 5.894920e-02 6.106045e-02 6.104809e-02
> [11,] 0.00  0.0  0.0 0.00 0.00000000 3.706805e-02 3.870698e-02 3.900862e-02
> [12,] 0.00  0.0  0.0 0.00 0.00000000 2.330888e-02 2.453684e-02 2.492580e-02
> [13,] 0.00  0.0  0.0 0.00 0.00000000 1.465694e-02 1.555421e-02 1.592713e-02
> [14,] 0.00  0.0  0.0 0.00 0.00000000 9.216478e-03 9.860008e-03 1.017715e-02
> [15,] 0.00  0.0  0.0 0.00 0.00000000 5.795444e-03 6.250383e-03 6.503013e-03
> [16,] 0.00  0.0  0.0 0.00 0.00000000 3.644253e-03 3.962196e-03 4.155307e-03
> [17,] 0.00  0.0  0.0 0.00 0.00000000 2.291555e-03 2.511685e-03 2.655166e-03
> [18,] 0.00  0.0  0.0 0.00 0.00000000 1.440960e-03 1.592189e-03 1.696603e-03
>               [,9]        [,10]        [,11]        [,12]        [,13]
>  [1,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
>  [2,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
>  [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,] 2.104752e-02 1.561422e-02 1.151723e-02 8.466546e-03 6.211279e-03
>  [7,] 2.230452e-01 2.089844e-01 1.931975e-01 1.765921e-01 1.597915e-01
>  [8,] 1.436456e-01 1.356365e-01 1.263515e-01 1.163640e-01 1.060774e-01
>  [9,] 9.251071e-02 8.803180e-02 8.263404e-02 7.667716e-02 7.041935e-02
> [10,] 5.957878e-02 5.713503e-02 5.404278e-02 5.052582e-02 4.674780e-02
> [11,] 3.836995e-02 3.708219e-02 3.534406e-02 3.329360e-02 3.103347e-02
> [12,] 2.471103e-02 2.406735e-02 2.311507e-02 2.193856e-02 2.060153e-02
> [13,] 1.591440e-02 1.562037e-02 1.511729e-02 1.445625e-02 1.367630e-02
> [14,] 1.024920e-02 1.013804e-02 9.886728e-03 9.525831e-03 9.078996e-03
> [15,] 6.600694e-03 6.579865e-03 6.465935e-03 6.276973e-03 6.027080e-03
> [16,] 4.250981e-03 4.270512e-03 4.228731e-03 4.136163e-03 4.001070e-03
> [17,] 2.737719e-03 2.771678e-03 2.765596e-03 2.725492e-03 2.656105e-03
> [18,] 1.763147e-03 1.798895e-03 1.808704e-03 1.795942e-03 1.763252e-03
>              [,14]        [,15]        [,16]        [,17]        [,18]
>  [1,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
>  [2,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
>  [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,] 4.551060e-03 3.331976e-03 2.438204e-03 1.783579e-03 1.304415e-03
>  [7,] 1.432476e-01 1.272995e-01 1.122010e-01 9.813457e-02 8.521888e-02
>  [8,] 9.579259e-02 8.574330e-02 7.611208e-02 6.703731e-02 5.861693e-02
>  [9,] 6.405844e-02 5.775289e-02 5.163098e-02 4.579427e-02 4.031905e-02
> [10,] 4.283717e-02 3.889979e-02 3.502411e-02 3.128281e-02 2.773304e-02
> [11,] 2.864608e-02 2.620118e-02 2.375877e-02 2.136979e-02 1.907589e-02
> [12,] 1.915621e-02 1.764795e-02 1.611687e-02 1.459805e-02 1.312115e-02
> [13,] 1.281015e-02 1.188688e-02 1.093295e-02 9.972165e-03 9.025248e-03
> [14,] 8.566404e-03 8.006476e-03 7.416420e-03 6.812147e-03 6.207923e-03
> [15,] 5.728527e-03 5.392808e-03 5.030963e-03 4.653488e-03 4.270056e-03
> [16,] 3.830782e-03 3.632356e-03 3.412777e-03 3.178873e-03 2.937114e-03
> [17,] 2.561721e-03 2.446595e-03 2.315073e-03 2.171539e-03 2.020263e-03
> [18,] 1.713075e-03 1.647918e-03 1.570441e-03 1.483414e-03 1.389617e-03
>              [,19]        [,20]        [,21]        [,22]        [,23]
>  [1,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
>  [2,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
>  [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,] 9.538316e-04 6.973975e-04 5.098658e-04 3.727414e-04 2.724849e-04
>  [7,] 7.351599e-02 6.303915e-02 5.376166e-02 4.562633e-02 3.855452e-02
>  [8,] 5.091162e-02 4.394899e-02 3.772861e-02 3.222771e-02 2.740702e-02
>  [9,] 3.525754e-02 3.063991e-02 2.647701e-02 2.276373e-02 1.948267e-02
> [10,] 2.441671e-02 2.136122e-02 1.858091e-02 1.607894e-02 1.384953e-02
> [11,] 1.690917e-02 1.489240e-02 1.303962e-02 1.135720e-02 9.845132e-03
> [12,] 1.171001e-02 1.038253e-02 9.150886e-03 8.022048e-03 6.998551e-03
> [13,] 8.109470e-03 7.238390e-03 6.421866e-03 5.666295e-03 4.975019e-03
> [14,] 5.616006e-03 5.046387e-03 4.506707e-03 4.002333e-03 3.536562e-03
> [15,] 3.889221e-03 3.518189e-03 3.162696e-03 2.827009e-03 2.514015e-03
> [16,] 2.693381e-03 2.452775e-03 2.219502e-03 1.996830e-03 1.787123e-03
> [17,] 1.865232e-03 1.710001e-03 1.557592e-03 1.410441e-03 1.270401e-03
> [18,] 1.291719e-03 1.192161e-03 1.093079e-03 9.962515e-04 9.030828e-04
>              [,24]        [,25]        [,26]        [,27]        [,28]
>  [1,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
>  [2,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
>  [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,] 1.991889e-04 1.456059e-04 1.064355e-04 7.780171e-05 5.687069e-05
>  [7,] 3.245447e-02 2.722820e-02 2.277695e-02 1.900522e-02 1.582343e-02
>  [8,] 2.321618e-02 1.959849e-02 1.649470e-02 1.384604e-02 1.159620e-02
>  [9,] 1.660761e-02 1.410672e-02 1.194520e-02 1.008737e-02 8.498270e-03
> [10,] 1.188019e-02 1.015383e-02 8.650523e-03 7.349041e-03 6.227956e-03
> [11,] 8.498453e-03 7.308589e-03 6.264571e-03 5.354060e-03 4.564156e-03
> [12,] 6.079337e-03 5.260623e-03 4.536702e-03 3.900640e-03 3.344841e-03
> [13,] 4.348831e-03 3.786525e-03 3.285407e-03 2.841767e-03 2.451266e-03
> [14,] 3.110920e-03 2.725490e-03 2.379240e-03 2.070337e-03 1.796410e-03
> [15,] 2.225386e-03 1.961770e-03 1.723007e-03 1.508320e-03 1.316499e-03
> [16,] 1.591921e-03 1.412056e-03 1.247774e-03 1.098870e-03 9.647960e-04
> [17,] 1.138775e-03 1.016379e-03 9.036183e-04 8.005690e-04 7.070506e-04
> [18,] 8.146185e-04 7.315756e-04 6.543859e-04 5.832454e-04 5.181619e-04
>              [,29]        [,30]        [,31]        [,32]        [,33]
>  [1,] 0.000000e+00 0.000000e+00 2.307477e+03 1.770474e+03 8.237340e+02
>  [2,] 0.000000e+00 0.000000e+00 2.307477e+03 1.770474e+03 8.237340e+02
>  [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,] 4.157052e-05 3.038649e-05 0.000000e+00 0.000000e+00 0.000000e+00
>  [7,] 1.314953e-02 1.090980e-02 0.000000e+00 0.000000e+00 0.000000e+00
>  [8,] 9.692744e-03 8.087878e-03 0.000000e+00 0.000000e+00 0.000000e+00
>  [9,] 7.144689e-03 5.995873e-03 0.000000e+00 0.000000e+00 0.000000e+00
> [10,] 5.266474e-03 4.444984e-03 0.000000e+00 0.000000e+00 0.000000e+00
> [11,] 3.882009e-03 3.295247e-03 0.000000e+00 0.000000e+00 0.000000e+00
> [12,] 2.861496e-03 2.442901e-03 0.000000e+00 0.000000e+00 0.000000e+00
> [13,] 2.109258e-03 1.811021e-03 0.000000e+00 0.000000e+00 0.000000e+00
> [14,] 1.554771e-03 1.342584e-03 0.000000e+00 0.000000e+00 0.000000e+00
> [15,] 1.146048e-03 9.953118e-04 0.000000e+00 0.000000e+00 0.000000e+00
> [16,] 8.447719e-04 7.378651e-04 0.000000e+00 0.000000e+00 0.000000e+00
> [17,] 6.226960e-04 5.470094e-04 0.000000e+00 0.000000e+00 0.000000e+00
> [18,] 4.590000e-04 4.055203e-04 0.000000e+00 0.000000e+00 0.000000e+00
>              [,34]        [,35]        [,36]        [,37]        [,38]
>  [1,] 2.262468e+02 4.937545e+01 1.015784e+01 2.063507e+00 4.181065e-01
>  [2,] 2.262468e+02 4.937545e+01 1.015784e+01 2.063507e+00 4.181065e-01
>  [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,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [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
>              [,39]        [,40]        [,41]        [,42]        [,43]
>  [1,] 8.467204e-02 1.714537e-02 3.471717e-03 7.029750e-04 1.423427e-04
>  [2,] 8.467204e-02 1.714537e-02 3.471717e-03 7.029750e-04 1.423427e-04
>  [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,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [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
>              [,44]        [,45]        [,46]        [,47]        [,48]
>  [1,] 2.882241e-05 5.836135e-06 1.181736e-06 2.392850e-07 4.845188e-08
>  [2,] 2.882241e-05 5.836135e-06 1.181736e-06 2.392850e-07 4.845188e-08
>  [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,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [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
>              [,49]        [,50]        [,51]        [,52]        [,53]
>  [1,] 9.810831e-09 1.986556e-09 4.022499e-10 8.145000e-11 1.649249e-11
>  [2,] 9.810831e-09 1.986556e-09 4.022499e-10 8.145000e-11 1.649249e-11
>  [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,] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
> [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
>              [,54]
>  [1,] 3.339499e-12
>  [2,] 3.339499e-12
>  [3,] 0.000000e+00
>  [4,] 0.000000e+00
>  [5,] 0.000000e+00
>  [6,] 0.000000e+00
>  [7,] 0.000000e+00
>  [8,] 0.000000e+00
>  [9,] 0.000000e+00
> [10,] 0.000000e+00
> [11,] 0.000000e+00
> [12,] 0.000000e+00
> [13,] 0.000000e+00
> [14,] 0.000000e+00
> [15,] 0.000000e+00
> [16,] 0.000000e+00
> [17,] 0.000000e+00
> [18,] 0.000000e+00
>  [ reached getOption("max.print") -- omitted 36 rows ]

Notice that a chunk of the entries are now zeros. We can also explore a summary of this MPM to see the reduction in estimated elements, particularly in comparison to the original ahistorical MPM.

summary(cypmatrix2fp_group)
> 
> This ahistorical lefkoMat object contains 15 matrices.
> 
> Each matrix is square with 54 rows and columns, and a total of 2916 elements.
> A total of 27165 survival transitions were estimated, with 1811 per matrix.
> A total of 720 fecundity transitions were estimated, with 48 per matrix.
> This lefkoMat object covers 1 population, 3 patches, and 5 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Observation estimated with 70 individuals and 303 individual transitions.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
> 
> 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.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050
> 1st Qu. 0.879 0.743 0.802 0.867 0.927 0.840 0.670 0.740 0.824 0.900 0.888 0.786
> Median  0.979 0.956 0.965 0.978 0.988 0.970 0.939 0.955 0.968 0.982 0.980 0.963
> Mean    0.865 0.819 0.838 0.858 0.887 0.850 0.796 0.818 0.843 0.877 0.872 0.832
> 3rd Qu. 0.999 1.000 1.000 1.000 0.999 0.998 1.000 0.999 0.999 0.999 0.997 0.999
> Max.    1.000 1.000 1.000 1.000 1.000 0.999 1.000 1.000 1.000 1.000 0.999 1.000
>         [,13] [,14] [,15]
> Min.    0.050 0.050 0.050
> 1st Qu. 0.838 0.885 0.936
> Median  0.970 0.980 0.988
> Mean    0.849 0.865 0.890
> 3rd Qu. 0.999 0.998 0.997
> Max.    1.000 0.999 0.999
summary(cypmatrix2fp)
> 
> This ahistorical lefkoMat object contains 15 matrices.
> 
> Each matrix is square with 54 rows and columns, and a total of 2916 elements.
> A total of 36165 survival transitions were estimated, with 2411 per matrix.
> A total of 720 fecundity transitions were estimated, with 48 per matrix.
> This lefkoMat object covers 1 population, 3 patches, and 5 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Observation estimated with 70 individuals and 303 individual transitions.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
> 
> 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.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050 0.050
> 1st Qu. 0.991 0.991 0.991 0.991 0.991 0.991 0.991 0.991 0.991 0.991 0.991 0.991
> Median  0.999 1.000 1.000 1.000 0.999 0.998 1.000 0.999 0.999 0.999 0.997 0.999
> Mean    0.916 0.918 0.917 0.917 0.914 0.915 0.917 0.917 0.916 0.914 0.913 0.916
> 3rd Qu. 1.000 1.000 1.000 1.000 1.000 0.999 1.000 1.000 1.000 0.999 0.999 1.000
> Max.    1.000 1.000 1.000 1.000 1.000 0.999 1.000 1.000 1.000 1.000 0.999 1.000
>         [,13] [,14] [,15]
> Min.    0.050 0.050 0.050
> 1st Qu. 0.991 0.991 0.991
> Median  0.999 0.998 0.997
> Mean    0.916 0.915 0.912
> 3rd Qu. 0.999 0.999 0.999
> Max.    1.000 0.999 0.999

We can see in the output that there are about 600 fewer non-zero survival terms per matrix in our new ahistorical MPM. This is because we have eliminated the return transitions from the reproductive to the non-reproductive stages.

Using stage groups provides a great deal of power in structuring your MPMs. Now let’s move on to the use of manually estimated vital rate models in function-based MPMs.

5.7 Environmental and individual covariates

Users may at times wish to develop MPMs conditioned on an extra variable, including environmental variables such as Spring precipitation, or a further individual state variable such as a frailty index. Package lefko3 can handle this easily with the modelsearch() function combined with the various function-based matrix estimators. Let’s see an example of how this might work.

For our example, let’s continue using the Cypripedium dataset with the life history model diagrammed in figure 5.1. We will see how total annual precipitation affects our vital rates, and use the resulting vital rates to create a new function-based MPM. First, we will need to add the proper variables to our original 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

summary(cypdata_env)
>     plantid       patch        X               Y              censor 
>  Min.   : 164.0   A:23   Min.   : 46.5   Min.   :-28.00   Min.   :1  
>  1st Qu.: 265.0   B:35   1st Qu.: 60.1   1st Qu.: 22.50   1st Qu.:1  
>  Median : 455.0   C:19   Median : 91.4   Median : 75.60   Median :1  
>  Mean   : 669.1          Mean   : 92.9   Mean   : 55.54   Mean   :1  
>  3rd Qu.: 829.0          3rd Qu.:142.2   3rd Qu.: 80.10   3rd Qu.:1  
>  Max.   :1560.0          Max.   :173.0   Max.   :142.40   Max.   :1  
>                                                                      
>     Inf2.04       Inf.04           Veg.04           Pod.04       Inf2.05       
>  Min.   :0    Min.   :0.0000   Min.   : 0.000   Min.   :0.0   Min.   :0.00000  
>  1st Qu.:0    1st Qu.:0.0000   1st Qu.: 1.000   1st Qu.:0.0   1st Qu.:0.00000  
>  Median :0    Median :0.0000   Median : 2.000   Median :0.0   Median :0.00000  
>  Mean   :0    Mean   :0.6923   Mean   : 2.923   Mean   :0.2   Mean   :0.04478  
>  3rd Qu.:0    3rd Qu.:1.0000   3rd Qu.: 4.000   3rd Qu.:0.0   3rd Qu.:0.00000  
>  Max.   :0    Max.   :8.0000   Max.   :12.000   Max.   :3.0   Max.   :1.00000  
>  NA's   :12   NA's   :12       NA's   :12       NA's   :12    NA's   :10       
>      Inf.05           Veg.05          Pod.05          Inf2.06  
>  Min.   : 0.000   Min.   :0.000   Min.   :0.0000   Min.   :0   
>  1st Qu.: 0.000   1st Qu.:1.000   1st Qu.:0.0000   1st Qu.:0   
>  Median : 0.000   Median :1.000   Median :0.0000   Median :0   
>  Mean   : 1.537   Mean   :2.134   Mean   :0.6567   Mean   :0   
>  3rd Qu.: 2.000   3rd Qu.:3.000   3rd Qu.:1.0000   3rd Qu.:0   
>  Max.   :18.000   Max.   :9.000   Max.   :7.0000   Max.   :0   
>  NA's   :10       NA's   :10      NA's   :10       NA's   :16  
>      Inf.06            Veg.06           Pod.06       Inf2.07       
>  Min.   : 0.0000   Min.   : 0.000   Min.   :0.0000   Mode:logical  
>  1st Qu.: 0.0000   1st Qu.: 1.000   1st Qu.:0.0000   NA's:77       
>  Median : 0.0000   Median : 2.000   Median :0.0000                 
>  Mean   : 0.9016   Mean   : 2.213   Mean   :0.3934                 
>  3rd Qu.: 1.0000   3rd Qu.: 3.000   3rd Qu.:0.0000                 
>  Max.   :18.0000   Max.   :13.000   Max.   :4.0000                 
>  NA's   :16        NA's   :16       NA's   :16                     
>      Inf.07           Veg.07           Pod.07          Inf2.08  
>  Min.   :0.0000   Min.   : 0.000   Min.   :0.0000   Min.   :0   
>  1st Qu.:0.0000   1st Qu.: 1.000   1st Qu.:0.0000   1st Qu.:0   
>  Median :0.0000   Median : 2.000   Median :0.0000   Median :0   
>  Mean   :0.6271   Mean   : 2.627   Mean   :0.0678   Mean   :0   
>  3rd Qu.:1.0000   3rd Qu.: 4.000   3rd Qu.:0.0000   3rd Qu.:0   
>  Max.   :7.0000   Max.   :13.000   Max.   :1.0000   Max.   :0   
>  NA's   :18       NA's   :18       NA's   :18       NA's   :24  
>      Inf.08            Veg.08          Pod.08          Inf2.09  
>  Min.   : 0.0000   Min.   : 0.00   Min.   :0.0000   Min.   :0   
>  1st Qu.: 0.0000   1st Qu.: 1.00   1st Qu.:0.0000   1st Qu.:0   
>  Median : 0.0000   Median : 2.00   Median :0.0000   Median :0   
>  Mean   : 0.8868   Mean   : 2.83   Mean   :0.1509   Mean   :0   
>  3rd Qu.: 1.0000   3rd Qu.: 4.00   3rd Qu.:0.0000   3rd Qu.:0   
>  Max.   :11.0000   Max.   :13.00   Max.   :2.0000   Max.   :0   
>  NA's   :24        NA's   :24      NA's   :24       NA's   :17  
>      Inf.09           Veg.09           Pod.09         prec.04    
>  Min.   : 0.000   Min.   : 0.000   Min.   :0.000   Min.   :92.2  
>  1st Qu.: 0.000   1st Qu.: 1.000   1st Qu.:0.000   1st Qu.:92.2  
>  Median : 1.000   Median : 1.000   Median :1.000   Median :92.2  
>  Mean   : 1.833   Mean   : 2.233   Mean   :1.133   Mean   :92.2  
>  3rd Qu.: 2.000   3rd Qu.: 3.000   3rd Qu.:1.000   3rd Qu.:92.2  
>  Max.   :11.000   Max.   :10.000   Max.   :8.000   Max.   :92.2  
>  NA's   :17       NA's   :17       NA's   :17                    
>     prec.05        prec.06      prec.07         prec.08         prec.09     
>  Min.   :57.6   Min.   :96   Min.   :109.8   Min.   :111.9   Min.   :106.8  
>  1st Qu.:57.6   1st Qu.:96   1st Qu.:109.8   1st Qu.:111.9   1st Qu.:106.8  
>  Median :57.6   Median :96   Median :109.8   Median :111.9   Median :106.8  
>  Mean   :57.6   Mean   :96   Mean   :109.8   Mean   :111.9   Mean   :106.8  
>  3rd Qu.:57.6   3rd Qu.:96   3rd Qu.:109.8   3rd Qu.:111.9   3rd Qu.:106.8  
>  Max.   :57.6   Max.   :96   Max.   :109.8   Max.   :111.9   Max.   :106.8  
> 

The dataset is still formatted for the most part using blocks of variables for each year, except for the precipitation variables, which are at the very end of the dataset. We will standardize the dataset using verticalize3() using blocks, but inputting all precipitation variables by name (users wishing to use the block feature for precipitation must reorder the columns in the data frame).

cypfb_env <- verticalize3(data = cypdata_env, noyears = 6, firstyear = 2004, 
  patchidcol = "patch", individcol = "plantid", blocksize = 4, 
  sizeacol = "Inf2.04", sizebcol = "Inf.04", sizeccol = "Veg.04", 
  repstracol = "Inf.04", repstrbcol = "Inf2.04", fecacol = "Pod.04",
  indcovacol = c("prec.04", "prec.05", "prec.06", "prec.07", "prec.08",
    "prec.09"), stageassign = cypframe_fb, stagesize = "sizeadded",
  NAas0 = TRUE, age_offset = 4)

summary(cypfb_env)
>      rowid       popid  patchid   individ              year2     
>  Min.   : 1.00   :320   A: 93   Length:320         Min.   :2004  
>  1st Qu.:21.00          B:154   Class :character   1st Qu.:2005  
>  Median :37.50          C: 73   Mode  :character   Median :2006  
>  Mean   :38.45                                     Mean   :2006  
>  3rd Qu.:56.00                                     3rd Qu.:2007  
>  Max.   :77.00                                     Max.   :2008  
>    firstseen       lastseen        obsage       obslifespan   
>  Min.   :2004   Min.   :2004   Min.   :5.000   Min.   :0.000  
>  1st Qu.:2004   1st Qu.:2009   1st Qu.:6.000   1st Qu.:5.000  
>  Median :2004   Median :2009   Median :7.000   Median :5.000  
>  Mean   :2004   Mean   :2009   Mean   :6.853   Mean   :4.556  
>  3rd Qu.:2004   3rd Qu.:2009   3rd Qu.:8.000   3rd Qu.:5.000  
>  Max.   :2008   Max.   :2009   Max.   :9.000   Max.   :5.000  
>      sizea1             sizeb1            sizec1       size1added    
>  Min.   :0.000000   Min.   : 0.0000   Min.   : 0.0   Min.   : 0.000  
>  1st Qu.:0.000000   1st Qu.: 0.0000   1st Qu.: 0.0   1st Qu.: 0.000  
>  Median :0.000000   Median : 0.0000   Median : 1.0   Median : 2.000  
>  Mean   :0.009375   Mean   : 0.7469   Mean   : 1.9   Mean   : 2.656  
>  3rd Qu.:0.000000   3rd Qu.: 1.0000   3rd Qu.: 3.0   3rd Qu.: 4.000  
>  Max.   :1.000000   Max.   :18.0000   Max.   :13.0   Max.   :21.000  
>     repstra1          repstrb1            feca1           indcova1     
>  Min.   : 0.0000   Min.   :0.000000   Min.   :0.0000   Min.   :  0.00  
>  1st Qu.: 0.0000   1st Qu.:0.000000   1st Qu.:0.0000   1st Qu.: 57.60  
>  Median : 0.0000   Median :0.000000   Median :0.0000   Median : 92.20  
>  Mean   : 0.7469   Mean   :0.009375   Mean   :0.2656   Mean   : 70.64  
>  3rd Qu.: 1.0000   3rd Qu.:0.000000   3rd Qu.:0.0000   3rd Qu.: 96.00  
>  Max.   :18.0000   Max.   :1.000000   Max.   :7.0000   Max.   :109.80  
>    juvgiven1   obsstatus1       repstatus1       fecstatus1    
>  Min.   :0   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
>  1st Qu.:0   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
>  Median :0   Median :1.0000   Median :0.0000   Median :0.0000  
>  Mean   :0   Mean   :0.7469   Mean   :0.2875   Mean   :0.1344  
>  3rd Qu.:0   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:0.0000  
>  Max.   :0   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
>    matstatus1         alive1          stage1           stage1index   
>  Min.   :0.0000   Min.   :0.0000   Length:320         Min.   : 0.00  
>  1st Qu.:1.0000   1st Qu.:1.0000   Class :character   1st Qu.: 6.00  
>  Median :1.0000   Median :1.0000   Mode  :character   Median : 8.00  
>  Mean   :0.7688   Mean   :0.7688                      Mean   :14.17  
>  3rd Qu.:1.0000   3rd Qu.:1.0000                      3rd Qu.:31.00  
>  Max.   :1.0000   Max.   :1.0000                      Max.   :51.00  
>      sizea2             sizeb2            sizec2         size2added    
>  Min.   :0.000000   Min.   : 0.0000   Min.   : 0.000   Min.   : 0.000  
>  1st Qu.:0.000000   1st Qu.: 0.0000   1st Qu.: 1.000   1st Qu.: 1.000  
>  Median :0.000000   Median : 0.0000   Median : 2.000   Median : 2.000  
>  Mean   :0.009375   Mean   : 0.8969   Mean   : 2.416   Mean   : 3.322  
>  3rd Qu.:0.000000   3rd Qu.: 1.0000   3rd Qu.: 3.000   3rd Qu.: 4.000  
>  Max.   :1.000000   Max.   :18.0000   Max.   :13.000   Max.   :24.000  
>     repstra2          repstrb2            feca2           indcova2     
>  Min.   : 0.0000   Min.   :0.000000   Min.   :0.0000   Min.   : 57.60  
>  1st Qu.: 0.0000   1st Qu.:0.000000   1st Qu.:0.0000   1st Qu.: 92.20  
>  Median : 0.0000   Median :0.000000   Median :0.0000   Median : 96.00  
>  Mean   : 0.8969   Mean   :0.009375   Mean   :0.2906   Mean   : 92.77  
>  3rd Qu.: 1.0000   3rd Qu.:0.000000   3rd Qu.:0.0000   3rd Qu.:109.80  
>  Max.   :18.0000   Max.   :1.000000   Max.   :7.0000   Max.   :111.90  
>    juvgiven2   obsstatus2       repstatus2       fecstatus2       matstatus2
>  Min.   :0   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :1   
>  1st Qu.:0   1st Qu.:1.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:1   
>  Median :0   Median :1.0000   Median :0.0000   Median :0.0000   Median :1   
>  Mean   :0   Mean   :0.9531   Mean   :0.3688   Mean   :0.1562   Mean   :1   
>  3rd Qu.:0   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:1   
>  Max.   :0   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1   
>      alive2     stage2           stage2index        sizea3        
>  Min.   :1   Length:320         Min.   : 6.00   Min.   :0.000000  
>  1st Qu.:1   Class :character   1st Qu.: 7.00   1st Qu.:0.000000  
>  Median :1   Mode  :character   Median :10.00   Median :0.000000  
>  Mean   :1                      Mean   :18.17   Mean   :0.009375  
>  3rd Qu.:1                      3rd Qu.:32.00   3rd Qu.:0.000000  
>  Max.   :1                      Max.   :54.00   Max.   :1.000000  
>      sizeb3           sizec3         size3added        repstra3     
>  Min.   : 0.000   Min.   : 0.000   Min.   : 0.000   Min.   : 0.000  
>  1st Qu.: 0.000   1st Qu.: 1.000   1st Qu.: 1.000   1st Qu.: 0.000  
>  Median : 0.000   Median : 1.000   Median : 2.000   Median : 0.000  
>  Mean   : 1.069   Mean   : 2.209   Mean   : 3.288   Mean   : 1.069  
>  3rd Qu.: 1.000   3rd Qu.: 3.000   3rd Qu.: 4.000   3rd Qu.: 1.000  
>  Max.   :18.000   Max.   :13.000   Max.   :24.000   Max.   :18.000  
>     repstrb3            feca3           indcova3       juvgiven3   obsstatus3 
>  Min.   :0.000000   Min.   :0.0000   Min.   : 57.6   Min.   :0   Min.   :0.0  
>  1st Qu.:0.000000   1st Qu.:0.0000   1st Qu.: 96.0   1st Qu.:0   1st Qu.:1.0  
>  Median :0.000000   Median :0.0000   Median :106.8   Median :0   Median :1.0  
>  Mean   :0.009375   Mean   :0.4562   Mean   : 96.1   Mean   :0   Mean   :0.9  
>  3rd Qu.:0.000000   3rd Qu.:0.0000   3rd Qu.:109.8   3rd Qu.:0   3rd Qu.:1.0  
>  Max.   :1.000000   Max.   :8.0000   Max.   :111.9   Max.   :0   Max.   :1.0  
>    repstatus3    fecstatus3       matstatus3     alive3      
>  Min.   :0.0   Min.   :0.0000   Min.   :1    Min.   :0.0000  
>  1st Qu.:0.0   1st Qu.:0.0000   1st Qu.:1    1st Qu.:1.0000  
>  Median :0.0   Median :0.0000   Median :1    Median :1.0000  
>  Mean   :0.4   Mean   :0.2219   Mean   :1    Mean   :0.9469  
>  3rd Qu.:1.0   3rd Qu.:0.0000   3rd Qu.:1    3rd Qu.:1.0000  
>  Max.   :1.0   Max.   :1.0000   Max.   :1    Max.   :1.0000  
>     stage3           stage3index   
>  Length:320         Min.   : 0.00  
>  Class :character   1st Qu.: 7.00  
>  Mode  :character   Median :10.00  
>                     Mean   :18.57  
>                     3rd Qu.:33.00  
>                     Max.   :54.00

Now we will proceed to 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 covariates. 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 instead). To prevent this function from taking too long to run, we will only do the ahistorical set, we will drop the patch term, and we will set suite = "main" to prevent the inclusion of interaction terms.

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

And now the summary.

summary(cypmodels2_env)
> This LefkoMod object includes 5 linear models.
> Best-fit model criterion used: aicc&k
> 
> 
> 
> Survival model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
>   Approximation) [glmerMod]
>  Family: binomial  ( logit )
> Formula: alive3 ~ indcova2 + (1 | year2) + (1 | individ)
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
>  87.6075 102.6808 -39.8038  79.6075      316 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 552.9   
>  year2   (Intercept) 394.6   
> Number of obs: 320, groups:  individ, 74; year2, 5
> Fixed Effects:
> (Intercept)     indcova2  
>     476.174       -3.278  
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 
> 
> 
> 
> Observation model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
>   Approximation) [glmerMod]
>  Family: binomial  ( logit )
> Formula: obsstatus3 ~ size2added + (1 | year2) + (1 | individ)
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
> 118.2567 133.1117 -55.1284 110.2567      299 
> Random effects:
>  Groups  Name        Std.Dev. 
>  individ (Intercept) 1.078e-05
>  year2   (Intercept) 8.776e-01
> Number of obs: 303, groups:  individ, 70; year2, 5
> Fixed Effects:
> (Intercept)   size2added  
>      2.4904       0.3134  
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 
> 
> 
> 
> Size model:
> Formula:          
> size3added ~ indcova2 + size2added + (1 | year2) + (1 | individ)
> Data: subdata
>       AIC       BIC    logLik  df.resid 
> 1001.3778 1023.3555 -494.6889       282 
> Random-effects (co)variances:
> 
> Conditional model:
>  Groups  Name        Std.Dev.
>  year2   (Intercept) 0.06554 
>  individ (Intercept) 0.96060 
> 
> Number of obs: 288 / Conditional model: year2, 5; individ, 70
> 
> Dispersion parameter for truncated_nbinom2 family (): 1.87e+12 
> 
> Fixed Effects:
> 
> Conditional model:
> (Intercept)     indcova2   size2added  
>    0.070927     0.004918     0.023979  
> 
> 
> 
> Secondary size model:
> [1] 1
> 
> 
> 
> Tertiary size model:
> [1] 1
> 
> 
> 
> Reproductive status model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
>   Approximation) [glmerMod]
>  Family: binomial  ( logit )
> Formula: repstatus3 ~ indcova2 + repstatus2 + size2added + (1 | year2) +  
>     (1 | individ)
>    Data: subdata
>       AIC       BIC    logLik  deviance  df.resid 
>  330.0418  352.0196 -159.0209  318.0418       282 
> Random effects:
>  Groups  Name        Std.Dev. 
>  individ (Intercept) 0.0002695
>  year2   (Intercept) 0.2479120
> Number of obs: 288, groups:  individ, 70; year2, 5
> Fixed Effects:
> (Intercept)     indcova2   repstatus2   size2added  
>    -4.23766      0.02954      1.71091      0.17187  
> 
> 
> 
> Fecundity model:
> Formula:          feca2 ~ size2added + (1 | year2) + (1 | individ)
> Zero inflation:         ~size2added + (1 | year2) + (1 | individ)
> Data: subdata
>       AIC       BIC    logLik  df.resid 
>  248.8609  271.0264 -116.4305       110 
> Random-effects (co)variances:
> 
> Conditional model:
>  Groups  Name        Std.Dev.
>  year2   (Intercept) 0.5760  
>  individ (Intercept) 0.1639  
> 
> Zero-inflation model:
>  Groups  Name        Std.Dev. 
>  year2   (Intercept) 1.642e-06
>  individ (Intercept) 3.089e-04
> 
> Number of obs: 118 / Conditional model: year2, 5; individ, 51 / Zero-inflation model: year2, 5; individ, 51
> 
> Fixed Effects:
> 
> Conditional model:
> (Intercept)   size2added  
>    -0.54014      0.06174  
> 
> Zero-inflation model:
> (Intercept)   size2added  
>       3.865       -1.574  
> 
> 
> Juvenile survival model:
> [1] 1
> 
> 
> 
> Juvenile observation model:
> [1] 1
> 
> 
> 
> Juvenile size model:
> [1] 1
> 
> 
> 
> Juvenile secondary size model:
> [1] 1
> 
> 
> 
> Juvenile tertiary size model:
> [1] 1
> 
> 
> 
> Juvenile reproduction model:
> [1] 1
> 
> 
> 
> Juvenile maturity model:
> [1] 1
> 
> 
> 
> 
> 
> Number of models in survival table: 8
> 
> Number of models in observation table: 8
> 
> Number of models in size table: 8
> 
> Number of models in secondary size table: 1
> 
> Number of models in tertiary size table: 1
> 
> Number of models in reproduction status table: 8
> 
> Number of models in fecundity table: 61
> 
> Number of models in juvenile survival table: 1
> 
> Number of models in juvenile observation table: 1
> 
> Number of models in juvenile size table: 1
> 
> Number of models in juvenile secondary size table: 1
> 
> Number of models in juvenile tertiary size table: 1
> 
> Number of models in juvenile reproduction table: 1
> 
> Number of models in juvenile maturity table: 1
> 
> 
> 
> 
> 
> General model parameter names (column 1), and 
> specific names used in these models (column 2): 
>                       parameter_names mainparams
> 1                              time t      year2
> 2                          individual    individ
> 3                               patch      patch
> 4                   alive in time t+1      surv3
> 5                observed in time t+1       obs3
> 6                   sizea in time t+1      size3
> 7                   sizeb in time t+1     sizeb3
> 8                   sizec in time t+1     sizec3
> 9     reproductive status in time t+1     repst3
> 10              fecundity in time t+1       fec3
> 11                fecundity in time t       fec2
> 12                    sizea in time t      size2
> 13                  sizea in time t-1      size1
> 14                    sizeb in time t     sizeb2
> 15                  sizeb in time t-1     sizeb1
> 16                    sizec in time t     sizec2
> 17                  sizec in time t-1     sizec1
> 18      reproductive status in time t     repst2
> 19    reproductive status in time t-1     repst1
> 20        maturity status in time t+1     matst3
> 21          maturity status in time t     matst2
> 22                      age in time t        age
> 23                  density in time t    density
> 24   individual covariate a in time t   indcova2
> 25 individual covariate a in time t-1   indcova1
> 26   individual covariate b in time t   indcovb2
> 27 individual covariate b in time t-1   indcovb1
> 28   individual covariate c in time t   indcovc2
> 29 individual covariate c in time t-1   indcovc1
> 30              stage group in time t     group2
> 31            stage group in time t-1     group1
> 
> 
> 
> 
> 
> Quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Survival accuracy is 1.
> Observation estimated with 70 individuals and 303 individual transitions.
> Observation accuracy is 0.95.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Primary size pseudo R-squared is 0.011.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Reproductive status accuracy is 0.719.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Fecundity pseudo R-squared is 0.32.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity probability not estimated.

Our individual covariate, total annual precipitation, is included in several best-fit vital rate models. The resulting lefkoMod object needs to be used with a set value of the individual covariate, or a set of such values corresponding to each time. Let’s try building an MPM with a single covariate value, which will be equal to 100cm of annual precipitation.

cypmatrix2f_env1 <- flefko2(stageframe = cypframe_fb, data = cypfb_env,
  supplement = cypsupp2_fb, modelsuite = cypmodels2_env, inda = 100)

summary(cypmatrix2f_env1)
> 
> This ahistorical lefkoMat object contains 5 matrices.
> 
> Each matrix is square with 54 rows and columns, and a total of 2916 elements.
> A total of 12055 survival transitions were estimated, with 2411 per matrix.
> A total of 240 fecundity transitions were estimated, with 48 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 5 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Observation estimated with 70 individuals and 303 individual transitions.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
> 
> Survival probability sum check (each matrix represented by column in order):
>          [,1]     [,2]  [,3]  [,4]  [,5]
> Min.    0.050 2.27e-22 0.050 0.050 0.050
> 1st Qu. 0.999 2.27e-22 0.999 0.999 0.999
> Median  1.000 2.27e-22 1.000 1.000 1.000
> Mean    0.923 8.89e-03 0.924 0.925 0.923
> 3rd Qu. 1.000 2.27e-22 1.000 1.000 1.000
> Max.    1.000 1.80e-01 1.000 1.000 1.000

Now let’s build a version with a different, set precipitation value for each year (time t only, meaning that there are five values in a six year dataset) based on the actual values.

precip_values <- c(92.2, 57.6, 96, 109.8, 111.9)

cypmatrix2f_env2 <- flefko2(stageframe = cypframe_fb, data = cypfb_env,
  supplement = cypsupp2_fb, modelsuite = cypmodels2_env, inda = precip_values)

summary(cypmatrix2f_env2)
> 
> This ahistorical lefkoMat object contains 5 matrices.
> 
> Each matrix is square with 54 rows and columns, and a total of 2916 elements.
> A total of 12055 survival transitions were estimated, with 2411 per matrix.
> A total of 240 fecundity transitions were estimated, with 48 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 5 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 74 individuals and 320 individual transitions.
> Observation estimated with 70 individuals and 303 individual transitions.
> Primary size estimated with 70 individuals and 288 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 70 individuals and 288 individual transitions.
> Fecundity estimated with 51 individuals and 118 individual transitions.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
> 
> Survival probability sum check (each matrix represented by column in order):
>          [,1]  [,2]  [,3]  [,4]  [,5]
> Min.    0.050 0.050 0.050 0.050 0.050
> 1st Qu. 0.999 1.000 1.000 0.999 0.999
> Median  1.000 1.000 1.000 1.000 1.000
> Mean    0.924 0.926 0.925 0.924 0.922
> 3rd Qu. 1.000 1.000 1.000 1.000 1.000
> Max.    1.000 1.000 1.000 1.000 1.000

Voilà! We’re done! Now we can move on to analysis.

5.8 Creating function-based MPMs without modelsearch()

Advanced users of lefko3 may wish to create their own models without the package’s automated model selection function. In that case, lefko3’s function-based MPM creation functions flefko2(), flefko3(), and aflefko2() can accommodate single models developed using a series of function covering generalized linear models and generalized linear mixed models in base R and a series of other packages. Currently, users can develop these models using the following functions:

  1. lm() and glm() from the core R package stats, with the latter limited to the following response distribution families: binomial, gaussian, Gamma, and poisson
  2. glm.nb() from the core R package MASS
  3. vglm() from package VGAM (Yee & Wild 1996; Yee 2015), with response distributions limited to pospoisson and posnegbinomial
  4. zeroinfl() from package pscl (Zeileis et al. 2008), with response distributions limited to poisson and negbin
  5. lmer and glmer from package lme4 (Bates et al. 2015), with the latter limited to the following response distribution families: binomial, gaussian, Gamma, and poisson
  6. glmmTMB from package glmmTMB (Brooks et al. 2017), with response distribution families limited to nbinom2, truncated_poisson, and truncated_nbinom2

Generally, survival probability, observation probability, and reproduction probability must be modeled using either function glm() or glmer(), in either case using the binomial response. Size and fecundity may be modeled as Gaussian, gamma, Poisson, or negative binomial, with the Poisson and negative binomial capable of being zero-inflated or zero-truncated. The binomial response is not allowed for size and fecundity.

Developing models manually requires a few key steps. First, the core dataset must be created and serve as the main dataset for modeling. Second, data subsets must be developed for conditional rates and probabilities. Third, the models may be developed. Fourth, a paramnames object must be made that identifies the names of variables from the dataset used in modeling. Finally, we can create our MPM using our manual models as inputs.

Let’s try an example of a manual modeling exercise. Let’s try to create a historical MPM for the Cypripedium candidum dataset, but forcing each model to include size and reproductive status in times t and t-1, and occasion in time t and individual identity. So, these will be mixed models. Let’s start off by looking over our core hfv format dataset.

summary(cypfb_v1)
>      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          juvgiven1
>  Min.   : 0.0000   Min.   :0.000000   Min.   :0.0000   Min.   :0  
>  1st Qu.: 0.0000   1st Qu.:0.000000   1st Qu.:0.0000   1st Qu.:0  
>  Median : 0.0000   Median :0.000000   Median :0.0000   Median :0  
>  Mean   : 0.7469   Mean   :0.009375   Mean   :0.2656   Mean   :0  
>  3rd Qu.: 1.0000   3rd Qu.:0.000000   3rd Qu.:0.0000   3rd Qu.:0  
>  Max.   :18.0000   Max.   :1.000000   Max.   :7.0000   Max.   :0  
>    obsstatus1       repstatus1       fecstatus1       matstatus1    
>  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
>  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:1.0000  
>  Median :1.0000   Median :0.0000   Median :0.0000   Median :1.0000  
>  Mean   :0.7469   Mean   :0.2875   Mean   :0.1344   Mean   :0.7688  
>  3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:1.0000  
>  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
>      alive1          stage1           stage1index        sizea2        
>  Min.   :0.0000   Length:320         Min.   : 0.00   Min.   :0.000000  
>  1st Qu.:1.0000   Class :character   1st Qu.: 6.00   1st Qu.:0.000000  
>  Median :1.0000   Mode  :character   Median : 8.00   Median :0.000000  
>  Mean   :0.7688                      Mean   :14.17   Mean   :0.009375  
>  3rd Qu.:1.0000                      3rd Qu.:31.00   3rd Qu.:0.000000  
>  Max.   :1.0000                      Max.   :51.00   Max.   :1.000000  
>      sizeb2            sizec2         size2added        repstra2      
>  Min.   : 0.0000   Min.   : 0.000   Min.   : 0.000   Min.   : 0.0000  
>  1st Qu.: 0.0000   1st Qu.: 1.000   1st Qu.: 1.000   1st Qu.: 0.0000  
>  Median : 0.0000   Median : 2.000   Median : 2.000   Median : 0.0000  
>  Mean   : 0.8969   Mean   : 2.416   Mean   : 3.322   Mean   : 0.8969  
>  3rd Qu.: 1.0000   3rd Qu.: 3.000   3rd Qu.: 4.000   3rd Qu.: 1.0000  
>  Max.   :18.0000   Max.   :13.000   Max.   :24.000   Max.   :18.0000  
>     repstrb2            feca2          juvgiven2   obsstatus2    
>  Min.   :0.000000   Min.   :0.0000   Min.   :0   Min.   :0.0000  
>  1st Qu.:0.000000   1st Qu.:0.0000   1st Qu.:0   1st Qu.:1.0000  
>  Median :0.000000   Median :0.0000   Median :0   Median :1.0000  
>  Mean   :0.009375   Mean   :0.2906   Mean   :0   Mean   :0.9531  
>  3rd Qu.:0.000000   3rd Qu.:0.0000   3rd Qu.:0   3rd Qu.:1.0000  
>  Max.   :1.000000   Max.   :7.0000   Max.   :0   Max.   :1.0000  
>    repstatus2       fecstatus2       matstatus2     alive2     stage2         
>  Min.   :0.0000   Min.   :0.0000   Min.   :1    Min.   :1   Length:320        
>  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:1    1st Qu.:1   Class :character  
>  Median :0.0000   Median :0.0000   Median :1    Median :1   Mode  :character  
>  Mean   :0.3688   Mean   :0.1562   Mean   :1    Mean   :1                     
>  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:1    3rd Qu.:1                     
>  Max.   :1.0000   Max.   :1.0000   Max.   :1    Max.   :1                     
>   stage2index        sizea3             sizeb3           sizec3      
>  Min.   : 6.00   Min.   :0.000000   Min.   : 0.000   Min.   : 0.000  
>  1st Qu.: 7.00   1st Qu.:0.000000   1st Qu.: 0.000   1st Qu.: 1.000  
>  Median :10.00   Median :0.000000   Median : 0.000   Median : 1.000  
>  Mean   :18.17   Mean   :0.009375   Mean   : 1.069   Mean   : 2.209  
>  3rd Qu.:32.00   3rd Qu.:0.000000   3rd Qu.: 1.000   3rd Qu.: 3.000  
>  Max.   :54.00   Max.   :1.000000   Max.   :18.000   Max.   :13.000  
>    size3added        repstra3         repstrb3            feca3       
>  Min.   : 0.000   Min.   : 0.000   Min.   :0.000000   Min.   :0.0000  
>  1st Qu.: 1.000   1st Qu.: 0.000   1st Qu.:0.000000   1st Qu.:0.0000  
>  Median : 2.000   Median : 0.000   Median :0.000000   Median :0.0000  
>  Mean   : 3.288   Mean   : 1.069   Mean   :0.009375   Mean   :0.4562  
>  3rd Qu.: 4.000   3rd Qu.: 1.000   3rd Qu.:0.000000   3rd Qu.:0.0000  
>  Max.   :24.000   Max.   :18.000   Max.   :1.000000   Max.   :8.0000  
>    juvgiven3   obsstatus3    repstatus3    fecstatus3       matstatus3
>  Min.   :0   Min.   :0.0   Min.   :0.0   Min.   :0.0000   Min.   :1   
>  1st Qu.:0   1st Qu.:1.0   1st Qu.:0.0   1st Qu.:0.0000   1st Qu.:1   
>  Median :0   Median :1.0   Median :0.0   Median :0.0000   Median :1   
>  Mean   :0   Mean   :0.9   Mean   :0.4   Mean   :0.2219   Mean   :1   
>  3rd Qu.:0   3rd Qu.:1.0   3rd Qu.:1.0   3rd Qu.:0.0000   3rd Qu.:1   
>  Max.   :0   Max.   :1.0   Max.   :1.0   Max.   :1.0000   Max.   :1   
>      alive3          stage3           stage3index   
>  Min.   :0.0000   Length:320         Min.   : 0.00  
>  1st Qu.:1.0000   Class :character   1st Qu.: 7.00  
>  Median :1.0000   Mode  :character   Median :10.00  
>  Mean   :0.9469                      Mean   :18.57  
>  3rd Qu.:1.0000                      3rd Qu.:33.00  
>  Max.   :1.0000                      Max.   :54.00

We will be estimating survival, and for that we need to assess whether individuals alive in time t are still alive in time t+1. Here, variables with a 2 at the end of their names represent variables coding for status in time t, while 1 codes for time t-1 and 3 codes for time t+1. The summary states that alive2 equals 1 in all cases, meaning that all rows in this dataset are for individuals alive in time t (this is the default output from the verticalize3() and historicalize3() functions). So, we can use this dataset to test for survival. Let’s create our survival dataset by copying this dataset over to something named more appropriately.

surv_data <- cypfb_v1

Observation status is assessed in time t+1, and needs to be assessed using individuals that are alive in that time. So, let’s create a new subset of these data to assess observation status.

obs_data <- subset(surv_data, alive3 == 1)

Primary size also needs to be assessed using individuals alive in time t+1. However, size can only be estimated for individuals that are observed. Hence, we need to subset obs_data to only those individuals observable in time t+1, as below.

size_data <- subset(obs_data, obsstatus3 == 1)

Reproduction status also requires that individuals are observed. Since this is the only requirement, we can simply copy size_data to a more appropriately named data frame for that model estimation.

repr_data <- size_data

Finally, there is the issue of fecundity. This rate takes some thought to handle properly, because it may be estimated differently in different situations, and hence its relationships to the variables in the dataset may change. In our case, we have utilized a pre-breeding model and use the number of fruits as a proxy for fecundity (the real fecundity rate needs to be multiplied by the number of seeds, the survival of seeds to the next time, and the germination rate, but this product can be handled through the supplement table). Fruits are produced in time t for germination in time t+1, and reproductive status is determined by the presence of flowers, so we need a subset of individuals alive in time t that also have repstatus2 == 1, as below.

fec_data <- subset(surv_data, repstatus2 == 1)

We have now built all of our data frames, and can proceed to develop our models. Here, each model is listed in a single block of code. We will rely on the glmer() function from package lme4, and the glmmTMB() function from package glmmTMB. The response distributions for survival, observation status, and reproduction status will be the binomial distribution. The response distributions for size and fecundity will be the zero-truncated negative binomial and the zero-inflated poisson, respectively.

library(lme4)
library(glmmTMB)

surv_mod <- glmer(alive3 ~ size2added + size1added + repstatus2 + repstatus1 +
  (1 | year2) + (1 | individ), data = surv_data, family = binomial)
> boundary (singular) fit: see help('isSingular')
obs_mod <- glmer(obsstatus3 ~ size2added + size1added + repstatus2 +
  repstatus1 + (1 | year2) + (1 | individ), data = obs_data, family = binomial)
size_mod <- glmmTMB(size3added ~ size2added + size1added + repstatus2 +
  repstatus1 + (1 | year2) + (1 | individ), data = size_data,
  family = truncated_nbinom2)
repr_mod <- glmer(repstatus3 ~ size2added + size1added + repstatus2 +
  repstatus1 + (1 | year2) + (1 | individ), data = repr_data, family = binomial)
> boundary (singular) fit: see help('isSingular')
fec_mod <- glmmTMB(feca3 ~ size2added + size1added + repstatus1 + (1 | year2) +
  (1 | individ), zi = ~., data = fec_data, family = poisson)
> Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
> NaN function evaluation

> Warning in (function (start, objective, gradient = NULL, hessian = NULL, : NA/
> NaN function evaluation
> Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
> Hessian matrix. See vignette('troubleshooting')
> Warning in fitTMB(TMBStruc): Model convergence problem; false convergence (8).
> See vignette('troubleshooting')

summary(surv_mod)
> Generalized linear mixed model fit by maximum likelihood (Laplace
>   Approximation) [glmerMod]
>  Family: binomial  ( logit )
> Formula: alive3 ~ size2added + size1added + repstatus2 + repstatus1 +  
>     (1 | year2) + (1 | individ)
>    Data: surv_data
> 
>      AIC      BIC   logLik deviance df.resid 
>    133.2    159.5    -59.6    119.2      313 
> 
> Scaled residuals: 
>     Min      1Q  Median      3Q     Max 
> -8.6555  0.0883  0.1947  0.3066  0.4706 
> 
> Random effects:
>  Groups  Name        Variance  Std.Dev.
>  individ (Intercept) 1.149e-13 3.39e-07
>  year2   (Intercept) 0.000e+00 0.00e+00
> Number of obs: 320, groups:  individ, 74; year2, 5
> 
> Fixed effects:
>             Estimate Std. Error z value Pr(>|z|)    
> (Intercept)  1.50738    0.45233   3.332 0.000861 ***
> size2added   0.48153    0.24076   2.000 0.045493 *  
> size1added   0.18764    0.20271   0.926 0.354645    
> repstatus2   0.42624    0.68931   0.618 0.536337    
> repstatus1  -0.05451    0.71646  -0.076 0.939349    
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Correlation of Fixed Effects:
>            (Intr) sz2ddd sz1ddd rpstt2
> size2added -0.592                     
> size1added -0.358 -0.174              
> repstatus2 -0.080 -0.222  0.042       
> repstatus1 -0.079  0.046 -0.318 -0.211
> optimizer (Nelder_Mead) convergence code: 0 (OK)
> boundary (singular) fit: see help('isSingular')
summary(obs_mod)
> Generalized linear mixed model fit by maximum likelihood (Laplace
>   Approximation) [glmerMod]
>  Family: binomial  ( logit )
> Formula: obsstatus3 ~ size2added + size1added + repstatus2 + repstatus1 +  
>     (1 | year2) + (1 | individ)
>    Data: obs_data
> 
>      AIC      BIC   logLik deviance df.resid 
>    122.3    148.3    -54.2    108.3      296 
> 
> Scaled residuals: 
>     Min      1Q  Median      3Q     Max 
> -6.2622  0.1121  0.1716  0.2435  0.6150 
> 
> Random effects:
>  Groups  Name        Variance  Std.Dev. 
>  individ (Intercept) 4.153e-07 0.0006444
>  year2   (Intercept) 7.302e-01 0.8545190
> Number of obs: 303, groups:  individ, 70; year2, 5
> 
> Fixed effects:
>             Estimate Std. Error z value Pr(>|z|)    
> (Intercept)   2.4945     0.6046   4.126 3.69e-05 ***
> size2added    0.3837     0.1979   1.939   0.0525 .  
> size1added   -0.1712     0.1507  -1.136   0.2559    
> repstatus2    0.4567     0.7075   0.645   0.5186    
> repstatus1    0.4619     0.7343   0.629   0.5293    
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Correlation of Fixed Effects:
>            (Intr) sz2ddd sz1ddd rpstt2
> size2added -0.257                     
> size1added -0.149 -0.569              
> repstatus2 -0.090 -0.111 -0.051       
> repstatus1 -0.043  0.026 -0.199 -0.276
summary(size_mod)
>  Family: truncated_nbinom2  ( log )
> Formula:          
> size3added ~ size2added + size1added + repstatus2 + repstatus1 +  
>     (1 | year2) + (1 | individ)
> Data: size_data
> 
>      AIC      BIC   logLik deviance df.resid 
>   1013.0   1042.3   -498.5    997.0      280 
> 
> Random effects:
> 
> Conditional model:
>  Groups  Name        Variance Std.Dev.
>  year2   (Intercept) 0.01378  0.1174  
>  individ (Intercept) 0.93532  0.9671  
> Number of obs: 288, groups:  year2, 5; individ, 70
> 
> Dispersion parameter for truncated_nbinom2 family (): 1.68e+07 
> 
> Conditional model:
>             Estimate Std. Error z value Pr(>|z|)   
> (Intercept)  0.51079    0.15819   3.229  0.00124 **
> size2added   0.02195    0.01362   1.611  0.10709   
> size1added  -0.00599    0.01159  -0.517  0.60540   
> repstatus2   0.06095    0.10005   0.609  0.54240   
> repstatus1   0.04482    0.10744   0.417  0.67657   
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(repr_mod)
> Generalized linear mixed model fit by maximum likelihood (Laplace
>   Approximation) [glmerMod]
>  Family: binomial  ( logit )
> Formula: repstatus3 ~ size2added + size1added + repstatus2 + repstatus1 +  
>     (1 | year2) + (1 | individ)
>    Data: repr_data
> 
>      AIC      BIC   logLik deviance df.resid 
>    334.4    360.1   -160.2    320.4      281 
> 
> Scaled residuals: 
>     Min      1Q  Median      3Q     Max 
> -2.1990 -0.6565 -0.3669  0.6895  2.8967 
> 
> Random effects:
>  Groups  Name        Variance Std.Dev.
>  individ (Intercept) 0.0000   0.0000  
>  year2   (Intercept) 0.5097   0.7139  
> Number of obs: 288, groups:  individ, 70; year2, 5
> 
> Fixed effects:
>             Estimate Std. Error z value Pr(>|z|)    
> (Intercept) -1.60119    0.40839  -3.921 8.83e-05 ***
> size2added   0.12168    0.06222   1.956   0.0505 .  
> size1added   0.07880    0.07391   1.066   0.2863    
> repstatus2   1.55291    0.31202   4.977 6.46e-07 ***
> repstatus1   0.43554    0.36684   1.187   0.2351    
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Correlation of Fixed Effects:
>            (Intr) sz2ddd sz1ddd rpstt2
> size2added -0.200                     
> size1added -0.186 -0.520              
> repstatus2 -0.193 -0.131  0.066       
> repstatus1 -0.110 -0.011 -0.172 -0.272
> optimizer (Nelder_Mead) convergence code: 0 (OK)
> boundary (singular) fit: see help('isSingular')
summary(fec_mod)
> Warning in sqrt(diag(vcovs)): NaNs produced
> Warning in sqrt(diag(vcovs)): NaNs produced
>  Family: poisson  ( log )
> Formula:          feca3 ~ size2added + size1added + repstatus1 + (1 | year2) +  
>     (1 | individ)
> Zero inflation:         ~.
> Data: fec_data
> 
>      AIC      BIC   logLik deviance df.resid 
>       NA       NA       NA       NA      106 
> 
> Random effects:
> 
> Conditional model:
>  Groups  Name        Variance Std.Dev.
>  year2   (Intercept) 0.6514   0.8071  
>  individ (Intercept) 0.7884   0.8879  
> Number of obs: 118, groups:  year2, 5; individ, 51
> 
> Zero-inflation model:
>  Groups  Name        Variance Std.Dev.
>  year2   (Intercept) 0.6388   0.7993  
>  individ (Intercept) 7.3858   2.7177  
> Number of obs: 118, groups:  year2, 5; individ, 51
> 
> Conditional model:
>             Estimate Std. Error z value Pr(>|z|)   
> (Intercept) -0.79519    0.26878  -2.958  0.00309 **
> size2added   0.05039    0.03558   1.416  0.15674   
> size1added   0.01121    0.04225   0.265  0.79080   
> repstatus1  -0.40228        NaN     NaN      NaN   
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Zero-inflation model:
>             Estimate Std. Error z value Pr(>|z|)    
> (Intercept)  3.48455    0.34185   10.19   <2e-16 ***
> size2added  -2.10471    0.07893  -26.66   <2e-16 ***
> size1added   0.71915        NaN     NaN      NaN    
> repstatus1  -5.15997        NaN     NaN      NaN    
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The last thing that we need to worry about is to create a paramnames object, which is a data frame that contains all of the information necessary for lefko3 to translate these models into MPMs. This should be composed of three columns: a column describing each parameter, a column giving the name of each parameter from lefko3’s parameter estimating protocols, and a column giving the name of each parameter within the modeling exercise. We can get the skeleton for this with the create_pm() function, as below.

our_pm <- create_pm()
our_pm
>                       parameter_names mainparams modelparams
> 1                              time t      year2        none
> 2                          individual    individ        none
> 3                               patch      patch        none
> 4                   alive in time t+1      surv3        none
> 5                observed in time t+1       obs3        none
> 6                   sizea in time t+1      size3        none
> 7                   sizeb in time t+1     sizeb3        none
> 8                   sizec in time t+1     sizec3        none
> 9     reproductive status in time t+1     repst3        none
> 10              fecundity in time t+1       fec3        none
> 11                fecundity in time t       fec2        none
> 12                    sizea in time t      size2        none
> 13                  sizea in time t-1      size1        none
> 14                    sizeb in time t     sizeb2        none
> 15                  sizeb in time t-1     sizeb1        none
> 16                    sizec in time t     sizec2        none
> 17                  sizec in time t-1     sizec1        none
> 18      reproductive status in time t     repst2        none
> 19    reproductive status in time t-1     repst1        none
> 20        maturity status in time t+1     matst3        none
> 21          maturity status in time t     matst2        none
> 22                      age in time t        age        none
> 23                  density in time t    density        none
> 24   individual covariate a in time t   indcova2        none
> 25 individual covariate a in time t-1   indcova1        none
> 26   individual covariate b in time t   indcovb2        none
> 27 individual covariate b in time t-1   indcovb1        none
> 28   individual covariate c in time t   indcovc2        none
> 29 individual covariate c in time t-1   indcovc1        none
> 30              stage group in time t     group2        none
> 31            stage group in time t-1     group1        none

Now we will fill in the modelparams column with the names of the variables that we are using. Bear in mind that only the terms that we actually used in modeling should be changed from none. Giving the names of other variables will cause errors in MPM creation. Here, we change the necessary terms and take a look at the final object.

our_pm$modelparams[1] <- "year2"
our_pm$modelparams[2] <- "individ"
our_pm$modelparams[4] <- "alive3"
our_pm$modelparams[5] <- "obsstatus3"
our_pm$modelparams[6] <- "size3added"
our_pm$modelparams[9] <- "repstatus3"
our_pm$modelparams[10] <- "feca3"
our_pm$modelparams[11] <- "feca2"
our_pm$modelparams[12] <- "size2added"
our_pm$modelparams[13] <- "size1added"
our_pm$modelparams[18] <- "repstatus2"
our_pm$modelparams[19] <- "repstatus1"

our_pm
>                       parameter_names mainparams modelparams
> 1                              time t      year2       year2
> 2                          individual    individ     individ
> 3                               patch      patch        none
> 4                   alive in time t+1      surv3      alive3
> 5                observed in time t+1       obs3  obsstatus3
> 6                   sizea in time t+1      size3  size3added
> 7                   sizeb in time t+1     sizeb3        none
> 8                   sizec in time t+1     sizec3        none
> 9     reproductive status in time t+1     repst3  repstatus3
> 10              fecundity in time t+1       fec3       feca3
> 11                fecundity in time t       fec2       feca2
> 12                    sizea in time t      size2  size2added
> 13                  sizea in time t-1      size1  size1added
> 14                    sizeb in time t     sizeb2        none
> 15                  sizeb in time t-1     sizeb1        none
> 16                    sizec in time t     sizec2        none
> 17                  sizec in time t-1     sizec1        none
> 18      reproductive status in time t     repst2  repstatus2
> 19    reproductive status in time t-1     repst1  repstatus1
> 20        maturity status in time t+1     matst3        none
> 21          maturity status in time t     matst2        none
> 22                      age in time t        age        none
> 23                  density in time t    density        none
> 24   individual covariate a in time t   indcova2        none
> 25 individual covariate a in time t-1   indcova1        none
> 26   individual covariate b in time t   indcovb2        none
> 27 individual covariate b in time t-1   indcovb1        none
> 28   individual covariate c in time t   indcovc2        none
> 29 individual covariate c in time t-1   indcovc1        none
> 30              stage group in time t     group2        none
> 31            stage group in time t-1     group1        none

Finally, let’s create our new MPM. Note that we will not set modelsuite this time. Instead, we will use some other model setting options.

cypmatrix3f_manual <- flefko3(stageframe = cypframe_fb, data = cypfb_v1,
  supplement = cypsupp3_fb, surv_model = surv_mod, obs_model = obs_mod,
  size_model = size_mod, repst_model = repr_mod, fec_model = fec_mod,
  paramnames = our_pm)
> Warning: Some models have not been specified, and so will be
>       set to a constant value of 1

summary(cypmatrix3f_manual)
> 
> This historical lefkoMat object contains 5 matrices.
> 
> Each matrix is square with 2916 rows and columns, and a total of 8503056 elements.
> A total of 589540 survival transitions were estimated, with 117908 per matrix.
> A total of 11760 fecundity transitions were estimated, with 2352 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 5 time steps.
> 
> Survival probability sum check (each matrix represented by column in order):
>          [,1]  [,2]  [,3]  [,4]  [,5]
> Min.    0.000 0.000 0.000 0.000 0.000
> 1st Qu. 0.992 0.992 0.992 0.992 0.992
> Median  0.999 1.000 1.000 0.999 0.999
> Mean    0.823 0.824 0.824 0.824 0.823
> 3rd Qu. 1.000 1.000 1.000 1.000 1.000
> Max.    1.000 1.000 1.000 1.000 1.000

Our summary is a little simpler than in previous cases, because we did not take advantage of the extra power provided by modelsearch(). However, we can still see the information of greatest importance. Please note that defining vital rate models manually is a highly sensitive exercise - it is likely that the first attempts will yield fatal errors, often involving odd error messages. If at first you don’t succeed, try finessing the dataset, the vital rate models, and the paramnames object to align completely with each other.

5.9 Points to remember

  1. Ahistorical MPMs should never be built using models that hold historical terms.
  2. The choices of response distributions for size and fecundity strongly influence the relationships between independent terms and their responses. Response distributions strongly influence the structure of best-fit models, and can strongly alter MPMs and their analyses.
  3. The more complex the global model, the longer it will take to run the modeling exercise. In extreme cases, in which highly complex models are run on slower computers with low memory levels, the function might run for several days. In these circumstances, the user is advised to carefully consider which terms are really needed for testing, and, if the resulting model is large, then setting the function to run on a powerful computer with lots of memory is advised.
  4. All else being equal, mixed models are better than GLMs. This is because they account for the lack of independence among transitions originating from the same individual by treating individual as a random factor. However, the complexity of mixed models may at times make the modeling exercise run very long, or possibly fail. In these cases, switching to GLMs might help yield results on a reasonably time frame.

References

Bartoń, K.A. (2014). MuMIn: Multi-model inference.
Bates, D., Maechler, M., Bolker, B. & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67, 1–48.
Brooks, M.E., Kristensen, K., Benthem, K.J. van, Magnusson, A., Berg, C.W., Nielsen, A., et al. (2017). glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. The R Journal, 9, 378–400.
Burnham, K.P. & Anderson, D.R. (2002). Model selection and multimodel inference: A practical information-theoretic approach. Springer-Verlag New York, Inc., New York, New York, USA.
Shefferson, R.P., Sandercock, B.K., Proper, J. & Beissinger, S.R. (2001). Estimating dormancy and survival of a rare herbaceous perennial using mark-recapture models. Ecology, 82, 145–156.
Yee, T.W. (2015). Vector Generalized Linear and Additive Models: With an Implementation in R.
Yee, T.W. & Wild, C.J. (1996). Vector generalized additive models. Journal of the Royal Statistical Society, Series B, 58, 481–493.
Zeileis, A., Kleiber, C. & Jackman, S. (2008). Regression models for count data in R. Journal of Statistical Software, 27, 1–25.