Basic theory and concepts

Life history models and population projection matrices

Population projection matrices, also sometimes referred to as matrix projection models (MPMs), are representations of the dynamics of a population across all life history stages deemed relevant, across the most relevant time interval (typically one year, but sometimes differing). They require a complete model of the organism’s life history prior to construction, and then attempt to model population dynamics as linear models of all transitions involved within that life history. Each stage is represented in the matrix by a single column and a single row, and the elements comprising the matrix show either probability of transition for an individual in stage j at time t (along the columns), to stage k at time t+1 (along the rows), or the mean rate of production of new recruits into stage k at time t+1 (along the rows) by individuals in reproductive stage j in time t (along the columns).

Here is an example of a model and matrix for one organism. This is a highly simplified stage-classified model based on a terrestrial orchid species, in this case Cypripedium parviflorum (Shefferson et al. 2001). Here, each stage is shown as a node, and each transition is shown as an arrow (a). The rates and probabilities are shown as mathematical symbols in (b), with S_kj denoting survival-transition probability from stage j in time t to stage k in time t+1, and F_F denoting the fecundity of reproductive stage F into whatever the main recruit stage in this life history.

Raw vs. function-derived matrix models, including integral projection models

The example above denotes matrix elements as symbols only and does not show how those transitions are actually estimated. Two methodologically differentiated classes of matrix models exist: raw and function-derived. The most common is the raw matrix, in which case transition probabilities in the matrix are estimated by counting all of the number of individuals alive in a particular stage at time t, and then dividing that number into the number of individuals from that set that are still alive in each possible stage in time t+1. For example, if 100 individuals are alive in stage D in 2010, and 20 of these are alive in stage D in year 2011, as are a further 40 in stage V and a further 25 in stage F, then the associated transition probabilities are 0.20, 0.40, and 0.25, respectively, and the overall survival probability for stage D is the sum of these transitions, or 0.85. Methods for estimating fecundity in raw matrices vary from system to system. In this example, we simply count the number of fruits in one year and then multiply by the mean number of seeds per fruit and the mean germination probability. In other systems, counts of actual recruits in the next year may be possible, as for example in studies of nesting birds.

Function-derived matrices are a more recent invention, and although they have been much less utilized than raw matrices, they are becoming increasingly common in the literature. In a function-derived matrix, most matrix elements are populated by kernels that link together functions representing key demographic processes governing each transition. Typically, demographic datasets are analyzed via linear models to estimate demographic parameters such as survival probability, the probability of becoming a certain size assuming survival, and fecundity rate. Elements in these matrices are then estimated as products of these linear functions set to particular inputs, which must involve all parameters governing the construction of the matrix. For example, if a matrix is from a size-classified life history, then size must determine at least some of these demographic parameters. The most common examples of function-derived matrix models are integral projection projection models (IPMs), which break up a life history into many fine-scale size classes using a continuous measure of size, and then estimate survival-transitions and fecundity according to these fine-scale size classes (Ellner and Rees 2006).

In both raw and function-derived matrices, most of the matrix is composed of survival-transition probabilities, which give the probability of an individual in stage j at time t surviving and transition to stage k at time t+1, with only some elements devoted to fecundity. In the above example, fecundity is shown in the top-right elements, which show the mean production of seeds dormant and germinating in the next time interval. Ignoring the fecundity elements, we have a survival-transition matrix (symbolized as either U or T), and the column sums correspond to the expected survival probabilities of individuals in each stage from time t to time t+1. Ignoring the survival-transition terms, we have the fecundity matrix (symbolized as F), and the column sums correspond to the expected overall fecundity of individuals in these respective stages.

Ahistorical vs. historical matrix models

The example shown above is an example of an ahistorical matrix, meaning that individual history is not incorporated into the matrix. It may sound strange to say that individual history is not incorporated into the matrices given that we generally work with datasets composed of individual histories. However, when we estimate a typical matrix, these individual histories are broken up into pairs of consecutive states across time and treated as independent of one another. This reflects a central assumption in the analysis of matrix projections, which is that the state of an individual in the next time step is only influenced by its current state. Conceptually, if an organism’s state in the next time interval is entirely determined by its current state, then its previous states do not influence these transitions. Thus, standard matrices are two dimensional and reflect only the current and next immediate state of individuals, given by the columns and rows, respectively. This is ultimately an extension of the iid assumption - that the states of individuals are independent and identically-distributed random variables. Although population projection matrices have been estimated for approximately a full century now, and we have never done a meta-analysis of all of these studies, nonetheless it is safe to say that well over 99% of published population projection matrix analyses have been conducted on ahistorical matrices. In fact, at the time of writing, we can only think of four examples of studies using a historical approach, meaning that they incorporated some degree of individual history into matrix estimation and analysis (Ehrlén 2000; Shefferson, Warren II, and Pulliam 2014; Shefferson, Mizuta, and Hutchings 2017; Shefferson et al. 2018).

The historical matrix is an extension of the population projection matrix that incorporates information on one extra previous time interval into the determination of survival and fecundity. Thus, the expected survival probability of an individual in stage j at time t to stage k at time t+1 also depends on its stage in time t-1. Population ecologists considering this problem analytically might be inclined to add an extra dimension to the matrix to deal with this, thus creating a 3d array or cube. However, this approach is computationally intensive and makes analysis challenging. Instead, we utilize the approach developed by Ehrlén (2000), in which life history stages are paired in representation in the matrix. Thus, columns now represent the From pair of stages (state in times t and t-1), and rows now represent the To pair of stages (state in times t and t+1), as in the following diagram.

Historical matrices are more data-hungry than ahistorical matrices. In the simplest terms, the historical matrix can be vastly bigger than its ahistorical counterpart, meaning that there are many more elements to estimate. The example historical matrix illustrates this issue, but a more realistic example might illustrate the issue more completely. For example, a life history with 10 stages would yield ahistorical matrices that have 10 columns and 10 rows, and so 100 elements. A historical version of such a matrix might have 100 columns and 100 rows, consisting of 10,000 elements, because the column and row stages represent pairs of stages. Clearly, 10,000 is a much bigger number than 100 - in fact it is two orders of magnitude bigger! If the life history includes 25 stages, then an ahistorical matrix has 25 rows and 25 columns, yielding 25² or 625 elements, while the associated historical matrix has 25² = 625 rows and 625 columns, yielding 25⁴ or 390,625 elements, which is three orders of magnitude larger.

The vast increase in the number of matrix elements can lead to natural worry that historical matrices may be overparameterized (and, indeed, in a short dataset of a small population that may very well be the case!). However, it turns out that most elements in a historical matrix are structural 0’s. Astute readers may have noticed that each To pair of stages and each From pair of stages includes state in time t. Most transitions in the matrix thus become logically impossible, because transition elements in historical matrices are only estimable if the stage at time t is equal in the column and row stage pairs. For example, the transition probability between stage 1 in time t-1 and stage 2 in time t (column stage), to stage 3 in time t and stage 4 in time t+1 (row stage) equals 0, because an organism cannot be in both stages 2 and 3 in time t. In fact, if there are m stages in a life history, yielding m² elements in an ahistorical matrix, then although there will be m⁴ elements in the historical matrix, only m³ will be potentially estimable while (m-1)m³ are structural 0’s (we say potentially because some of the logically possible transitions may still be biologically impossible, and so may also be structural 0’s).

Step 1: Data characterization and reorganization

The basic workflow to analysis with package lefko3 starts with the development of a life history model that encapsulates all of the appropriate life stages relevant to population dynamics. We leave it to the reader to figure out how to do this properly, though many resources exist in the literature. Once this is done, we need to characterize the main life stages in a way that is relevant to the dataset. We do this with the sf_create() function.

Function sf_create() creates what we have termed a stageframe object. This term is a combination of the terms ‘life history stage’ and ‘data frame’, because it creates a data frame in R that describes all of the life history stages. This object should be created in a way that shows how stage relates to size, reproductive status, observation status, propagule status, maturity status, and presence in the dataset. It describes stages in ways that matrix-estimating functions utilize to conduct proper computations, and also allows the user to stipulate extra information to remind him or herself of what exactly the stages mean. If the user wishes to create an IPM, then this function also allows that to be specified.

Once the stageframe has been created, the data needs to be organized into a format that lefko3 can use. We have developed two functions for this purpose, depending on the format of the data. If the data is arranged horizontally, meaning that individual life histories are recorded in a spreadsheet with rows corresponding to unique individuals and columns corresponding to their states at different times, then the verticalize3() function can turn this dataset into the proper format. If the dataset is in vertical, ahistorical format, in which individuals states are recorded in across rows in a spreadsheet, with a single row corresponding to a single observation, then the historicalize3() function can organize the data properly.

Step 2: Develop supplemental information for matrix estimation

Matrices are often estimated only partially from available demographic datasets. Some information is sometimes provided from other studies to parameterize key transitions. For example, if fecundity is a function of seed production and germination probability, then germination probability might be estimated via a separate field germination study, and then used to develop the matrix. In lefko3, this information is provided in two ways.

First, the user may specify specific rates and probabilities to incorporate directly into the matrix via the overwrite() function. This function provides an opportunity to show the specific transition, in either historical or ahistorical format, along with the exact number to use. It also allows the user to specify that a specific transition should actually be equal to another transition to be estimated within the matrix estimation procedure. For example, if there are no juveniles in a study, then it may still be valid in some circumstances to set transitions from juvenile to adult stages as equal to transitions from the smallest adult class. Such proxy estimations can be specified with overwrite(). Additionally, overwrite() allows large groups of transitions to be altered with a single line, for example all reproductive transitions.

Second, the user may provide information on which stage yields what recruit stages via fecundity, and also any modifiers on those numbers, by creating a numeric matrix equal in both dimension and row/column designation to the ahistorical matrix associated with the study. Assuming the same order of stages as in the associated stageframe object for analysis, a zero matrix is created and fecundity elements can be changed appropriately. Changing all fecundity elements to 1 simply yields fecundity estimated as equal to whatever variable corresponds to fecundity in the dataset. Multiplying by, say, 0.25 will yield estimates 1/4 of that magnitude.

The next step in the process depends on whether we wish to build raw matrices or function-derived matrices. If the former, then we can skip ahead to Step 4, but if the latter, then we need to develop models of demographic parameters.

Step 3: Develop models of demographic parameters (for function-derived matrices only)

Function-derived matrices are generated using models of key demographic parameters. Package lefko3 can create models to estimate up to nine different vital rates:

1. Survival probability - This is the probability of surviving in time step t to time step t+1, given that the individual is in some stage j. In lefko3, this parameter may be modeled as a function of size, reproductive status, patch, year, and individual identity. This parameter is required in all function-derived matrices.

2. Observation probability - This is an optional parameter denoting the probability of observation in time step t+1 of an individual in stage k given survival from time t to time t+1. This parameter is only used when at least one stage is technically unobservable. For example, some plants are capable of vegetative dormancy, in which case they are alive do not necessarily sprout in all years. In these cases, the probability of sprouting can 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 0. In lefko3, this parameter may be modeled as a function of size, reproductive status, patch, year, and individual identity.

3. Size transition probability - This is the probability of becoming size k in time step t+1 assuming survival from time t to time step t+1 and observation in that time. In lefko3, this parameter may be modeled as a function of size, reproductive status, patch, year, and individual identity. This parameter is required in all function-derived size-classified matrices.

4. Reproduction probability - This is an optional parameter denoting the probability of reproducing in time step t+1 given survival from time t to time 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, then this parameter is not needed. In lefko3, this parameter may be modeled as a function of size, reproductive status, patch, year, and individual identity.

5. Fecundity rate - This is the rate of successful production of offspring in time t by individuals alive, observable, and reproductive in that time, and the survival of those offspring into time 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 seed. In lefko3, this parameter may be modeled as a function of size, reproductive status, patch, year, and individual identity. This parameter is required in all function-derived matrices.

6. Juvenile survival probability - This is the probability of surviving from juvenile stage j in time step t to a mature stage in time step t+1. In lefko3, this parameter may be modeled as a function of patch, year, and individual identity.

7. Juvenile observation probability - This is a parameter denoting the probability of observation in time step t+1 of an individual in mature stage k given survival from a juvenile stage in time t to time t+1. In lefko3, this parameter may be modeled as a function of patch, year, and individual identity, and all other caveats noted in (2) above apply.

8. Juvenile size transition probability - This is the probability of becoming mature size k in time step t+1 assuming survival from juvenile stage j in time t to time step t+1 and observation in that time. In lefko3, this parameter may be modeled as a function of patch, year, and individual identity.

9. Juvenile reproduction probability - This is a parameter denoting the probability of reproducing in mature stage k in time step t+1 given survival from juvenile stage j in time t to time t+1, and observation in that time. In lefko3, this parameter may be modeled as a function of patch, year, and individual identity, and all other caveats in (4) apply. Note that this parameter denotes an exit from the juvenile stages of life.

Of these nine vital rates, most users will estimate at least parameters (1) survival probability, (3) size transition probability, and (5) fecundity. Parameters (2) observation probability and (4) 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 (6) through (9) should only be added if the dataset contains juvenile individuals, making the estimation of these parameters possible.

Package lefko3 includes one main function to handle model estimation: modelsearch(). Function modelsearch() handles the entire modeling process, including the development of global models, exhaustive model building, and the selection of best-fit models. Users need to provide this function with information about the following.

1. Individual history - Are the matrices to be built historical or ahistorical? If the former, then the state of the individual in time t-1 will be included in modeling.

2. Modeling approach - Should the models be estimated as generalized linear models (GLMs) or mixed linear models? While most function-derived matrix models are estimated as the former, the latter approach can correct for pseudo-replication by treating individual identity as a random factor. They also allow time and patch to be treated as random variables.

3. Suite of factors - Should both size and reproductive status be tested as causal factors? Or just one? Or is just a constant required? Should two-way interactions be included?

4. Suite of vital rates - Which demographic parameters should be estimated? The defaults are (1) survival, (3) size, and (5) fecundity.

5. Juvenile stage estimation - Should parameters (6) through (9) also be modeled?

6. Best-fit criterion - If a model with fewer parameters exists within 2.0 AICc units of the model with the lowest AICc, then should this model be used as the best-fit model, or should the model with the lowest AICc always be chosen?

7. Size distribution - Should size be modeled as a continuous variable under a Gaussian distribution, or as a count variable under either a Poisson or negative binomial distribution?

8. Fecundity distribution - Should fecundity be modeled as a continuous variable under a Gaussian distribution, or as a count variable under either a Poisson or negative binomial distribution?

9. Timing of fecundity - Modelsearch() assumes that linear models of fecundity use a metric counted or measured in time t as the response. This applies well with plants, where flowers or seeds produced might be the fecundity response. However, users not wishing to follow this default behavior can use the fectime option to stipulate a fecundity metric measured in time t+1.

10. Censoring - Should data points marked as questionable be used or eliminated?

11. Variable names - The names of all relevant variables in the dataset need to be specified. Note that some defaults can be used, particularly if the dataset has been produced using the historicalize() or verticalize() functions.

Once all inputs are provided, modelsearch() goes to work. The result will be a lefkoMod object, which is essentially a list in which first several 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 an element detailing the best-fit criterion used, and ending on a data frame showing the number of individuals and number of unique transitions used in the estimation of each model. Model selection proceeds through the dredge() function in package MuMIn (Bartoń 2014).

Advanced users of lefko3 may wish to create their own models without the package’s automated model selection function. In that case, lefko3’s matrix functions can accommodate single models developed using base R, package lme4 (Bates et al. 2015), and package glmmTMB (Brooks et al. 2017). Plans are in the works to add nlme compatibility, as well as non-linear models (e.g. general additibe models).

Step 4: Estimate matrices

Now we may create the matrices themselves. This step can be accomplished with four different functions.

1. rlefko2() - This function creates a raw ahistorical matrix given a dataset, a stageframe, a reproductive matrix, and an overwrite matrix.

2. rlefko3() - This function creates a raw historical matrix given a dataset, a stageframe, a reproductive matrix, and an overwrite matrix.

3. flefko2() - This function creates a function-derived ahistorical matrix given a dataset, a set of models, a stageframe, a reproductive matrix, and an overwrite matrix.

4. flefko3() - This function creates a function-derived historical matrix given a dataset, a set of models, a stageframe, a reproductive matrix, and an overwrite matrix.

These functions are provided with 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 (30s or so on RPS’ MacBook Pro 2019 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 4.0.

Arithmetic vs. geometric mean matrices

The function lmean() estimates mean matrices using lefkoMat objects as the input. This function is quite powerful, and can estimate both arithmetic and geometric mean matrices under sparse and non-sparse settings even for extremely large matrices with speed. It separates these means as grand mean population and patch-level matrices. Among the most interesting aspects of this function is its ability to handle time geometrically.

Demographic processes typically have geometric relationships with time. A simple example might illustrate: Imagine a population of 5 mice, each of which survives this year with 100% annual survival and produces exactly 5 offspring. It should be obvious that in the next time step there will be 30 mice - the 5 parents and 25 new juvenile mice. Let’s say that in the next year, survival is still 100% but fecundity is now 2 instead of 5, and that our previous new recruits are reproductively mature. This leads to the addition of 60 new mice into the population, leading to a total population size of 90.

Now, let’s see that we wish to represent the mousey dynamics above in a mean matrix. We take the arithmetic mean of the two annual fecundity values - 2 and 5 - and end up with 3.5, which we put into the mean matrix. Survival is 100% in both years, so averages to 100% here. We now take the initial vector of individual numbers representing population size (originally just 5 mice), and we multiply that vector by our mean matrix twice to see what the population size is supposed to be after 2 years. Since this is a very simple example, we can just take the original number of mice - 5 - and multiply by 3.5, then add back the original number of mice (5), and repeat. However, doing this leads to a total expected number of mice equal to 22.5 after the first year, and 101.25 expected mice after the second year. In both cases, our expected numbers have overshot reality!

The reason that we have overshot reality with our estimates is that reproduction is actually a geometric process, meaning that population size multiplies over time. Because of this reality, variability in reproduction and in population growth rate leads to lower actual population growth - indeed, lower annual values of fecundity and population growth rate have greater impacts on reproduction and population growth over time. This is so true that an annual population growth rate of 0 leads to all population growth ending - population extinction, with no further growth possible.

One method of dealing with this particular problem is to estimate the geometric mean of matrix elements across time, instead of estimating the arithmetic mean. The geometric mean is the nth root of the product of n items. The lmean() function includes the capacity to generate geometric mean matrices, using a logarithm-based approach in C++ that it applies to average annual matrices (note: the arithmetic mean is ALWAYS applied across patches and populations, while across time the user may stipulate the use of either geometric or arithmetic mean). This makes the lmean() quick at generating geometric means (or arithmetic means) of even the largest matrices.

One further consideration with regards to mean matrices is the sparseness of the matrices themselves relative to the data. Particularly when raw matrices have been estimated, are there cases in which matrix elements are occasionally 0 simply because in some years no individual moves through that particular stage? Such situations may cause difficulties in mean matrix estimation because they result in 0 values for matrix elements when no data exists to properly estimate them. The impact is to drag the population growth rate down under both arithmetic and geometric mean scenarios, although the latter is impacted more strongly than the former. One method of dealing with this situation is to design the life history model appropriately enough that all life history stages are occupied in all years, but this may not be possible in some cases when datasets are particularly small. Thus, lmean() also includes the capacity to ignore 0 values unless they are structural (the sparse option), where structural zeroes are defined as those elements occurring as zeroes in ALL matrices. Note, however, that this may also have unintended consequences, such as resulting in artificially high population growth rate estimates under some scenarios.

One of us (RPS) personally believes that geometric means are severely (and unfortunately) underutilized in population ecology, and wonders how our population viability analyses in particular would change if they were used….

Stable stage distribution and reproductive value

In ahistorical analyses, the stable stage distribution and the reproductive value of stages are estimated as standardized right and left eigenvectors associated with the dominant eigenvalue of the matrix. Standardization in the stable stage distribution is handled by dividing each respective element of the right eigenvector by the sum of the elements in that eigenvector. Standardization in the reproductive value case is handled by dividing each element in the left eigenvector by the value of the first non-zero element in that eigenvector.

The methods mentioned above apply to historical matrices as well. However, as described, they only provide the stable stage distribution and reproductive values of stage pairs. We provide two functions, stablestage3() and repvalue3() to allow the estimation of these vectors, and also to estimate the associated stable stage distributions and reproductive values of the original stages themselves from the historical matrices. In the stable stage distribution case, this is handled simply as the sum of all stable stage proportions for stage j in time t across all stages in time t-1. In the reproductive value case, this is calculated as the sum of all reproductive values for stage j in time t across all stages in time t-1, weighted by the stable stage proportion of the respective stage pair divided by the stable stage proportion of the respective stage at time t (Ehrlén 2000). The influence of history can make these numbers differ quite dramatically from those produced by ahistorical matrices.

Further analyses

Users can take the matrices produced by matrix creation functions in package lefko3 and plug them into matrix analysis functions in other packages. Technically, matrices are produced and stored within lefkoMat objects, which are S3 objects structured as lists. These lists include a number of elements, but among the most important are elements $A, $U, and $F, which are 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. The $labels element is a dataset giving a description of each matrix, including its population, patch, and year designations. This allows the use of all functions that work with matrices, including function in base R such as eigen(), as well as in packages such as popbio (Stubben and Milligan 2007) and popdemo (Stott, Hodgson, and Townley 2012).

Now let’s turn to some case studies using two datasets included in the package.