Chapter 6 Matrix Models III: Age-by-Stage and Leslie MPMs

“The really frightening thing about middle age is that you know you’ll grow out of it.”

— Doris Day

The first MPMs were age-classified, and are also known as Leslie MPMs in honor of their creator (Leslie 1945). However, in many organisms age is seen as less important in determining demography than other variables. Stage-classified, or Lefkovitch, MPMs were developed with this in mind, as they allow the life history of an organism to be stratified by life history stages (Lefkovitch 1965). Life history stages can be defined as necessary by any status variable or combination of status variables (Caswell 2001). So stage-based MPMs provide a great deal of flexibility and power to analyze population dynamics.

For much of the history of the MPM, the dichotomy between age-classified and stage-classified MPMs was driven by the lack of computing power to analyze complex relationships in large matrices. With the advent of powerful home computers, more and more ecologists have addressed questions using age-by-stage MPMs, which are MPMs that incorporate both and age and stage (Caswell & Salguero-Gómez 2013). Here, we provide a short introduction to the core theory before building an age-by-stage MPM. We strongly advise readers interested in this style of MPM to read Caswell et al. (2018), which is the most comprehensive overview of the theory behind age-by-stage MPMs that we are aware of.

6.1 What is an age-by-stage MPM?

To answer the question of what an age-by-stage MPM is, let us first address the question of what an age-classified MPM is. Age-classified, or Leslie, MPMs treat the life history of an organism as a series of ages rather than stages. Each age is of equal duration, and the number of ages to include is decided through demographic studies that attempt to determine which groups of ages have unique demographic characteristics. The only transitions that are possible are one-way survival transitions from one age to the next, survival transitions between the final age and itself if the lifespan has no hard limit, and fecundity from any reproductive age to age 0 (if post-breeding) or 1 (if pre-breeding).

Let’s look at an example. Let us suppose that we are studying an organism that can live for potentially many years, such as a long-lived seabird. We conduct several years of monitoring, and find that the seabird at our monitored population seems to follow a basic pattern in which survival is relatively low in the first year, higher in the second, and from the third year on survival is typically quite high and stable. Although the maximum longevity is unknown, nonetheless birds seem to live not much longer than 20 years. Fecundity occurs from the third year on, with the first two years spent as a juvenile. In this circumstance, we might develop a Leslie matrix that looks like this:

\[\begin{equation} \tag{6.1} \left[\begin{array} {rrr} 0 & 0 & F _{1,3} \\ S _{2,1} & 0 & 0 \\ 0 & S _{3,2} & S _{3,3} \\ \end{array}\right] \end{equation}\]

This matrix includes only three rows and columns. Survival-transition probabilities are typically just below the diagonal, and represent the probability of surviving from age \(i\) (column \(i\)) to age \(i+1\) (row \(i+1\)). There is a notable exception in element \(S _{3,3}\), which is the survival-transition probability for an organism to stay alive within the final age. This allows us to model the tendency of organisms to live on without an explicit maximum longevity in a state in which the probability of survival is not expected to change. Fecundity is in the top row (\(F_{1,3}\)), since all offspring begin in the first age.

Age-by-stage MPMs differ from the above because they incorporate both age and stage, and so the matrix needs to be expanded to allow both characteristics. The standard approach is to create block matrices in which each age has potentially several stages, and survival-transition probabilities must describe the probability of transitioning from each stage at one age to each stage in the next.

As an example, let’s take the three stage life history above, which has a single newborn stage and two adult stages, only one of which is reproductive. Let’s further assume that age impacts survival and fecundity only for the first four years, beyond which survival and fecundity transitions remain essentially the same. This leads to the following age-by-stage matrix:

\[\begin{equation} \tag{6.2} \left[\begin{array} {rrrrrrr} 0 & 0 & F _{1A,2C} & 0 & F _{1A,3C} & 0 & F _{1A,4+C} \\ S _{2B,1A} & 0 & 0 & 0 & 0 & 0 & 0 \\ S _{2C,1A} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & S _{3B,2B} & S _{3B,2C} & 0 & 0 & 0 & 0 \\ 0 & S _{3C,2B} & S _{3C,2C} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & S _{4B,3B} & S _{4B,3C} & S _{5+B,4B} & S _{5+B,4+C} \\ 0 & 0 & 0 & S _{4C,3B} & S _{4C,3C} & S _{5+C,4B} & S _{5+C,4+C} \\ \end{array}\right] \end{equation}\]

Here, we have only one possible initial stage, followed by two possible stages in all later ages. Looking only at the columns, column 1 refers to age 1, which only includes the seedling stage; columns 2 and 3 refer to age 2, which includes non-reproductive and reproductive adult stages; columns 4 and 5 refer to age 3, which also includes non-reproductive and reproductive adult stages; and columns 6 and 7 refer to ages 4 and higher, once again also including non-reproductive and reproductive adult stages. While the columns refer to the From ages and stages, the rows refer to the To ages and stages and are in the same order as the columns. If we ignore stage and consider only the placement of age-related transitions, then we see that we are still generally following the Leslie MPM pattern of placing survival-transition probabilities below the diagonal for all but the final transitions. We also have more fecundity terms because although only one stage is reproductive, multiple reproductive ages occur, and all of these stay within the first row (more entry rows are possible only if the life history of the organism can start from multiple different stages).

The above matrix is ahistorical. However, unlike a pure Leslie MPM, the possibility of different stages within each age allows for the estimation of historical age-by-stage MPMs. Package lefko3 does not currently support historical age-by-stage matrices, but this capability will eventually be added.

6.2 Developing function-based age-by-stage MPMs in lefko3

Package lefko3 can produce both raw (empirical) and function-based age-by-stage MPMs. Let’s start with the latter. For our example, we will use the lathyrus dataset that comes with lefko3 to illustrate the estimation of age-by-stage function-based MPMs. First, let’s load the dataset, and then look at its dimensions.

data(lathyrus)
dim(lathyrus)
> [1] 1119   38

This dataset includes information on 1,119 individuals, so there are 1,119 rows with data (not counting the header). There are 38 columns. The first two columns are variables identifying each individual (SUBPLOT refers to the patch, and GENET refers to individual identity), with each individual’s data entirely restricted to one row. This is followed by four sets of nine columns, each named VolumeXX, lnVolXX, FCODEXX, FlowXX, IntactseedXX, Dead19XX, DormantXX, Missing19XX, and SeedlingXX, where XX corresponds to the year of observation and with years organized consecutively. Thus, columns 3-11 refer to year 1988, columns 12-20 refer to year 1989, etc.

To begin, we will create a stageframe that describes the organism’s life cycle for this dataset. In this case, the life history model is a life cycle graph (figure 6.1).

Figure 6.1: Life history model of Lathyrus vernus using log leaf volume as the size classification metric

This model is based on the life history model provided in Ehrlén (2000), but we utilize a different size classification based on the log leaf volume to make the size distribution more closely match a symmetric and somewhat normal distribution.

Our stageframe needs to include complete descriptions of all stages that occur in the dataset, with each stage defined uniquely, and also needs to describe the ages for each stage portrayed in our life history model. Since this object can be used for automated classification of individuals, all sizes, reproductive states, and other characteristics defining each proper stage in the dataset need to be accounted for explicitly. This can be difficult if a few data points exist outside the range of sizes specified in the stageframe. Such points can cause problems, because rare stages can cause an overestimation of survival for existing stages under some circumstances, and can also yield spurious values for survival-transition probabilities and fecundity rates. The final description of each stage occurring in the dataset must also avoid complete overlap with any other stage also found in the dataset, although partial overlap is allowed and expected.

Before creating the stage frame, let’s explore the possible size variables. We will particularly look at summaries of the distribution of original and log sizes.

summary(c(lathyrus$Volume88, lathyrus$Volume89, lathyrus$Volume90,
    lathyrus$Volume91))
>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
>     1.8    14.7   123.0   484.2   732.5  7032.0    1248
summary(c(lathyrus$lnVol88, lathyrus$lnVol89, lathyrus$lnVol90,
    lathyrus$lnVol91))
>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
>   0.600   2.700   4.800   4.777   6.600   8.900    1248

The upper summary shows the original size, while the lower line shows the size given in logarithmic terms. We should note the size minima and maxima, because we have been using 0 as the size of vegetatively dormant individuals. The lowest uncorrected size is 1.8, with a maximum of 7032.0. The minimum corrected size is 0.600, and the maximum corrected size is 8.900. Since the minimum corrected size is above 0 (i.e. all log sizes should be positive), and since the number of NAs has not increased (increased NAs would suggest that some unusable log sizes occur in the dataset), we are still able to use the log size value 0 as an indicator of vegetative dormancy. Note, however, that vegetative dormancy is also currently included in the many NAs that occur in size variables in this dataset.

It can also help to take a look at plots of these distributions. We will plot raw and log volume (figure 6.2).

par(mfrow=c(1,2))
plot(density(c(lathyrus$Volume88, lathyrus$Volume89, lathyrus$Volume90,
    lathyrus$Volume91), na.rm = TRUE), main = "", xlab = "Volume", bty = "n")
title("a)", adj = 0)
plot(density(c(lathyrus$lnVol88, lathyrus$lnVol89, lathyrus$lnVol90,
    lathyrus$lnVol91), na.rm = TRUE), main = "", xlab = "Log volume", bty = "n")
title("b)", adj = 0)
Density plot of aboveground plant volume (a) and log volume (b).

Figure 6.2: Density plot of aboveground plant volume (a) and log volume (b).

Note how highly skewed the raw volume distribution is. This might cause some difficulty if we use raw size untransformed and with a Gaussian distribution. Certainly, a gamma distribution would be perfectly justified, and users are urged to try that approach. We will use the log volume here, which looks ‘better’ than the raw volume distribution in the sense that it is closer to some semblance of a Gaussian distribution, mostly through an increased level of symmetry. We can then assume that log volume is Gaussian-distributed and that the mean bears no relationship to the variance.

We will now develop a stageframe that incorporates the log volume of size. We will build this by creating vectors of the values describing each stage, always in the same order. Because we wish to build an age-by-stage MPM, we will also incorporate age information for each stage. Here, we include minimum and maximum ages for each stage via the vectors minima and maxima (NA in the maximum age vector is interpreted as meaning that there is no maximum).

sizevector <- c(0, 4.6, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9)
stagevector <- c("Sd", "Sdl", "Dorm", "Sz1nr", "Sz2nr", "Sz3nr", "Sz4nr",
  "Sz5nr","Sz6nr", "Sz7nr", "Sz8nr", "Sz9nr", "Sz1r", "Sz2r", "Sz3r", "Sz4r",
  "Sz5r", "Sz6r", "Sz7r", "Sz8r", "Sz9r")
repvector <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1)
obsvector <- c(0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
matvector <- c(0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
immvector <- c(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
indataset <- c(0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
binvec <- c(0, 4.6, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 
            0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5)
minima <- c(1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)
maxima <- c(NA, 1, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
  NA, NA, NA, NA, NA)
comments <- c("Dormant seed", "Seedling", "Dormant", "Size 1 Veg", "Size 2 Veg",
  "Size 3 Veg", "Size 4 Veg", "Size 5 Veg", "Size 6 Veg", "Size 7 Veg",
  "Size 8 Veg", "Size 9 Veg", "Size 1 Flo", "Size 2 Flo", "Size 3 Flo",
  "Size 4 Flo", "Size 5 Flo", "Size 6 Flo", "Size 7 Flo", "Size 8 Flo",
  "Size 9 Flo")
lathframeln <- sf_create(sizes = sizevector, stagenames = stagevector,
  repstatus = repvector, obsstatus = obsvector, propstatus = propvector, 
  immstatus = immvector, matstatus = matvector, indataset = indataset,
  binhalfwidth = binvec, minage = minima, maxage = maxima, comments = comments)

lathframeln
>    stage size size_b size_c min_age max_age repstatus obsstatus propstatus
> 1     Sd  0.0     NA     NA       1      NA         0         0          1
> 2    Sdl  4.6     NA     NA       1       1         0         1          0
> 3   Dorm  0.0     NA     NA       2      NA         0         0          0
> 4  Sz1nr  1.0     NA     NA       2      NA         0         1          0
> 5  Sz2nr  2.0     NA     NA       2      NA         0         1          0
> 6  Sz3nr  3.0     NA     NA       2      NA         0         1          0
> 7  Sz4nr  4.0     NA     NA       2      NA         0         1          0
> 8  Sz5nr  5.0     NA     NA       2      NA         0         1          0
> 9  Sz6nr  6.0     NA     NA       2      NA         0         1          0
> 10 Sz7nr  7.0     NA     NA       2      NA         0         1          0
> 11 Sz8nr  8.0     NA     NA       2      NA         0         1          0
> 12 Sz9nr  9.0     NA     NA       2      NA         0         1          0
> 13  Sz1r  1.0     NA     NA       2      NA         1         1          0
> 14  Sz2r  2.0     NA     NA       2      NA         1         1          0
> 15  Sz3r  3.0     NA     NA       2      NA         1         1          0
> 16  Sz4r  4.0     NA     NA       2      NA         1         1          0
> 17  Sz5r  5.0     NA     NA       2      NA         1         1          0
> 18  Sz6r  6.0     NA     NA       2      NA         1         1          0
> 19  Sz7r  7.0     NA     NA       2      NA         1         1          0
> 20  Sz8r  8.0     NA     NA       2      NA         1         1          0
> 21  Sz9r  9.0     NA     NA       2      NA         1         1          0
>    immstatus matstatus indataset binhalfwidth_raw sizebin_min sizebin_max
> 1          1         0         0              0.0         0.0         0.0
> 2          1         0         1              4.6         0.0         9.2
> 3          0         1         1              0.5        -0.5         0.5
> 4          0         1         1              0.5         0.5         1.5
> 5          0         1         1              0.5         1.5         2.5
> 6          0         1         1              0.5         2.5         3.5
> 7          0         1         1              0.5         3.5         4.5
> 8          0         1         1              0.5         4.5         5.5
> 9          0         1         1              0.5         5.5         6.5
> 10         0         1         1              0.5         6.5         7.5
> 11         0         1         1              0.5         7.5         8.5
> 12         0         1         1              0.5         8.5         9.5
> 13         0         1         1              0.5         0.5         1.5
> 14         0         1         1              0.5         1.5         2.5
> 15         0         1         1              0.5         2.5         3.5
> 16         0         1         1              0.5         3.5         4.5
> 17         0         1         1              0.5         4.5         5.5
> 18         0         1         1              0.5         5.5         6.5
> 19         0         1         1              0.5         6.5         7.5
> 20         0         1         1              0.5         7.5         8.5
> 21         0         1         1              0.5         8.5         9.5
>    sizebin_center sizebin_width binhalfwidthb_raw sizebinb_min sizebinb_max
> 1             0.0           0.0                NA           NA           NA
> 2             4.6           9.2                NA           NA           NA
> 3             0.0           1.0                NA           NA           NA
> 4             1.0           1.0                NA           NA           NA
> 5             2.0           1.0                NA           NA           NA
> 6             3.0           1.0                NA           NA           NA
> 7             4.0           1.0                NA           NA           NA
> 8             5.0           1.0                NA           NA           NA
> 9             6.0           1.0                NA           NA           NA
> 10            7.0           1.0                NA           NA           NA
> 11            8.0           1.0                NA           NA           NA
> 12            9.0           1.0                NA           NA           NA
> 13            1.0           1.0                NA           NA           NA
> 14            2.0           1.0                NA           NA           NA
> 15            3.0           1.0                NA           NA           NA
> 16            4.0           1.0                NA           NA           NA
> 17            5.0           1.0                NA           NA           NA
> 18            6.0           1.0                NA           NA           NA
> 19            7.0           1.0                NA           NA           NA
> 20            8.0           1.0                NA           NA           NA
> 21            9.0           1.0                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
>    sizebinc_center sizebinc_width group     comments
> 1               NA             NA     0 Dormant seed
> 2               NA             NA     0     Seedling
> 3               NA             NA     0      Dormant
> 4               NA             NA     0   Size 1 Veg
> 5               NA             NA     0   Size 2 Veg
> 6               NA             NA     0   Size 3 Veg
> 7               NA             NA     0   Size 4 Veg
> 8               NA             NA     0   Size 5 Veg
> 9               NA             NA     0   Size 6 Veg
> 10              NA             NA     0   Size 7 Veg
> 11              NA             NA     0   Size 8 Veg
> 12              NA             NA     0   Size 9 Veg
> 13              NA             NA     0   Size 1 Flo
> 14              NA             NA     0   Size 2 Flo
> 15              NA             NA     0   Size 3 Flo
> 16              NA             NA     0   Size 4 Flo
> 17              NA             NA     0   Size 5 Flo
> 18              NA             NA     0   Size 6 Flo
> 19              NA             NA     0   Size 7 Flo
> 20              NA             NA     0   Size 8 Flo
> 21              NA             NA     0   Size 9 Flo

Once the stageframe is created, we can reorganize the dataset into historically-formatted vertical (hfv) format. To handle this, we will use the verticalize3() function, which creates historically-formatted vertical datasets, as below. We will also replace NAs in size and fecundity variables with zeros for modelsearch to work properly when we build models of vital rates, so we will now set NAas0 = TRUE. Some care needs to be taken with this last step, since some authors give missing values extra meaning not present in a value of zero. In our case, a missing value indicates that a plant was dead (both size and fecundity are missing), was alive but not sprouting (size was missing), or was alive but did not produce seed (fecundity was missing). In all cases, these NA values may be replaced by zero, because other variables indicate those conditions.

We also have two choices for use as our reproductive status and fecundity variables. The first choice, FCODE88 indicates whether a plant flowered. The second choice, Intactseed88, indicates the number of seed produced. The choice of which to use depends strongly on the aims of the study. In our case, we would like to treat all plants that flowered as reproductive, but treat fecundity in terms of real seed produced. The reason is that we believe that flowering plants have a different demography than non-flowering plants, either reflecting reproductive costs, or, conversely, because flowering plants might have more resources and hence higher survival than non-flowering plants, and so we wish to separate transitions among these two groups. So, let’s use FCODE88 to indicate reproductive status, and Intactseed88 to indicate fecundity. Once complete, we will look at a summary.

Finally, note that in the input to the following function, we utilize a strictly repeating pattern of variable names arranged in the same order for each monitoring occasion. This arrangement allows us to enter only the first variable in each set, as long as noyears and blocksize are set properly and no gaps or shuffles appear in the dataset. The data management functions that we have created for lefko3 do not require such repeating patterns, but they do make the required input in the function much shorter and more succinct.

lathvertln <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
  patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
  juvcol = "Seedling1988", sizeacol = "lnVol88", repstracol = "FCODE88",
  fecacol = "Intactseed88", deadacol = "Dead1988", nonobsacol = "Dormant1988",
  stageassign = lathframeln, stagesize = "sizea", censorcol = "Missing1988",
  censorkeep = NA, NAas0 = TRUE, censor = TRUE)

dim(lathvertln)
> [1] 2552   42
summary(lathvertln)
>      rowid        popid   patchid   individ              year2     
>  Min.   :   1.0   :2552   1:674   Length:2552        Min.   :1988  
>  1st Qu.: 239.8           2:206   Class :character   1st Qu.:1988  
>  Median : 523.0           3:560   Mode  :character   Median :1989  
>  Mean   : 538.7           4:478                      Mean   :1989  
>  3rd Qu.: 823.2           5:312                      3rd Qu.:1990  
>  Max.   :1119.0           6:322                      Max.   :1990  
>    firstseen       lastseen        obsage       obslifespan       sizea1     
>  Min.   :1988   Min.   :1988   Min.   :1.000   Min.   :0.00   Min.   :0.000  
>  1st Qu.:1988   1st Qu.:1991   1st Qu.:1.000   1st Qu.:2.00   1st Qu.:0.000  
>  Median :1988   Median :1991   Median :2.000   Median :3.00   Median :2.200  
>  Mean   :1988   Mean   :1991   Mean   :1.829   Mean   :2.45   Mean   :2.947  
>  3rd Qu.:1988   3rd Qu.:1991   3rd Qu.:2.000   3rd Qu.:3.00   3rd Qu.:6.400  
>  Max.   :1990   Max.   :1991   Max.   :3.000   Max.   :3.00   Max.   :8.900  
>     repstra1          feca1           juvgiven1         obsstatus1    
>  Min.   :0.0000   Min.   : 0.0000   Min.   :0.00000   Min.   :0.0000  
>  1st Qu.:0.0000   1st Qu.: 0.0000   1st Qu.:0.00000   1st Qu.:0.0000  
>  Median :0.0000   Median : 0.0000   Median :0.00000   Median :1.0000  
>  Mean   :0.1795   Mean   : 0.9793   Mean   :0.06309   Mean   :0.5537  
>  3rd Qu.:0.0000   3rd Qu.: 0.0000   3rd Qu.:0.00000   3rd Qu.:1.0000  
>  Max.   :1.0000   Max.   :66.0000   Max.   :1.00000   Max.   :1.0000  
>    repstatus1       fecstatus1        matstatus1         alive1      
>  Min.   :0.0000   Min.   :0.00000   Min.   :0.0000   Min.   :0.0000  
>  1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.0000  
>  Median :0.0000   Median :0.00000   Median :1.0000   Median :1.0000  
>  Mean   :0.1795   Mean   :0.08895   Mean   :0.5192   Mean   :0.5823  
>  3rd Qu.:0.0000   3rd Qu.:0.00000   3rd Qu.:1.0000   3rd Qu.:1.0000  
>  Max.   :1.0000   Max.   :1.00000   Max.   :1.0000   Max.   :1.0000  
>     stage1           stage1index         sizea2         repstra2     
>  Length:2552        Min.   : 0.000   Min.   :0.000   Min.   :0.0000  
>  Class :character   1st Qu.: 0.000   1st Qu.:2.500   1st Qu.:0.0000  
>  Mode  :character   Median : 3.000   Median :4.700   Median :0.0000  
>                     Mean   : 6.115   Mean   :4.562   Mean   :0.2355  
>                     3rd Qu.:10.000   3rd Qu.:6.600   3rd Qu.:0.0000  
>                     Max.   :21.000   Max.   :8.900   Max.   :1.0000  
>      feca2          juvgiven2        obsstatus2       repstatus2    
>  Min.   : 0.000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
>  1st Qu.: 0.000   1st Qu.:0.0000   1st Qu.:1.0000   1st Qu.:0.0000  
>  Median : 0.000   Median :0.0000   Median :1.0000   Median :0.0000  
>  Mean   : 1.136   Mean   :0.1117   Mean   :0.9459   Mean   :0.2355  
>  3rd Qu.: 0.000   3rd Qu.:0.0000   3rd Qu.:1.0000   3rd Qu.:0.0000  
>  Max.   :66.000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
>    fecstatus2       matstatus2         alive2     stage2         
>  Min.   :0.0000   Min.   :0.0000   Min.   :1   Length:2552       
>  1st Qu.:0.0000   1st Qu.:1.0000   1st Qu.:1   Class :character  
>  Median :0.0000   Median :1.0000   Median :1   Mode  :character  
>  Mean   :0.1046   Mean   :0.8883   Mean   :1                     
>  3rd Qu.:0.0000   3rd Qu.:1.0000   3rd Qu.:1                     
>  Max.   :1.0000   Max.   :1.0000   Max.   :1                     
>   stage2index         sizea3         repstra3          feca3          juvgiven3
>  Min.   : 2.000   Min.   :0.000   Min.   :0.0000   Min.   : 0.000   Min.   :0  
>  1st Qu.: 5.000   1st Qu.:2.300   1st Qu.:0.0000   1st Qu.: 0.000   1st Qu.:0  
>  Median : 8.000   Median :4.300   Median :0.0000   Median : 0.000   Median :0  
>  Mean   : 9.305   Mean   :4.003   Mean   :0.2159   Mean   : 1.277   Mean   :0  
>  3rd Qu.:10.000   3rd Qu.:6.300   3rd Qu.:0.0000   3rd Qu.: 0.000   3rd Qu.:0  
>  Max.   :21.000   Max.   :8.800   Max.   :1.0000   Max.   :66.000   Max.   :0  
>    obsstatus3       repstatus3       fecstatus3       matstatus3
>  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :1   
>  1st Qu.:1.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:1   
>  Median :1.0000   Median :0.0000   Median :0.0000   Median :1   
>  Mean   :0.8264   Mean   :0.2159   Mean   :0.1105   Mean   :1   
>  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:1   
>  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1   
>      alive3          stage3           stage3index    
>  Min.   :0.0000   Length:2552        Min.   : 0.000  
>  1st Qu.:1.0000   Class :character   1st Qu.: 5.000  
>  Median :1.0000   Mode  :character   Median : 7.000  
>  Mean   :0.9138                      Mean   : 8.632  
>  3rd Qu.:1.0000                      3rd Qu.:10.000  
>  Max.   :1.0000                      Max.   :21.000

This dataset has 2,552 rows representing our original dataset of over 1100 individuals. Ordinarily, we would now go on to produce either vital rate models to create function-based MPMs. However, the fact that we are incorporating age in our analysis leads to the problem that there are many individuals in our dataset that are of unknown age. The hfv dataset that we have created includes an estimated age for each individual in each year, but this is estimated from the time in which the individual is first seen. Many individuals are first seen in the first year of the study, by which time many could have already been alive for years. So, we will subset our data to eliminate individuals first seen in the very first year of the study.

lathvertln_small <- subset(lathvertln, firstseen > 1988)
dim(lathvertln_small)
> [1] 542  42
summary(lathvertln_small)
>      rowid        popid  patchid   individ              year2     
>  Min.   :  30.0   :542   1:131   Length:542         Min.   :1989  
>  1st Qu.: 302.5          2: 72   Class :character   1st Qu.:1989  
>  Median : 583.0          3:101   Mode  :character   Median :1990  
>  Mean   : 586.8          4:132                      Mean   :1990  
>  3rd Qu.: 801.0          5: 70                      3rd Qu.:1990  
>  Max.   :1097.0          6: 36                      Max.   :1990  
>    firstseen       lastseen        obsage       obslifespan       sizea1     
>  Min.   :1989   Min.   :1989   Min.   :1.000   Min.   :0.00   Min.   :0.000  
>  1st Qu.:1989   1st Qu.:1990   1st Qu.:1.000   1st Qu.:1.00   1st Qu.:0.000  
>  Median :1989   Median :1991   Median :1.000   Median :2.00   Median :0.000  
>  Mean   :1989   Mean   :1991   Mean   :1.352   Mean   :1.36   Mean   :1.186  
>  3rd Qu.:1989   3rd Qu.:1991   3rd Qu.:2.000   3rd Qu.:2.00   3rd Qu.:2.200  
>  Max.   :1990   Max.   :1991   Max.   :2.000   Max.   :2.00   Max.   :8.400  
>     repstra1          feca1           juvgiven1        obsstatus1    
>  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.:0.0000  
>  Median :0.0000   Median : 0.0000   Median :0.0000   Median :0.0000  
>  Mean   :0.0203   Mean   : 0.1218   Mean   :0.1624   Mean   :0.3524  
>  3rd Qu.:0.0000   3rd Qu.: 0.0000   3rd Qu.:0.0000   3rd Qu.:1.0000  
>  Max.   :1.0000   Max.   :34.0000   Max.   :1.0000   Max.   :1.0000  
>    repstatus1       fecstatus1        matstatus1       alive1      
>  Min.   :0.0000   Min.   :0.00000   Min.   :0.00   Min.   :0.0000  
>  1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.00   1st Qu.:0.0000  
>  Median :0.0000   Median :0.00000   Median :0.00   Median :0.0000  
>  Mean   :0.0203   Mean   :0.01292   Mean   :0.19   Mean   :0.3524  
>  3rd Qu.:0.0000   3rd Qu.:0.00000   3rd Qu.:0.00   3rd Qu.:1.0000  
>  Max.   :1.0000   Max.   :1.00000   Max.   :1.00   Max.   :1.0000  
>     stage1           stage1index         sizea2         repstra2      
>  Length:542         Min.   : 0.000   Min.   :0.000   Min.   :0.00000  
>  Class :character   1st Qu.: 0.000   1st Qu.:2.200   1st Qu.:0.00000  
>  Mode  :character   Median : 0.000   Median :2.500   Median :0.00000  
>                     Mean   : 1.893   Mean   :3.059   Mean   :0.03506  
>                     3rd Qu.: 2.000   3rd Qu.:3.600   3rd Qu.:0.00000  
>                     Max.   :20.000   Max.   :8.400   Max.   :1.00000  
>      feca2           juvgiven2        obsstatus2       repstatus2     
>  Min.   : 0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00000  
>  1st Qu.: 0.0000   1st Qu.:0.0000   1st Qu.:1.0000   1st Qu.:0.00000  
>  Median : 0.0000   Median :0.0000   Median :1.0000   Median :0.00000  
>  Mean   : 0.1421   Mean   :0.3469   Mean   :0.9815   Mean   :0.03506  
>  3rd Qu.: 0.0000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:0.00000  
>  Max.   :34.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00000  
>    fecstatus2        matstatus2         alive2     stage2         
>  Min.   :0.00000   Min.   :0.0000   Min.   :1   Length:542        
>  1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:1   Class :character  
>  Median :0.00000   Median :1.0000   Median :1   Mode  :character  
>  Mean   :0.01661   Mean   :0.6531   Mean   :1                     
>  3rd Qu.:0.00000   3rd Qu.:1.0000   3rd Qu.:1                     
>  Max.   :1.00000   Max.   :1.0000   Max.   :1                     
>   stage2index        sizea3         repstra3           feca3        
>  Min.   : 2.00   Min.   :0.000   Min.   :0.00000   Min.   : 0.0000  
>  1st Qu.: 2.00   1st Qu.:0.000   1st Qu.:0.00000   1st Qu.: 0.0000  
>  Median : 5.00   Median :2.400   Median :0.00000   Median : 0.0000  
>  Mean   : 5.19   Mean   :2.303   Mean   :0.04059   Mean   : 0.2454  
>  3rd Qu.: 7.00   3rd Qu.:3.200   3rd Qu.:0.00000   3rd Qu.: 0.0000  
>  Max.   :20.00   Max.   :8.800   Max.   :1.00000   Max.   :48.0000  
>    juvgiven3   obsstatus3       repstatus3        fecstatus3        matstatus3
>  Min.   :0   Min.   :0.0000   Min.   :0.00000   Min.   :0.00000   Min.   :1   
>  1st Qu.:0   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:1   
>  Median :0   Median :1.0000   Median :0.00000   Median :0.00000   Median :1   
>  Mean   :0   Mean   :0.7306   Mean   :0.04059   Mean   :0.01292   Mean   :1   
>  3rd Qu.:0   3rd Qu.:1.0000   3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:1   
>  Max.   :0   Max.   :1.0000   Max.   :1.00000   Max.   :1.00000   Max.   :1   
>      alive3          stage3           stage3index    
>  Min.   :0.0000   Length:542         Min.   : 0.000  
>  1st Qu.:1.0000   Class :character   1st Qu.: 3.000  
>  Median :1.0000   Mode  :character   Median : 5.000  
>  Mean   :0.7989                      Mean   : 4.969  
>  3rd Qu.:1.0000                      3rd Qu.: 6.000  
>  Max.   :1.0000                      Max.   :21.000

We are now down to only 542 rows, so we have lost close to 80% of the transition data here. That is a steep price to pay for accurate inference, but it is necessary in this case.

Let’s look at fecundity. Fecundity is integer-based, suggesting that it can be treated as a count variable. This package currently allows eight choices of count distributions: Gaussian, gamma, Poisson, negative binomial, zero-inflated Poisson, zero-inflated negative binomial, zero-truncated Poisson, and zero-truncated negative binomial. To assess which to use, we should first assess whether the mean and variance of the count are equal using a dispersion test. This test allows us to test whether the variance is greater than that expected under our mean value for fecundity using a chi-squared test. If it is not significantly different, then we may use some variant of the Poisson distribution. If the data are overdispersed, then we should use the negative binomial distribution. We should also test whether the number of zeroes is more than expected under these distributions, and make the distribution zero-inflated if so. Note that, because we have not excluded 0s from fecundity using reproductive status, we should not use a zero-truncated distribution.

Let’s start off by looking at a bar plot of the distribution of fecundity (figure 6.3).

hist(subset(lathvertln_small, repstatus2 == 1)$feca2, main = "Fecundity",
  xlab = "Intact seeds produced in occasion t")
Histogram of fecundity in occasion t

Figure 6.3: Histogram of fecundity in occasion t

We see that the distribution conforms to a classic count variable with a very low mean value. The first bar suggests that there may be too many zeroes, as well. Let’s now formally test our assumptions, that the mean equals the variance and that the number of zeroes meets expectation. Both tests use chi-squared distribution-based approaches, with the zero-inflation test in particular based on Broek (1995).

sf_distrib(lathvertln_small, sizea = c("sizea3", "sizea2"),
  fec = c("feca3", "feca2"), obs3 = "obsstatus3",
  repst = c("repstatus3", "repstatus2"))
> Non-integer values detected, so will not test for overdispersion and zero-inflation in sizea
> Mean fec is 4.6
> 
> The variance in fec is 73.54
> 
> The probability of this dispersion level by chance assuming that
> the true mean fec = variance in fec,
> and an alternative hypothesis of overdispersion, is 3.787e-145
> 
> Fecundity is significantly overdispersed.
> 
> 
> Mean lambda in fec is 0.01005
> The actual number of 0s in fec is 7
> The expected number of 0s in fec under the null hypothesis is 0.1508
> The probability of this deviation in 0s from expectation by chance is 1.122e-73
> 
> Fecundity is significantly zero-inflated.

We see that fecundity is significantly overdispersed, and has a significant excess of zeros. So, we should use a zero-inflated negative binomial distribution here.

Next, let’s create an ahistorical supplement table organizing the extra data that we need to incorporate into our matrices. Each row refers to a specific transition. The first two of these transitions are set to specific probabilities, which are the probabilities of germination and seed dormancy, estimated from a separate study. The final two terms are fecundity multipliers, which mark which transitions correspond to fecundity and provide information on what multiple of fecundity estimated via linear modeling applies to each case. Note that we can also include proxy transitions, in which we define a specific transition as being equal to another in the matrix. The latter approach is useful when some transitions cannot be estimated given a particular dataset, and so need to be set to other, proxy values that are estimable.

lathsupp2 <- supplemental(stage3 = c("Sd", "Sdl", "Sd", "Sdl"), 
  stage2 = c("Sd", "Sd", "rep", "rep"), givenrate = c(0.345, 0.054, NA, NA),
  multiplier = c(NA, NA, 0.345, 0.054), type = c(1, 1, 3, 3),
  stageframe = lathframeln, historical = FALSE)

lathsupp2
>   stage3 stage2 stage1 eststage3 eststage2 eststage1 givenrate multiplier
> 1     Sd     Sd   <NA>      <NA>      <NA>      <NA>     0.345      1.000
> 2    Sdl     Sd   <NA>      <NA>      <NA>      <NA>     0.054      1.000
> 3     Sd    rep   <NA>      <NA>      <NA>      <NA>        NA      0.345
> 4    Sdl    rep   <NA>      <NA>      <NA>      <NA>        NA      0.054
>   convtype convtype_t12
> 1        1            1
> 2        1            1
> 3        3            1
> 4        3            1

Next we will run the modelsearch function with the new vertical dataset. This function will develop our best-fit vital rate models for us. This function looks simple, but it automates several crucial and complex tasks in MPM analysis. Specifically, it automates 1) the building of global models for each vital rate requested, 2) the exhaustive construction of all reduced models, and 3) the selection of the best-fit models. In relation to our previous uses of this function in chapter 5, the most noteworthy difference is the inclusion of an age term (age = "obsage", which we know from looking at the summary of the vertical dataset). We will only develop the full effects models here, which include main effects and all two-way interactions between age and other terms, including size and reproductive status. Let’s start off by generating a set of vital rate models that covers the entire population.

lathmodelsln2 <- modelsearch(lathvertln_small, historical = FALSE,
  approach = "mixed", suite = "full", bestfit = "AICc&k", juvestimate = "Sdl",
  vitalrates = c("surv", "obs", "size", "repst", "fec"), sizedist = "gaussian",
  fecdist = "negbin", indiv = "individ", year = "year2", age = "obsage",
  year.as.random = TRUE, patch.as.random = TRUE, show.model.tables = TRUE,
  fec.zero = TRUE, quiet = TRUE)

Let’s see a summary of the lefkoMod object that we have created.

summary(lathmodelsln2)
> This LefkoMod object includes 8 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 ~ obsage + repstatus2 + (1 | year2) + (1 | individ)
>    Data: subdata
>       AIC       BIC    logLik  deviance  df.resid 
>  314.3559  333.7024 -152.1779  304.3559       349 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 0       
>  year2   (Intercept) 0       
> Number of obs: 354, groups:  individ, 190; year2, 2
> Fixed Effects:
> (Intercept)       obsage   repstatus2  
>      2.5636      -0.6252      29.6185  
> 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 ~ obsage + sizea2 + (1 | year2) + (1 | individ) +  
>     obsage:sizea2
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
> 150.6174 172.7595 -69.3087 138.6174      290 
> Random effects:
>  Groups  Name        Std.Dev. 
>  individ (Intercept) 1.741e+01
>  year2   (Intercept) 1.408e-04
> Number of obs: 296, groups:  individ, 162; year2, 2
> Fixed Effects:
>   (Intercept)         obsage         sizea2  obsage:sizea2  
>        -9.892         15.058          6.489         -4.448  
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 2 lme4 warnings 
> 
> 
> 
> Size model:
> Linear mixed model fit by REML ['lmerMod']
> Formula: sizea3 ~ sizea2 + (1 | year2) + (1 | individ)
>    Data: subdata
> REML criterion at convergence: 700.2483
> Random effects:
>  Groups   Name        Std.Dev.
>  individ  (Intercept) 0.0000  
>  year2    (Intercept) 0.4826  
>  Residual             0.8850  
> Number of obs: 266, groups:  individ, 158; year2, 2
> Fixed Effects:
> (Intercept)       sizea2  
>      0.7088       0.7777  
> optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 
> 
> 
> 
> 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 + sizea2 + (1 | year2) + (1 | individ)
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
>  74.2007  92.1181 -32.1003  64.2007      261 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 101.78  
>  year2   (Intercept)  52.35  
> Number of obs: 266, groups:  individ, 158; year2, 2
> Fixed Effects:
> (Intercept)   repstatus2       sizea2  
>     -184.21       -13.28        19.93  
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 
> 
> 
> 
> Fecundity model:
> Formula:          feca2 ~ obsage + sizea2 + (1 | year2) + (1 | individ)
> Zero inflation:         ~(1 | year2) + (1 | individ)
> Data: subdata
>       AIC       BIC    logLik  df.resid 
>  80.13600  88.63595 -31.06800        10 
> Random-effects (co)variances:
> 
> Conditional model:
>  Groups  Name        Std.Dev. 
>  year2   (Intercept) 2.487e-04
>  individ (Intercept) 2.897e-05
> 
> Zero-inflation model:
>  Groups  Name        Std.Dev.
>  year2   (Intercept) 0.001679
>  individ (Intercept) 2.021180
> 
> Number of obs: 19 / Conditional model: year2, 2; individ, 16 / Zero-inflation model: year2, 2; individ, 16
> 
> Dispersion parameter for nbinom2 family (): 2.45e+08 
> 
> Fixed Effects:
> 
> Conditional model:
> (Intercept)       obsage       sizea2  
>     -7.1322       0.4786       1.2036  
> 
> Zero-inflation model:
> (Intercept)  
>     -0.1414  
> 
> 
> Juvenile survival model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
>   Approximation) [glmerMod]
>  Family: binomial  ( logit )
> Formula: alive3 ~ (1 | year2) + (1 | individ)
>    Data: subdata
>       AIC       BIC    logLik  deviance  df.resid 
>  225.7193  235.4286 -109.8596  219.7193       185 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 0.3755  
>  year2   (Intercept) 0.0000  
> Number of obs: 188, groups:  individ, 155; year2, 2
> Fixed Effects:
> (Intercept)  
>       1.024  
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 
> 
> 
> 
> Juvenile observation model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
>   Approximation) [glmerMod]
>  Family: binomial  ( logit )
> Formula: obsstatus3 ~ (1 | year2) + (1 | individ)
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
>  28.1909  36.9509 -11.0955  22.1909      134 
> Random effects:
>  Groups  Name        Std.Dev. 
>  individ (Intercept) 5.638e+01
>  year2   (Intercept) 1.312e-05
> Number of obs: 137, groups:  individ, 122; year2, 2
> Fixed Effects:
> (Intercept)  
>        13.8  
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 
> 
> 
> 
> Juvenile size model:
> Linear mixed model fit by REML ['lmerMod']
> Formula: sizea3 ~ (1 | year2) + (1 | individ)
>    Data: subdata
> REML criterion at convergence: 145.4426
> Random effects:
>  Groups   Name        Std.Dev. 
>  individ  (Intercept) 4.504e-06
>  year2    (Intercept) 1.430e-01
>  Residual             4.139e-01
> Number of obs: 130, groups:  individ, 115; year2, 2
> Fixed Effects:
> (Intercept)  
>       2.288  
> optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 
> 
> 
> 
> Juvenile secondary size model:
> [1] 1
> 
> 
> 
> Juvenile tertiary size model:
> [1] 1
> 
> 
> 
> Juvenile reproduction model:
> [1] 0
> 
> 
> 
> Juvenile maturity model:
> [1] 1
> 
> 
> 
> 
> 
> Number of models in survival table: 18
> 
> Number of models in observation table: 18
> 
> Number of models in size table: 18
> 
> Number of models in secondary size table: 1
> 
> Number of models in tertiary size table: 1
> 
> Number of models in reproduction status table: 17
> 
> Number of models in fecundity table: 285
> 
> 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 190 individuals and 354 individual transitions.
> Survival accuracy is 0.836.
> Observation estimated with 162 individuals and 296 individual transitions.
> Observation accuracy is 0.973.
> Primary size estimated with 158 individuals and 266 individual transitions.
> Primary size pseudo R-squared is 0.692.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 158 individuals and 266 individual transitions.
> Reproductive status accuracy is 1.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Fecundity pseudo R-squared is 0.863.
> Juvenile survival estimated with 155 individuals and 188 individual transitions.
> Juvenile survival accuracy is 0.729.
> Juvenile observation estimated with 122 individuals and 137 individual transitions.
> Juvenile observation accuracy is 1.
> Juvenile primary size estimated with 115 individuals and 130 individual transitions.
> Juvenile primary size pseudo R-squared is 0.107.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity probability not estimated.

We see that age is influential in survival, observation, and fecundity, though not in other vital rates. Accuracy in our binomial models is high though not fabulous.

Next, we will create a second set that includes patch as a random factor. This model set will allow us to explore patch dynamics in addition to the population dynamics of the previous set.

lathmodelsln2p <- modelsearch(lathvertln_small, historical = FALSE,
  approach = "mixed", suite = "full", bestfit = "AICc&k", juvestimate = "Sdl",
  vitalrates = c("surv", "obs", "size", "repst", "fec"), sizedist = "gaussian",
  fecdist = "negbin", indiv = "individ", patch = "patchid", year = "year2",
  age = "obsage", year.as.random = TRUE, patch.as.random = TRUE,
  show.model.tables = TRUE, fec.zero = TRUE, quiet = TRUE)

And a summary.

summary(lathmodelsln2p)
> This LefkoMod object includes 8 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 ~ obsage + repstatus2 + (1 | year2) + (1 | patchid) +  
>     (1 | individ)
>    Data: subdata
>       AIC       BIC    logLik  deviance  df.resid 
>  307.2136  330.4294 -147.6068  295.2136       348 
> Random effects:
>  Groups  Name        Std.Dev. 
>  individ (Intercept) 7.285e-07
>  patchid (Intercept) 5.828e-01
>  year2   (Intercept) 0.000e+00
> Number of obs: 354, groups:  individ, 190; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept)       obsage   repstatus2  
>      2.7389      -0.6484      21.7793  
> 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 ~ obsage + sizea2 + (1 | year2) + (1 | patchid) +  
>     (1 | individ) + obsage:sizea2
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
> 153.0167 178.8493 -69.5084 139.0167      289 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 17.8608 
>  patchid (Intercept)  0.1586 
>  year2   (Intercept)  0.5055 
> Number of obs: 296, groups:  individ, 162; patchid, 6; year2, 2
> Fixed Effects:
>   (Intercept)         obsage         sizea2  obsage:sizea2  
>       -12.146         16.752          7.071         -4.800  
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 2 lme4 warnings 
> 
> 
> 
> Size model:
> Linear mixed model fit by REML ['lmerMod']
> Formula: sizea3 ~ sizea2 + (1 | year2) + (1 | patchid) + (1 | individ)
>    Data: subdata
> REML criterion at convergence: 700.2483
> Random effects:
>  Groups   Name        Std.Dev.
>  individ  (Intercept) 0.0000  
>  patchid  (Intercept) 0.0000  
>  year2    (Intercept) 0.4826  
>  Residual             0.8850  
> Number of obs: 266, groups:  individ, 158; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept)       sizea2  
>      0.7088       0.7777  
> optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 
> 
> 
> 
> 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 ~ obsage + sizea2 + (1 | year2) + (1 | patchid) +  
>     (1 | individ)
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
>  75.5686  97.0696 -31.7843  63.5686      260 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 105.936 
>  patchid (Intercept)   7.812 
>  year2   (Intercept)   1.106 
> Number of obs: 266, groups:  individ, 158; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept)       obsage       sizea2  
>     -236.07        66.23        18.97  
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 2 lme4 warnings 
> 
> 
> 
> Fecundity model:
> Formula:          feca2 ~ sizea2 + (1 | year2) + (1 | patchid) + (1 | individ)
> Zero inflation:         ~(1 | year2) + (1 | patchid) + (1 | individ)
> Data: subdata
>       AIC       BIC    logLik  df.resid 
>  83.29320  92.73759 -31.64660         9 
> Random-effects (co)variances:
> 
> Conditional model:
>  Groups  Name        Std.Dev. 
>  year2   (Intercept) 1.936e-04
>  patchid (Intercept) 6.955e-05
>  individ (Intercept) 2.645e-04
> 
> Zero-inflation model:
>  Groups  Name        Std.Dev. 
>  year2   (Intercept) 6.761e-03
>  patchid (Intercept) 8.938e-06
>  individ (Intercept) 1.777e+00
> 
> Number of obs: 19 / Conditional model: year2, 2; patchid, 6; individ, 16 / Zero-inflation model: year2, 2; patchid, 6; individ, 16
> 
> Dispersion parameter for nbinom2 family (): 2.61e+09 
> 
> Fixed Effects:
> 
> Conditional model:
> (Intercept)       sizea2  
>      -6.265        1.157  
> 
> Zero-inflation model:
> (Intercept)  
>     -0.1917  
> 
> 
> Juvenile survival model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
>   Approximation) [glmerMod]
>  Family: binomial  ( logit )
> Formula: alive3 ~ (1 | year2) + (1 | patchid) + (1 | individ)
>    Data: subdata
>       AIC       BIC    logLik  deviance  df.resid 
>  227.5465  240.4922 -109.7732  219.5465       184 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 0.3386  
>  patchid (Intercept) 0.2282  
>  year2   (Intercept) 0.0000  
> Number of obs: 188, groups:  individ, 155; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept)  
>       1.015  
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 
> 
> 
> 
> Juvenile observation model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
>   Approximation) [glmerMod]
>  Family: binomial  ( logit )
> Formula: obsstatus3 ~ (1 | year2) + (1 | patchid) + (1 | individ)
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
>  30.1909  41.8709 -11.0955  22.1909      133 
> Random effects:
>  Groups  Name        Std.Dev. 
>  individ (Intercept) 5.638e+01
>  patchid (Intercept) 0.000e+00
>  year2   (Intercept) 1.494e-04
> Number of obs: 137, groups:  individ, 122; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept)  
>        13.8  
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 
> 
> 
> 
> Juvenile size model:
> Linear mixed model fit by REML ['lmerMod']
> Formula: sizea3 ~ (1 | year2) + (1 | patchid) + (1 | individ)
>    Data: subdata
> REML criterion at convergence: 145.4426
> Random effects:
>  Groups   Name        Std.Dev.
>  individ  (Intercept) 0.0000  
>  patchid  (Intercept) 0.0000  
>  year2    (Intercept) 0.1430  
>  Residual             0.4139  
> Number of obs: 130, groups:  individ, 115; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept)  
>       2.288  
> optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 
> 
> 
> 
> Juvenile secondary size model:
> [1] 1
> 
> 
> 
> Juvenile tertiary size model:
> [1] 1
> 
> 
> 
> Juvenile reproduction model:
> [1] 0
> 
> 
> 
> Juvenile maturity model:
> [1] 1
> 
> 
> 
> 
> 
> Number of models in survival table: 18
> 
> Number of models in observation table: 18
> 
> Number of models in size table: 18
> 
> Number of models in secondary size table: 1
> 
> Number of models in tertiary size table: 1
> 
> Number of models in reproduction status table: 7
> 
> Number of models in fecundity table: 265
> 
> 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 190 individuals and 354 individual transitions.
> Survival accuracy is 0.836.
> Observation estimated with 162 individuals and 296 individual transitions.
> Observation accuracy is 0.973.
> Primary size estimated with 158 individuals and 266 individual transitions.
> Primary size pseudo R-squared is 0.692.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 158 individuals and 266 individual transitions.
> Reproductive status accuracy is 1.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Fecundity pseudo R-squared is 0.877.
> Juvenile survival estimated with 155 individuals and 188 individual transitions.
> Juvenile survival accuracy is 0.729.
> Juvenile observation estimated with 122 individuals and 137 individual transitions.
> Juvenile observation accuracy is 1.
> Juvenile primary size estimated with 115 individuals and 130 individual transitions.
> Juvenile primary size pseudo R-squared is 0.107.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity probability not estimated.

Incorporating patch has led to some changes in the models, most notably that now reproductive status is a function of age while fecundity is not. Interesting patterns worth exploring….

Next, we will estimate the ahistorical sets of matrices. We will match the ahistorical age-by-stage matrix estimation function, aflefko2(), with the appropriate ahistorical input, including the ahistorical lefkoMod objects lathmodelsln2 and lathmodelsln2p. Model sets that include historical terms should not be used to create ahistorical matrices, since the coefficients in the best-fit models are estimated assuming a specific model structure that either relies on historical terms or does not. Historical vital rate models may yield biased results if used to construct ahistorical matrices. Also note that lefko3 does not currently allow the construction of historical age-by-stage MPMs. Let’s start off by developing the population-only MPM.

lathmat2age <- aflefko2(year = "all", stageframe = lathframeln,
  modelsuite = lathmodelsln2, data = lathvertln_small, supplement = lathsupp2,
  continue = TRUE, reduce = FALSE)

summary(lathmat2age)
> 
> This ahistorical lefkoMat object contains 2 matrices.
> 
> Each matrix is square with 42 rows and columns, and a total of 1764 elements.
> A total of 682 survival transitions were estimated, with 341 per matrix.
> A total of 36 fecundity transitions were estimated, with 18 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 2 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 190 individuals and 354 individual transitions.
> Observation estimated with 162 individuals and 296 individual transitions.
> Primary size estimated with 158 individuals and 266 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 158 individuals and 266 individual transitions.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival estimated with 155 individuals and 188 individual transitions.
> Juvenile observation estimated with 122 individuals and 137 individual transitions.
> Juvenile primary size estimated with 115 individuals and 130 individual transitions.
> 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]
> Min.    0.000 0.000
> 1st Qu. 0.000 0.000
> Median  0.172 0.172
> Mean    0.395 0.409
> 3rd Qu. 0.788 0.788
> Max.    1.000 1.000

This first model set led to the development of two matrices, because although there are four years of data, we have limited our use of the data to individuals first seen from the second year on. Function aflefko2() has assessed patterns in the data set and has found that the maximum age is 2 years (it would have been 3 years if we allowed the first year). We have informed the function that this age is not terminal and that the demography of that age should continue onward (continue = TRUE). This has resulted in a block matrix with two ages and 21 stages, and so 42 age-stage combinations and 42 rows and columns. Of course, this particular plant species is actually long-lived, and so there may very well be further vital rate variability in across the lifespan of the plant. However, as our dataset only includes four years of study and we do not have absolute ages for any plant, we can only incur four years of relative age at best. Thus, our settings are actually the most parsimonious possible under the circumstances.

The quality control section gives us a sense of the amount of data used to model each vital rate, and also shows us that the survival-transition (U) matrices are composed entirely of proper probabilities yielding stage survival probability falling between 0.0 and 1.0. Matrix estimation can sometimes create spurious values, such as stage survival greater than 1.0. Such values can occur for a variety of reasons, but the most common is the inclusion through a supplemental table or overwrite table of externally-determined survival probabilities that are too high. Make sure to check your matrix column sums each time you estimate MPMs to prevent this problem. Survival greater than 1.0 can lead to strange effects on metrics of population dynamics.

Let’s now develop the patch-level MPMs.

lathmat2agep <- aflefko2(year = "all", patch = "all", stageframe = lathframeln,
  modelsuite = lathmodelsln2p, data = lathvertln_small, supplement = lathsupp2,
  continue = TRUE, reduce = FALSE)

summary(lathmat2agep)
> 
> This ahistorical lefkoMat object contains 12 matrices.
> 
> Each matrix is square with 42 rows and columns, and a total of 1764 elements.
> A total of 3876 survival transitions were estimated, with 323 per matrix.
> A total of 216 fecundity transitions were estimated, with 18 per matrix.
> This lefkoMat object covers 1 population, 6 patches, and 2 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 190 individuals and 354 individual transitions.
> Observation estimated with 162 individuals and 296 individual transitions.
> Primary size estimated with 158 individuals and 266 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 158 individuals and 266 individual transitions.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival estimated with 155 individuals and 188 individual transitions.
> Juvenile observation estimated with 122 individuals and 137 individual transitions.
> Juvenile primary size estimated with 115 individuals and 130 individual transitions.
> 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.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
> Median  0.172 0.172 0.168 0.172 0.172 0.172 0.148 0.172 0.166 0.172 0.172 0.172
> Mean    0.415 0.430 0.389 0.402 0.412 0.427 0.369 0.382 0.387 0.401 0.409 0.424
> 3rd Qu. 0.879 0.879 0.757 0.757 0.865 0.865 0.667 0.667 0.750 0.750 0.851 0.851
> Max.    1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000

This second model set led to the development of 12 matrices, reflecting our subset number of three years and six patches. The rest of the output looks quite similar and even the survival-transition matrix column sum summaries look extremely similar, suggesting little impact of patch.

We can get a sense of what these matrices look like by visualizing them. Let’s use the image3() function to look at just one (figure 6.4).

image3(lathmat2age, used = 1)
Visualization of 1st A Matrix. Red area corresponds to non-zero elements

Figure 6.4: Visualization of 1st A Matrix. Red area corresponds to non-zero elements

> [[1]]
> NULL

The clear squares refer to zero elements, and the red elements refer to non-zero values corresponding to survival transitions and fecundity. The vast number of zeros may be surprising, but this matrix is a supermatrix organized by age first, with stage organizing within-age blocks. The first age is age 1, which cannot be adult, and so we find zeros in adult stages at age 1. The adult block occurs from age 2, and this block can perpetuate indefinitely. The number of elements estimated is greater now than in the typical ahistorical MPM, because now we have added age as a major factor for analysis. This matrix is overwhelmingly composed of elements that must be zero, and so it is a rather sparse matrix ((341 + 18) / 1764 = 20.4% of elements).

Given the amount of white space, we might prefer to remove impossible age-stage combinations and reduce the size of the matrices. We can do this by recreating our matrices with reduce = TRUE. Let’s try that with the main population set.

lathmat2age_red <- aflefko2(year = "all", stageframe = lathframeln,
  modelsuite = lathmodelsln2, data = lathvertln_small, supplement = lathsupp2,
  continue = TRUE, reduce = TRUE)

summary(lathmat2age_red)
> 
> This ahistorical lefkoMat object contains 2 matrices.
> 
> Each matrix is square with 22 rows and columns, and a total of 484 elements.
> A total of 682 survival transitions were estimated, with 341 per matrix.
> A total of 36 fecundity transitions were estimated, with 18 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 2 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 190 individuals and 354 individual transitions.
> Observation estimated with 162 individuals and 296 individual transitions.
> Primary size estimated with 158 individuals and 266 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 158 individuals and 266 individual transitions.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival estimated with 155 individuals and 188 individual transitions.
> Juvenile observation estimated with 122 individuals and 137 individual transitions.
> Juvenile primary size estimated with 115 individuals and 130 individual transitions.
> 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]
> Min.    0.000 0.000
> 1st Qu. 0.752 0.781
> Median  0.788 0.788
> Mean    0.755 0.782
> 3rd Qu. 0.997 0.995
> Max.    1.000 1.000

This exercise has eliminated 20 rows and columns, yielding a matrix with 22 rows and 22 columns. The total number of estimated elements has not changed, meaning that our matrices are now much denser ((341 + 18) / 484 = 74.2%). Let’s take a look at an image of the first matrix (figure 6.5).

image3(lathmat2age_red, used = 1)
Visualization of 1st A Matrix. Red area corresponds to non-zero elements

Figure 6.5: Visualization of 1st A Matrix. Red area corresponds to non-zero elements

> [[1]]
> NULL

Some of the white space has been reduced, and we see greater coverage of the matrix by non-zero elements.

We can see the order of ages and stages using the agestages element of the lefkoMat object we produced, as below. Note that our matrix is 22 rows by 22 columns, and this object gives us the exact order used.

lathmat2age_red$agestages
>      stage_id stage age
> 1           1    Sd   1
> 2           2   Sdl   1
> 1.1         1    Sd   2
> 3.1         3  Dorm   2
> 4.1         4 Sz1nr   2
> 5.1         5 Sz2nr   2
> 6.1         6 Sz3nr   2
> 7.1         7 Sz4nr   2
> 8.1         8 Sz5nr   2
> 9.1         9 Sz6nr   2
> 10.1       10 Sz7nr   2
> 11.1       11 Sz8nr   2
> 12.1       12 Sz9nr   2
> 13.1       13  Sz1r   2
> 14.1       14  Sz2r   2
> 15.1       15  Sz3r   2
> 16.1       16  Sz4r   2
> 17.1       17  Sz5r   2
> 18.1       18  Sz6r   2
> 19.1       19  Sz7r   2
> 20.1       20  Sz8r   2
> 21.1       21  Sz9r   2

Now let’s estimate the element-wise arithmetic mean matrices. The first lefkoMat object created will include a single mean matrix, while the second will include six mean matrices for the patches followed by a grand mean matrix, yielding a total of seven matrices.

lathmat2mean <- lmean(lathmat2age)
lathmat2pmean <- lmean(lathmat2agep)

summary(lathmat2mean)
> 
> This ahistorical lefkoMat object contains 1 matrix.
> 
> Each matrix is square with 42 rows and columns, and a total of 1764 elements.
> A total of 361 survival transitions were estimated, with 361 per matrix.
> A total of 18 fecundity transitions were estimated, with 18 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 0 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 190 individuals and 354 individual transitions.
> Observation estimated with 162 individuals and 296 individual transitions.
> Primary size estimated with 158 individuals and 266 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 158 individuals and 266 individual transitions.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival estimated with 155 individuals and 188 individual transitions.
> Juvenile observation estimated with 122 individuals and 137 individual transitions.
> Juvenile primary size estimated with 115 individuals and 130 individual transitions.
> 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]
> Min.    0.000
> 1st Qu. 0.000
> Median  0.172
> Mean    0.402
> 3rd Qu. 0.788
> Max.    1.000
summary(lathmat2pmean)
> 
> This ahistorical lefkoMat object contains 7 matrices.
> 
> Each matrix is square with 42 rows and columns, and a total of 1764 elements.
> A total of 2275 survival transitions were estimated, with 325 per matrix.
> A total of 126 fecundity transitions were estimated, with 18 per matrix.
> This lefkoMat object covers 1 population, 7 patches, and 0 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 190 individuals and 354 individual transitions.
> Observation estimated with 162 individuals and 296 individual transitions.
> Primary size estimated with 158 individuals and 266 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 158 individuals and 266 individual transitions.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival estimated with 155 individuals and 188 individual transitions.
> Juvenile observation estimated with 122 individuals and 137 individual transitions.
> Juvenile primary size estimated with 115 individuals and 130 individual transitions.
> 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]
> Min.    0.000 0.000 0.000 0.000 0.000 0.000 0.000
> 1st Qu. 0.000 0.000 0.000 0.000 0.000 0.000 0.000
> Median  0.172 0.172 0.172 0.172 0.172 0.172 0.172
> Mean    0.423 0.395 0.420 0.375 0.394 0.417 0.404
> 3rd Qu. 0.879 0.757 0.865 0.667 0.750 0.851 0.795
> Max.    1.000 1.000 1.000 1.000 1.000 1.000 1.000

We see one overall population mean in the first case, and a set of six patch-level means and one population mean in the second case. Note that the population mean in each case should be a little different, because the population mean in the second set weights each patch equally, while the overall population mean in the first case weights individual transitions equally regardless of patch of origin.

6.3 Developing raw (empirical) age-by-stage MPMs in lefko3

Now let’s focus our attention on raw (empirical) age-by-stage MPMs. For this, we will use the function arlefko2(). Raw MPMs generally require a great deal more data to properly parameterize. We will also need a stageframe with fewer stages, in order to limit the number of artificial zeros entering our matrices. Let’s start with the following life history model, which assumes a grand total of seven stages.

Figure 6.6: Life history model of Lathyrus vernus using log leaf volume as the size classification metric

Let’s now look at our stageframe for this model. The actual stage definitions in terms of size bins are no longer based on the log leaf volume. Instead, they are based on the raw leaf volume. In order to deal with the fact that small plants are far more common than large plants, the bins have increasing widths with increasing leaf size. We will also lump all reproduction into a single flowering class that can include individuals of any size, provided that they have sprouted.

sizevector <- c(0, 100, 13, 127, 3730, 3800, 0)
stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm")
repvector <- c(0, 0, 0, 0, 0, 1, 0)
obsvector <- c(0, 1, 1, 1, 1, 1, 0)
matvector <- c(0, 0, 1, 1, 1, 1, 1)
immvector <- c(1, 1, 0, 0, 0, 0, 0)
propvector <- c(1, 0, 0, 0, 0, 0, 0)
indataset <- c(0, 1, 1, 1, 1, 1, 1)
binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5)

lathframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector, 
  repstatus = repvector, obsstatus = obsvector, matstatus = matvector, 
  immstatus = immvector, indataset = indataset, binhalfwidth = binvec, 
  propstatus = propvector)
lathframe_raw
>   stage size size_b size_c min_age max_age repstatus obsstatus propstatus
> 1    Sd    0     NA     NA      NA      NA         0         0          1
> 2   Sdl  100     NA     NA      NA      NA         0         1          0
> 3   VSm   13     NA     NA      NA      NA         0         1          0
> 4    Sm  127     NA     NA      NA      NA         0         1          0
> 5   VLa 3730     NA     NA      NA      NA         0         1          0
> 6   Flo 3800     NA     NA      NA      NA         1         1          0
> 7  Dorm    0     NA     NA      NA      NA         0         0          0
>   immstatus matstatus indataset binhalfwidth_raw sizebin_min sizebin_max
> 1         1         0         0              0.0         0.0         0.0
> 2         1         0         1            100.0         0.0       200.0
> 3         0         1         1             11.0         2.0        24.0
> 4         0         1         1            103.0        24.0       230.0
> 5         0         1         1           3500.0       230.0      7230.0
> 6         0         1         1           3800.0         0.0      7600.0
> 7         0         1         1              0.5        -0.5         0.5
>   sizebin_center sizebin_width binhalfwidthb_raw sizebinb_min sizebinb_max
> 1              0             0                NA           NA           NA
> 2            100           200                NA           NA           NA
> 3             13            22                NA           NA           NA
> 4            127           206                NA           NA           NA
> 5           3730          7000                NA           NA           NA
> 6           3800          7600                NA           NA           NA
> 7              0             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
>   sizebinc_center sizebinc_width group       comments
> 1              NA             NA     0 No description
> 2              NA             NA     0 No description
> 3              NA             NA     0 No description
> 4              NA             NA     0 No description
> 5              NA             NA     0 No description
> 6              NA             NA     0 No description
> 7              NA             NA     0 No description

Let’s now restandardize our core dataset, using the new stageframe for stage assignment and the original size variable for size classification. We will also cut sightings of individuals in the first year, as before, to result in a dataset in which individuals are ages that we can be reasonably sure about.

lathvert_raw <- verticalize3(lathyrus, noyears = 4, firstyear = 1988, 
  patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9, 
  juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88", 
  fecacol = "Intactseed88", deadacol = "Dead1988", nonobsacol = "Dormant1988", 
  stageassign = lathframe_raw, stagesize = "sizea", censorcol = "Missing1988", 
  censorkeep = NA, censor = TRUE)

lathvert_raw_small <- subset(lathvert_raw, firstseen > 1988)
dim(lathvert_raw_small)
> [1] 542  45
summary(lathvert_raw_small)
>      rowid        popid  patchid   individ              year2     
>  Min.   :  30.0   :542   1:131   Length:542         Min.   :1989  
>  1st Qu.: 302.5          2: 72   Class :character   1st Qu.:1989  
>  Median : 583.0          3:101   Mode  :character   Median :1990  
>  Mean   : 586.8          4:132                      Mean   :1990  
>  3rd Qu.: 801.0          5: 70                      3rd Qu.:1990  
>  Max.   :1097.0          6: 36                      Max.   :1990  
>                                                                   
>    firstseen       lastseen        obsage       obslifespan       sizea1       
>  Min.   :1989   Min.   :1989   Min.   :1.000   Min.   :0.00   Min.   :   0.00  
>  1st Qu.:1989   1st Qu.:1990   1st Qu.:1.000   1st Qu.:1.00   1st Qu.:   0.00  
>  Median :1989   Median :1991   Median :1.000   Median :2.00   Median :   0.00  
>  Mean   :1989   Mean   :1991   Mean   :1.352   Mean   :1.36   Mean   :  59.41  
>  3rd Qu.:1989   3rd Qu.:1991   3rd Qu.:2.000   3rd Qu.:2.00   3rd Qu.:   9.00  
>  Max.   :1990   Max.   :1991   Max.   :2.000   Max.   :2.00   Max.   :4394.20  
>                                                                                
>     repstra1          feca1       fec1added         juvgiven1     
>  Min.   :0.0000   Min.   : 0    Min.   : 0.0000   Min.   :0.0000  
>  1st Qu.:0.0000   1st Qu.: 0    1st Qu.: 0.0000   1st Qu.:0.0000  
>  Median :0.0000   Median : 4    Median : 0.0000   Median :0.0000  
>  Mean   :0.0576   Mean   : 6    Mean   : 0.1218   Mean   :0.1624  
>  3rd Qu.:0.0000   3rd Qu.: 6    3rd Qu.: 0.0000   3rd Qu.:0.0000  
>  Max.   :1.0000   Max.   :34    Max.   :34.0000   Max.   :1.0000  
>  NA's   :351      NA's   :531                                     
>    obsstatus1       repstatus1       fecstatus1        matstatus1  
>  Min.   :0.0000   Min.   :0.0000   Min.   :0.00000   Min.   :0.00  
>  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.00  
>  Median :0.0000   Median :0.0000   Median :0.00000   Median :0.00  
>  Mean   :0.3524   Mean   :0.0203   Mean   :0.01292   Mean   :0.19  
>  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:0.00000   3rd Qu.:0.00  
>  Max.   :1.0000   Max.   :1.0000   Max.   :1.00000   Max.   :1.00  
>                                                                    
>      alive1          stage1           stage1index      sizea2       
>  Min.   :0.0000   Length:542         Min.   :0.0   Min.   :   0.00  
>  1st Qu.:0.0000   Class :character   1st Qu.:0.0   1st Qu.:   9.00  
>  Median :0.0000   Mode  :character   Median :0.0   Median :  12.60  
>  Mean   :0.3524                      Mean   :1.1   Mean   :  97.24  
>  3rd Qu.:1.0000                      3rd Qu.:2.0   3rd Qu.:  37.40  
>  Max.   :1.0000                      Max.   :6.0   Max.   :4394.20  
>                                                                     
>     repstra2           feca2          fec2added         juvgiven2     
>  Min.   :0.00000   Min.   : 0.000   Min.   : 0.0000   Min.   :0.0000  
>  1st Qu.:0.00000   1st Qu.: 0.000   1st Qu.: 0.0000   1st Qu.:0.0000  
>  Median :0.00000   Median : 0.000   Median : 0.0000   Median :0.0000  
>  Mean   :0.03571   Mean   : 4.053   Mean   : 0.1421   Mean   :0.3469  
>  3rd Qu.:0.00000   3rd Qu.: 6.000   3rd Qu.: 0.0000   3rd Qu.:1.0000  
>  Max.   :1.00000   Max.   :34.000   Max.   :34.0000   Max.   :1.0000  
>  NA's   :10        NA's   :523                                        
>    obsstatus2       repstatus2        fecstatus2        matstatus2    
>  Min.   :0.0000   Min.   :0.00000   Min.   :0.00000   Min.   :0.0000  
>  1st Qu.:1.0000   1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.0000  
>  Median :1.0000   Median :0.00000   Median :0.00000   Median :1.0000  
>  Mean   :0.9815   Mean   :0.03506   Mean   :0.01661   Mean   :0.6531  
>  3rd Qu.:1.0000   3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:1.0000  
>  Max.   :1.0000   Max.   :1.00000   Max.   :1.00000   Max.   :1.0000  
>                                                                       
>      alive2     stage2           stage2index        sizea3       
>  Min.   :1   Length:542         Min.   :2.000   Min.   :   0.00  
>  1st Qu.:1   Class :character   1st Qu.:2.000   1st Qu.:   0.00  
>  Median :1   Mode  :character   Median :3.000   Median :  10.50  
>  Mean   :1                      Mean   :3.159   Mean   :  76.21  
>  3rd Qu.:1                      3rd Qu.:4.000   3rd Qu.:  24.50  
>  Max.   :1                      Max.   :7.000   Max.   :6645.80  
>                                                                  
>     repstra3           feca3          fec3added         juvgiven3
>  Min.   :0.00000   Min.   : 0.000   Min.   : 0.0000   Min.   :0  
>  1st Qu.:0.00000   1st Qu.: 0.000   1st Qu.: 0.0000   1st Qu.:0  
>  Median :0.00000   Median : 0.000   Median : 0.0000   Median :0  
>  Mean   :0.05556   Mean   : 6.045   Mean   : 0.2454   Mean   :0  
>  3rd Qu.:0.00000   3rd Qu.: 7.500   3rd Qu.: 0.0000   3rd Qu.:0  
>  Max.   :1.00000   Max.   :48.000   Max.   :48.0000   Max.   :0  
>  NA's   :146       NA's   :520                                   
>    obsstatus3       repstatus3        fecstatus3        matstatus3
>  Min.   :0.0000   Min.   :0.00000   Min.   :0.00000   Min.   :1   
>  1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:1   
>  Median :1.0000   Median :0.00000   Median :0.00000   Median :1   
>  Mean   :0.7306   Mean   :0.04059   Mean   :0.01292   Mean   :1   
>  3rd Qu.:1.0000   3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:1   
>  Max.   :1.0000   Max.   :1.00000   Max.   :1.00000   Max.   :1   
>                                                                   
>      alive3          stage3           stage3index   
>  Min.   :0.0000   Length:542         Min.   :0.000  
>  1st Qu.:1.0000   Class :character   1st Qu.:3.000  
>  Median :1.0000   Mode  :character   Median :3.000  
>  Mean   :0.7989                      Mean   :3.037  
>  3rd Qu.:1.0000                      3rd Qu.:4.000  
>  Max.   :1.0000                      Max.   :7.000  
> 

Once again we have a dataset with 542 rows.

Now let’s build our matrices. We will use the same supplement table as before in ahistorical analysis, since there are no stages used in that supplement table unique to the function-based analysis.

lathmat2p_raw <- arlefko2(data = lathvert_raw_small, stageframe = lathframe_raw,
  supplement = lathsupp2, stages = c("stage3", "stage2", "stage1"),
  patch = "all", patchcol = "patchid", yearcol = "year2", agecol = "obsage",
  indivcol = "individ")

summary(lathmat2p_raw)
> 
> This ahistorical lefkoMat object contains 12 matrices.
> 
> Each matrix is square with 14 rows and columns, and a total of 196 elements.
> A total of 138 survival transitions were estimated, with 11.5 per matrix.
> A total of 12 fecundity transitions were estimated, with 1 per matrix.
> This lefkoMat object covers 1 population, 6 patches, and 2 time steps.
> 
> The dataset contains a total of 232 unique individuals and 542 unique transitions.
> 
> Survival probability sum check (each matrix represented by column in order):
>          [,1]  [,2]  [,3]   [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11]
> Min.    0.000 0.000 0.000 0.0000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
> 1st Qu. 0.000 0.000 0.000 0.0000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
> Median  0.199 0.000 0.199 0.0625 0.199 0.000 0.199 0.000 0.000 0.000 0.000
> Mean    0.393 0.220 0.350 0.2802 0.383 0.320 0.374 0.293 0.306 0.211 0.319
> 3rd Qu. 0.893 0.399 0.688 0.4748 0.842 0.725 0.781 0.577 0.635 0.383 0.600
> Max.    1.000 1.000 1.000 1.0000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
>         [,12]
> Min.    0.000
> 1st Qu. 0.000
> Median  0.000
> Mean    0.184
> 3rd Qu. 0.299
> Max.    1.000

Here we have 12 matrices, each with 14 rows and columns. There are a total of 196 elements per matrix, but only 12.5 are estimated per matrix due to the sparsity of data relative to the size of the matrix. Let’s take a look at the first matrix.

lathmat2p_raw$A[[1]]
>        [,1]      [,2]      [,3]   [,4]  [,5] [,6]  [,7]  [,8] [,9] [,10] [,11]
>  [1,] 0.000 0.0000000 0.0000000 0.0000 0.000    0 1.035 0.000    0     0     0
>  [2,] 0.000 0.0000000 0.0000000 0.0000 0.000    0 0.162 0.000    0     0     0
>  [3,] 0.000 0.0000000 0.0000000 0.0000 0.000    0 0.000 0.000    0     0     0
>  [4,] 0.000 0.0000000 0.0000000 0.0000 0.000    0 0.000 0.000    0     0     0
>  [5,] 0.000 0.0000000 0.0000000 0.0000 0.000    0 0.000 0.000    0     0     0
>  [6,] 0.000 0.0000000 0.0000000 0.0000 0.000    0 0.000 0.000    0     0     0
>  [7,] 0.000 0.0000000 0.0000000 0.0000 0.000    0 0.000 0.000    0     0     0
>  [8,] 0.345 0.0000000 0.0000000 0.0000 0.000    0 0.000 0.345    0     0     0
>  [9,] 0.054 0.0000000 0.0000000 0.0000 0.000    0 0.000 0.054    0     0     0
> [10,] 0.000 0.9047619 0.7142857 0.1875 0.250    0 0.000 0.000    0     0     0
> [11,] 0.000 0.0000000 0.1428571 0.6875 0.500    0 1.000 0.000    0     0     0
> [12,] 0.000 0.0000000 0.0000000 0.0000 0.000    0 0.000 0.000    0     0     0
> [13,] 0.000 0.0000000 0.0000000 0.0625 0.125    0 0.000 0.000    0     0     0
> [14,] 0.000 0.0000000 0.0000000 0.0000 0.125    0 0.000 0.000    0     0     0
>       [,12] [,13] [,14]
>  [1,]     0     0     0
>  [2,]     0     0     0
>  [3,]     0     0     0
>  [4,]     0     0     0
>  [5,]     0     0     0
>  [6,]     0     0     0
>  [7,]     0     0     0
>  [8,]     0     0     0
>  [9,]     0     0     0
> [10,]     0     0     0
> [11,]     0     0     0
> [12,]     0     0     0
> [13,]     0     0     0
> [14,]     0     0     0

One of the more interesting problems that we can see here is the fact that since this matrix reflects a very small dataset, both in terms of numbers of observed transitions and in terms of years of observation, we cannot really estimate transitions staying within the final age. To yield such transitions, we would typically need a large and long enough dataset that we could deliberately reduce the number of ages to a number less than the longest seen, perhaps reflecting stable demographic patterns beyond some age. Clearly, here is one area in which function-based MPMs have a clear advantage. To help increase the numbers of estimated elements, let’s develop an MPM that does not discriminate patches.

lathmat2_raw <- arlefko2(data = lathvert_raw_small, stageframe = lathframe_raw,
  supplement = lathsupp2, stages = c("stage3", "stage2", "stage1"),
  patchcol = "patchid", yearcol = "year2", agecol = "obsage",
  indivcol = "individ")

summary(lathmat2_raw)
> 
> This ahistorical lefkoMat object contains 2 matrices.
> 
> Each matrix is square with 14 rows and columns, and a total of 196 elements.
> A total of 40 survival transitions were estimated, with 20 per matrix.
> A total of 6 fecundity transitions were estimated, with 3 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 2 time steps.
> 
> The dataset contains a total of 232 unique individuals and 542 unique transitions.
> 
> Survival probability sum check (each matrix represented by column in order):
>          [,1]  [,2]
> Min.    0.000 0.000
> 1st Qu. 0.000 0.000
> Median  0.199 0.199
> Mean    0.374 0.356
> 3rd Qu. 0.833 0.758
> Max.    1.000 1.000

Note the increase in estimated elements - roughly twice the number per matrix as last time. Let’s take a look at the first matrix here.

lathmat2_raw$A[[1]]
>        [,1]       [,2]       [,3]       [,4]      [,5] [,6]       [,7]  [,8]
>  [1,] 0.000 0.00000000 0.00000000 0.00000000 0.0000000    0 2.07000000 0.000
>  [2,] 0.000 0.00000000 0.00000000 0.00000000 0.0000000    0 0.32400000 0.000
>  [3,] 0.000 0.00000000 0.00000000 0.00000000 0.0000000    0 0.00000000 0.000
>  [4,] 0.000 0.00000000 0.00000000 0.00000000 0.0000000    0 0.00000000 0.000
>  [5,] 0.000 0.00000000 0.00000000 0.00000000 0.0000000    0 0.00000000 0.000
>  [6,] 0.000 0.00000000 0.00000000 0.00000000 0.0000000    0 0.00000000 0.000
>  [7,] 0.000 0.00000000 0.00000000 0.00000000 0.0000000    0 0.00000000 0.000
>  [8,] 0.345 0.00000000 0.00000000 0.00000000 0.0000000    0 0.00000000 0.345
>  [9,] 0.054 0.00000000 0.00000000 0.00000000 0.0000000    0 0.00000000 0.054
> [10,] 0.000 0.69491525 0.78787879 0.25490196 0.1052632    0 0.00000000 0.000
> [11,] 0.000 0.02542373 0.12121212 0.52941176 0.3157895    0 0.36363636 0.000
> [12,] 0.000 0.00000000 0.00000000 0.00000000 0.2105263    0 0.27272727 0.000
> [13,] 0.000 0.02542373 0.03030303 0.05882353 0.1052632    0 0.09090909 0.000
> [14,] 0.000 0.00000000 0.00000000 0.01960784 0.1578947    0 0.27272727 0.000
>       [,9] [,10] [,11] [,12] [,13] [,14]
>  [1,]    0     0     0     0     0     0
>  [2,]    0     0     0     0     0     0
>  [3,]    0     0     0     0     0     0
>  [4,]    0     0     0     0     0     0
>  [5,]    0     0     0     0     0     0
>  [6,]    0     0     0     0     0     0
>  [7,]    0     0     0     0     0     0
>  [8,]    0     0     0     0     0     0
>  [9,]    0     0     0     0     0     0
> [10,]    0     0     0     0     0     0
> [11,]    0     0     0     0     0     0
> [12,]    0     0     0     0     0     0
> [13,]    0     0     0     0     0     0
> [14,]    0     0     0     0     0     0

Let’s now move on to Leslie MPMs.

6.4 Age-classified (Leslie) MPMs

Package lefko3 can also be used to estimate Leslie MPMs, which are purely age-based and so include no individual history. Here, we will illustrate how to create both raw and function-based Leslie MPMs using the rleslie() and fleslie() functions, respectively. We will ignore the dormant seed stage, as including dormant seeds would require using an age-by-stage approach. Our final matrices will take the following form.

\[\begin{equation} \tag{6.3} \left[\begin{array} {rrrr} F _{0,0} & F _{0,1} & F _{0,2} & F _{0,3} \\ S _{1,0} & 0 & 0 & 0 \\ 0 & S _{2,1} & 0 & 0 \\ 0 & 0 & S _{3,2} & S_{3,3} \\ \end{array}\right] \end{equation}\]

Normally, we would need to develop a stageframe. However, in the purely age-based case, a stageframe is unnecessary and instead all we need is the standardized hfv dataset. So, let’s start by creating that. Most of the settings will be as before. However, we need to include both NRasRep = TRUE and NOasObs = TRUE to make sure that stage classification ignores whether the individual was actually reproductive and actually observed at each time (the latter has to do with the fact that individuals can be alive with a size of zero if they are vegetatively dormant). We will also subset our data to only those individuals whose age we are reasonably sure of, by eliminating those individuals first observed in the first monitoring session.

lathvert_base <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
  patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
  sizeacol = "Volume88", repstracol = "FCODE88", fecacol = "Intactseed88",
  deadacol = "Dead1988", censorcol = "Missing1988", censorkeep = NA,
  censor = TRUE, NRasRep = TRUE, NOasObs = TRUE)
lathvert_age <- subset(lathvert_base, firstseen > 1988)

summary(lathvert_age)
>      rowid        popid  patchid   individ              year2     
>  Min.   :  30.0   :542   1:131   Length:542         Min.   :1989  
>  1st Qu.: 302.5          2: 72   Class :character   1st Qu.:1989  
>  Median : 583.0          3:101   Mode  :character   Median :1990  
>  Mean   : 586.8          4:132                      Mean   :1990  
>  3rd Qu.: 801.0          5: 70                      3rd Qu.:1990  
>  Max.   :1097.0          6: 36                      Max.   :1990  
>                                                                   
>    firstseen       lastseen        obsage       obslifespan       sizea1       
>  Min.   :1989   Min.   :1989   Min.   :1.000   Min.   :0.00   Min.   :   0.00  
>  1st Qu.:1989   1st Qu.:1990   1st Qu.:1.000   1st Qu.:1.00   1st Qu.:   0.00  
>  Median :1989   Median :1991   Median :1.000   Median :2.00   Median :   0.00  
>  Mean   :1989   Mean   :1991   Mean   :1.352   Mean   :1.36   Mean   :  59.41  
>  3rd Qu.:1989   3rd Qu.:1991   3rd Qu.:2.000   3rd Qu.:2.00   3rd Qu.:   9.00  
>  Max.   :1990   Max.   :1991   Max.   :2.000   Max.   :2.00   Max.   :4394.20  
>                                                                                
>     repstra1          feca1       fec1added         juvgiven1   obsstatus1    
>  Min.   :0.0000   Min.   : 0    Min.   : 0.0000   Min.   :0   Min.   :0.0000  
>  1st Qu.:0.0000   1st Qu.: 0    1st Qu.: 0.0000   1st Qu.:0   1st Qu.:0.0000  
>  Median :0.0000   Median : 4    Median : 0.0000   Median :0   Median :0.0000  
>  Mean   :0.0576   Mean   : 6    Mean   : 0.1218   Mean   :0   Mean   :0.3524  
>  3rd Qu.:0.0000   3rd Qu.: 6    3rd Qu.: 0.0000   3rd Qu.:0   3rd Qu.:1.0000  
>  Max.   :1.0000   Max.   :34    Max.   :34.0000   Max.   :0   Max.   :1.0000  
>  NA's   :351      NA's   :531                                                 
>    repstatus1       fecstatus1        matstatus1     alive1      
>  Min.   :0.0000   Min.   :0.00000   Min.   :1    Min.   :0.0000  
>  1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:1    1st Qu.:0.0000  
>  Median :0.0000   Median :0.00000   Median :1    Median :0.0000  
>  Mean   :0.0203   Mean   :0.01292   Mean   :1    Mean   :0.3524  
>  3rd Qu.:0.0000   3rd Qu.:0.00000   3rd Qu.:1    3rd Qu.:1.0000  
>  Max.   :1.0000   Max.   :1.00000   Max.   :1    Max.   :1.0000  
>                                                                  
>     stage1           stage1index     sizea2           repstra2      
>  Length:542         Min.   :0    Min.   :   0.00   Min.   :0.00000  
>  Class :character   1st Qu.:0    1st Qu.:   9.00   1st Qu.:0.00000  
>  Mode  :character   Median :0    Median :  12.60   Median :0.00000  
>                     Mean   :0    Mean   :  97.24   Mean   :0.03571  
>                     3rd Qu.:0    3rd Qu.:  37.40   3rd Qu.:0.00000  
>                     Max.   :0    Max.   :4394.20   Max.   :1.00000  
>                                                    NA's   :10       
>      feca2          fec2added         juvgiven2   obsstatus2    
>  Min.   : 0.000   Min.   : 0.0000   Min.   :0   Min.   :0.0000  
>  1st Qu.: 0.000   1st Qu.: 0.0000   1st Qu.:0   1st Qu.:1.0000  
>  Median : 0.000   Median : 0.0000   Median :0   Median :1.0000  
>  Mean   : 4.053   Mean   : 0.1421   Mean   :0   Mean   :0.9815  
>  3rd Qu.: 6.000   3rd Qu.: 0.0000   3rd Qu.:0   3rd Qu.:1.0000  
>  Max.   :34.000   Max.   :34.0000   Max.   :0   Max.   :1.0000  
>  NA's   :523                                                    
>    repstatus2        fecstatus2        matstatus2     alive2 
>  Min.   :0.00000   Min.   :0.00000   Min.   :1    Min.   :1  
>  1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:1    1st Qu.:1  
>  Median :0.00000   Median :0.00000   Median :1    Median :1  
>  Mean   :0.03506   Mean   :0.01661   Mean   :1    Mean   :1  
>  3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:1    3rd Qu.:1  
>  Max.   :1.00000   Max.   :1.00000   Max.   :1    Max.   :1  
>                                                              
>     stage2           stage2index     sizea3           repstra3      
>  Length:542         Min.   :0    Min.   :   0.00   Min.   :0.00000  
>  Class :character   1st Qu.:0    1st Qu.:   0.00   1st Qu.:0.00000  
>  Mode  :character   Median :0    Median :  10.50   Median :0.00000  
>                     Mean   :0    Mean   :  76.21   Mean   :0.05556  
>                     3rd Qu.:0    3rd Qu.:  24.50   3rd Qu.:0.00000  
>                     Max.   :0    Max.   :6645.80   Max.   :1.00000  
>                                                    NA's   :146      
>      feca3          fec3added         juvgiven3   obsstatus3    
>  Min.   : 0.000   Min.   : 0.0000   Min.   :0   Min.   :0.0000  
>  1st Qu.: 0.000   1st Qu.: 0.0000   1st Qu.:0   1st Qu.:0.0000  
>  Median : 0.000   Median : 0.0000   Median :0   Median :1.0000  
>  Mean   : 6.045   Mean   : 0.2454   Mean   :0   Mean   :0.7306  
>  3rd Qu.: 7.500   3rd Qu.: 0.0000   3rd Qu.:0   3rd Qu.:1.0000  
>  Max.   :48.000   Max.   :48.0000   Max.   :0   Max.   :1.0000  
>  NA's   :520                                                    
>    repstatus3        fecstatus3        matstatus3     alive3      
>  Min.   :0.00000   Min.   :0.00000   Min.   :1    Min.   :0.0000  
>  1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:1    1st Qu.:1.0000  
>  Median :0.00000   Median :0.00000   Median :1    Median :1.0000  
>  Mean   :0.04059   Mean   :0.01292   Mean   :1    Mean   :0.7528  
>  3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:1    3rd Qu.:1.0000  
>  Max.   :1.00000   Max.   :1.00000   Max.   :1    Max.   :1.0000  
>                                                                   
>     stage3           stage3index
>  Length:542         Min.   :0   
>  Class :character   1st Qu.:0   
>  Mode  :character   Median :0   
>                     Mean   :0   
>                     3rd Qu.:0   
>                     Max.   :0   
> 

So far so good! Now let’s create the vital rate models. This will go much more quickly than last time, because we no longer care about size and reproductive status as factors determining vital rates. So, we will set suite = "cons" to prevent these factors from being tested, and set age = "obsage" to incorporate our age at time t variable into all models. We will also use global.only = TRUE to prevent the age term from being dropped (normally we would not set this, but in our case the dataset is so small that age will drop out in the best-fit models, so we will force the models to include age in this case).

lathmodels2_age <- modelsearch(lathvert_age, historical = FALSE,
  approach = "mixed", suite = "cons", bestfit = "AICc&k", age = "obsage",
  vitalrates = c("surv", "fec"), fecdist = "negbin", indiv = "individ",
  year = "year2", year.as.random = TRUE, patch.as.random = TRUE,
  show.model.tables = TRUE, fec.zero = TRUE, global.only = TRUE, quiet = TRUE)

Let’s see a summary of the models.

summary(lathmodels2_age)
> This LefkoMod object includes 2 linear models.
> Best-fit model criterion used: global model only
> 
> 
> 
> Survival model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
>   Approximation) [glmerMod]
>  Family: binomial  ( logit )
> Formula: alive3 ~ obsage + (1 | year2) + (1 | individ)
>    Data: subdata
>       AIC       BIC    logLik  deviance  df.resid 
>  608.1500  625.3311 -300.0750  600.1500       538 
> Random effects:
>  Groups  Name        Std.Dev. 
>  individ (Intercept) 0.0002353
>  year2   (Intercept) 0.3738235
> Number of obs: 542, groups:  individ, 232; year2, 2
> Fixed Effects:
> (Intercept)       obsage  
>      0.9445       0.1817  
> 
> 
> 
> Observation model:
> [1] 1
> 
> 
> 
> Size model:
> [1] 1
> 
> 
> 
> Secondary size model:
> [1] 1
> 
> 
> 
> Tertiary size model:
> [1] 1
> 
> 
> 
> Reproductive status model:
> [1] 1
> 
> 
> 
> Fecundity model:
> Formula:          feca2 ~ obsage + (1 | year2) + (1 | individ)
> Zero inflation:         ~.
> Data: subdata
>      AIC      BIC   logLik df.resid 
>       NA       NA       NA       10 
> Random-effects (co)variances:
> 
> Conditional model:
>  Groups  Name        Std.Dev.
>  year2   (Intercept) 0.02271 
>  individ (Intercept) 1.19287 
> 
> Zero-inflation model:
>  Groups  Name        Std.Dev.
>  year2   (Intercept) 0.09092 
>  individ (Intercept) 0.17023 
> 
> Number of obs: 19 / Conditional model: year2, 2; individ, 16 / Zero-inflation model: year2, 2; individ, 16
> 
> Dispersion parameter for nbinom2 family (): 88.6 
> 
> Fixed Effects:
> 
> Conditional model:
> (Intercept)       obsage  
>      1.3220       0.3938  
> 
> Zero-inflation model:
> (Intercept)       obsage  
>      -3.842        2.635  
> 
> 
> 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: 1
> 
> Number of models in observation table: 1
> 
> Number of models in size table: 1
> 
> Number of models in secondary size table: 1
> 
> Number of models in tertiary size table: 1
> 
> Number of models in reproduction status table: 1
> 
> Number of models in fecundity table: 1
> 
> 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 232 individuals and 542 individual transitions.
> Survival accuracy is 0.753.
> Observation probability not estimated.
> Primary size transition not estimated.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproduction probability not estimated.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Fecundity pseudo R-squared is NA.
> 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 two models all include observed age, as well as our random terms. Now we will create the raw matrices.

lathmat2ageonly_raw <- rleslie(data = lathvert_age, year = "all",
  yearcol = "year2", indivcol = "individ")

summary(lathmat2ageonly_raw)
> 
> This ahistorical lefkoMat object contains 1 matrix.
> 
> Each matrix is square with 3 rows and columns, and a total of 9 elements.
> A total of 2 survival transitions were estimated, with 2 per matrix.
> A total of 2 fecundity transitions were estimated, with 2 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 1 time step.
> 
> The dataset contains a total of 232 unique individuals and 542 unique transitions.
> 
> Survival probability sum check (each matrix represented by column in order):
>          [,1]
> Min.    0.000
> 1st Qu. 0.332
> Median  0.664
> Mean    0.462
> 3rd Qu. 0.693
> Max.    0.723

We have just created the raw Leslie MPM. We have a single matrix, because we needed to eliminate the first year of data to estimate age properly, leaving only three observation periods. The quality control looks reasonable in terms of the survival-transition probabilities falling within the range of 0.0 to 1.0. Let’s also take a look at our one matrix, to get a handle on its structure.

lathmat2ageonly_raw$A[[1]]
>            [,1]       [,2] [,3]
> [1,] 0.02521008 0.04188482    0
> [2,] 0.66386555 0.00000000    0
> [3,] 0.00000000 0.72251309    0

Our matrix certainly looks like a valid Leslie matrix, with survival transitions in the subdiagonal and fecundity at the top. The last column is entirely composed of zeros, and that is likely due to a combination of a small dataset and a short dataset in which we cannot find individuals whom we know are older than three years old.

Let’s now build the function-based versions and see what happens.

lathmat2ageonly_func <- fleslie(year = "all", data = lathvert_age,
  modelsuite = lathmodels2_age)

summary(lathmat2ageonly_func)
> 
> This ahistorical lefkoMat object contains 2 matrices.
> 
> Each matrix is square with 3 rows and columns, and a total of 9 elements.
> A total of 6 survival transitions were estimated, with 3 per matrix.
> A total of 6 fecundity transitions were estimated, with 3 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 2 time steps.
> 
> Vital rate modeling quality control:
> 
> Survival estimated with 232 individuals and 542 individual transitions.
> Observation probability not estimated.
> Primary size transition not estimated.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproduction probability not estimated.
> Fecundity estimated with 16 individuals and 19 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]
> Min.    0.813 0.685
> 1st Qu. 0.826 0.704
> Median  0.839 0.723
> Mean    0.838 0.722
> 3rd Qu. 0.850 0.740
> Max.    0.862 0.758
lathmat2ageonly_func$A[[1]]
>            [,1]      [,2]      [,3]
> [1,] 0.02118612 0.2318637 0.8080455
> [2,] 0.81273691 0.0000000 0.0000000
> [3,] 0.00000000 0.8388330 0.8619104

Notice that we have one more matrix here (reflecting the two estimable years), and we have estimates survival transition probabilities staying within the last age. Clearly, function-based approaches have an advantage in this small, short dataset.

Everything looks quite good! Congratulations on creating your Leslie MPMs!

6.5 Points to remember

  1. Age-based (Leslie) MPMs can be be built in raw or function-based form using functions rleslie() and fleslie(), respectively. Age-by-stage MPMs can be built ion raw or function-based forms using functions arlefko2() and aflefko2(), respectively.
  2. Age-by-stage MPMs are larger than either age-based or stage-based MPMs, and so generally require more data to parameterize properly.
  3. Dot all possible ages need to be modeled in age-based and age-by-stage MPMs. Often, the final age in these MPMs represents a long stretch of the adult span of life during which vital rates are not expected to change dramatically, and so is capable of self-transition while other ages are not.

References

Broek, J. van den. (1995). A score test for zero inflation in a Poisson distribution. Biometrics, 51, 738–743.
Caswell, H. (2001). Matrix population models: Construction, analysis, and interpretation. Second edition. Sinauer Associates, Inc., Sunderland, Massachusetts, USA.
Caswell, H. & Salguero-Gómez, R. (2013). Age, stage and senescence in plants. Journal of Ecology, 101, 585–595.
Caswell, H., Vries, C. de, Hartemink, N., Roth, G. & Daalen, S.F. van. (2018). Age × stage-classified demographic analysis: A comprehensive approach. Ecological Monographs, 88, 560–584.
Ehrlén, J. (2000). The dynamics of plant populations: Does the history of individuals matter? Ecology, 81, 1675–1684.
Lefkovitch, L.P. (1965). The study of population growth in organisms grouped by stages. Biometrics, 21, 1–18.
Leslie, P.H. (1945). On the use of matrices in certain population mathematics. Biometrika, 33, 183–212.