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 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 in population dynamics using age-by-stage MPMs, which are MPMs that incorporate both 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, and also showing how to build a standard Leslie 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 like an albatross. 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} \left[\begin{array} {rrr} 0 & 0 & F _{1,3} \\ S _{2,1} & 0 & 0 \\ 0 & S _{3,2} & S _{3,3} \\ \end{array}\right] \tag{6.1} \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 include 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. These conditions lead to the following age-by-stage matrix.

\[\begin{equation} \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] \tag{6.2} \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 be added in the future.

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 and a summary.

data(lathyrus)

dim(lathyrus)
> [1] 1119   38
summary(lathyrus)
>     SUBPLOT          GENET          Volume88         lnVol88     
>  Min.   :1.000   Min.   :  1.0   Min.   :   3.4   Min.   :1.200  
>  1st Qu.:2.000   1st Qu.: 48.0   1st Qu.:  63.0   1st Qu.:4.100  
>  Median :3.000   Median : 97.0   Median : 732.5   Median :6.600  
>  Mean   :3.223   Mean   :110.2   Mean   : 749.4   Mean   :5.538  
>  3rd Qu.:4.000   3rd Qu.:167.5   3rd Qu.:1025.5   3rd Qu.:6.900  
>  Max.   :6.000   Max.   :284.0   Max.   :7032.0   Max.   :8.900  
>                                  NA's   :404      NA's   :404    
>     FCODE88           Flow88       Intactseed88 Dead1988       Dormant1988   
>  Min.   :0.0000   Min.   : 1.00   Min.   : 0    Mode:logical   Mode:logical  
>  1st Qu.:0.0000   1st Qu.: 4.00   1st Qu.: 0    NA's:1119      NA's:1119     
>  Median :0.0000   Median : 8.00   Median : 0                                 
>  Mean   :0.3399   Mean   :11.86   Mean   : 3                                 
>  3rd Qu.:1.0000   3rd Qu.:15.00   3rd Qu.: 4                                 
>  Max.   :1.0000   Max.   :66.00   Max.   :34                                 
>  NA's   :404      NA's   :910     NA's   :875                                
>  Missing1988     Seedling1988      Volume89         lnVol89     
>  Mode:logical   Min.   :1.000   Min.   :   1.8   Min.   :0.600  
>  NA's:1119      1st Qu.:2.000   1st Qu.:  15.6   1st Qu.:2.700  
>                 Median :2.000   Median : 118.8   Median :4.800  
>                 Mean   :2.144   Mean   : 573.3   Mean   :4.855  
>                 3rd Qu.:3.000   3rd Qu.: 968.8   3rd Qu.:6.900  
>                 Max.   :3.000   Max.   :6539.4   Max.   :8.800  
>                 NA's   :1022    NA's   :294      NA's   :294    
>     FCODE89           Flow89       Intactseed89       Dead1989   
>  Min.   :0.0000   Min.   : 1.00   Min.   : 0.000   Min.   :1     
>  1st Qu.:0.0000   1st Qu.: 5.00   1st Qu.: 0.000   1st Qu.:1     
>  Median :0.0000   Median :11.00   Median : 5.000   Median :1     
>  Mean   :0.2667   Mean   :14.88   Mean   : 8.273   Mean   :1     
>  3rd Qu.:1.0000   3rd Qu.:20.00   3rd Qu.:13.000   3rd Qu.:1     
>  Max.   :1.0000   Max.   :97.00   Max.   :66.000   Max.   :1     
>  NA's   :294      NA's   :906     NA's   :899      NA's   :1077  
>   Dormant1989    Missing1989    Seedling1989      Volume90         lnVol90     
>  Min.   :1      Min.   :1      Min.   :1.000   Min.   :   2.1   Min.   :0.700  
>  1st Qu.:1      1st Qu.:1      1st Qu.:2.000   1st Qu.:  12.6   1st Qu.:2.500  
>  Median :1      Median :1      Median :2.000   Median :  61.0   Median :4.100  
>  Mean   :1      Mean   :1      Mean   :2.136   Mean   : 244.1   Mean   :4.207  
>  3rd Qu.:1      3rd Qu.:1      3rd Qu.:2.000   3rd Qu.: 295.2   3rd Qu.:5.700  
>  Max.   :1      Max.   :1      Max.   :3.000   Max.   :4242.8   Max.   :8.400  
>  NA's   :1046   NA's   :1112   NA's   :1001    NA's   :245      NA's   :245    
>     FCODE90           Flow90        Intactseed90       Dead1990   
>  Min.   :0.0000   Min.   : 1.000   Min.   : 0.000   Min.   :1     
>  1st Qu.:0.0000   1st Qu.: 3.000   1st Qu.: 0.000   1st Qu.:1     
>  Median :0.0000   Median : 6.000   Median : 0.000   Median :1     
>  Mean   :0.1581   Mean   : 8.104   Mean   : 2.514   Mean   :1     
>  3rd Qu.:0.0000   3rd Qu.:10.750   3rd Qu.: 1.000   3rd Qu.:1     
>  Max.   :1.0000   Max.   :54.000   Max.   :37.000   Max.   :1     
>  NA's   :246      NA's   :985      NA's   :981      NA's   :1007  
>   Dormant1990    Missing1990    Seedling1990      Volume91         lnVol91     
>  Min.   :1      Min.   :1      Min.   :1.000   Min.   :   4.0   Min.   :1.400  
>  1st Qu.:1      1st Qu.:1      1st Qu.:2.000   1st Qu.:  12.0   1st Qu.:2.500  
>  Median :1      Median :1      Median :2.000   Median : 118.5   Median :4.800  
>  Mean   :1      Mean   :1      Mean   :2.186   Mean   : 418.7   Mean   :4.642  
>  3rd Qu.:1      3rd Qu.:1      3rd Qu.:2.000   3rd Qu.: 689.7   3rd Qu.:6.500  
>  Max.   :1      Max.   :1      Max.   :3.000   Max.   :6645.8   Max.   :8.800  
>  NA's   :1054   NA's   :1105   NA's   :1049    NA's   :305      NA's   :305    
>     FCODE91           Flow91       Intactseed91       Dead1991    Dormant1991  
>  Min.   :0.0000   Min.   : 1.00   Min.   : 0.000   Min.   :1     Min.   :1     
>  1st Qu.:0.0000   1st Qu.: 4.00   1st Qu.: 0.000   1st Qu.:1     1st Qu.:1     
>  Median :0.0000   Median : 8.00   Median : 3.500   Median :1     Median :1     
>  Mean   :0.2525   Mean   :11.12   Mean   : 5.805   Mean   :1     Mean   :1     
>  3rd Qu.:1.0000   3rd Qu.:15.00   3rd Qu.:10.000   3rd Qu.:1     3rd Qu.:1     
>  Max.   :1.0000   Max.   :48.00   Max.   :48.000   Max.   :1     Max.   :1     
>  NA's   :307      NA's   :954     NA's   :919      NA's   :925   NA's   :1034  
>   Missing1991    Seedling1991  
>  Min.   :1      Min.   :1.000  
>  1st Qu.:1      1st Qu.:2.000  
>  Median :1      Median :2.000  
>  Mean   :1      Mean   :1.973  
>  3rd Qu.:1      3rd Qu.:2.000  
>  Max.   :1      Max.   :3.000  
>  NA's   :1095   NA's   :1082

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. Before doing this, we need to alter the dataset slightly. Currently, the variable GENET lists the individual number, but only within each subpopulation. We wish to identify each individual within each subpopulation uniquely, and this requires us to develop a new ID variable. For this purpose, we will create a new variable that concatenates the subpopulation and genet number together into one string, as below.

lathyrus$indiv_id <- paste(lathyrus$SUBPLOT, lathyrus$GENET)

Now let’s 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 = "indiv_id", blocksize = 9,
  juvcol = "Seedling1988", sizeacol = "lnVol88", repstracol = "FCODE88",
  fecacol = "Intactseed88", deadacol = "Dead1988", nonobsacol = "Dormant1988",
  stageassign = lathframeln, stagesize = "sizea", censorcol = "Missing1988",
  censorkeep = NA, censorRepeat = TRUE, NAas0 = TRUE, censor = TRUE)

summary_hfv(lathvertln)
> 
> This hfv dataset contains 2527 rows, 42 variables, 1 population, 
> 6 patches, 1053 individuals, and 3 time steps.
>      rowid        popid   patchid   individ              year2     
>  Min.   :   1.0   :2527   1:673   Length:2527        Min.   :1988  
>  1st Qu.: 237.0           2:202   Class :character   1st Qu.:1988  
>  Median : 522.0           3:556   Mode  :character   Median :1989  
>  Mean   : 537.3           4:469                      Mean   :1989  
>  3rd Qu.: 820.5           5:308                      3rd Qu.:1990  
>  Max.   :1118.0           6:319                      Max.   :1990  
>    firstseen       lastseen        obsage       obslifespan        sizea1     
>  Min.   :   0   Min.   :   0   Min.   :1.000   Min.   :0.000   Min.   :0.000  
>  1st Qu.:1988   1st Qu.:1991   1st Qu.:1.000   1st Qu.:2.000   1st Qu.:0.000  
>  Median :1988   Median :1991   Median :2.000   Median :3.000   Median :2.200  
>  Mean   :1979   Mean   :1981   Mean   :1.822   Mean   :2.437   Mean   :2.957  
>  3rd Qu.:1988   3rd Qu.:1991   3rd Qu.:2.000   3rd Qu.:3.000   3rd Qu.:6.400  
>  Max.   :1990   Max.   :1991   Max.   :3.000   Max.   :3.000   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.1805   Mean   : 0.9889   Mean   :0.06292   Mean   :0.5548  
>  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.1805   Mean   :0.08983   Mean   :0.5204   Mean   :0.5833  
>  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:2527        Min.   : 0.000   Min.   :0.000   Min.   :0.000  
>  Class :character   1st Qu.: 0.000   1st Qu.:2.500   1st Qu.:0.000  
>  Mode  :character   Median : 3.000   Median :4.700   Median :0.000  
>                     Mean   : 6.137   Mean   :4.565   Mean   :0.237  
>                     3rd Qu.:10.000   3rd Qu.:6.600   3rd Qu.:0.000  
>                     Max.   :21.000   Max.   :8.900   Max.   :1.000  
>      feca2         juvgiven2        obsstatus2       repstatus2   
>  Min.   : 0.00   Min.   :0.0000   Min.   :0.0000   Min.   :0.000  
>  1st Qu.: 0.00   1st Qu.:0.0000   1st Qu.:1.0000   1st Qu.:0.000  
>  Median : 0.00   Median :0.0000   Median :1.0000   Median :0.000  
>  Mean   : 1.14   Mean   :0.1112   Mean   :0.9454   Mean   :0.237  
>  3rd Qu.: 0.00   3rd Qu.:0.0000   3rd Qu.:1.0000   3rd Qu.:0.000  
>  Max.   :66.00   Max.   :1.0000   Max.   :1.0000   Max.   :1.000  
>    fecstatus2       matstatus2         alive2     stage2         
>  Min.   :0.0000   Min.   :0.0000   Min.   :1   Length:2527       
>  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.1053   Mean   :0.8888   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.000   Min.   : 0.00   Min.   :0  
>  1st Qu.: 5.000   1st Qu.:2.300   1st Qu.:0.000   1st Qu.: 0.00   1st Qu.:0  
>  Median : 8.000   Median :4.300   Median :0.000   Median : 0.00   Median :0  
>  Mean   : 9.323   Mean   :4.043   Mean   :0.218   Mean   : 1.29   Mean   :0  
>  3rd Qu.:10.000   3rd Qu.:6.300   3rd Qu.:0.000   3rd Qu.: 0.00   3rd Qu.:0  
>  Max.   :21.000   Max.   :8.800   Max.   :1.000   Max.   :66.00   Max.   :0  
>    obsstatus3       repstatus3      fecstatus3       matstatus3
>  Min.   :0.0000   Min.   :0.000   Min.   :0.0000   Min.   :1   
>  1st Qu.:1.0000   1st Qu.:0.000   1st Qu.:0.0000   1st Qu.:1   
>  Median :1.0000   Median :0.000   Median :0.0000   Median :1   
>  Mean   :0.8346   Mean   :0.218   Mean   :0.1116   Mean   :1   
>  3rd Qu.:1.0000   3rd Qu.:0.000   3rd Qu.:0.0000   3rd Qu.:1   
>  Max.   :1.0000   Max.   :1.000   Max.   :1.0000   Max.   :1   
>      alive3          stage3           stage3index    
>  Min.   :0.0000   Length:2527        Min.   : 0.000  
>  1st Qu.:1.0000   Class :character   1st Qu.: 5.000  
>  Median :1.0000   Mode  :character   Median : 7.000  
>  Mean   :0.9224                      Mean   : 8.717  
>  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 1053 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. Lathyrus is a long-lived plant, and this moniroting study lasted only four years, leading to a poor understanding of the ages of most plants. 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)

summary_hfv(lathvertln_small)
> 
> This hfv dataset contains 531 rows, 42 variables, 1 population, 
> 6 patches, 345 individuals, and 2 time steps.
>      rowid        popid  patchid   individ              year2     
>  Min.   :  30.0   :531   1:130   Length:531         Min.   :1989  
>  1st Qu.: 301.0          2: 69   Class :character   1st Qu.:1989  
>  Median : 582.0          3:101   Mode  :character   Median :1990  
>  Mean   : 585.3          4:127                      Mean   :1990  
>  3rd Qu.: 800.5          5: 69                      3rd Qu.:1990  
>  Max.   :1097.0          6: 35                      Max.   :1990  
>    firstseen       lastseen        obsage      obslifespan        sizea1     
>  Min.   :1989   Min.   :1989   Min.   :1.00   Min.   :0.000   Min.   :0.000  
>  1st Qu.:1989   1st Qu.:1990   1st Qu.:1.00   1st Qu.:1.000   1st Qu.:0.000  
>  Median :1989   Median :1991   Median :1.00   Median :2.000   Median :0.000  
>  Mean   :1989   Mean   :1991   Mean   :1.35   Mean   :1.367   Mean   :1.175  
>  3rd Qu.:1989   3rd Qu.:1991   3rd Qu.:2.00   3rd Qu.:2.000   3rd Qu.:2.200  
>  Max.   :1990   Max.   :1991   Max.   :2.00   Max.   :2.000   Max.   :8.400  
>     repstra1           feca1           juvgiven1       obsstatus1    
>  Min.   :0.00000   Min.   : 0.0000   Min.   :0.000   Min.   :0.0000  
>  1st Qu.:0.00000   1st Qu.: 0.0000   1st Qu.:0.000   1st Qu.:0.0000  
>  Median :0.00000   Median : 0.0000   Median :0.000   Median :0.0000  
>  Mean   :0.02072   Mean   : 0.1243   Mean   :0.162   Mean   :0.3503  
>  3rd Qu.:0.00000   3rd Qu.: 0.0000   3rd Qu.:0.000   3rd Qu.:1.0000  
>  Max.   :1.00000   Max.   :34.0000   Max.   :1.000   Max.   :1.0000  
>    repstatus1        fecstatus1        matstatus1         alive1      
>  Min.   :0.00000   Min.   :0.00000   Min.   :0.0000   Min.   :0.0000  
>  1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.0000  
>  Median :0.00000   Median :0.00000   Median :0.0000   Median :0.0000  
>  Mean   :0.02072   Mean   :0.01318   Mean   :0.1883   Mean   :0.3503  
>  3rd Qu.:0.00000   3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:1.0000  
>  Max.   :1.00000   Max.   :1.00000   Max.   :1.0000   Max.   :1.0000  
>     stage1           stage1index         sizea2         repstra2      
>  Length:531         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.881   Mean   :3.072   Mean   :0.03578  
>                     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.000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00000  
>  1st Qu.: 0.000   1st Qu.:0.0000   1st Qu.:1.0000   1st Qu.:0.00000  
>  Median : 0.000   Median :0.0000   Median :1.0000   Median :0.00000  
>  Mean   : 0.145   Mean   :0.3465   Mean   :0.9849   Mean   :0.03578  
>  3rd Qu.: 0.000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:0.00000  
>  Max.   :34.000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00000  
>    fecstatus2        matstatus2         alive2     stage2         
>  Min.   :0.00000   Min.   :0.0000   Min.   :1   Length:531        
>  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.01695   Mean   :0.6535   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.000   Min.   :0.000   Min.   :0.00000   Min.   : 0.0000  
>  1st Qu.: 2.000   1st Qu.:0.000   1st Qu.:0.00000   1st Qu.: 0.0000  
>  Median : 5.000   Median :2.400   Median :0.00000   Median : 0.0000  
>  Mean   : 5.209   Mean   :2.351   Mean   :0.04143   Mean   : 0.2505  
>  3rd Qu.: 7.000   3rd Qu.:3.200   3rd Qu.:0.00000   3rd Qu.: 0.0000  
>  Max.   :20.000   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.7458   Mean   :0.04143   Mean   :0.01318   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:531         Min.   : 0.00  
>  1st Qu.:1.0000   Class :character   1st Qu.: 3.00  
>  Median :1.0000   Mode  :character   Median : 5.00  
>  Mean   :0.8117                      Mean   : 5.06  
>  3rd Qu.:1.0000                      3rd Qu.: 6.00  
>  Max.   :1.0000                      Max.   :21.00

We are now down to only 531 rows and 345 individuals, 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 zeros 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 zeros, as well. Let’s now formally test our assumptions, that the mean equals the variance and that the number of zeros meets expectation. Both tests use chi-squared distribution-based approaches, with the zero-inflation test in particular based on Broek (1995). We will do this through a quality control assessment of the entire dataset for modeling, using hfv_qc(), into which we need to specify which vital rates we intend to model later, and the names of some key variables..

hfv_qc(lathvertln_small, indiv = "individ", year = "year2", age = "obsage",
  vitalrates = c("surv", "obs", "size", "repst", "fec"))
> Survival:
> 
>   Data subset has 43 variables and 531 transitions.
> 
>   Variable alive3 has 0 missing values.
>   Variable alive3 is a binomial variable.
> 
> 
> Observation status:
> 
>   Data subset has 43 variables and 431 transitions.
> 
>   Variable obsstatus3 has 0 missing values.
>   Variable obsstatus3 is a binomial variable.
> 
> 
> Primary size:
> 
>   Data subset has 43 variables and 396 transitions.
> 
>   Variable sizea3 has 0 missing values.
>   Variable sizea3 appears to be a floating point variable.
>   379 elements are not integers.
>   The minimum value of sizea3 is 0.7 and the maximum is 8.8.
>   The mean value of sizea3 is 3.153 and the variance is 1.818.
>   The value of the Shapiro-Wilk test of normality is 0.8367 with P = 8.608e-20.
>   Variable sizea3 differs significantly from a Gaussian distribution.
> 
>   Variable sizea3 is fully positive, lacking even 0s.
> 
> 
> Reproductive status:
> 
>   Data subset has 43 variables and 396 transitions.
> 
>   Variable repstatus3 has 0 missing values.
>   Variable repstatus3 is a binomial variable.
> 
> 
> Fecundity:
> 
>   Data subset has 43 variables and 19 transitions.
> 
>   Variable feca2 has 0 missing values.
>   Variable feca2 appears to be an integer variable.
> 
>   Variable feca2 is fully non-negative.
> 
>   Overdispersion test:
>     Mean feca2 is 4.053
>     The variance in feca2 is 61.05
>     The probability of this dispersion level by chance assuming that
>     the true mean feca2 = variance in feca2,
>     and an alternative hypothesis of overdispersion, is 8.249e-152
>     Variable feca2 is significantly overdispersed.
> 
>   Zero-inflation and truncation tests:
>     Mean lambda in feca2 is 0.01738
>     The actual number of 0s in feca2 is 10
>     The expected number of 0s in feca2 under the null hypothesis is 0.3302
>     The probability of this deviation in 0s from expectation by chance is 1.719e-69
>     Variable feca2 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. Note that our output also suggests problems with using the Gaussian distribution for size, but we will ignore that this time.

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 + (1 | year2) + (1 | individ)
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
> 182.9812 198.3785 -87.4906 174.9812      343 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 62.23   
>  year2   (Intercept)  0.00   
> Number of obs: 347, groups:  individ, 247; year2, 2
> Fixed Effects:
> (Intercept)       obsage  
>       39.46       -13.65  
> 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 + repstatus2 + sizea2 + (1 | year2) + (1 |  
>     individ) + obsage:repstatus2 + obsage:sizea2
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
> 128.8981 158.3667 -56.4490 112.8981      286 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 63.3445 
>  year2   (Intercept)  0.1541 
> Number of obs: 294, groups:  individ, 203; year2, 2
> Fixed Effects:
>       (Intercept)             obsage         repstatus2             sizea2  
>            35.382              3.959             23.704              2.048  
> obsage:repstatus2      obsage:sizea2  
>           -18.983             -4.057  
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 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) 5.346e-05
>  year2    (Intercept) 4.826e-01
>  Residual             8.850e-01
> Number of obs: 266, groups:  individ, 191; 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 + (1 | year2) + (1 | individ)
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
>  71.3773  85.7112 -31.6886  63.3773      262 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 83.68658
>  year2   (Intercept)  0.06152
> Number of obs: 266, groups:  individ, 191; year2, 2
> Fixed Effects:
> (Intercept)       obsage  
>      -44.60        15.32  
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 2 lme4 warnings 
> 
> 
> 
> Fecundity model:
> Formula:          
> feca2 ~ obsage + repstatus2 + sizea2 + (1 | year2) + (1 | individ) +  
>     obsage:repstatus2
> Zero inflation:         ~obsage + sizea2 + (1 | year2) + (1 | individ)
> Data: subdata
>       AIC       BIC    logLik  df.resid 
> -2765.266 -2752.988  1395.633         6 
> Random-effects (co)variances:
> 
> Conditional model:
>  Groups  Name        Std.Dev. 
>  year2   (Intercept) 1.907e-10
>  individ (Intercept) 1.258e-06
> 
> Zero-inflation model:
>  Groups  Name        Std.Dev. 
>  year2   (Intercept) 4.337e-12
>  individ (Intercept) 2.905e+01
> 
> Number of obs: 19 / Conditional model: year2, 2; individ, 16 / Zero-inflation model: year2, 2; individ, 16
> 
> Dispersion parameter for nbinom2 family (): 3.87e+16 
> 
> Fixed Effects:
> 
> Conditional model:
>       (Intercept)             obsage         repstatus2             sizea2  
>           -3.3689             0.2283            -3.3689             1.1553  
> obsage:repstatus2  
>            0.2283  
> 
> Zero-inflation model:
> (Intercept)       obsage       sizea2  
>      -9.468       15.147       -1.786  
> 
> 
> 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 
>  215.1077  224.7525 -104.5539  209.1077       181 
> Random effects:
>  Groups  Name        Std.Dev. 
>  individ (Intercept) 6.015e-05
>  year2   (Intercept) 0.000e+00
> Number of obs: 184, groups:  individ, 184; year2, 2
> Fixed Effects:
> (Intercept)  
>        1.07  
> 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.1925  36.9524 -11.0962  22.1925      134 
> Random effects:
>  Groups  Name        Std.Dev. 
>  individ (Intercept) 5.587e+01
>  year2   (Intercept) 1.151e-05
> Number of obs: 137, groups:  individ, 137; year2, 2
> Fixed Effects:
> (Intercept)  
>       13.77  
> 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)
>    Data: subdata
> REML criterion at convergence: 145.4426
> Random effects:
>  Groups   Name        Std.Dev.
>  year2    (Intercept) 0.1430  
>  Residual             0.4139  
> Number of obs: 130, groups:  year2, 2
> Fixed Effects:
> (Intercept)  
>       2.288  
> 
> 
> 
> 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: 17
> 
> Number of models in observation table: 16
> 
> 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: 16
> 
> Number of models in fecundity table: 283
> 
> Number of models in juvenile survival table: 1
> 
> Number of models in juvenile observation table: 1
> 
> Number of models in juvenile size table: 1
> 
> Number of models in juvenile secondary size table: 1
> 
> Number of models in juvenile tertiary size table: 1
> 
> Number of models in juvenile reproduction table: 1
> 
> Number of models in juvenile maturity table: 1
> 
> 
> 
> 
> 
> General model parameter names (column 1), and 
> specific names used in these models (column 2): 
>                       parameter_names mainparams
> 1                              time t      year2
> 2                          individual    individ
> 3                               patch      patch
> 4                   alive in time t+1      surv3
> 5                observed in time t+1       obs3
> 6                   sizea in time t+1      size3
> 7                   sizeb in time t+1     sizeb3
> 8                   sizec in time t+1     sizec3
> 9     reproductive status in time t+1     repst3
> 10              fecundity in time t+1       fec3
> 11                fecundity in time t       fec2
> 12                    sizea in time t      size2
> 13                  sizea in time t-1      size1
> 14                    sizeb in time t     sizeb2
> 15                  sizeb in time t-1     sizeb1
> 16                    sizec in time t     sizec2
> 17                  sizec in time t-1     sizec1
> 18      reproductive status in time t     repst2
> 19    reproductive status in time t-1     repst1
> 20        maturity status in time t+1     matst3
> 21          maturity status in time t     matst2
> 22                      age in time t        age
> 23                  density in time t    density
> 24   individual covariate a in time t   indcova2
> 25 individual covariate a in time t-1   indcova1
> 26   individual covariate b in time t   indcovb2
> 27 individual covariate b in time t-1   indcovb1
> 28   individual covariate c in time t   indcovc2
> 29 individual covariate c in time t-1   indcovc1
> 30              stage group in time t     group2
> 31            stage group in time t-1     group1
> 
> 
> 
> 
> 
> Quality control:
> 
> Survival model estimated with 247 individuals and 347 individual transitions.
> Survival model accuracy is 1.
> Observation status model estimated with 203 individuals and 294 individual transitions.
> Observation status model accuracy is 1.
> Primary size model estimated with 191 individuals and 266 individual transitions.
> Primary size model R-squared is 0.632.
> Secondary size model not estimated.
> Tertiary size model not estimated.
> Reproductive status model estimated with 191 individuals and 266 individual transitions.
> Reproductive status model accuracy is 1.
> Fecundity model estimated with 16 individuals and 19 individual transitions.
> Fecundity model R-squared is 0.965.
> Juvenile survival model estimated with 184 individuals and 184 individual transitions.
> Juvenile survival model accuracy is 0.745.
> Juvenile observation status model estimated with 137 individuals and 137 individual transitions.
> Juvenile observation status model accuracy is 1.
> Juvenile primary size model estimated with 130 individuals and 130 individual transitions.
> Juvenile primary size model R-squared is 0.058.
> Juvenile secondary size model not estimated.
> Juvenile tertiary size model not estimated.
> Juvenile reproductive status model not estimated.
> Juvenile maturity status model not estimated.

We see that age is influential in survival, observation, reproductive status, and fecundity, though not in other vital rates. Accuracy in our adult binomial models is high, but primary size for adults and juveniles is explained by models with low R2, suggesting problems.

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 + sizea2 + (1 | year2) + (1 | patchid) +  
>     (1 | individ) + obsage:sizea2
>    Data: subdata
>       AIC       BIC    logLik  deviance  df.resid 
>  288.3454  319.1400 -136.1727  272.3454       339 
> Random effects:
>  Groups  Name        Std.Dev. 
>  individ (Intercept) 7.285e-01
>  patchid (Intercept) 5.452e-01
>  year2   (Intercept) 3.494e-05
> Number of obs: 347, groups:  individ, 247; patchid, 6; year2, 2
> Fixed Effects:
>   (Intercept)         obsage     repstatus2         sizea2  obsage:sizea2  
>        6.0406        -2.8058        20.6037        -0.8776         0.6609  
> 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 + repstatus2 + sizea2 + (1 | year2) + (1 |  
>     patchid) + (1 | individ) + obsage:sizea2 + repstatus2:sizea2
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
> 131.4644 164.6166 -56.7322 113.4644      285 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 43.8384 
>  patchid (Intercept)  0.9045 
>  year2   (Intercept)  0.0000 
> Number of obs: 294, groups:  individ, 203; patchid, 6; year2, 2
> Fixed Effects:
>       (Intercept)             obsage         repstatus2             sizea2  
>            21.103              8.302             41.223              4.327  
>     obsage:sizea2  repstatus2:sizea2  
>            -4.843             -6.055  
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 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, 191; 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 + (1 | year2) + (1 | patchid) + (1 | individ)
>    Data: subdata
>      AIC      BIC   logLik deviance df.resid 
>  73.3774  91.2949 -31.6887  63.3774      261 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 83.45090
>  patchid (Intercept)  0.06441
>  year2   (Intercept)  0.15574
> Number of obs: 266, groups:  individ, 191; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept)       obsage  
>      -44.58        15.32  
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 2 lme4 warnings 
> 
> 
> 
> Fecundity model:
> Formula:          
> feca2 ~ obsage + repstatus2 + sizea2 + (1 | year2) + (1 | patchid) +  
>     (1 | individ) + obsage:sizea2
> Zero inflation:         ~obsage + (1 | year2) + (1 | patchid) + (1 | individ)
> Data: subdata
>        AIC        BIC     logLik   df.resid 
> -1337.1622 -1323.9401   682.5811          5 
> Random-effects (co)variances:
> 
> Conditional model:
>  Groups  Name        Std.Dev. 
>  year2   (Intercept) 2.013e-09
>  patchid (Intercept) 4.076e-10
>  individ (Intercept) 4.174e-08
> 
> Zero-inflation model:
>  Groups  Name        Std.Dev. 
>  year2   (Intercept) 2.659e-07
>  patchid (Intercept) 7.703e-07
>  individ (Intercept) 2.159e+01
> 
> 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 (): 1.69e+16 
> 
> Fixed Effects:
> 
> Conditional model:
>   (Intercept)         obsage     repstatus2         sizea2  obsage:sizea2  
>        2.2006       -10.4697         2.2006        -0.4917         1.6159  
> 
> Zero-inflation model:
> (Intercept)       obsage  
>      -25.42        18.00  
> 
> 
> 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 
>  217.1077  229.9675 -104.5539  209.1077       180 
> Random effects:
>  Groups  Name        Std.Dev. 
>  individ (Intercept) 2.663e-07
>  patchid (Intercept) 0.000e+00
>  year2   (Intercept) 0.000e+00
> Number of obs: 184, groups:  individ, 184; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept)  
>        1.07  
> 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.1925  41.8724 -11.0962  22.1925      133 
> Random effects:
>  Groups  Name        Std.Dev. 
>  individ (Intercept) 5.587e+01
>  patchid (Intercept) 3.555e-05
>  year2   (Intercept) 0.000e+00
> Number of obs: 137, groups:  individ, 137; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept)  
>       13.77  
> 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)
>    Data: subdata
> REML criterion at convergence: 145.4426
> Random effects:
>  Groups   Name        Std.Dev.
>  patchid  (Intercept) 0.0000  
>  year2    (Intercept) 0.1430  
>  Residual             0.4139  
> Number of obs: 130, groups:  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: 17
> 
> Number of models in observation table: 16
> 
> 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: 18
> 
> Number of models in fecundity table: 268
> 
> Number of models in juvenile survival table: 1
> 
> Number of models in juvenile observation table: 1
> 
> Number of models in juvenile size table: 1
> 
> Number of models in juvenile secondary size table: 1
> 
> Number of models in juvenile tertiary size table: 1
> 
> Number of models in juvenile reproduction table: 1
> 
> Number of models in juvenile maturity table: 1
> 
> 
> 
> 
> 
> General model parameter names (column 1), and 
> specific names used in these models (column 2): 
>                       parameter_names mainparams
> 1                              time t      year2
> 2                          individual    individ
> 3                               patch      patch
> 4                   alive in time t+1      surv3
> 5                observed in time t+1       obs3
> 6                   sizea in time t+1      size3
> 7                   sizeb in time t+1     sizeb3
> 8                   sizec in time t+1     sizec3
> 9     reproductive status in time t+1     repst3
> 10              fecundity in time t+1       fec3
> 11                fecundity in time t       fec2
> 12                    sizea in time t      size2
> 13                  sizea in time t-1      size1
> 14                    sizeb in time t     sizeb2
> 15                  sizeb in time t-1     sizeb1
> 16                    sizec in time t     sizec2
> 17                  sizec in time t-1     sizec1
> 18      reproductive status in time t     repst2
> 19    reproductive status in time t-1     repst1
> 20        maturity status in time t+1     matst3
> 21          maturity status in time t     matst2
> 22                      age in time t        age
> 23                  density in time t    density
> 24   individual covariate a in time t   indcova2
> 25 individual covariate a in time t-1   indcova1
> 26   individual covariate b in time t   indcovb2
> 27 individual covariate b in time t-1   indcovb1
> 28   individual covariate c in time t   indcovc2
> 29 individual covariate c in time t-1   indcovc1
> 30              stage group in time t     group2
> 31            stage group in time t-1     group1
> 
> 
> 
> 
> 
> Quality control:
> 
> Survival model estimated with 247 individuals and 347 individual transitions.
> Survival model accuracy is 0.847.
> Observation status model estimated with 203 individuals and 294 individual transitions.
> Observation status model accuracy is 1.
> Primary size model estimated with 191 individuals and 266 individual transitions.
> Primary size model R-squared is 0.632.
> Secondary size model not estimated.
> Tertiary size model not estimated.
> Reproductive status model estimated with 191 individuals and 266 individual transitions.
> Reproductive status model accuracy is 1.
> Fecundity model estimated with 16 individuals and 19 individual transitions.
> Fecundity model R-squared is 0.95.
> Juvenile survival model estimated with 184 individuals and 184 individual transitions.
> Juvenile survival model accuracy is 0.745.
> Juvenile observation status model estimated with 137 individuals and 137 individual transitions.
> Juvenile observation status model accuracy is 1.
> Juvenile primary size model estimated with 130 individuals and 130 individual transitions.
> Juvenile primary size model R-squared is 0.058.
> Juvenile secondary size model not estimated.
> Juvenile tertiary size model not estimated.
> Juvenile reproductive status model not estimated.
> Juvenile maturity status model not estimated.

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 716 survival transitions were estimated, with 358 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 247 individuals and 347 individual transitions.
> Observation estimated with 203 individuals and 294 individual transitions.
> Primary size estimated with 191 individuals and 266 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 191 individuals and 266 individual transitions.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival estimated with 184 individuals and 184 individual transitions.
> Juvenile observation estimated with 137 individuals and 137 individual transitions.
> Juvenile primary size estimated with 130 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.442 0.459
> 3rd Qu. 1.000 1.000
> 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 two years (it would have been three 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 4260 survival transitions were estimated, with 355 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 247 individuals and 347 individual transitions.
> Observation estimated with 203 individuals and 294 individual transitions.
> Primary size estimated with 191 individuals and 266 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 191 individuals and 266 individual transitions.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival estimated with 184 individuals and 184 individual transitions.
> Juvenile observation estimated with 137 individuals and 137 individual transitions.
> Juvenile primary size estimated with 130 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]
> Min.    0.000 0.000 0.000 0.000 0.000 0.000 0.0000 0.000 0.000 0.000 0.000
> 1st Qu. 0.000 0.000 0.000 0.000 0.000 0.000 0.0000 0.000 0.000 0.000 0.000
> Median  0.152 0.172 0.114 0.172 0.148 0.172 0.0903 0.149 0.115 0.172 0.147
> Mean    0.423 0.437 0.410 0.422 0.422 0.436 0.3985 0.409 0.410 0.422 0.422
> 3rd Qu. 0.965 0.977 0.945 0.955 0.962 0.975 0.9315 0.937 0.945 0.955 0.961
> Max.    1.000 1.000 1.000 1.000 1.000 1.000 1.0000 1.000 1.000 1.000 1.000
>         [,12]
> Min.    0.000
> 1st Qu. 0.000
> Median  0.172
> Mean    0.436
> 3rd Qu. 0.975
> Max.    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 ((358 + 18) / 1764 = 21.3% 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 716 survival transitions were estimated, with 358 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 247 individuals and 347 individual transitions.
> Observation estimated with 203 individuals and 294 individual transitions.
> Primary size estimated with 191 individuals and 266 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 191 individuals and 266 individual transitions.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival estimated with 184 individuals and 184 individual transitions.
> Juvenile observation estimated with 137 individuals and 137 individual transitions.
> Juvenile primary size estimated with 130 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.813 0.947
> Median  1.000 1.000
> Mean    0.843 0.876
> 3rd Qu. 1.000 1.000
> 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 ((358 + 18) / 484 = 77.7%). 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 360 survival transitions were estimated, with 360 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 247 individuals and 347 individual transitions.
> Observation estimated with 203 individuals and 294 individual transitions.
> Primary size estimated with 191 individuals and 266 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 191 individuals and 266 individual transitions.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival estimated with 184 individuals and 184 individual transitions.
> Juvenile observation estimated with 137 individuals and 137 individual transitions.
> Juvenile primary size estimated with 130 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.450
> 3rd Qu. 1.000
> 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 2499 survival transitions were estimated, with 357 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 247 individuals and 347 individual transitions.
> Observation estimated with 203 individuals and 294 individual transitions.
> Primary size estimated with 191 individuals and 266 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 191 individuals and 266 individual transitions.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival estimated with 184 individuals and 184 individual transitions.
> Juvenile observation estimated with 137 individuals and 137 individual transitions.
> Juvenile primary size estimated with 130 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.152 0.172 0.120 0.153 0.172 0.169
> Mean    0.430 0.416 0.429 0.404 0.416 0.429 0.421
> 3rd Qu. 0.969 0.955 0.968 0.931 0.955 0.968 0.961
> 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 = "indiv_id", blocksize = 9, 
  juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88", 
  fecacol = "Intactseed88", deadacol = "Dead1988", nonobsacol = "Dormant1988", 
  stageassign = lathframe_raw, stagesize = "sizea", censorcol = "Missing1988", 
  censorkeep = NA, censorRepeat = TRUE, censor = TRUE)

lathvert_raw_small <- subset(lathvert_raw, firstseen > 1988)
dim(lathvert_raw_small)
> [1] 531  48
summary_hfv(lathvert_raw_small)
> 
> This hfv dataset contains 531 rows, 48 variables, 1 population, 
> 6 patches, 345 individuals, and 2 time steps.
>      rowid        popid  patchid   individ              year2     
>  Min.   :  30.0   :531   1:130   Length:531         Min.   :1989  
>  1st Qu.: 301.0          2: 69   Class :character   1st Qu.:1989  
>  Median : 582.0          3:101   Mode  :character   Median :1990  
>  Mean   : 585.3          4:127                      Mean   :1990  
>  3rd Qu.: 800.5          5: 69                      3rd Qu.:1990  
>  Max.   :1097.0          6: 35                      Max.   :1990  
>                                                                   
>    firstseen       lastseen        obsage      obslifespan        sizea1       
>  Min.   :1989   Min.   :1989   Min.   :1.00   Min.   :0.000   Min.   :   3.60  
>  1st Qu.:1989   1st Qu.:1990   1st Qu.:1.00   1st Qu.:1.000   1st Qu.:   9.00  
>  Median :1989   Median :1991   Median :1.00   Median :2.000   Median :  12.60  
>  Mean   :1989   Mean   :1991   Mean   :1.35   Mean   :1.367   Mean   : 168.22  
>  3rd Qu.:1989   3rd Qu.:1991   3rd Qu.:2.00   3rd Qu.:2.000   3rd Qu.:  65.22  
>  Max.   :1990   Max.   :1991   Max.   :2.00   Max.   :2.000   Max.   :4394.20  
>                                                               NA's   :345      
>    size1added         repstra1          feca1       fec1added      
>  Min.   :   0.00   Min.   :0.0000   Min.   : 0    Min.   : 0.0000  
>  1st Qu.:   0.00   1st Qu.:0.0000   1st Qu.: 0    1st Qu.: 0.0000  
>  Median :   0.00   Median :0.0000   Median : 4    Median : 0.0000  
>  Mean   :  58.92   Mean   :0.0591   Mean   : 6    Mean   : 0.1243  
>  3rd Qu.:   9.00   3rd Qu.:0.0000   3rd Qu.: 6    3rd Qu.: 0.0000  
>  Max.   :4394.20   Max.   :1.0000   Max.   :34    Max.   :34.0000  
>                    NA's   :345      NA's   :520                    
>    juvgiven1       obsstatus1       repstatus1        fecstatus1     
>  Min.   :0.000   Min.   :0.0000   Min.   :0.00000   Min.   :0.00000  
>  1st Qu.:0.000   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.00000  
>  Median :0.000   Median :0.0000   Median :0.00000   Median :0.00000  
>  Mean   :0.162   Mean   :0.3503   Mean   :0.02072   Mean   :0.01318  
>  3rd Qu.:0.000   3rd Qu.:1.0000   3rd Qu.:0.00000   3rd Qu.:0.00000  
>  Max.   :1.000   Max.   :1.0000   Max.   :1.00000   Max.   :1.00000  
>                                                                      
>    matstatus1         alive1          stage1           stage1index  
>  Min.   :0.0000   Min.   :0.0000   Length:531         Min.   :0.00  
>  1st Qu.:0.0000   1st Qu.:0.0000   Class :character   1st Qu.:0.00  
>  Median :0.0000   Median :0.0000   Mode  :character   Median :0.00  
>  Mean   :0.1883   Mean   :0.3503                      Mean   :1.09  
>  3rd Qu.:0.0000   3rd Qu.:1.0000                      3rd Qu.:2.00  
>  Max.   :1.0000   Max.   :1.0000                      Max.   :6.00  
>                                                                     
>      sizea2          size2added         repstra2           feca2       
>  Min.   :   1.80   Min.   :   0.00   Min.   :0.00000   Min.   : 0.000  
>  1st Qu.:   9.00   1st Qu.:   9.00   1st Qu.:0.00000   1st Qu.: 0.000  
>  Median :  12.60   Median :  12.60   Median :0.00000   Median : 0.000  
>  Mean   :  99.38   Mean   :  97.89   Mean   :0.03633   Mean   : 4.053  
>  3rd Qu.:  40.80   3rd Qu.:  37.40   3rd Qu.:0.00000   3rd Qu.: 6.000  
>  Max.   :4394.20   Max.   :4394.20   Max.   :1.00000   Max.   :34.000  
>  NA's   :8                           NA's   :8         NA's   :512     
>    fec2added        juvgiven2        obsstatus2       repstatus2     
>  Min.   : 0.000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00000  
>  1st Qu.: 0.000   1st Qu.:0.0000   1st Qu.:1.0000   1st Qu.:0.00000  
>  Median : 0.000   Median :0.0000   Median :1.0000   Median :0.00000  
>  Mean   : 0.145   Mean   :0.3465   Mean   :0.9849   Mean   :0.03578  
>  3rd Qu.: 0.000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:0.00000  
>  Max.   :34.000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00000  
>                                                                      
>    fecstatus2        matstatus2         alive2     stage2         
>  Min.   :0.00000   Min.   :0.0000   Min.   :1   Length:531        
>  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.01695   Mean   :0.6535   Mean   :1                     
>  3rd Qu.:0.00000   3rd Qu.:1.0000   3rd Qu.:1                     
>  Max.   :1.00000   Max.   :1.0000   Max.   :1                     
>                                                                   
>   stage2index        sizea3         size3added         repstra3      
>  Min.   :2.000   Min.   :   2.1   Min.   :   0.00   Min.   :0.00000  
>  1st Qu.:2.000   1st Qu.:  10.0   1st Qu.:   0.00   1st Qu.:0.00000  
>  Median :3.000   Median :  12.6   Median :  10.50   Median :0.00000  
>  Mean   :3.149   Mean   : 104.3   Mean   :  77.79   Mean   :0.05556  
>  3rd Qu.:4.000   3rd Qu.:  41.1   3rd Qu.:  24.50   3rd Qu.:0.00000  
>  Max.   :7.000   Max.   :6645.8   Max.   :6645.80   Max.   :1.00000  
>                  NA's   :135                        NA's   :135      
>      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.2505   Mean   :0   Mean   :0.7458  
>  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   :509                                                    
>    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.04143   Mean   :0.01318   Mean   :1    Mean   :0.8117  
>  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:531         Min.   :0.000  
>  Class :character   1st Qu.:3.000  
>  Mode  :character   Median :3.000  
>                     Mean   :3.073  
>                     3rd Qu.:4.000  
>                     Max.   :7.000  
> 

Once again we have a dataset with 531 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 345 unique individuals and 531 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.0714 0.199 0.000 0.199 0.000 0.000 0.000 0.000
> Mean    0.393 0.220 0.362 0.2815 0.383 0.320 0.374 0.303 0.306 0.223 0.323
> 3rd Qu. 0.893 0.399 0.729 0.4748 0.842 0.725 0.781 0.600 0.635 0.399 0.645
> 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 345 unique individuals and 531 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.378 0.361
> 3rd Qu. 0.835 0.766
> 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.70085470 0.78787879 0.25490196 0.1111111    0 0.00000000 0.000
> [11,] 0.000 0.02564103 0.12121212 0.52941176 0.3333333    0 0.36363636 0.000
> [12,] 0.000 0.00000000 0.00000000 0.00000000 0.2222222    0 0.27272727 0.000
> [13,] 0.000 0.02564103 0.03030303 0.05882353 0.1111111    0 0.09090909 0.000
> [14,] 0.000 0.00000000 0.00000000 0.01960784 0.1666667    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} \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] \tag{6.3} \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 = "indiv_id", blocksize = 9,
  sizeacol = "Volume88", repstracol = "FCODE88", fecacol = "Intactseed88",
  deadacol = "Dead1988", censorcol = "Missing1988", censorkeep = NA,
  censorRepeat= TRUE, censor = TRUE, NRasRep = TRUE, NOasObs = TRUE)
lathvert_age <- subset(lathvert_base, firstseen > 1988)

summary_hfv(lathvert_age)
> 
> This hfv dataset contains 531 rows, 48 variables, 1 population, 
> 6 patches, 345 individuals, and 2 time steps.
>      rowid        popid  patchid   individ              year2     
>  Min.   :  30.0   :531   1:130   Length:531         Min.   :1989  
>  1st Qu.: 301.0          2: 69   Class :character   1st Qu.:1989  
>  Median : 582.0          3:101   Mode  :character   Median :1990  
>  Mean   : 585.3          4:127                      Mean   :1990  
>  3rd Qu.: 800.5          5: 69                      3rd Qu.:1990  
>  Max.   :1097.0          6: 35                      Max.   :1990  
>                                                                   
>    firstseen       lastseen        obsage      obslifespan        sizea1       
>  Min.   :1989   Min.   :1989   Min.   :1.00   Min.   :0.000   Min.   :   3.60  
>  1st Qu.:1989   1st Qu.:1990   1st Qu.:1.00   1st Qu.:1.000   1st Qu.:   9.00  
>  Median :1989   Median :1991   Median :1.00   Median :2.000   Median :  12.60  
>  Mean   :1989   Mean   :1991   Mean   :1.35   Mean   :1.367   Mean   : 168.22  
>  3rd Qu.:1989   3rd Qu.:1991   3rd Qu.:2.00   3rd Qu.:2.000   3rd Qu.:  65.22  
>  Max.   :1990   Max.   :1991   Max.   :2.00   Max.   :2.000   Max.   :4394.20  
>                                                               NA's   :345      
>    size1added         repstra1          feca1       fec1added         juvgiven1
>  Min.   :   0.00   Min.   :0.0000   Min.   : 0    Min.   : 0.0000   Min.   :0  
>  1st Qu.:   0.00   1st Qu.:0.0000   1st Qu.: 0    1st Qu.: 0.0000   1st Qu.:0  
>  Median :   0.00   Median :0.0000   Median : 4    Median : 0.0000   Median :0  
>  Mean   :  58.92   Mean   :0.0591   Mean   : 6    Mean   : 0.1243   Mean   :0  
>  3rd Qu.:   9.00   3rd Qu.:0.0000   3rd Qu.: 6    3rd Qu.: 0.0000   3rd Qu.:0  
>  Max.   :4394.20   Max.   :1.0000   Max.   :34    Max.   :34.0000   Max.   :0  
>                    NA's   :345      NA's   :520                                
>    obsstatus1       repstatus1        fecstatus1        matstatus1
>  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 :0.0000   Median :0.00000   Median :0.00000   Median :1   
>  Mean   :0.3503   Mean   :0.02072   Mean   :0.01318   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   
>                                                                   
>      alive1          stage1           stage1index     sizea2       
>  Min.   :0.0000   Length:531         Min.   :0    Min.   :   1.80  
>  1st Qu.:0.0000   Class :character   1st Qu.:0    1st Qu.:   9.00  
>  Median :0.0000   Mode  :character   Median :0    Median :  12.60  
>  Mean   :0.3503                      Mean   :0    Mean   :  99.38  
>  3rd Qu.:1.0000                      3rd Qu.:0    3rd Qu.:  40.80  
>  Max.   :1.0000                      Max.   :0    Max.   :4394.20  
>                                                   NA's   :8        
>    size2added         repstra2           feca2          fec2added     
>  Min.   :   0.00   Min.   :0.00000   Min.   : 0.000   Min.   : 0.000  
>  1st Qu.:   9.00   1st Qu.:0.00000   1st Qu.: 0.000   1st Qu.: 0.000  
>  Median :  12.60   Median :0.00000   Median : 0.000   Median : 0.000  
>  Mean   :  97.89   Mean   :0.03633   Mean   : 4.053   Mean   : 0.145  
>  3rd Qu.:  37.40   3rd Qu.:0.00000   3rd Qu.: 6.000   3rd Qu.: 0.000  
>  Max.   :4394.20   Max.   :1.00000   Max.   :34.000   Max.   :34.000  
>                    NA's   :8         NA's   :512                      
>    juvgiven2   obsstatus2       repstatus2        fecstatus2        matstatus2
>  Min.   :0   Min.   :0.0000   Min.   :0.00000   Min.   :0.00000   Min.   :1   
>  1st Qu.:0   1st Qu.:1.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.9849   Mean   :0.03578   Mean   :0.01695   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   
>                                                                               
>      alive2     stage2           stage2index     sizea3         size3added     
>  Min.   :1   Length:531         Min.   :0    Min.   :   2.1   Min.   :   0.00  
>  1st Qu.:1   Class :character   1st Qu.:0    1st Qu.:  10.0   1st Qu.:   0.00  
>  Median :1   Mode  :character   Median :0    Median :  12.6   Median :  10.50  
>  Mean   :1                      Mean   :0    Mean   : 104.3   Mean   :  77.79  
>  3rd Qu.:1                      3rd Qu.:0    3rd Qu.:  41.1   3rd Qu.:  24.50  
>  Max.   :1                      Max.   :0    Max.   :6645.8   Max.   :6645.80  
>                                              NA's   :135                       
>     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.2505   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   :135       NA's   :509                                   
>    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.7458   Mean   :0.04143   Mean   :0.01318   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:531         Min.   :0   
>  1st Qu.:1.0000   Class :character   1st Qu.:0   
>  Median :1.0000   Mode  :character   Median :0   
>  Mean   :0.7608                      Mean   :0   
>  3rd Qu.:1.0000                      3rd Qu.:0   
>  Max.   :1.0000                      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 
>  588.5991  605.6982 -290.2996  580.5991       527 
> Random effects:
>  Groups  Name        Std.Dev.
>  individ (Intercept) 0.0000  
>  year2   (Intercept) 0.2988  
> Number of obs: 531, groups:  individ, 345; year2, 2
> Fixed Effects:
> (Intercept)       obsage  
>     1.08983      0.09083  
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings 
> 
> 
> 
> 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 model estimated with 345 individuals and 531 individual transitions.
> Survival model accuracy is 0.761.
> Observation status model not estimated.
> Primary size model not estimated.
> Secondary size model not estimated.
> Tertiary size model not estimated.
> Reproductive status model not estimated.
> Fecundity model estimated with 16 individuals and 19 individual transitions.
> Fecundity model R-squared is 0.834.
> Juvenile survival model not estimated.
> Juvenile observation status model not estimated.
> Juvenile primary size model not estimated.
> Juvenile secondary size model not estimated.
> Juvenile tertiary size model not estimated.
> Juvenile reproductive status model not estimated.
> Juvenile maturity status model not estimated.

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 345 unique individuals and 531 unique transitions.
> 
> Survival probability sum check (each matrix represented by column in order):
>          [,1]
> Min.    0.000
> 1st Qu. 0.343
> Median  0.687
> Mean    0.473
> 3rd Qu. 0.709
> Max.    0.731

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.02608696 0.04301075    0
> [2,] 0.68695652 0.00000000    0
> [3,] 0.00000000 0.73118280    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 345 individuals and 531 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.809 0.713
> 1st Qu. 0.816 0.722
> Median  0.823 0.731
> Mean    0.822 0.731
> 3rd Qu. 0.829 0.740
> Max.    0.835 0.749
lathmat2ageonly_func$A[[1]]
>            [,1]      [,2]      [,3]
> [1,] 0.02118613 0.2318637 0.8080455
> [2,] 0.80895594 0.0000000 0.0000000
> [3,] 0.00000000 0.8226009 0.8354695

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.