# 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.”*

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:

**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*.**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.**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*.**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.**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.**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.**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.**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.**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.**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.**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.**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.**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.**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:

**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.**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.**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?**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?**Juvenile vital rate estimation**- Should juvenile parameters (8) through (14) also be modeled?**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?**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?**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?**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.**Age**- Is an age-classified (Leslie) or age-by-stage MPM the main goal?**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?**Stage groups**- If stage groups are noted in the stageframe, then should they also be used as fixed, independent categorical predictors in modeling?**Censoring**- Should data points marked as questionable be used in or excluded from modeling?**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.**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?**Accuracy**- Should the accuracy of binomial models and pseudo-*R*^{2}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-*R*^{2} 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).

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)
<- c("SD", "P1", "P2", "P3", "SL", "D", "V1", "V2", "V3", "V4",
stagevector "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")
<- c(0, 0, 0, 0, 0, rep(1, 49))
indataset <- c(0, 0, 0, 0, 0, seq(from = 0, t = 24), seq(from = 1, to = 24))
sizevector <- c(0, 0, 0, 0, 0, rep(0, 25), rep(1, 24))
repvector <- c(0, 0, 0, 0, 0, 0, rep(1, 48))
obsvector <- c(0, 0, 0, 0, 0, rep(1, 49))
matvector <- c(0, 1, 1, 1, 1, rep(0, 49))
immvector <- c(1, rep(0, 53))
propvector <- c("Dormant seed", "Yr1 protocorm", "Yr2 protocorm", "Yr3 protocorm",
comments "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")
<- sf_create(sizes = sizevector, stagenames = stagevector,
cypframe_fb repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
propstatus = propvector, immstatus = immvector, indataset = indataset,
comments = comments)
<- verticalize3(data = cypdata, noyears = 6, firstyear = 2004,
cypfb_v1 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.

```
<- 5000
seeds_per_fruit <- 0.7
sl_mult
<- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "SL", "D",
cypsupp2_fb "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)
<- supplemental(stage3 = c("SD", "SD", "P1", "P1", "P2", "P3", "SL",
cypsupp3_fb "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,
1, 1, 1, 1,
sl_mult, sl_mult, 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", "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…).

```
<- modelsearch(cypfb_v1, historical = TRUE, approach = "mixed",
cypmodels3p 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-R^{2} 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-R^{2} 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.**

```
<- modelsearch(cypfb_v1, historical = FALSE, approach = "mixed",
cypmodels2p 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))`

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:

`rlefko2()`

- Creates raw ahistorical MPMs given a dataset, a stageframe, and either a supplement table or a reproductive matrix and an overwrite table.`rlefko3()`

- Creates raw historical MPMs given a dataset, a stageframe, and either a supplement table or a reproductive matrix and an overwrite table.`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.`rleslie()`

- Creates raw age-based (Leslie) MPMs given a dataset.`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.`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.`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.`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.

```
<- flefko2(stageframe = cypframe_fb, supplement = cypsupp2_fb,
cypmatrix2fp 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.

```
<- flefko3(stageframe = cypframe_fb, supplement = cypsupp3_fb,
cypmatrix3fp 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.

```
<- flefko3(stageframe = cypframe_fb, supplement = cypsupp3_fb,
cypmatrix3fp_deV 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.

```
<- lmean(cypmatrix2fp)
cyp2fp_mean <- lmean(cypmatrix3fp)
cyp3fp_mean
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.

```
<- modelsearch(cypfb_v1, historical = TRUE, approach = "mixed",
cypmodels3 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)
<- modelsearch(cypfb_v1, historical = FALSE, approach = "mixed",
cypmodels2 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.

```
<- flefko2(stageframe = cypframe_fb, supplement = cypsupp2_fb,
cypmatrix2f modelsuite = cypmodels2, data = cypfb_v1)
<- flefko3(stageframe = cypframe_fb, supplement = cypsupp3_fb,
cypmatrix3f 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).

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`

.

```
<- c(rep(0, 5), rep(1, 25), rep(2, 24))
groupvec
<- sf_create(sizes = sizevector, stagenames = stagevector,
cypframe_fb_group 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.

```
<- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "SL",
cypsupp2_fb_group "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.

```
<- flefko2(stageframe = cypframe_fb_group, data = cypfb_v1,
cypmatrix2fp_group supplement = cypsupp2_fb_group, modelsuite = cypmodels2p)
$A[[1]]
cypmatrix2fp_group> [,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
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
cypdata_env
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).

```
<- verticalize3(data = cypdata_env, noyears = 6, firstyear = 2004,
cypfb_env 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.

```
<- modelsearch(cypfb_env, historical = FALSE, approach = "mixed",
cypmodels2_env 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.

```
<- flefko2(stageframe = cypframe_fb, data = cypfb_env,
cypmatrix2f_env1 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.

```
<- c(92.2, 57.6, 96, 109.8, 111.9)
precip_values
<- flefko2(stageframe = cypframe_fb, data = cypfb_env,
cypmatrix2f_env2 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:

`lm()`

and`glm()`

from the core R package`stats`

, with the latter limited to the following response distribution families:`binomial`

,`gaussian`

,`Gamma`

, and`poisson`

`glm.nb()`

from the core R package`MASS`

`vglm()`

from package`VGAM`

(Yee & Wild 1996; Yee 2015), with response distributions limited to`pospoisson`

and`posnegbinomial`

`zeroinfl()`

from package`pscl`

(Zeileis*et al.*2008), with response distributions limited to`poisson`

and`negbin`

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

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

`<- cypfb_v1 surv_data `

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.

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

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.

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

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.

`<- size_data repr_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.

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

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)
<- glmer(alive3 ~ size2added + size1added + repstatus2 + repstatus1 +
surv_mod 1 | year2) + (1 | individ), data = surv_data, family = binomial)
(> boundary (singular) fit: see help('isSingular')
<- glmer(obsstatus3 ~ size2added + size1added + repstatus2 +
obs_mod + (1 | year2) + (1 | individ), data = obs_data, family = binomial)
repstatus1 <- glmmTMB(size3added ~ size2added + size1added + repstatus2 +
size_mod + (1 | year2) + (1 | individ), data = size_data,
repstatus1 family = truncated_nbinom2)
<- glmer(repstatus3 ~ size2added + size1added + repstatus2 +
repr_mod + (1 | year2) + (1 | individ), data = repr_data, family = binomial)
repstatus1 > boundary (singular) fit: see help('isSingular')
<- glmmTMB(feca3 ~ size2added + size1added + repstatus1 + (1 | year2) +
fec_mod 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.

```
<- create_pm()
our_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.

```
$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
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.

```
<- flefko3(stageframe = cypframe_fb, data = cypfb_v1,
cypmatrix3f_manual 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

- Ahistorical MPMs should never be built using models that hold historical terms.
- 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.
- 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.
- 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

*lme4*.

*Journal of Statistical Software*, 67, 1–48.

*et al.*(2017). glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling.

*The R Journal*, 9, 378–400.

*Model selection and multimodel inference: A practical information-theoretic approach*. Springer-Verlag New York, Inc., New York, New York, USA.

*Ecology*, 82, 145–156.

*Journal of the Royal Statistical Society, Series B*, 58, 481–493.

*Journal of Statistical Software*, 27, 1–25.