# Chapter 3 Cypripedium candidum

#### Richard P. Shefferson

This document was built in Markdown in R 4.0.2.

## CASE STUDIES OF AMERICAN Cypripedium candidum POPULATION

### ORGANISM AND POPULATION

In this tutorial, we will focus on dataset `cypdata`, which is a demographic dataset for an American population of the white lady’s slipper, Cypripedium candidum. This population was monitored annually from 2004 to 2009. This species is of conservation concern, and the population is located within a state nature preserve located in northwestern Illinois, USA. More about this population and its characteristics is given in Shefferson et al. (2001) and Shefferson et al. (2017).

Population matrix projection modeling requires an appropriate life history model showing how stages are related. The figure below shows a very general life history model detailing how stages are related to one another in Cypripedium candidum. The first stage of life is a dormant seed stage, although an individual may germinate in the year following seed production. The first germinated stage is a protocorm, which is an underground, mycoheterotrophic stage unique to the families Orchidaceae and Pyrolaceae. There are three years of protocorm stages, followed by a seedling stage, and finally a set of stages that comprise the size-classified adult portion of life. The figure shows 49 such stages, each for a different number of stems and one of two reproductive statuses. These stages may be compressed for different circumstances (more on this later).

We can see a variety of transitions within this figure. The juvenile stages have fairly simple transitions. New recruits may enter the population directly from germination of a seed produced the previous year, in which case they start in the protocorm 1 stage, or they may begin as dormant seed. Dormant seed may remain dormant, die, or germinate into the protocorm 1 stage. Protocorms exist for up to 3 years, yielding the protocorm 1, 2, and 3 stages, without any possibility of staying within each of these stages for more than a single year. This leads to a seedling stage, in which the plant may persist for many years, before becoming mature. Here, maturity does not really refer to reproduction per se, but rather to a morphology indistinguishable from a reproductive plant except for the lack of a flower. The first mature stage is usually either vegetative dormancy (dorm), during which time the plant does not sprout, or a small, non-flowering adult (1V). Once in this portion of the life history, the plant may move among 49 mature stages, including vegetative dormancy, 1-24 shoots without flowers, or 1-24 shoots with at least one flower.

The dataset `cypdata`, and the equivalent dataset `cypvert`, which is simply structured differently, both include only data for the adult stages, and so later we will need to set juvenile transitions to constants.

### ANALYSES WITH CYPRIPEDIUM DATA

We will analyze these data in two different ways to illustrate the power and flexibility of package `lefko3`:

1. through the estimation of raw matrices using a simplified life history; and

2. through the estimation of function-derived matrices using a count-based size metric and the general life history model shown above.

### Analysis 1. Raw matrix estimation

In this portion of the tutorial, we will create raw matrices with these data. Here, we use the term ‘raw’ to refer to the fact that we will estimate matrix elements as exact proportions of individuals surviving and transitioning to different stages. This requires us to develop a life history model that is both biologically realistic, and numerically parsimonious. The former requirement means that stages need to be defined in biologically meaningful ways, while the latter requirement means that we need to design our life stages in such a way that most years include some individuals in each stage. We also need to consider the fact that very low numbers of stages result in biased matrix analyses, so we want to make sure that we have at least 7 stages in the final model (Salguero-Gómez and Plotkin 2010).

#### Step 1a. Data characterization and reorganization (horizontal file)

First let’s wipe the memory, load `lefko3`, and then load the data.

    rm(list=ls(all=TRUE))

    library(lefko3)
    data(cypdata)

    summary(cypdata)
    #>     plantid       patch      censor     Inf2.04       Inf.04           Veg.04           Pod.04       Inf2.05            Inf.05
    #>  Min.   : 164.0   A:23   Min.   :1   Min.   :0    Min.   :0.0000   Min.   : 0.000   Min.   :0.0   Min.   :0.00000   Min.   : 0.000
    #>  1st Qu.: 265.0   B:35   1st Qu.:1   1st Qu.:0    1st Qu.:0.0000   1st Qu.: 1.000   1st Qu.:0.0   1st Qu.:0.00000   1st Qu.: 0.000
    #>  Median : 455.0   C:19   Median :1   Median :0    Median :0.0000   Median : 2.000   Median :0.0   Median :0.00000   Median : 0.000
    #>  Mean   : 669.1          Mean   :1   Mean   :0    Mean   :0.6923   Mean   : 2.923   Mean   :0.2   Mean   :0.04478   Mean   : 1.537
    #>  3rd Qu.: 829.0          3rd Qu.:1   3rd Qu.:0    3rd Qu.:1.0000   3rd Qu.: 4.000   3rd Qu.:0.0   3rd Qu.:0.00000   3rd Qu.: 2.000
    #>  Max.   :1560.0          Max.   :1   Max.   :0    Max.   :8.0000   Max.   :12.000   Max.   :3.0   Max.   :1.00000   Max.   :18.000
    #>                                      NA's   :12   NA's   :12       NA's   :12       NA's   :12    NA's   :10        NA's   :10
    #>      Veg.05          Pod.05          Inf2.06       Inf.06            Veg.06           Pod.06          Inf2.07      Inf.07
    #>  Min.   :0.000   Min.   :0.0000   Min.   :0    Min.   : 0.0000   Min.   : 0.000   Min.   :0.0000   Min.   :0   Min.   :0.0000
    #>  1st Qu.:1.000   1st Qu.:0.0000   1st Qu.:0    1st Qu.: 0.0000   1st Qu.: 1.000   1st Qu.:0.0000   1st Qu.:0   1st Qu.:0.0000
    #>  Median :1.000   Median :0.0000   Median :0    Median : 0.0000   Median : 2.000   Median :0.0000   Median :0   Median :0.0000
    #>  Mean   :2.134   Mean   :0.6567   Mean   :0    Mean   : 0.9016   Mean   : 2.213   Mean   :0.3934   Mean   :0   Mean   :0.6271
    #>  3rd Qu.:3.000   3rd Qu.:1.0000   3rd Qu.:0    3rd Qu.: 1.0000   3rd Qu.: 3.000   3rd Qu.:0.0000   3rd Qu.:0   3rd Qu.:1.0000
    #>  Max.   :9.000   Max.   :7.0000   Max.   :0    Max.   :18.0000   Max.   :13.000   Max.   :4.0000   Max.   :0   Max.   :7.0000
    #>  NA's   :10      NA's   :10       NA's   :16   NA's   :16        NA's   :16       NA's   :16                   NA's   :18
    #>      Veg.07           Pod.07          Inf2.08       Inf.08            Veg.08          Pod.08          Inf2.09       Inf.09
    #>  Min.   : 0.000   Min.   :0.0000   Min.   :0    Min.   : 0.0000   Min.   : 0.00   Min.   :0.0000   Min.   :0    Min.   : 0.000
    #>  1st Qu.: 1.000   1st Qu.:0.0000   1st Qu.:0    1st Qu.: 0.0000   1st Qu.: 1.00   1st Qu.:0.0000   1st Qu.:0    1st Qu.: 0.000
    #>  Median : 2.000   Median :0.0000   Median :0    Median : 0.0000   Median : 2.00   Median :0.0000   Median :0    Median : 1.000
    #>  Mean   : 2.627   Mean   :0.0678   Mean   :0    Mean   : 0.8868   Mean   : 2.83   Mean   :0.1509   Mean   :0    Mean   : 1.833
    #>  3rd Qu.: 4.000   3rd Qu.:0.0000   3rd Qu.:0    3rd Qu.: 1.0000   3rd Qu.: 4.00   3rd Qu.:0.0000   3rd Qu.:0    3rd Qu.: 2.000
    #>  Max.   :13.000   Max.   :1.0000   Max.   :0    Max.   :11.0000   Max.   :13.00   Max.   :2.0000   Max.   :0    Max.   :11.000
    #>  NA's   :18       NA's   :18       NA's   :24   NA's   :24        NA's   :24      NA's   :24       NA's   :17   NA's   :17
    #>      Veg.09           Pod.09
    #>  Min.   : 0.000   Min.   :0.000
    #>  1st Qu.: 1.000   1st Qu.:0.000
    #>  Median : 1.000   Median :1.000
    #>  Mean   : 2.233   Mean   :1.133
    #>  3rd Qu.: 3.000   3rd Qu.:1.000
    #>  Max.   :10.000   Max.   :8.000
    #>  NA's   :17       NA's   :17

The dataset that we have provided is organized in horizontal format, meaning that rows correspond to unique individuals and columns correspond to state in particular years. Looking at the original Excel spreadsheet (below), you will note a repeating pattern in the names of the columns. Package `lefko3` includes functions to handle data in horizontal format, as well as functions to handle vertically formatted data (i.e. data for individuals is broken up across rows, where each row is a unique combination of individual and year in time t).

In this dataset, there are 77 individuals, so there are 77 rows with data (not counting the header). There are 27 columns. Note the the first 3 columns are variables giving specific information about each individual, with each individual’s data entirely restricted to one row. This is followed by a number of sets of 4 columns, each named `Inf2.XX`, `Inf.XX`, `Veg.XX`, and `Pod.XX`. The XX in each case corresponds to a specific year, which are organized consecutively. Thus, columns 4-7 refer to year 04 (short for 2004), columns 8-11 refer to year 05, columns 12-15 refer to year 06, etc. To properly conduct this exercise, we need to know the exact number of years used, which is six years here (includes all years from 2004 to 2009). Note that each year MUST utilize exactly the same pattern of columns. Also note that this package also includes dataset `cypvert`, which is the same dataset but set in ahistorical vertical format.

Now on to the assessment of size. The full sizes are actually the sums of columns (representing sprouts) governing different years. We will take these sums, and then assess the distribution of sizes across years. We will look at all years and look for general patterns and abnormalities.

    size.04 <- cypdata$Inf2.04 + cypdata$Inf.04 + cypdata$Veg.04
    size.05 <- cypdata$Inf2.05 + cypdata$Inf.05 + cypdata$Veg.05
    size.06 <- cypdata$Inf2.06 + cypdata$Inf.06 + cypdata$Veg.06
    size.07 <- cypdata$Inf2.07 + cypdata$Inf.07 + cypdata$Veg.07
    size.08 <- cypdata$Inf2.08 + cypdata$Inf.08 + cypdata$Veg.08
    size.09 <- cypdata$Inf2.09 + cypdata$Inf.09 + cypdata$Veg.09

    summary(c(size.04, size.05, size.06, size.07, size.08, size.09))
    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
    #>   1.000   1.000   2.000   3.581   5.000  24.000      97

Now let’s quickly plot this.

    plot(density(c(size.04, size.05, size.06, size.07, size.08, size.09),
                 na.rm = TRUE))

plot of chunk Ch3an1st1aln79pl

This exercise gives us a reasonable idea of size classes to use for adults. We will have a dormant class (size = 0 shoots), extra small class (1 shoot), small class (2-3 shoots), medium class (4-5 shoots), large class (6-10 shoots), and extra large class (>10 shoots). Let’s define a `stageframe` that shows this.

    sizevector <- c(0, 0, 0, 0, 0, 0, 1, 2.5, 4.5, 8, 17.5)
    stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "XSm", "Sm", "Md", "Lg", "XLg")
    repvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
    obsvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
    matvector <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)
    immvector <- c(0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0)
    propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
    indataset <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)
    binvec <- c(0, 0, 0, 0, 0, 0.5, 0.5, 1, 1, 2.5, 7)

    cypframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector,
                          repstatus = repvector, obsstatus = obsvector,
                          matstatus = matvector, propstatus = propvector,
                          immstatus = immvector, indataset = indataset,
                          binhalfwidth = binvec)

Now we will add some comments to the stageframe for our later use in interpretation.

    cypframe_raw$comments[(cypframe_raw$stagenames == "SD")] <- "Dormant seed"
    cypframe_raw$comments[(cypframe_raw$stagenames == "P1")] <- "1st yr protocorm"
    cypframe_raw$comments[(cypframe_raw$stagenames == "P2")] <- "2nd yr protocorm"
    cypframe_raw$comments[(cypframe_raw$stagenames == "P3")] <- "3rd yr protocorm"
    cypframe_raw$comments[(cypframe_raw$stagenames == "SL")] <- "Seedling"
    cypframe_raw$comments[(cypframe_raw$stagenames == "D")] <- "Dormant adult"
    cypframe_raw$comments[(cypframe_raw$stagenames == "XSm")] <- "Extra small adult (1 shoot)"
    cypframe_raw$comments[(cypframe_raw$stagenames == "Sm")] <- "Small adult (2-3 shoots)"
    cypframe_raw$comments[(cypframe_raw$stagenames == "Md")] <- "Medium adult (4-5 shoots)"
    cypframe_raw$comments[(cypframe_raw$stagenames == "Lg")] <- "Large adult (6-10 shoots)"
    cypframe_raw$comments[(cypframe_raw$stagenames == "XLg")] <- "Extra large adult (>10 shoots)"

    cypframe_raw
    #>    stagenames size repstatus obsstatus propstatus immstatus matstatus indataset binhalfwidth_raw sizebin_min sizebin_max sizebin_center
    #> 1          SD  0.0         0         0          1         0         0         0              0.0         0.0         0.0            0.0
    #> 2          P1  0.0         0         0          0         1         0         0              0.0         0.0         0.0            0.0
    #> 3          P2  0.0         0         0          0         1         0         0              0.0         0.0         0.0            0.0
    #> 4          P3  0.0         0         0          0         1         0         0              0.0         0.0         0.0            0.0
    #> 5          SL  0.0         0         0          0         1         0         0              0.0         0.0         0.0            0.0
    #> 6           D  0.0         0         0          0         0         1         1              0.5        -0.5         0.5            0.0
    #> 7         XSm  1.0         1         1          0         0         1         1              0.5         0.5         1.5            1.0
    #> 8          Sm  2.5         1         1          0         0         1         1              1.0         1.5         3.5            2.5
    #> 9          Md  4.5         1         1          0         0         1         1              1.0         3.5         5.5            4.5
    #> 10         Lg  8.0         1         1          0         0         1         1              2.5         5.5        10.5            8.0
    #> 11        XLg 17.5         1         1          0         0         1         1              7.0        10.5        24.5           17.5
    #>    sizebin_width                       comments
    #> 1              0                   Dormant seed
    #> 2              0               1st yr protocorm
    #> 3              0               2nd yr protocorm
    #> 4              0               3rd yr protocorm
    #> 5              0                       Seedling
    #> 6              1                  Dormant adult
    #> 7              1    Extra small adult (1 shoot)
    #> 8              2       Small adult (2-3 shoots)
    #> 9              2      Medium adult (4-5 shoots)
    #> 10             5      Large adult (6-10 shoots)
    #> 11            14 Extra large adult (>10 shoots)

Next we will create the vertical dataset. Note that because we are lumping reproductive and non-reproductive individuals into the non-dormant adult classes, we need to set `NRasRP = TRUE`. Otherwise, `verticalize3()` will attempt to use the reproductive status of individuals in classification, and will fail due to the presence of non-reproductive adults.

    cypraw_v1 <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004,
                             patchidcol = "patch", individcol = "plantid",
                             blocksize = 4, size1col = "Inf2.04", size2col = "Inf.04",
                             size3col = "Veg.04", repstr1col = "Inf.04",
                             repstr2col = "Inf2.04", fec1col = "Pod.04",
                             stageassign = cypframe_raw, stagesize = "sizeadded",
                             NAas0 = TRUE, NRasRep = TRUE)

    summary(cypraw_v1)
    #>      rowid        popid         patchid    individ           year2          sizea1             sizea2             sizea3
    #>  Min.   : 1.00   Mode:logical   A: 93   Min.   : 164.0   Min.   :2004   Min.   :0.000000   Min.   :0.000000   Min.   :0.000000
    #>  1st Qu.:21.00   NA's:320       B:154   1st Qu.: 391.0   1st Qu.:2005   1st Qu.:0.000000   1st Qu.:0.000000   1st Qu.:0.000000
    #>  Median :37.50                  C: 73   Median : 453.0   Median :2006   Median :0.000000   Median :0.000000   Median :0.000000
    #>  Mean   :38.45                          Mean   : 651.5   Mean   :2006   Mean   :0.009375   Mean   :0.009375   Mean   :0.009375
    #>  3rd Qu.:56.00                          3rd Qu.: 476.0   3rd Qu.:2007   3rd Qu.:0.000000   3rd Qu.:0.000000   3rd Qu.:0.000000
    #>  Max.   :77.00                          Max.   :1560.0   Max.   :2008   Max.   :1.000000   Max.   :1.000000   Max.   :1.000000
    #>     repstra1          repstra2          repstra3          feca1            feca2            feca3            sizeb1
    #>  Min.   : 0.0000   Min.   : 0.0000   Min.   : 0.000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   : 0.0000
    #>  1st Qu.: 0.0000   1st Qu.: 0.0000   1st Qu.: 0.000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.: 0.0000
    #>  Median : 0.0000   Median : 0.0000   Median : 0.000   Median :0.0000   Median :0.0000   Median :0.0000   Median : 0.0000
    #>  Mean   : 0.7469   Mean   : 0.8969   Mean   : 1.069   Mean   :0.2656   Mean   :0.2906   Mean   :0.4562   Mean   : 0.7469
    #>  3rd Qu.: 1.0000   3rd Qu.: 1.0000   3rd Qu.: 1.000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.: 1.0000
    #>  Max.   :18.0000   Max.   :18.0000   Max.   :18.000   Max.   :7.0000   Max.   :7.0000   Max.   :8.0000   Max.   :18.0000
    #>      sizeb2            sizeb3          repstrb1           repstrb2           repstrb3            sizec1         sizec2
    #>  Min.   : 0.0000   Min.   : 0.000   Min.   :0.000000   Min.   :0.000000   Min.   :0.000000   Min.   : 0.0   Min.   : 0.000
    #>  1st Qu.: 0.0000   1st Qu.: 0.000   1st Qu.:0.000000   1st Qu.:0.000000   1st Qu.:0.000000   1st Qu.: 0.0   1st Qu.: 1.000
    #>  Median : 0.0000   Median : 0.000   Median :0.000000   Median :0.000000   Median :0.000000   Median : 1.0   Median : 2.000
    #>  Mean   : 0.8969   Mean   : 1.069   Mean   :0.009375   Mean   :0.009375   Mean   :0.009375   Mean   : 1.9   Mean   : 2.416
    #>  3rd Qu.: 1.0000   3rd Qu.: 1.000   3rd Qu.:0.000000   3rd Qu.:0.000000   3rd Qu.:0.000000   3rd Qu.: 3.0   3rd Qu.: 3.000
    #>  Max.   :18.0000   Max.   :18.000   Max.   :1.000000   Max.   :1.000000   Max.   :1.000000   Max.   :13.0   Max.   :13.000
    #>      sizec3         size1added       size2added       size3added       obsstatus1       obsstatus2       obsstatus3    repstatus1
    #>  Min.   : 0.000   Min.   : 0.000   Min.   : 0.000   Min.   : 0.000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0   Min.   :0.0000
    #>  1st Qu.: 1.000   1st Qu.: 0.000   1st Qu.: 1.000   1st Qu.: 1.000   1st Qu.:0.0000   1st Qu.:1.0000   1st Qu.:1.0   1st Qu.:0.0000
    #>  Median : 1.000   Median : 2.000   Median : 2.000   Median : 2.000   Median :1.0000   Median :1.0000   Median :1.0   Median :0.0000
    #>  Mean   : 2.209   Mean   : 2.656   Mean   : 3.322   Mean   : 3.288   Mean   :0.7469   Mean   :0.9531   Mean   :0.9   Mean   :0.2875
    #>  3rd Qu.: 3.000   3rd Qu.: 4.000   3rd Qu.: 4.000   3rd Qu.: 4.000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0   3rd Qu.:1.0000
    #>  Max.   :13.000   Max.   :21.000   Max.   :24.000   Max.   :24.000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0   Max.   :1.0000
    #>    repstatus2       repstatus3    fecstatus1       fecstatus2       fecstatus3       firstseen       lastseen        alive1
    #>  Min.   :0.0000   Min.   :0.0   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :2004   Min.   :2004   Min.   :0.0000
    #>  1st Qu.:0.0000   1st Qu.:0.0   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:2004   1st Qu.:2009   1st Qu.:1.0000
    #>  Median :0.0000   Median :0.0   Median :0.0000   Median :0.0000   Median :0.0000   Median :2004   Median :2009   Median :1.0000
    #>  Mean   :0.3688   Mean   :0.4   Mean   :0.1344   Mean   :0.1562   Mean   :0.2219   Mean   :2004   Mean   :2009   Mean   :0.7688
    #>  3rd Qu.:1.0000   3rd Qu.:1.0   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:2004   3rd Qu.:2009   3rd Qu.:1.0000
    #>  Max.   :1.0000   Max.   :1.0   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :2008   Max.   :2009   Max.   :1.0000
    #>      alive2      alive3           obsage       obslifespan       stage1             stage2             stage3            matstatus1
    #>  Min.   :1   Min.   :0.0000   Min.   :0.000   Min.   :0.000   Length:320         Length:320         Length:320         Min.   :1
    #>  1st Qu.:1   1st Qu.:1.0000   1st Qu.:1.000   1st Qu.:5.000   Class :character   Class :character   Class :character   1st Qu.:1
    #>  Median :1   Median :1.0000   Median :2.000   Median :5.000   Mode  :character   Mode  :character   Mode  :character   Median :1
    #>  Mean   :1   Mean   :0.9469   Mean   :1.853   Mean   :4.556                                                            Mean   :1
    #>  3rd Qu.:1   3rd Qu.:1.0000   3rd Qu.:3.000   3rd Qu.:5.000                                                            3rd Qu.:1
    #>  Max.   :1   Max.   :1.0000   Max.   :4.000   Max.   :5.000                                                            Max.   :1
    #>    matstatus2   matstatus3
    #>  Min.   :1    Min.   :1
    #>  1st Qu.:1    1st Qu.:1
    #>  Median :1    Median :1
    #>  Mean   :1    Mean   :1
    #>  3rd Qu.:1    3rd Qu.:1
    #>  Max.   :1    Max.   :1
    dim(cypraw_v1)
    #> [1] 320  48

#### Step 1b. Data characterization and reorganization (vertical file)

We may also wish to see how to proceed if our original dataset is already in vertical, but ahistorical, format. Here, we see such a situation and use the `historicalize3()` function. First, let’s load the ahistorical vertical raw data file.

    data(cypvert)

    summary(cypvert)
    #>     plantid          patch               censor      year2          Inf2.2             Inf.2            Veg.2            Pod.2
    #>  Min.   : 164.0   Length:331         Min.   :1   Min.   :2004   Min.   :0.000000   Min.   : 0.000   Min.   : 0.000   Min.   :0.0000
    #>  1st Qu.: 391.0   Class :character   1st Qu.:1   1st Qu.:2005   1st Qu.:0.000000   1st Qu.: 0.000   1st Qu.: 1.000   1st Qu.:0.0000
    #>  Median : 454.0   Mode  :character   Median :1   Median :2006   Median :0.000000   Median : 0.000   Median : 2.000   Median :0.0000
    #>  Mean   : 662.8                      Mean   :1   Mean   :2006   Mean   :0.009646   Mean   : 0.941   Mean   : 2.534   Mean   :0.3049
    #>  3rd Qu.: 664.0                      3rd Qu.:1   3rd Qu.:2007   3rd Qu.:0.000000   3rd Qu.: 1.000   3rd Qu.: 3.000   3rd Qu.:0.0000
    #>  Max.   :1560.0                      Max.   :1   Max.   :2008   Max.   :1.000000   Max.   :18.000   Max.   :13.000   Max.   :7.0000
    #>                                                                 NA's   :20         NA's   :26       NA's   :26       NA's   :26
    #>      Inf2.3             Inf.3            Veg.3           Pod.3
    #>  Min.   :0.000000   Min.   : 0.000   Min.   : 0.00   Min.   :0.0000
    #>  1st Qu.:0.000000   1st Qu.: 0.000   1st Qu.: 1.00   1st Qu.:0.0000
    #>  Median :0.000000   Median : 0.000   Median : 1.00   Median :0.0000
    #>  Mean   :0.009804   Mean   : 1.173   Mean   : 2.39   Mean   :0.4933
    #>  3rd Qu.:0.000000   3rd Qu.: 1.000   3rd Qu.: 3.00   3rd Qu.:0.0000
    #>  Max.   :1.000000   Max.   :18.000   Max.   :13.00   Max.   :8.0000
    #>  NA's   :25         NA's   :31       NA's   :31      NA's   :31

And let’s also look at its dimensions.

    dim(cypvert)
    #> [1] 331  12

This dataset appears to be longer and narrower, with more rows and fewer columns. This is because we now split data for each individual across multiple columns. A single column designates time in year t, given as `year2`. This dataset then includes columns for only pairs of consecutive years corresponding to times t and t+1. State in time t-1 is not presented because this is an ahistorical dataset. Fortunately, this dataset includes the `plantid` variable, which is an individual identity term and must be supplied for conversion. The `historicalize3()` function can use individual identity variables to reorganize datasets into historical vertical format.

    cypraw_v2 <- historicalize3(data = cypvert, patchidcol = "patch", individcol = "plantid",
                               year2col = "year2", sizea2col = "Inf2.2", sizea3col = "Inf2.3",
                               sizeb2col = "Inf.2", sizeb3col = "Inf.3", sizec2col = "Veg.2",
                               sizec3col = "Veg.3", repstra2col = "Inf2.2", repstra3col = "Inf2.3",
                               repstrb2col = "Inf.2", repstrb3col = "Inf.3", feca2col = "Pod.2",
                               feca3col = "Pod.3", repstrrel = 2, stageassign = cypframe_raw,
                               stagesize = "sizeadded", censorcol = "censor", censor = FALSE,
                               NAas0 = TRUE, NRasRep = TRUE, reduce = TRUE)
    summary(cypraw_v2)
    #>   popid           patchid              rowid           sizea1             sizeb1            sizec1          repstra1
    #>  Mode:logical   Length:320         Min.   :  1.0   Min.   :0.000000   Min.   : 0.0000   Min.   : 0.000   Min.   :0.000000
    #>  NA's:320       Class :character   1st Qu.: 41.0   1st Qu.:0.000000   1st Qu.: 0.0000   1st Qu.: 0.000   1st Qu.:0.000000
    #>                 Mode  :character   Median : 95.5   Median :0.000000   Median : 0.0000   Median : 1.000   Median :0.000000
    #>                                    Mean   :111.3   Mean   :0.009375   Mean   : 0.7375   Mean   : 1.881   Mean   :0.009375
    #>                                    3rd Qu.:178.2   3rd Qu.:0.000000   3rd Qu.: 1.0000   3rd Qu.: 3.000   3rd Qu.:0.000000
    #>                                    Max.   :311.0   Max.   :1.000000   Max.   :18.0000   Max.   :13.000   Max.   :1.000000
    #>     repstrb1           feca1           stage1              year2          sizea2             sizeb2            sizec2
    #>  Min.   : 0.0000   Min.   :0.0000   Length:320         Min.   :2004   Min.   :0.000000   Min.   : 0.0000   Min.   : 0.000
    #>  1st Qu.: 0.0000   1st Qu.:0.0000   Class :character   1st Qu.:2005   1st Qu.:0.000000   1st Qu.: 0.0000   1st Qu.: 1.000
    #>  Median : 0.0000   Median :0.0000   Mode  :character   Median :2006   Median :0.000000   Median : 0.0000   Median : 1.500
    #>  Mean   : 0.7375   Mean   :0.2625                      Mean   :2006   Mean   :0.009375   Mean   : 0.8812   Mean   : 2.378
    #>  3rd Qu.: 1.0000   3rd Qu.:0.0000                      3rd Qu.:2007   3rd Qu.:0.000000   3rd Qu.: 1.0000   3rd Qu.: 3.000
    #>  Max.   :18.0000   Max.   :7.0000                      Max.   :2008   Max.   :1.000000   Max.   :18.0000   Max.   :13.000
    #>     repstra2           repstrb2           feca2           stage2              sizea3             sizeb3           sizec3
    #>  Min.   :0.000000   Min.   : 0.0000   Min.   :0.0000   Length:320         Min.   :0.000000   Min.   : 0.000   Min.   : 0.000
    #>  1st Qu.:0.000000   1st Qu.: 0.0000   1st Qu.:0.0000   Class :character   1st Qu.:0.000000   1st Qu.: 0.000   1st Qu.: 1.000
    #>  Median :0.000000   Median : 0.0000   Median :0.0000   Mode  :character   Median :0.000000   Median : 0.000   Median : 1.000
    #>  Mean   :0.009375   Mean   : 0.8812   Mean   :0.2875                      Mean   :0.009375   Mean   : 1.059   Mean   : 2.209
    #>  3rd Qu.:0.000000   3rd Qu.: 1.0000   3rd Qu.:0.0000                      3rd Qu.:0.000000   3rd Qu.: 1.000   3rd Qu.: 3.000
    #>  Max.   :1.000000   Max.   :18.0000   Max.   :7.0000                      Max.   :1.000000   Max.   :18.000   Max.   :13.000
    #>     repstra3           repstrb3          feca3           stage3            firstseen       lastseen        alive1           alive2
    #>  Min.   :0.000000   Min.   : 0.000   Min.   :0.0000   Length:320         Min.   :2004   Min.   :2005   Min.   :0.0000   Min.   :1
    #>  1st Qu.:0.000000   1st Qu.: 0.000   1st Qu.:0.0000   Class :character   1st Qu.:2004   1st Qu.:2009   1st Qu.:1.0000   1st Qu.:1
    #>  Median :0.000000   Median : 0.000   Median :0.0000   Mode  :character   Median :2004   Median :2009   Median :1.0000   Median :1
    #>  Mean   :0.009375   Mean   : 1.059   Mean   :0.4437                      Mean   :2004   Mean   :2009   Mean   :0.7688   Mean   :1
    #>  3rd Qu.:0.000000   3rd Qu.: 1.000   3rd Qu.:0.0000                      3rd Qu.:2004   3rd Qu.:2009   3rd Qu.:1.0000   3rd Qu.:1
    #>  Max.   :1.000000   Max.   :18.000   Max.   :8.0000                      Max.   :2008   Max.   :2009   Max.   :1.0000   Max.   :1
    #>      alive3    obsstatus1       obsstatus2       obsstatus3       matstatus1   matstatus2   matstatus3   repstatus1       repstatus2
    #>  Min.   :1   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :1    Min.   :1    Min.   :1    Min.   :0.0000   Min.   :0.0000
    #>  1st Qu.:1   1st Qu.:0.0000   1st Qu.:1.0000   1st Qu.:1.0000   1st Qu.:1    1st Qu.:1    1st Qu.:1    1st Qu.:0.0000   1st Qu.:0.0000
    #>  Median :1   Median :1.0000   Median :1.0000   Median :1.0000   Median :1    Median :1    Median :1    Median :0.0000   Median :0.0000
    #>  Mean   :1   Mean   :0.7344   Mean   :0.9281   Mean   :0.9125   Mean   :1    Mean   :1    Mean   :1    Mean   :0.2812   Mean   :0.3563
    #>  3rd Qu.:1   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1    3rd Qu.:1    3rd Qu.:1    3rd Qu.:1.0000   3rd Qu.:1.0000
    #>  Max.   :1   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1    Max.   :1    Max.   :1    Max.   :1.0000   Max.   :1.0000
    #>    repstatus3       fecstatus1       fecstatus2       fecstatus3        individ           obsage       obslifespan      size1added
    #>  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   : 164.0   Min.   :0.000   Min.   :1.000   Min.   : 0.000
    #>  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.: 391.0   1st Qu.:1.000   1st Qu.:5.000   1st Qu.: 0.000
    #>  Median :0.0000   Median :0.0000   Median :0.0000   Median :0.0000   Median : 453.0   Median :2.000   Median :5.000   Median : 1.000
    #>  Mean   :0.4031   Mean   :0.1313   Mean   :0.1531   Mean   :0.2125   Mean   : 651.5   Mean   :1.853   Mean   :4.706   Mean   : 2.628
    #>  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.: 476.0   3rd Qu.:3.000   3rd Qu.:5.000   3rd Qu.: 4.000
    #>  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1560.0   Max.   :4.000   Max.   :5.000   Max.   :21.000
    #>    size2added       size3added
    #>  Min.   : 0.000   Min.   : 0.000
    #>  1st Qu.: 1.000   1st Qu.: 1.000
    #>  Median : 2.000   Median : 2.000
    #>  Mean   : 3.269   Mean   : 3.278
    #>  3rd Qu.: 4.000   3rd Qu.: 4.000
    #>  Max.   :24.000   Max.   :24.000

We can compare the dimensions of these datasets.

    dim(cypraw_v1)
    #> [1] 320  48
    dim(cypraw_v2)
    #> [1] 320  48

The lengths of the datasets are the same in terms of rows and columns, and the variables and data are the same although the order of the rows might not exactly match.

#### Step 2. Provide supplemental information for matrix estimation

For our next step, we need to create a reproductive matrix, which essentially tells R not only which stages are reproductive, but which stages they lead to the reproduction of, and at what level. This matrix is mostly composed of 0s, but fecundity is noted as non-zero entries equal to a modifier to the full fecundity estimated by R. This matrix has as many rows and columns as the number of stages described in the stageframe for this dataset, and the rows and columns refer to these stages in the same order as in the stageframe. In many ways, it looks like a nearly empty population matrix, but notes the per-individual mean modifiers on fecundity for each stage that actually reproduces. Here, we first create a 0 matrix with dimensionality equal to the number of rows in `cypframe_raw`. Then we modify elements corresponding to fecundity by dividing fecundity evenly between dormant seeds (row 1) and germinating seeds (row 2).

    rep_cyp_raw <- matrix(0, 11, 11)
    rep_cyp_raw[1:2,7:11] <- 0.5

    rep_cyp_raw
    #>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
    #>  [1,]    0    0    0    0    0    0  0.5  0.5  0.5   0.5   0.5
    #>  [2,]    0    0    0    0    0    0  0.5  0.5  0.5   0.5   0.5
    #>  [3,]    0    0    0    0    0    0  0.0  0.0  0.0   0.0   0.0
    #>  [4,]    0    0    0    0    0    0  0.0  0.0  0.0   0.0   0.0
    #>  [5,]    0    0    0    0    0    0  0.0  0.0  0.0   0.0   0.0
    #>  [6,]    0    0    0    0    0    0  0.0  0.0  0.0   0.0   0.0
    #>  [7,]    0    0    0    0    0    0  0.0  0.0  0.0   0.0   0.0
    #>  [8,]    0    0    0    0    0    0  0.0  0.0  0.0   0.0   0.0
    #>  [9,]    0    0    0    0    0    0  0.0  0.0  0.0   0.0   0.0
    #> [10,]    0    0    0    0    0    0  0.0  0.0  0.0   0.0   0.0
    #> [11,]    0    0    0    0    0    0  0.0  0.0  0.0   0.0   0.0

Next we will provide some given transitions. In this case, we are providing the seed dormancy probability and germination rate, which in this case are provided as transitions from the dormant seed stage to another year of seed dormancy or to a first-year protocorm, respectively. We are also providing the survival-transition probabilities between different year protocorm stages (P1, P2, and P3), to the seedling stage (SL), and from the seedling stage to some of the adult stages (XSm, Sm, and D). Let’s start with the ahistorical case.

    cypover2r <- overwrite(stage3 = c("SD", "P1", "P2", "P3", "D", "XSm", "Sm", "SL", "SL"),
                           stage2 = c("SD", "SD", "P1", "P2", "P3", "P3", "P3", "P3", "SL"),
                           eststage3 = c(NA, NA, NA, NA, "D", "XSm", "Sm", NA, NA),
                           eststage2 = c(NA, NA, NA, NA, "D", "D", "D", NA, NA),
                           givenrate = c(0.1, 0.2, 0.2, 0.2, NA, NA, NA, 0.25, 0.4),
                           type = c("S", "S", "S", "S", "S", "S", "S", "S", "S"))
    cypover2r
    #>   stage3 stage2 stage1 eststage3 eststage2 eststage1 givenrate convtype
    #> 1     SD     SD     NA      <NA>      <NA>        NA      0.10        1
    #> 2     P1     SD     NA      <NA>      <NA>        NA      0.20        1
    #> 3     P2     P1     NA      <NA>      <NA>        NA      0.20        1
    #> 4     P3     P2     NA      <NA>      <NA>        NA      0.20        1
    #> 5      D     P3     NA         D         D        NA        NA        1
    #> 6    XSm     P3     NA       XSm         D        NA        NA        1
    #> 7     Sm     P3     NA        Sm         D        NA        NA        1
    #> 8     SL     P3     NA      <NA>      <NA>        NA      0.25        1
    #> 9     SL     SL     NA      <NA>      <NA>        NA      0.40        1

This overwrite table shows us that we have survival-transition probabilities (`convtype = 1`), that the given transitions are ahistorical, and particularly outlines probabilities for transitions that we cannot estimate with our dataset, which in this case refers to the immature stages of life. While 6 of these survival-transitions are given in the `givenrate` column, we also mark 3 of them as survival-transitions that we wish to use other estimates as proxies for. This is indicated via the `eststageX` columns, which have entries corresponding the stages to use as proxies (note that the `givenrate` entries are blank for these cases).

And now the historical case. Here we need to show the states in time step t-1 for this to work properly. Because of the extra time step, we will use the short-hand term `"rep"` to code for reproductive stages leading to the seeds and first-year protocorms that must survive to the next year.< /p>

    cypover3r <- overwrite(stage3 = c("SD", "SD", "P1", "P1", "P2", "P3", "D", "XSm", "Sm",
                           "SL", "SL", "SL"), stage2 = c("SD", "SD", "SD", "SD", "P1", "P2",
                           "P3", "P3", "P3", "P3", "SL", "SL"), stage1 = c("SD", "rep", "SD",
                           "rep", "SD", "P1", "P2", "P2", "P2", "P2", "P3", "SL"),
                           eststage3 = c(NA, NA, NA, NA, NA, NA, "D", "XSm", "Sm", NA, NA, NA),
                           eststage2 = c(NA, NA, NA, NA, NA, NA, "D", "D", "D", NA, NA, NA),
                           eststage1 = c(NA, NA, NA, NA, NA, NA, "D", "D", "D", NA, NA, NA),
                           givenrate = c(0.1, 0.1, 0.2, 0.2, 0.2, 0.2, NA, NA, NA, 0.25, 0.4, 0.4),
                           type = c("S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S"))
    cypover3r
    #>    stage3 stage2 stage1 eststage3 eststage2 eststage1 givenrate convtype
    #> 1      SD     SD     SD      <NA>      <NA>      <NA>      0.10        1
    #> 2      SD     SD    rep      <NA>      <NA>      <NA>      0.10        1
    #> 3      P1     SD     SD      <NA>      <NA>      <NA>      0.20        1
    #> 4      P1     SD    rep      <NA>      <NA>      <NA>      0.20        1
    #> 5      P2     P1     SD      <NA>      <NA>      <NA>      0.20        1
    #> 6      P3     P2     P1      <NA>      <NA>      <NA>      0.20        1
    #> 7       D     P3     P2         D         D         D        NA        1
    #> 8     XSm     P3     P2       XSm         D         D        NA        1
    #> 9      Sm     P3     P2        Sm         D         D        NA        1
    #> 10     SL     P3     P2      <NA>      <NA>      <NA>      0.25        1
    #> 11     SL     SL     P3      <NA>      <NA>      <NA>      0.40        1
    #> 12     SL     SL     SL      <NA>      <NA>      <NA>      0.40        1

Now we are read to create some matrices!

#### Step 3. Estimate matrices

We will begin with the creation of a set of ahistorical matrices for the Cypripedium candidum dataset. This should not take terribly long. The `rlefko2` function was created to deal with the construction of standard ahistorical matrices using raw data. Matrices may strongly differ, particularly if the demographic dataset is somewhat sparse.

    cypmatrix2r <- rlefko2(data = cypraw_v1, stageframe = cypframe_raw, year = "all",
                           patch = "all", stages = c("stage3", "stage2", "stage1"),
                           size = c("size3added", "size2added"),
                           repmatrix = rep_cyp_raw, overwrite = cypover2r,
                           yearcol = "year2", patchcol = "patchid",
                           indivcol = "individ")

The input for the `rlefko2()` function includes `year = "all"`, but can be set to focus on any set of years included within the data. Package `lefko3` actually includes a great deal of flexibility here, and can estimate many matrices covering all of the populations, patches, and years occurring in a specific dataset. The function-based matrix approach in the next section will showcase some more of this flexibility.

The output from this analysis includes is a `lefkoMat` object, which is a list object with the following elements:

A: a list of full population projection matrices, in order of population, patch, and year

U: a list of matrices showing only survival-transition elements, in the same order as A

F: a list of matrices showing only fecundity elements, in the same order as A

hstages: a data frame showing the order of paired stages (given if matrices are historical, otherwise NA)

ahstages: this is the stageframe used in analysis, with stages reordered and edited as they occur in the matrix

labels: a table showing the order of matrices, according to population, patch, and year

matrixqc: a short vector used in `summary` statements to describe the overall quality of each matrix

dataqc: a short vector used in `summary` statements to describe key sampling aspects of the dataset

Objects of class `lefkoMat` have their own `summary` statements, and we can use this to understand more about this particular object.

    summary(cypmatrix2r)
    #>
    #> This lefkoMat object contains 15 matrices.
    #>
    #> Each matrix is a square matrix with 11 rows and columns, and a total of 121 elements.
    #> A total of 258 survival-transitions were estimated, with 17.2 per matrix.
    #> A total of 66 fecundity transitions were estimated, with 4.4 per matrix.
    #>
    #> The dataset contains a total of 74 unique individuals and 320 unique transitions.
    #> NULL

We start off learning that 15 matrices were estimated, and we learn the dimensionality of those matrices. Of note here is the output telling us how many elements were actually estimated, both overall and per matrix, and the number of individuals and transitions the matrices are based on. Matrices are often overparameterized in population ecology, meaning that the number of elements estimated is quite high given the length of the data. It is typical for population ecologists to consider the total number of transitions in a dataset as a measure of the statistical power of a matrix, but the number of individuals used is just as important because each transition that an individual experiences is dependent on the other transitions that it also experiences. Indeed, this is the fundamental point that led to the development of historical matrices and of this package - the assumption that the status of an individual in the next time step is dependent only on its current state is too simplistic and leads to pseudoreplication. So, this output can be very helpful to understand the degree to which estimated matrices might be overparameterized.

Let’s look at a sample matrix here.

    cypmatrix2r$A[[1]]
    #>       [,1] [,2] [,3] [,4] [,5] [,6]      [,7] [,8] [,9] [,10] [,11]
    #>  [1,]  0.1  0.0  0.0 0.00  0.0    0 0.0000000  0.1  0.5   0.0     0
    #>  [2,]  0.2  0.0  0.0 0.00  0.0    0 0.0000000  0.1  0.5   0.0     0
    #>  [3,]  0.0  0.2  0.0 0.00  0.0    0 0.0000000  0.0  0.0   0.0     0
    #>  [4,]  0.0  0.0  0.2 0.00  0.0    0 0.0000000  0.0  0.0   0.0     0
    #>  [5,]  0.0  0.0  0.0 0.25  0.4    0 0.0000000  0.0  0.0   0.0     0
    #>  [6,]  0.0  0.0  0.0 0.00  0.0    0 0.0000000  0.0  0.0   0.0     0
    #>  [7,]  0.0  0.0  0.0 0.00  0.0    0 0.6363636  0.2  0.0   0.0     0
    #>  [8,]  0.0  0.0  0.0 0.00  0.0    0 0.2727273  0.6  0.5   0.0     0
    #>  [9,]  0.0  0.0  0.0 0.00  0.0    0 0.0000000  0.0  0.5   0.5     0
    #> [10,]  0.0  0.0  0.0 0.00  0.0    0 0.0000000  0.2  0.0   0.5     0
    #> [11,]  0.0  0.0  0.0 0.00  0.0    0 0.0000000  0.0  0.0   0.0     0

The reader will note that although this is an ahistorical matrix, it is predominantly composed of 0 elements. This is a result of the sparseness of the data, and will likely lead to different elements shifting between 0 and positive elements across time.

Now we will create some historical matrices. Historical matrix construction parses the data much more finely across many more stages than ahistorical matrix construction, so historical matrices are even more likely to differ strongly across time, particularly as the number of individuals in a dataset decreases. Let’s see what these matrices look like here.

    cypmatrix3r <- rlefko3(data = cypraw_v1, stageframe = cypframe_raw, year = "all",
                           patch = "all", stages = c("stage3", "stage2", "stage1"),
                           size = c("size3added", "size2added", "size1added"),
                           repmatrix = rep_cyp_raw, overwrite = cypover3r,
                           yearcol = "year2", patchcol = "patchid", indivcol = "individ")

    summary(cypmatrix3r)
    #>
    #> This lefkoMat object contains 12 matrices.
    #>
    #> Each matrix is a square matrix with 121 rows and columns, and a total of 14641 elements.
    #> A total of 388 survival-transitions were estimated, with 32.3333333333333 per matrix.
    #> A total of 72 fecundity transitions were estimated, with 6 per matrix.
    #>
    #> The dataset contains a total of 74 unique individuals and 320 unique transitions.
    #> NULL

There are at least two things to note here. First, there are 3 fewer matrices here than in the ahistorical case. There are 3 patches that we are estimating matrices for, and 6 years of data for each patch, leading to 5 possible ahistorical time transitions and 15 possible ahistorical matrices. Since historical matrices require 3 years of transitions, that means that only 4 historical transitions are possible per patch, leading to 12 total matrices. Second, the dimensionality of the matrices is the square of the dimensionality of the ahistorical matrices. This leads to vastly more matrix elements within each matrix, although it turns out that most of these matrix elements are structural 0s because they reflect impossible transitions. Indeed, in this case, although there are 14641 elements in each matrix, on average only 38.333 are actually estimated, with the rest set to 0.

Let’s look at the first matrix, corresponding to the transition from 2004 and 2005 to 2006 in the first patch. Because this is a huge matrix, we will only look at the top corner, and a middle section.

    cypmatrix3r$A[[1]][1:25,1:10]
    #>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
    #>  [1,]  0.1  0.0    0    0    0    0    0    0    0     0
    #>  [2,]  0.2  0.0    0    0    0    0    0    0    0     0
    #>  [3,]  0.0  0.0    0    0    0    0    0    0    0     0
    #>  [4,]  0.0  0.0    0    0    0    0    0    0    0     0
    #>  [5,]  0.0  0.0    0    0    0    0    0    0    0     0
    #>  [6,]  0.0  0.0    0    0    0    0    0    0    0     0
    #>  [7,]  0.0  0.0    0    0    0    0    0    0    0     0
    #>  [8,]  0.0  0.0    0    0    0    0    0    0    0     0
    #>  [9,]  0.0  0.0    0    0    0    0    0    0    0     0
    #> [10,]  0.0  0.0    0    0    0    0    0    0    0     0
    #> [11,]  0.0  0.0    0    0    0    0    0    0    0     0
    #> [12,]  0.0  0.0    0    0    0    0    0    0    0     0
    #> [13,]  0.0  0.0    0    0    0    0    0    0    0     0
    #> [14,]  0.0  0.2    0    0    0    0    0    0    0     0
    #> [15,]  0.0  0.0    0    0    0    0    0    0    0     0
    #> [16,]  0.0  0.0    0    0    0    0    0    0    0     0
    #> [17,]  0.0  0.0    0    0    0    0    0    0    0     0
    #> [18,]  0.0  0.0    0    0    0    0    0    0    0     0
    #> [19,]  0.0  0.0    0    0    0    0    0    0    0     0
    #> [20,]  0.0  0.0    0    0    0    0    0    0    0     0
    #> [21,]  0.0  0.0    0    0    0    0    0    0    0     0
    #> [22,]  0.0  0.0    0    0    0    0    0    0    0     0
    #> [23,]  0.0  0.0    0    0    0    0    0    0    0     0
    #> [24,]  0.0  0.0    0    0    0    0    0    0    0     0
    #> [25,]  0.0  0.0    0    0    0    0    0    0    0     0
    cypmatrix3r$A[[1]][66:90,73:81]
    #>             [,1]      [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
    #>  [1,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #>  [2,] 0.07142857 0.0000000    0    0    0    0    0    0    0
    #>  [3,] 0.07142857 0.0000000    0    0    0    0    0    0    0
    #>  [4,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #>  [5,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #>  [6,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #>  [7,] 0.14285714 0.0000000    0    0    0    0    0    0    0
    #>  [8,] 0.71428571 0.0000000    0    0    0    0    0    0    0
    #>  [9,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #> [10,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #> [11,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #> [12,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #> [13,] 0.00000000 0.3333333    0    0    0    0    0    0    0
    #> [14,] 0.00000000 0.3333333    0    0    0    0    0    0    0
    #> [15,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #> [16,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #> [17,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #> [18,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #> [19,] 0.00000000 0.6666667    0    0    0    0    0    0    0
    #> [20,] 0.00000000 0.3333333    0    0    0    0    0    0    0
    #> [21,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #> [22,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #> [23,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #> [24,] 0.00000000 0.0000000    0    0    0    0    0    0    0
    #> [25,] 0.00000000 0.0000000    0    0    0    0    0    0    0

The full matrix is actually not shown here, but we can focus on portions of it if we wish. These matrices may also be exported to Excel or another spreadsheet program to look over in detail.

Next we will create a mean ahistorical matrix. Note that we will treat time arithmetically in this mean, meaning that matrix elements in the mean matrix are arithmetic means across time. With regards to which should be estimated, arithmetic or geometric mean: Our own feeling is that the geometric mean is a better choice in principle when taking a mean matrix elements across time, but the literature is generally full of the former case. The reason that the geometric mean may be better is that survival and reproduction are multiplicative over time, and the overall impact of variability in multiplicative processes is not symmetrical with respect to outcomes. For example, if I take the arithmetic mean matrix of all 4 matrices resulting from a 5 year study, and multiply the matrix four times over by the original vector of individuals across their respective stages, then I will end up projecting more individuals alive in the 5th year than actually exist! In multiplicative (or geometric) processes, the impact of lower values on the outcome is greater than the impact of higher values. Thus, in principle, the geometric mean is preferable. However, in practice, occasional lack of individuals moving through some stages may lead to sparse instances of elements dropping to 0. These situations cause the geometric mean to drop to 0, as well. In raw matrix projection, such situations may end up leading to geometric mean matrices overwhelmingly composed of 0s if the stages are not designed to absorb as many individuals as possible in each year. This can lead to a situation in the geometric mean projects a far lower number of individuals alive after some period.

The `lmean()` function in `lefko3` incorporates a number of interesting options to help deal with these situations. Let’s do a comparison. First the pure arithmetic grand mean matrix, which is the last A matrix in this object.

    cypr2mean <- lmean(cypmatrix2r, time = "arithmetic")
    cypr2mean$A[[4]]
    #>       [,1] [,2] [,3]       [,4] [,5]       [,6]        [,7]       [,8]       [,9]      [,10]      [,11]
    #>  [1,]  0.1  0.0  0.0 0.00000000  0.0 0.00000000 0.007870370 0.07174603 0.18898148 0.29238095 0.46666667
    #>  [2,]  0.2  0.0  0.0 0.00000000  0.0 0.00000000 0.007870370 0.07174603 0.18898148 0.29238095 0.46666667
    #>  [3,]  0.0  0.2  0.0 0.00000000  0.0 0.00000000 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000
    #>  [4,]  0.0  0.0  0.2 0.00000000  0.0 0.00000000 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000
    #>  [5,]  0.0  0.0  0.0 0.25000000  0.4 0.00000000 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000
    #>  [6,]  0.0  0.0  0.0 0.02222222  0.0 0.02222222 0.067255892 0.02814815 0.01944444 0.00000000 0.00000000
    #>  [7,]  0.0  0.0  0.0 0.21111111  0.0 0.21111111 0.519745070 0.22907407 0.06944444 0.06666667 0.00000000
    #>  [8,]  0.0  0.0  0.0 0.15555556  0.0 0.15555556 0.254737855 0.40314815 0.19611111 0.18222222 0.02222222
    #>  [9,]  0.0  0.0  0.0 0.00000000  0.0 0.05555556 0.007407407 0.21351852 0.48814815 0.16793651 0.00000000
    #> [10,]  0.0  0.0  0.0 0.00000000  0.0 0.00000000 0.037037037 0.06500000 0.22685185 0.54349206 0.08888889
    #> [11,]  0.0  0.0  0.0 0.00000000  0.0 0.02222222 0.000000000 0.00000000 0.00000000 0.03968254 0.48888889
    writeLines("\nColumn sums (survival probabilities) for grand arithmetic mean matrix")
    #>
    #> Column sums (survival probabilities) for grand arithmetic mean matrix
    summary(colSums(cypr2mean$U[[4]]))
    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    #>  0.2000  0.3500  0.6000  0.6028  0.9125  1.0000

For comparison, let’s look first at a geometric mean matrix.

    lmean(cypmatrix2r, time = "geometric")$A[[4]]
    #>       [,1] [,2] [,3] [,4] [,5] [,6]      [,7]       [,8]      [,9]     [,10]     [,11]
    #>  [1,]  0.1  0.0  0.0 0.00  0.0    0 0.0000000 0.00000000 0.0000000 0.1123006 0.0000000
    #>  [2,]  0.2  0.0  0.0 0.00  0.0    0 0.0000000 0.00000000 0.0000000 0.1123006 0.0000000
    #>  [3,]  0.0  0.2  0.0 0.00  0.0    0 0.0000000 0.00000000 0.0000000 0.0000000 0.0000000
    #>  [4,]  0.0  0.0  0.2 0.00  0.0    0 0.0000000 0.00000000 0.0000000 0.0000000 0.0000000
    #>  [5,]  0.0  0.0  0.0 0.25  0.4    0 0.0000000 0.00000000 0.0000000 0.0000000 0.0000000
    #>  [6,]  0.0  0.0  0.0 0.00  0.0    0 0.0000000 0.00000000 0.0000000 0.0000000 0.0000000
    #>  [7,]  0.0  0.0  0.0 0.00  0.0    0 0.4952254 0.08359622 0.0000000 0.0000000 0.0000000
    #>  [8,]  0.0  0.0  0.0 0.00  0.0    0 0.0000000 0.31081295 0.0000000 0.0000000 0.0000000
    #>  [9,]  0.0  0.0  0.0 0.00  0.0    0 0.0000000 0.00000000 0.1469767 0.0000000 0.0000000
    #> [10,]  0.0  0.0  0.0 0.00  0.0    0 0.0000000 0.00000000 0.0000000 0.3689368 0.0000000
    #> [11,]  0.0  0.0  0.0 0.00  0.0    0 0.0000000 0.00000000 0.0000000 0.0000000 0.2901835

    writeLines("\nColumn sums (survival probabilities) for grand geometric mean matrix")
    #>
    #> Column sums (survival probabilities) for grand geometric mean matrix
    summary(colSums(lmean(cypmatrix2r, time = "geometric")$U[[4]]))
    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    #>  0.0000  0.2000  0.2902  0.2769  0.3817  0.4952

Note that this matrix suggests extremely low survival probabilities relative to the arithmetic case. This is due to the fact that our data is relatively sparse. So, a single year in which no individuals start in a particular stage will result in a zero entry in transition probabilities starting from that stage, leading the geometric mean matrix to have zero entries in that column. In comparison, we can look at a geometric mean that scraps 0s unless they are structural.

    lmean(cypmatrix2r, time = "geometric", sparse = TRUE)$A[[4]]
    #>       [,1] [,2] [,3]      [,4] [,5]      [,6]       [,7]      [,8]       [,9]      [,10]     [,11]
    #>  [1,]  0.1  0.0  0.0 0.0000000  0.0 0.0000000 0.03935185 0.1793293 0.28972162 0.38437657 0.5206737
    #>  [2,]  0.2  0.0  0.0 0.0000000  0.0 0.0000000 0.03935185 0.1793293 0.28972162 0.38437657 0.5206737
    #>  [3,]  0.0  0.2  0.0 0.0000000  0.0 0.0000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
    #>  [4,]  0.0  0.0  0.2 0.0000000  0.0 0.0000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
    #>  [5,]  0.0  0.0  0.0 0.2500000  0.4 0.0000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
    #>  [6,]  0.0  0.0  0.0 0.1111111  0.0 0.1111111 0.20004424 0.1407407 0.09722222 0.00000000 0.0000000
    #>  [7,]  0.0  0.0  0.0 0.5088982  0.0 0.5088982 0.49522542 0.2967068 0.34722222 0.33333333 0.0000000
    #>  [8,]  0.0  0.0  0.0 0.5257834  0.0 0.5257834 0.36473518 0.4117641 0.43129808 0.51356794 0.1111111
    #>  [9,]  0.0  0.0  0.0 0.0000000  0.0 0.2777778 0.03703704 0.2888765 0.60017161 0.39258730 0.0000000
    #> [10,]  0.0  0.0  0.0 0.0000000  0.0 0.0000000 0.14814815 0.1933207 0.39360704 0.56387133 0.3027494
    #> [11,]  0.0  0.0  0.0 0.0000000  0.0 0.1111111 0.00000000 0.0000000 0.00000000 0.08908708 0.5031646

    writeLines("\nColumn sums (survival probabilities) for grand geometric mean matrix")
    #>
    #> Column sums (survival probabilities) for grand geometric mean matrix
    summary(colSums(lmean(cypmatrix2r, time = "geometric", sparse = TRUE)$U[[4]]))
    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    #>   0.200   0.350   1.245   1.026   1.465   1.892

Clearly, no stage can have a probability of survival over 1.0, but eliminating 0 entries from the estimation of mean matrices can sometimes result in this effect under both arithmetic and geometric scenarios. So, here we scrap this approach in favor of the arithmetic mean.

And now the historical grand mean matrix, with a peek at a middle portion with some non-zero values.

    cypr3mean <- lmean(cypmatrix3r, time = "arithmetic")
    cypr3mean$A[[4]][66:90,73:80]
    #>              [,1]       [,2] [,3] [,4] [,5] [,6] [,7] [,8]
    #>  [1,] 0.000000000 0.00000000    0    0    0    0    0    0
    #>  [2,] 0.005952381 0.00000000    0    0    0    0    0    0
    #>  [3,] 0.005952381 0.00000000    0    0    0    0    0    0
    #>  [4,] 0.000000000 0.00000000    0    0    0    0    0    0
    #>  [5,] 0.000000000 0.00000000    0    0    0    0    0    0
    #>  [6,] 0.000000000 0.00000000    0    0    0    0    0    0
    #>  [7,] 0.136904762 0.00000000    0    0    0    0    0    0
    #>  [8,] 0.537301587 0.00000000    0    0    0    0    0    0
    #>  [9,] 0.225000000 0.00000000    0    0    0    0    0    0
    #> [10,] 0.000000000 0.00000000    0    0    0    0    0    0
    #> [11,] 0.000000000 0.00000000    0    0    0    0    0    0
    #> [12,] 0.000000000 0.00000000    0    0    0    0    0    0
    #> [13,] 0.000000000 0.02777778    0    0    0    0    0    0
    #> [14,] 0.000000000 0.02777778    0    0    0    0    0    0
    #> [15,] 0.000000000 0.00000000    0    0    0    0    0    0
    #> [16,] 0.000000000 0.00000000    0    0    0    0    0    0
    #> [17,] 0.000000000 0.00000000    0    0    0    0    0    0
    #> [18,] 0.000000000 0.00000000    0    0    0    0    0    0
    #> [19,] 0.000000000 0.30138889    0    0    0    0    0    0
    #> [20,] 0.000000000 0.09861111    0    0    0    0    0    0
    #> [21,] 0.000000000 0.10000000    0    0    0    0    0    0
    #> [22,] 0.000000000 0.02777778    0    0    0    0    0    0
    #> [23,] 0.000000000 0.00000000    0    0    0    0    0    0
    #> [24,] 0.000000000 0.00000000    0    0    0    0    0    0
    #> [25,] 0.000000000 0.00000000    0    0    0    0    0    0

Do not fear the prevalence of 0s in this matrix - this is normal, both because many elements are structural 0s and so cannot equal anything else, and because this is a raw matrix, meaning that transitions that do not actually occur in the dataset cannot equal anything other than 0.

To understand the dominance of structural 0s in the historical case, let’s take a look at the `hstages` object associated with this mean matrix.

    cypr3mean$hstages
    #>     stcod3 stcod2 stage3 stage2
    #> 1       SD     SD      1      1
    #> 2       P1     SD      2      1
    #> 3       P2     SD      3      1
    #> 4       P3     SD      4      1
    #> 5       SL     SD      5      1
    #> 6        D     SD      6      1
    #> 7      XSm     SD      7      1
    #> 8       Sm     SD      8      1
    #> 9       Md     SD      9      1
    #> 10      Lg     SD     10      1
    #> 11     XLg     SD     11      1
    #> 13      SD     P1      1      2
    #> 14      P1     P1      2      2
    #> 15      P2     P1      3      2
    #> 16      P3     P1      4      2
    #> 17      SL     P1      5      2
    #> 18       D     P1      6      2
    #> 19     XSm     P1      7      2
    #> 20      Sm     P1      8      2
    #> 21      Md     P1      9      2
    #> 22      Lg     P1     10      2
    #> 23     XLg     P1     11      2
    #> 25      SD     P2      1      3
    #> 26      P1     P2      2      3
    #> 27      P2     P2      3      3
    #> 28      P3     P2      4      3
    #> 29      SL     P2      5      3
    #> 30       D     P2      6      3
    #> 31     XSm     P2      7      3
    #> 32      Sm     P2      8      3
    #> 33      Md     P2      9      3
    #> 34      Lg     P2     10      3
    #> 35     XLg     P2     11      3
    #> 37      SD     P3      1      4
    #> 38      P1     P3      2      4
    #> 39      P2     P3      3      4
    #> 40      P3     P3      4      4
    #> 41      SL     P3      5      4
    #> 42       D     P3      6      4
    #> 43     XSm     P3      7      4
    #> 44      Sm     P3      8      4
    #> 45      Md     P3      9      4
    #> 46      Lg     P3     10      4
    #> 47     XLg     P3     11      4
    #> 49      SD     SL      1      5
    #> 50      P1     SL      2      5
    #> 51      P2     SL      3      5
    #> 52      P3     SL      4      5
    #> 53      SL     SL      5      5
    #> 54       D     SL      6      5
    #> 55     XSm     SL      7      5
    #> 56      Sm     SL      8      5
    #> 57      Md     SL      9      5
    #> 58      Lg     SL     10      5
    #> 59     XLg     SL     11      5
    #> 61      SD      D      1      6
    #> 62      P1      D      2      6
    #> 63      P2      D      3      6
    #> 64      P3      D      4      6
    #> 65      SL      D      5      6
    #> 66       D      D      6      6
    #> 67     XSm      D      7      6
    #> 68      Sm      D      8      6
    #> 69      Md      D      9      6
    #> 70      Lg      D     10      6
    #> 71     XLg      D     11      6
    #> 73      SD    XSm      1      7
    #> 74      P1    XSm      2      7
    #> 75      P2    XSm      3      7
    #> 76      P3    XSm      4      7
    #> 77      SL    XSm      5      7
    #> 78       D    XSm      6      7
    #> 79     XSm    XSm      7      7
    #> 80      Sm    XSm      8      7
    #> 81      Md    XSm      9      7
    #> 82      Lg    XSm     10      7
    #> 83     XLg    XSm     11      7
    #> 85      SD     Sm      1      8
    #> 86      P1     Sm      2      8
    #> 87      P2     Sm      3      8
    #> 88      P3     Sm      4      8
    #> 89      SL     Sm      5      8
    #> 90       D     Sm      6      8
    #> 91     XSm     Sm      7      8
    #> 92      Sm     Sm      8      8
    #> 93      Md     Sm      9      8
    #> 94      Lg     Sm     10      8
    #> 95     XLg     Sm     11      8
    #> 97      SD     Md      1      9
    #> 98      P1     Md      2      9
    #> 99      P2     Md      3      9
    #> 100     P3     Md      4      9
    #> 101     SL     Md      5      9
    #> 102      D     Md      6      9
    #> 103    XSm     Md      7      9
    #> 104     Sm     Md      8      9
    #> 105     Md     Md      9      9
    #> 106     Lg     Md     10      9
    #> 107    XLg     Md     11      9
    #> 109     SD     Lg      1     10
    #> 110     P1     Lg      2     10
    #> 111     P2     Lg      3     10
    #> 112     P3     Lg      4     10
    #> 113     SL     Lg      5     10
    #> 114      D     Lg      6     10
    #> 115    XSm     Lg      7     10
    #> 116     Sm     Lg      8     10
    #> 117     Md     Lg      9     10
    #> 118     Lg     Lg     10     10
    #> 119    XLg     Lg     11     10
    #> 121     SD    XLg      1     11
    #> 122     P1    XLg      2     11
    #> 123     P2    XLg      3     11
    #> 124     P3    XLg      4     11
    #> 125     SL    XLg      5     11
    #> 126      D    XLg      6     11
    #> 127    XSm    XLg      7     11
    #> 128     Sm    XLg      8     11
    #> 129     Md    XLg      9     11
    #> 130     Lg    XLg     10     11
    #> 131    XLg    XLg     11     11

There are 121 pairs of ahistorical stages, and these pairs correspond to the rows and columns of the historical matrices output by `rlefko3` in this case. The pairs are interpreted so that matrix columns represent the states of the individual in times t-1 and t, and the rows represent states in times t and t+1. For an element in the matrix to contain a number other than 0, it must represent the same stage at time t in both the column stage pairs and the row stage pairs. The element [1, 1], for example, represents the transition probability from dormant seed at times t-1 and t (column pair) to times t and t+1 (row pair) - the time t stages match, and so this element is possible. However, element [1, 2] represents the transition probability from seed in time t-1 and protocorm 1 in time t (column pair), to dormant seed in time t and in time t+1 (row pair). Clearly [1, 2] is a structural 0 because it is impossible for an individual to be in two stages at once in time t!

Error-checking is more difficult with historical matrices because they are typically one or two orders of magnitude bigger than their ahistorical counterparts. Take a look at the column sums here to see the difficulty.

    colSums(cypr3mean$U[[4]])
    #>   [1] 0.30000000 0.20000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
    #>  [13] 0.00000000 0.20000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
    #>  [25] 0.00000000 0.33333333 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
    #>  [37] 0.00000000 0.40000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
    #>  [49] 0.40000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
    #>  [61] 0.08333333 0.16666667 0.16666667 0.00000000 0.00000000 0.08333333 0.30000000 0.00000000 0.00000000 0.00000000 0.00000000 0.50000000
    #>  [73] 0.89920635 0.52777778 0.08333333 0.25000000 0.00000000 0.30000000 0.00000000 0.00000000 0.00000000 0.00000000 0.16666667 0.77500000
    #>  [85] 0.88888889 0.50000000 0.41666667 0.00000000 0.30000000 0.00000000 0.00000000 0.00000000 0.00000000 0.16666667 0.25000000 0.50000000
    #>  [97] 0.75000000 0.50000000 0.00000000 0.30000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.08333333 0.41666667 0.50000000
    #> [109] 0.83333333 0.08333333 0.30000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.08333333 0.00000000 0.16666667
    #> [121] 0.58333333

While not too bad here, other historical matrices often have more than 100 columns, sometimes many, many more (some historical matrices used in Shefferson et al. (2014) had dimensions of over 2500 x 2500!). In these cases we can assess the distribution of survival estimates for historical stages, which is given as the set of column sums in the survival-transition matrix, as below.

    summary(colSums(cypr3mean$U[[4]]))
    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    #>  0.0000  0.0000  0.0000  0.1137  0.1667  0.8992

As long as all of the numbers above are between 0 and 1, then all is probably well. Fine-scale error-checking would require outputting the matrix into a spreadsheet and assessing it using the `hstages` output as a guide to what the elements refer to.

#### Step 4. Analyses

Now let’s estimate the deterministic population growth rate in each case. We will start by looking at the annual population growth rate estimated from the ahistorical analyses, followed by the population growth rate associated with the mean matrix from that analysis. First the annual ahistorical matrices.

    lambda3(cypmatrix2r)
    #>    pop patch year2    lambda
    #> 1   NA     A  2004 0.9782070
    #> 2   NA     A  2005 1.0000000
    #> 3   NA     A  2006 1.0000000
    #> 4   NA     A  2007 1.0000000
    #> 5   NA     A  2008 1.0000000
    #> 6   NA     B  2004 0.8964607
    #> 7   NA     B  2005 0.8597354
    #> 8   NA     B  2006 1.0000000
    #> 9   NA     B  2007 1.0065871
    #> 10  NA     B  2008 1.0055269
    #> 11  NA     C  2004 1.0000000
    #> 12  NA     C  2005 1.0000000
    #> 13  NA     C  2006 1.0000000
    #> 14  NA     C  2007 1.0000000
    #> 15  NA     C  2008 1.0039623

Now the ahistorical mean matrix.

    lambda3(cypr2mean)
    #>   pop patch    lambda
    #> 1   1     A 0.8992573
    #> 2   1     B 0.9321193
    #> 3   1     C 0.9577990
    #> 4   1   All 0.9332170

And we will finish by looking at the same numbers for the historical analyses. First the annual matrices.

    lambda3(cypmatrix3r)
    #>    pop patch year2    lambda
    #> 1   NA     A  2005 0.7142857
    #> 2   NA     A  2006 1.0000000
    #> 3   NA     A  2007 1.0000000
    #> 4   NA     A  2008 1.0000000
    #> 5   NA     B  2005 0.6666667
    #> 6   NA     B  2006 1.0000000
    #> 7   NA     B  2007 1.0000000
    #> 8   NA     B  2008 1.0000000
    #> 9   NA     C  2005 1.0000000
    #> 10  NA     C  2006 1.0000000
    #> 11  NA     C  2007 1.0000000
    #> 12  NA     C  2008 1.0000000

Now the mean.

    lambda3(cypr3mean)
    #>   pop patch    lambda
    #> 1   1     A 0.7802134
    #> 2   1     B 0.7601081
    #> 3   1     C 1.0000000
    #> 4   1   All 0.6905696

Readers will likely observe both that there are fewer lambda estimates in the historical case, and that the mean lambda is typically lower. Both are expected. First, because there are 6 years of data, there are three ahistorical transitions possible for estimation: year 1 to 2, year 2 to 3, and year 3 to 4. However, in the historical case, only two are possible: from years 1 and 2 to years 2 and 3, and from years 2 and 3 to years 3 and 4. Second, historical matrices integrate temporal autocorrelation in vital rates in ways that ahistorical matrices generally cannot, and these autocorrelations are likely to be most strongly impacted by trade-offs operating across years (Shefferson and Roach 2010). One particularly common such trade-off is the cost of growth: an individual that grows a great deal in one time step due to great environmental conditions in that year might pay a large cost of survival, growth, or reproduction in the next if those environmental conditions deteriorate (Shefferson, Warren II, and Pulliam 2014). While we do not argue that the drop in lambda is due to this specific trade-off, we do argue that this lambda is likely to be more realistic than the higher lambda estimated in the ahistorical case.

We can also take a peek at the stable stage distributions, as follows for the ahistorical case:

    stablestage3(cypr2mean)
    #>    matrix new_stage_id orig_stage_id original_size     ss_prop
    #> 1       1            1            SD           0.0 0.090037090
    #> 2       1            2            P1           0.0 0.100049436
    #> 3       1            3            P2           0.0 0.022251576
    #> 4       1            4            P3           0.0 0.004948866
    #> 5       1            5            SL           0.0 0.002478141
    #> 6       1            6             D           0.0 0.030513881
    #> 7       1            7           XSm           1.0 0.249352147
    #> 8       1            8            Sm           2.5 0.290121522
    #> 9       1            9            Md           4.5 0.101212097
    #> 10      1           10            Lg           8.0 0.109035244
    #> 11      1           11           XLg          17.5 0.000000000
    #> 12      2            1            SD           0.0 0.147677868
    #> 13      2            2            P1           0.0 0.163521115
    #> 14      2            3            P2           0.0 0.035085868
    #> 15      2            4            P3           0.0 0.007528201
    #> 16      2            5            SL           0.0 0.003536880
    #> 17      2            6             D           0.0 0.027764530
    #> 18      2            7           XSm           1.0 0.172012760
    #> 19      2            8            Sm           2.5 0.199776799
    #> 20      2            9            Md           4.5 0.093236400
    #> 21      2           10            Lg           8.0 0.109216888
    #> 22      2           11           XLg          17.5 0.040642691
    #> 23      3            1            SD           0.0 0.085711562
    #> 24      3            2            P1           0.0 0.094660376
    #> 25      3            3            P2           0.0 0.019766226
    #> 26      3            4            P3           0.0 0.004127438
    #> 27      3            5            SL           0.0 0.001849871
    #> 28      3            6             D           0.0 0.020105859
    #> 29      3            7           XSm           1.0 0.178308076
    #> 30      3            8            Sm           2.5 0.182244061
    #> 31      3            9            Md           4.5 0.221104786
    #> 32      3           10            Lg           8.0 0.192121744
    #> 33      3           11           XLg          17.5 0.000000000
    #> 34      4            1            SD           0.0 0.114482899
    #> 35      4            2            P1           0.0 0.126750486
    #> 36      4            3            P2           0.0 0.027164217
    #> 37      4            4            P3           0.0 0.005821612
    #> 38      4            5            SL           0.0 0.002729480
    #> 39      4            6             D           0.0 0.022921317
    #> 40      4            7           XSm           1.0 0.177582482
    #> 41      4            8            Sm           2.5 0.203366188
    #> 42      4            9            Md           4.5 0.158605194
    #> 43      4           10            Lg           8.0 0.146358608
    #> 44      4           11           XLg          17.5 0.014217515

We can do this for the historical case as well, but need to bear in mind that the distribution is in paired stages, and so can take some effort to interpret properly. We can look first at the paired stage distribution. Let’s see the summary and first 6 entries.

    cypr3ss <- stablestage3(cypr3mean)
    summary(cypr3ss$hist)
    #>      matrix     orig_stage_id_2    orig_stage_id_1    new_stage_id_2 new_stage_id_1    ss_prop
    #>  Min.   :1.00   Length:484         Length:484         Min.   : 1     Min.   : 1     Min.   :0.0000000
    #>  1st Qu.:1.75   Class :character   Class :character   1st Qu.: 3     1st Qu.: 3     1st Qu.:0.0000000
    #>  Median :2.50   Mode  :character   Mode  :character   Median : 6     Median : 6     Median :0.0000000
    #>  Mean   :2.50                                         Mean   : 6     Mean   : 6     Mean   :0.0082645
    #>  3rd Qu.:3.25                                         3rd Qu.: 9     3rd Qu.: 9     3rd Qu.:0.0000079
    #>  Max.   :4.00                                         Max.   :11     Max.   :11     Max.   :0.5273180
    head(cypr3ss$hist)
    #>   matrix orig_stage_id_2 orig_stage_id_1 new_stage_id_2 new_stage_id_1     ss_prop
    #> 1      1              SD              SD              1              1 0.003083154
    #> 2      1              P1              SD              2              1 0.006166308
    #> 3      1              P2              SD              3              1 0.000000000
    #> 4      1              P3              SD              4              1 0.000000000
    #> 5      1              SL              SD              5              1 0.000000000
    #> 6      1               D              SD              6              1 0.000000000

And we can look at the resulting life history stable stage distribution associated with the historical case.

    cypr3ss$ahist
    #>    matrix new_stage_id      ss_prop
    #> 1       1            1 0.0240552215
    #> 2       1            2 0.0271383756
    #> 3       1            3 0.0015806534
    #> 4       1            4 0.0004052074
    #> 5       1            5 0.0002664555
    #> 6       1            6 0.0578366170
    #> 7       1            7 0.4955790205
    #> 8       1            8 0.3584457196
    #> 9       1            9 0.0317567933
    #> 10      1           10 0.0029359364
    #> 11      1           11 0.0000000000
    #> 12      2            1 0.0931549377
    #> 13      2            2 0.1054104247
    #> 14      2            3 0.0064493456
    #> 15      2            4 0.0016969480
    #> 16      2            5 0.0011780848
    #> 17      2            6 0.0807390019
    #> 18      2            7 0.2526776338
    #> 19      2            8 0.3393131233
    #> 20      2            9 0.0370767622
    #> 21      2           10 0.0518189154
    #> 22      2           11 0.0304848226
    #> 23      3            1 0.2197158418
    #> 24      3            2 0.2416874140
    #> 25      3            3 0.0087886169
    #> 26      3            4 0.0017577354
    #> 27      3            5 0.0007323598
    #> 28      3            6 0.0000000000
    #> 29      3            7 0.0000000000
    #> 30      3            8 0.0000000000
    #> 31      3            9 0.0000000000
    #> 32      3           10 0.0000000000
    #> 33      3           11 0.5273180322
    #> 34      4            1 0.0662983940
    #> 35      4            2 0.0758989111
    #> 36      4            3 0.0055609200
    #> 37      4            4 0.0016105376
    #> 38      4            5 0.0013856524
    #> 39      4            6 0.0674429087
    #> 40      4            7 0.4016556123
    #> 41      4            8 0.2403823061
    #> 42      4            9 0.0930652020
    #> 43      4           10 0.0453846335
    #> 44      4           11 0.0013149224

These values are considerably different from those of the purely ahistorical stable stage distribution. Indeed, here seeds and protocorms are a much smaller part of the stable stage structure, while dormant, extra small, and small adults are a much larger component. It seems that individual history has a strong impact here.

Let’s look at the reproductive values next. First, the ahistorical case.

    repvalue3(cypr2mean)
    #>    matrix new_stage_id orig_stage_id original_size left_vector   rep_value
    #> 1       1            1            SD           0.0   0.0014030    1.000000
    #> 2       1            2            P1           0.0   0.0056066    3.996151
    #> 3       1            3            P2           0.0   0.0252089   17.967855
    #> 4       1            4            P3           0.0   0.1133465   80.788667
    #> 5       1            5            SL           0.0   0.0000000    0.000000
    #> 6       1            6             D           0.0   0.1533835  109.325374
    #> 7       1            7           XSm           1.0   0.2536043  180.758589
    #> 8       1            8            Sm           2.5   0.3188277  227.247113
    #> 9       1            9            Md           4.5   0.5400528  384.927156
    #> 10      1           10            Lg           8.0   0.7108684  506.677406
    #> 11      1           11           XLg          17.5   0.0000000    0.000000
    #> 12      2            1            SD           0.0  -0.0030046    1.000000
    #> 13      2            2            P1           0.0  -0.0125010    4.160620
    #> 14      2            3            P2           0.0  -0.0582622   19.391000
    #> 15      2            4            P3           0.0  -0.2715368   90.373694
    #> 16      2            5            SL           0.0   0.0000000    0.000000
    #> 17      2            6             D           0.0  -0.2952139   98.253977
    #> 18      2            7           XSm           1.0  -0.3418877  113.788092
    #> 19      2            8            Sm           2.5  -0.3693806  122.938361
    #> 20      2            9            Md           4.5  -0.4824601  160.573820
    #> 21      2           10            Lg           8.0  -0.4897617  163.003961
    #> 22      2           11           XLg          17.5  -0.3310488  110.180656
    #> 23      3            1            SD           0.0  -0.0003103    1.000000
    #> 24      3            2            P1           0.0  -0.0013308    4.288753
    #> 25      3            3            P2           0.0  -0.0063730   20.538189
    #> 26      3            4            P3           0.0  -0.0305204   98.357718
    #> 27      3            5            SL           0.0   0.0000000    0.000000
    #> 28      3            6             D           0.0  -0.0713462  229.926523
    #> 29      3            7           XSm           1.0  -0.2923241  942.069288
    #> 30      3            8            Sm           2.5  -0.3880920 1250.699323
    #> 31      3            9            Md           4.5  -0.3910292 1260.165002
    #> 32      3           10            Lg           8.0  -0.3836712 1236.452465
    #> 33      3           11           XLg          17.5  -0.6765793 2180.403803
    #> 34      4            1            SD           0.0  -0.0018065    1.000000
    #> 35      4            2            P1           0.0  -0.0075259    4.166012
    #> 36      4            3            P2           0.0  -0.0351166   19.439026
    #> 37      4            4            P3           0.0  -0.1638568   90.704013
    #> 38      4            5            SL           0.0   0.0000000    0.000000
    #> 39      4            6             D           0.0  -0.1992797  110.312593
    #> 40      4            7           XSm           1.0  -0.3691083  204.322336
    #> 41      4            8            Sm           2.5  -0.4536172  251.102795
    #> 42      4            9            Md           4.5  -0.5397128  298.761583
    #> 43      4           10            Lg           8.0  -0.5288866  292.768669
    #> 44      4           11           XLg          17.5  -0.1382934   76.553224

Now the historical case, beginning with the stage pair distribution. We will first look at the summary.

    cypr3rv <- repvalue3(cypr3mean)
    summary(cypr3rv$hist)
    #>      matrix     orig_stage_id_2    orig_stage_id_1    new_stage_id_2 new_stage_id_1  left_vector          rep_value
    #>  Min.   :1.00   Length:484         Length:484         Min.   : 1     Min.   : 1     Min.   :-0.561342   Min.   : 0.0000
    #>  1st Qu.:1.75   Class :character   Class :character   1st Qu.: 3     1st Qu.: 3     1st Qu.: 0.000000   1st Qu.: 0.0000
    #>  Median :2.50   Mode  :character   Mode  :character   Median : 6     Median : 6     Median : 0.000000   Median : 0.0000
    #>  Mean   :2.50                                         Mean   : 6     Mean   : 6     Mean   :-0.006077   Mean   : 0.1994
    #>  3rd Qu.:3.25                                         3rd Qu.: 9     3rd Qu.: 9     3rd Qu.: 0.000000   3rd Qu.: 0.0000
    #>  Max.   :4.00                                         Max.   :11     Max.   :11     Max.   : 1.000000   Max.   :27.5077

Now the first 6 rows.

    head(cypr3rv$hist)
    #>   matrix orig_stage_id_2 orig_stage_id_1 new_stage_id_2 new_stage_id_1 left_vector rep_value
    #> 1      1              SD              SD              1              1           0         0
    #> 2      1              P1              SD              2              1           0         0
    #> 3      1              P2              SD              3              1           0         0
    #> 4      1              P3              SD              4              1           0         0
    #> 5      1              SL              SD              5              1           0         0
    #> 6      1               D              SD              6              1           0         0

And now the life history stage distribution associated with the historical case.

    cypr3rv$ahist
    #>    matrix new_stage_id  rep_value
    #> 1       1            1 0.00000000
    #> 2       1            2 0.00000000
    #> 3       1            3 0.00000000
    #> 4       1            4 0.00000000
    #> 5       1            5 0.00000000
    #> 6       1            6 1.00000000
    #> 7       1            7 4.50823457
    #> 8       1            8 9.00497739
    #> 9       1            9 0.62246947
    #> 10      1           10 0.06189564
    #> 11      1           11 0.00000000
    #> 12      2            1 0.00000000
    #> 13      2            2 0.00000000
    #> 14      2            3 0.00000000
    #> 15      2            4 0.00000000
    #> 16      2            5 0.00000000
    #> 17      2            6 1.00000000
    #> 18      2            7 3.00926000
    #> 19      2            8 3.88439088
    #> 20      2            9 0.88071161
    #> 21      2           10 1.12439876
    #> 22      2           11 0.02567988
    #> 23      3            1 0.00000000
    #> 24      3            2 0.00000000
    #> 25      3            3 0.00000000
    #> 26      3            4 0.00000000
    #> 27      3            5 0.00000000
    #> 28      3            6 0.00000000
    #> 29      3            7 0.00000000
    #> 30      3            8 0.00000000
    #> 31      3            9 0.00000000
    #> 32      3           10 0.00000000
    #> 33      3           11 1.00000000
    #> 34      4            1 0.00000000
    #> 35      4            2 0.00000000
    #> 36      4            3 0.00000000
    #> 37      4            4 0.00000000
    #> 38      4            5 0.00000000
    #> 39      4            6 1.00000000
    #> 40      4            7 6.64583147
    #> 41      4            8 4.53914861
    #> 42      4            9 1.35352344
    #> 43      4           10 0.77936516
    #> 44      4           11 0.01444315

The reproductive values of juveniles stages have essentially dropped to 0, leaving only the adult stages, starting from the dormant stage. So, indeed, we see dramatically different reproductive values here.

Now on to function-based matrices!

### Analysis 2. Function-based matrix estimation

### Step 1. Data characterization and reorganization

Before we go further, we need to describe the life history characterizing dataset properly with a `stageframe` for our Cypripedium candidum dataset, which will be different from the one that we created for the raw matrix example. Since this analysis will be function-based, we will include all possible size classes here. If constructing raw matrices, as in the previous example, all sizes that occur in the dataset need to be accounted for in a way that is both natural and parsimonious with respect to transition estimation. If constructing function-based matrices, such as IPMs, then representative sizes at systematic increments will be satisfactory. Since size is count-based in the Cypripedium candidum case, we will use all numbers of stems that might occur from 0 to the maximum in the dataset, representing the life history diagram shown in the beginning of this chapter.

    sizevector <- c(0, 0, 0, 0, 0, seq(from = 0, t = 24), seq(from = 1, to = 24))
    stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "V1", "V2", "V3", "V4", "V5",
                     "V6", "V7", "V8", "V9", "V10", "V11", "V12", "V13", "V14", "V15",
                     "V16", "V17", "V18", "V19", "V20", "V21", "V22", "V23", "V24",
                     "F1", "F2", "F3", "F4", "F5", "F6", "F7", "F8", "F9", "F10",
                     "F11", "F12", "F13", "F14", "F15", "F16", "F17", "F18", "F19",
                     "F20", "F21", "F22", "F23", "F24")
    repvector <- c(0, 0, 0, 0, 0, rep(0, 25), rep(1, 24))
    obsvector <- c(0, 0, 0, 0, 0, 0, rep(1, 48))
    matvector <- c(0, 0, 0, 0, 0, rep(1, 49))
    immvector <- c(0, 1, 1, 1, 1, rep(0, 49))
    propvector <- c(1, rep(0, 53))
    indataset <- c(0, 0, 0, 0, 0, rep(1, 49))

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

A close look at the output object, `cypframe`, shows a data frame that includes in order for each stage: the stage’s name, the associated size, its reproductive status, its status as an observable stage, its status as a propagule stage, its status as an immature stage, its status as a mature stage, whether it occurs in the dataset, the half-width of a size class bin, the minima and maxima of size class bins, the centroid of the size class bin, the full size class bin width, and comments. Stage names and combinations of characteristics must be unique to prevent estimation errors, and the comments field may be edited to include any information deemed pertinent. We may edit the comments field as below.

    cypframe$comments[(cypframe$stagenames == "SD")] <- "Dormant seed"
    cypframe$comments[(cypframe$stagenames == "P1")] <- "1st yr protocorm"
    cypframe$comments[(cypframe$stagenames == "P2")] <- "2nd yr protocorm"
    cypframe$comments[(cypframe$stagenames == "P3")] <- "3rd yr protocorm"
    cypframe$comments[(cypframe$stagenames == "SL")] <- "Seedling"
    cypframe$comments[(cypframe$stagenames == "D")] <- "Dormant mature"
    cypframe$comments[(cypframe$stagenames == "V1")] <- "Non-reproductive mature individual with 1 stem"
    cypframe$comments[(cypframe$stagenames == "V2")] <- "Non-reproductive mature individual with 2 stems"
    cypframe$comments[(cypframe$stagenames == "V3")] <- "Non-reproductive mature individual with 3 stems"
    cypframe$comments[(cypframe$stagenames == "V4")] <- "Non-reproductive mature individual with 4 stems"
    cypframe$comments[(cypframe$stagenames == "V5")] <- "Non-reproductive mature individual with 5 stems"
    cypframe$comments[(cypframe$stagenames == "V6")] <- "Non-reproductive mature individual with 6 stems"
    cypframe$comments[(cypframe$stagenames == "V7")] <- "Non-reproductive mature individual with 7 stems"
    cypframe$comments[(cypframe$stagenames == "V8")] <- "Non-reproductive mature individual with 8 stems"
    cypframe$comments[(cypframe$stagenames == "V9")] <- "Non-reproductive mature individual with 9 stems"
    cypframe$comments[(cypframe$stagenames == "V10")] <- "Non-reproductive mature individual with 10 stems"
    cypframe$comments[(cypframe$stagenames == "V11")] <- "Non-reproductive mature individual with 11 stems"
    cypframe$comments[(cypframe$stagenames == "V12")] <- "Non-reproductive mature individual with 12 stems"
    cypframe$comments[(cypframe$stagenames == "V13")] <- "Non-reproductive mature individual with 13 stems"
    cypframe$comments[(cypframe$stagenames == "V14")] <- "Non-reproductive mature individual with 14 stems"
    cypframe$comments[(cypframe$stagenames == "V15")] <- "Non-reproductive mature individual with 15 stems"
    cypframe$comments[(cypframe$stagenames == "V16")] <- "Non-reproductive mature individual with 16 stems"
    cypframe$comments[(cypframe$stagenames == "V17")] <- "Non-reproductive mature individual with 17 stems"
    cypframe$comments[(cypframe$stagenames == "V18")] <- "Non-reproductive mature individual with 18 stems"
    cypframe$comments[(cypframe$stagenames == "V19")] <- "Non-reproductive mature individual with 19 stems"
    cypframe$comments[(cypframe$stagenames == "V20")] <- "Non-reproductive mature individual with 20 stems"
    cypframe$comments[(cypframe$stagenames == "V21")] <- "Non-reproductive mature individual with 21 stems"
    cypframe$comments[(cypframe$stagenames == "V22")] <- "Non-reproductive mature individual with 22 stems"
    cypframe$comments[(cypframe$stagenames == "V23")] <- "Non-reproductive mature individual with 23 stems"
    cypframe$comments[(cypframe$stagenames == "V24")] <- "Non-reproductive mature individual with 24 stems"
    cypframe$comments[(cypframe$stagenames == "F1")] <- "Flowering mature individual with 1 stem"
    cypframe$comments[(cypframe$stagenames == "F2")] <- "Flowering mature individual with 2 stems"
    cypframe$comments[(cypframe$stagenames == "F3")] <- "Flowering mature individual with 3 stems"
    cypframe$comments[(cypframe$stagenames == "F4")] <- "Flowering mature individual with 4 stems"
    cypframe$comments[(cypframe$stagenames == "F5")] <- "Flowering mature individual with 5 stems"
    cypframe$comments[(cypframe$stagenames == "F6")] <- "Flowering mature individual with 6 stems"
    cypframe$comments[(cypframe$stagenames == "F7")] <- "Flowering mature individual with 7 stems"
    cypframe$comments[(cypframe$stagenames == "F8")] <- "Flowering mature individual with 8 stems"
    cypframe$comments[(cypframe$stagenames == "F9")] <- "Flowering mature individual with 9 stems"
    cypframe$comments[(cypframe$stagenames == "F10")] <- "Flowering mature individual with 10 stems"
    cypframe$comments[(cypframe$stagenames == "F11")] <- "Flowering mature individual with 11 stems"
    cypframe$comments[(cypframe$stagenames == "F12")] <- "Flowering mature individual with 12 stems"
    cypframe$comments[(cypframe$stagenames == "F13")] <- "Flowering mature individual with 13 stems"
    cypframe$comments[(cypframe$stagenames == "F14")] <- "Flowering mature individual with 14 stems"
    cypframe$comments[(cypframe$stagenames == "F15")] <- "Flowering mature individual with 15 stems"
    cypframe$comments[(cypframe$stagenames == "F16")] <- "Flowering mature individual with 16 stems"
    cypframe$comments[(cypframe$stagenames == "F17")] <- "Flowering mature individual with 17 stems"
    cypframe$comments[(cypframe$stagenames == "F18")] <- "Flowering mature individual with 18 stems"
    cypframe$comments[(cypframe$stagenames == "F19")] <- "Flowering mature individual with 19 stems"
    cypframe$comments[(cypframe$stagenames == "F20")] <- "Flowering mature individual with 20 stems"
    cypframe$comments[(cypframe$stagenames == "F21")] <- "Flowering mature individual with 21 stems"
    cypframe$comments[(cypframe$stagenames == "F22")] <- "Flowering mature individual with 22 stems"
    cypframe$comments[(cypframe$stagenames == "F23")] <- "Flowering mature individual with 23 stems"
    cypframe$comments[(cypframe$stagenames == "F24")] <- "Flowering mature individual with 24 stems"

Now we may transform our vertical dataset into a historically-formatted vertical file. The resulting dataset will have each individual’s observed life history broken up into states corresponding to three consecutive years per row, with plant identity marked in each row. To handle this, we use the `verticalize3()` function, as below.

    vertdata <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004,
                             patchidcol = "patch", individcol = "plantid",
                             blocksize = 4, size1col = "Inf2.04", size2col = "Inf.04",
                             size3col = "Veg.04", repstr1col = "Inf.04",
                             fec1col = "Pod.04", stageassign = cypframe,
                             stagesize = "sizeadded", NAas0 = TRUE)

    summary(vertdata)
    #>      rowid        popid         patchid    individ           year2          sizea1             sizea2             sizea3
    #>  Min.   : 1.00   Mode:logical   A: 93   Min.   : 164.0   Min.   :2004   Min.   :0.000000   Min.   :0.000000   Min.   :0.000000
    #>  1st Qu.:21.00   NA's:320       B:154   1st Qu.: 391.0   1st Qu.:2005   1st Qu.:0.000000   1st Qu.:0.000000   1st Qu.:0.000000
    #>  Median :37.50                  C: 73   Median : 453.0   Median :2006   Median :0.000000   Median :0.000000   Median :0.000000
    #>  Mean   :38.45                          Mean   : 651.5   Mean   :2006   Mean   :0.009375   Mean   :0.009375   Mean   :0.009375
    #>  3rd Qu.:56.00                          3rd Qu.: 476.0   3rd Qu.:2007   3rd Qu.:0.000000   3rd Qu.:0.000000   3rd Qu.:0.000000
    #>  Max.   :77.00                          Max.   :1560.0   Max.   :2008   Max.   :1.000000   Max.   :1.000000   Max.   :1.000000
    #>     repstra1          repstra2          repstra3          feca1            feca2            feca3            sizeb1
    #>  Min.   : 0.0000   Min.   : 0.0000   Min.   : 0.000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   : 0.0000
    #>  1st Qu.: 0.0000   1st Qu.: 0.0000   1st Qu.: 0.000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.: 0.0000
    #>  Median : 0.0000   Median : 0.0000   Median : 0.000   Median :0.0000   Median :0.0000   Median :0.0000   Median : 0.0000
    #>  Mean   : 0.7469   Mean   : 0.8969   Mean   : 1.069   Mean   :0.2656   Mean   :0.2906   Mean   :0.4562   Mean   : 0.7469
    #>  3rd Qu.: 1.0000   3rd Qu.: 1.0000   3rd Qu.: 1.000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.: 1.0000
    #>  Max.   :18.0000   Max.   :18.0000   Max.   :18.000   Max.   :7.0000   Max.   :7.0000   Max.   :8.0000   Max.   :18.0000
    #>      sizeb2            sizeb3           sizec1         sizec2           sizec3         size1added       size2added       size3added
    #>  Min.   : 0.0000   Min.   : 0.000   Min.   : 0.0   Min.   : 0.000   Min.   : 0.000   Min.   : 0.000   Min.   : 0.000   Min.   : 0.000
    #>  1st Qu.: 0.0000   1st Qu.: 0.000   1st Qu.: 0.0   1st Qu.: 1.000   1st Qu.: 1.000   1st Qu.: 0.000   1st Qu.: 1.000   1st Qu.: 1.000
    #>  Median : 0.0000   Median : 0.000   Median : 1.0   Median : 2.000   Median : 1.000   Median : 2.000   Median : 2.000   Median : 2.000
    #>  Mean   : 0.8969   Mean   : 1.069   Mean   : 1.9   Mean   : 2.416   Mean   : 2.209   Mean   : 2.656   Mean   : 3.322   Mean   : 3.288
    #>  3rd Qu.: 1.0000   3rd Qu.: 1.000   3rd Qu.: 3.0   3rd Qu.: 3.000   3rd Qu.: 3.000   3rd Qu.: 4.000   3rd Qu.: 4.000   3rd Qu.: 4.000
    #>  Max.   :18.0000   Max.   :18.000   Max.   :13.0   Max.   :13.000   Max.   :13.000   Max.   :21.000   Max.   :24.000   Max.   :24.000
    #>    obsstatus1       obsstatus2       obsstatus3    repstatus1       repstatus2       repstatus3    fecstatus1       fecstatus2
    #>  Min.   :0.0000   Min.   :0.0000   Min.   :0.0   Min.   :0.0000   Min.   :0.0000   Min.   :0.0   Min.   :0.0000   Min.   :0.0000
    #>  1st Qu.:0.0000   1st Qu.:1.0000   1st Qu.:1.0   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0   1st Qu.:0.0000   1st Qu.:0.0000
    #>  Median :1.0000   Median :1.0000   Median :1.0   Median :0.0000   Median :0.0000   Median :0.0   Median :0.0000   Median :0.0000
    #>  Mean   :0.7469   Mean   :0.9531   Mean   :0.9   Mean   :0.2875   Mean   :0.3688   Mean   :0.4   Mean   :0.1344   Mean   :0.1562
    #>  3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0   3rd Qu.:0.0000   3rd Qu.:0.0000
    #>  Max.   :1.0000   Max.   :1.0000   Max.   :1.0   Max.   :1.0000   Max.   :1.0000   Max.   :1.0   Max.   :1.0000   Max.   :1.0000
    #>    fecstatus3       firstseen       lastseen        alive1           alive2      alive3           obsage       obslifespan
    #>  Min.   :0.0000   Min.   :2004   Min.   :2004   Min.   :0.0000   Min.   :1   Min.   :0.0000   Min.   :0.000   Min.   :0.000
    #>  1st Qu.:0.0000   1st Qu.:2004   1st Qu.:2009   1st Qu.:1.0000   1st Qu.:1   1st Qu.:1.0000   1st Qu.:1.000   1st Qu.:5.000
    #>  Median :0.0000   Median :2004   Median :2009   Median :1.0000   Median :1   Median :1.0000   Median :2.000   Median :5.000
    #>  Mean   :0.2219   Mean   :2004   Mean   :2009   Mean   :0.7688   Mean   :1   Mean   :0.9469   Mean   :1.853   Mean   :4.556
    #>  3rd Qu.:0.0000   3rd Qu.:2004   3rd Qu.:2009   3rd Qu.:1.0000   3rd Qu.:1   3rd Qu.:1.0000   3rd Qu.:3.000   3rd Qu.:5.000
    #>  Max.   :1.0000   Max.   :2008   Max.   :2009   Max.   :1.0000   Max.   :1   Max.   :1.0000   Max.   :4.000   Max.   :5.000
    #>     stage1             stage2             stage3            matstatus1   matstatus2   matstatus3
    #>  Length:320         Length:320         Length:320         Min.   :1    Min.   :1    Min.   :1
    #>  Class :character   Class :character   Class :character   1st Qu.:1    1st Qu.:1    1st Qu.:1
    #>  Mode  :character   Mode  :character   Mode  :character   Median :1    Median :1    Median :1
    #>                                                           Mean   :1    Mean   :1    Mean   :1
    #>                                                           3rd Qu.:1    3rd Qu.:1    3rd Qu.:1
    #>                                                           Max.   :1    Max.   :1    Max.   :1
    dim(vertdata)
    #> [1] 320  45

In the above code, we described the input dataset in a way that allows R to reorganize it appropriately. For this to work properly, the input dataset needs to be arranged in blocks of columns for each year, with variables in the same order every year. The output dataset includes a number of summary variables, but the data is essentially broken down into groups of three censecutive monitoring occasions each (time t+1, t, and t-1, corresponding to `year3`, `year2`, and `year1` in the output, respectively), with individuals spread across multiple rows. The output dataset is further limited to those entries in which the individual is alive in time t (`year2`), meaning that all rows in which an individual is dead or not yet recruited in time t are dropped. Thus, we have 320 rows of data and 45 variables.

This reorganized dataset includes a set of interesting terms, the `sizeadded` group of three variables. These are sums of the size variables for each time. So, `size1added` is calculated as `sizea1 + sizeb1 + sizec1`. This may or may not be sensical depending on the dataset. In this particular dataset, the full size of the individual in each time step is this sum, because `Veg` gives the number of non-reproductive stems, `Inf` gives the number of single-flowered inflorescences, and `Inf2` gives the number of double-flowered inflorescences per plant per time-step (an inflorescence takes a single stem, and no inflorescence has more than two flowers). So, we will use the `sizeadded` group of variables to code individual size in our analyses.

### Step 2. Provide supplemental information for matrix estimation

The next steps involve the creation of a reproduction matrix and an `overwrite` data frame. These are optional, and only to be set if the life history of the organism calls for it. The reproduction matrix tells where where fecundity rates need to be set, and at what level. Cypripedium candidum can create seed that germinates by the following growing season (stage P1, or a first year protocorm), or that remains dormant for the next year (stage SD). In the following matrix, we detail that the fecundity of each reproductive stage needs to be split into two between each of these output stages. The actual split places 50% of the fecundity of a stage into each category of recruit, where the full fecundity is estimated by the linear models that we have already created.

    rep.assumptions <- matrix(0, 54, 54)
    rep.assumptions[1:2,31:54] <- 0.5

Next we create a data frame that outlines transitions that cannot be estimated from the data set and need to be set by other means. For this task, we will use the `overwrite` function. The function will handle two kinds of given transitions:

1. transitions that will be set to specific probabilities or rates that we specify, and

2. transitions that will be set to the values of other transitions that are to be estimated.

Here is an example for the Cypripedium candidum analysis. Note that each row refers to a specific transition, and there are codes for 13 given transitions. Most of these transitions are set to specific probabilities, but four are transitions that will be set to other, estimated transitions (these are the non-NA transitions in `eststage` parameter set below).

    cypover <- overwrite(stage3 = c("SD", "SD", "P1", "P1", "P2", "P3", "D", "V1",
                         "V2", "V3", "SL", "SL", "SL"), stage2 = c("SD", "SD", "SD",
                         "SD", "P1", "P2", "P3", "P3", "P3", "P3", "P3", "SL", "SL"),
                         stage1 = c("SD", "rep", "SD", "rep", "SD", "P1", "P2", "P2",
                         "P2", "P2", "P2", "P3", "SL"), eststage3 = c(NA, NA, NA, NA,
                         NA, NA, "D", "V1", "V2", "V3", NA, NA, NA), eststage2 = c(NA,
                         NA, NA, NA, NA, NA, "D", "D", "D", "D", NA, NA, NA),
                         eststage1 = c(NA, NA, NA, NA, NA, NA, "D", "D", "D", "D", NA,
                         NA, NA), givenrate = c(0.1, 0.1, 0.2, 0.2, 0.2, 0.2, NA, NA,
                         NA, NA, 0.25, 0.4, 0.4), type = c("S", "S", "S", "S", "S",
                         "S", "S", "S", "S", "S", "S", "S", "S"))
    cypover
    #>    stage3 stage2 stage1 eststage3 eststage2 eststage1 givenrate convtype
    #> 1      SD     SD     SD      <NA>      <NA>      <NA>      0.10        1
    #> 2      SD     SD    rep      <NA>      <NA>      <NA>      0.10        1
    #> 3      P1     SD     SD      <NA>      <NA>      <NA>      0.20        1
    #> 4      P1     SD    rep      <NA>      <NA>      <NA>      0.20        1
    #> 5      P2     P1     SD      <NA>      <NA>      <NA>      0.20        1
    #> 6      P3     P2     P1      <NA>      <NA>      <NA>      0.20        1
    #> 7       D     P3     P2         D         D         D        NA        1
    #> 8      V1     P3     P2        V1         D         D        NA        1
    #> 9      V2     P3     P2        V2         D         D        NA        1
    #> 10     V3     P3     P2        V3         D         D        NA        1
    #> 11     SL     P3     P2      <NA>      <NA>      <NA>      0.25        1
    #> 12     SL     SL     P3      <NA>      <NA>      <NA>      0.40        1
    #> 13     SL     SL     SL      <NA>      <NA>      <NA>      0.40        1

We can now proceed with matrix estimation!

### Step 3. Develop models of demographic parameters

Matrix creation can proceed either as raw matrix creation, as initially outlined in Ehrlén (2000), or via the creation of function-based matrices, in many ways equivalent to complex integral projection models per Ellner and Rees (2006) and as further described in the non-Gaussian case in Shefferson et al. (2014). In the raw matrix case, no vital rate estimation is conducted prior to the construction of a projection matrix. In function-based matrix case, vital rates are first analyzed via linear or non-linear models of the raw demographic dataset, and functions are created that estimate these vital rates according the inputs given. Matrices are then estimated using these functions, as opposed to the raw data.

Package `lefko3` estimates up to nine different sets of vital rates:

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

2. Observation probability, which is an optional parameter denoting the probability of observation in time step t+1 of an individual in stage k given survival from time t to time t+1. This parameter is only used when at least one stage is technically unobservable. For example, some plants are capable of vegetative dormancy, in which case they are alive but do not necessarily sprout in all years. In these cases, the probability of sprouting can be estimated as the observation probability. Note that this probability does NOT refer to observer effort, and so should ONLY be used when some stages are completely unobservable or when size is a count variable and a size of 0 is possible. In `lefko3`, this parameter may be modeled as a function of size, reproductive status, patch, year, and individual identity.

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

4. Reproduction probability, which is an optional parameter denoting the probability of reproducing in time step t+1 given survival from time t to time t+1, and observation in that time. Note that this should be used to separate breeding from non-breeding mature stages. If all adult stages are potentially reproductive, then this parameter is not needed. In `lefko3`, this parameter may be modeled as a function of size, reproductive status, patch, year, and individual identity.

5. Fecundity rate, which is the rate of successful production of offspring in time t by individuals alive, observable, and reproductive in that time, and the survival of those offspring into time t+1 in whatever juvenile class is possible. Thus, fecundity rate would be in terms of seedlings produced if seeds are not capable of dormancy. In `lefko3`, this parameter may be modeled as a function of size, reproductive status, patch, year, and individual identity. This parameter is required in all function-based matrices.

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

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

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

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

All models require data, and so analyses of datasets that do not include the fates of juveniles should not attempt to estimate juvenile-associated rates and probabilities. Please see chapter 1 for further information.

Prior to vital rate estimation, a number of key decisions need to be made regarding the assumptions underlying the vital rates, and their relationships with the factors under investigation. These decisions include the general modeling strategy, and the size and fecundity distributions.

#### Step 3a. General modeling strategy

Most function-based matrices, whether integral projection models or otherwise, use either a generalized linear modeling (GLM) strategy, or a generalized linear mixed modeling (GLMM) strategy. The former is more common because of its simplicity, but the latter is theoretically more correct because it allows a correction for pseudoreplication of individuals and a more realistic approach handling time. The difference between the two with regards to vital rate modeling is strongly related to assumptions regarding the individual and the nature of temporal variation in vital rates.

In both cases, the underlying dataset utilized is a vertical dataset. Here, each row of data gives the state of the individual in either two consecutive time steps (the ahistorical case), or three consecutive time steps (the historical case). Under a GLM framework, time step is a fixed categorical variable, and individual identity is ignored. Using time step as a fixed categorical variable implies that the monitoring occasions worked with are the only time steps for which inference is wanted. Thus, it would not be correct to infer future population dynamics after 2009 for a dataset collected between 2004 and 2009, if year is treated as fixed. Ignoring individual identity treats all data points as independent, even though data originating from the same sampled individual is clearly not independent. This leads to pseudoreplication, which is observable in demographic modeling as higher significance of terms, and more significant terms, in modeling.

Under a GLMM framework, both time step and individual identity can be treated as random categorical terms. This has two major implications. First, both time step and individuals can be assumed to be random samples from a broader population of time steps and individuals for which we want to make inferences. Thus, sampled years represent a greater universe of years for which inference can be made, and so their associated coefficients are assumed to come from a normal distribution with mean = 0. Second, treating individual as a random categorical term eliminates the pseudoreplication that is inherent in the GLM approach to vital rate estimation when individuals are monitored potentially many times. We encourage researchers to use the GLMM approach in their demographic work, but we have also included easy-to-use GLM functionality, since many will find the GLM approach particularly useful in cases where mixed modeling breaks down.

#### Step 3b. Size and fecundity distributions

Once a general approach is decided upon, the next step is to decide the underlying distributions. The probabilities of survival, observation, and reproductive status are automatically set to the binomial distribution, and this cannot be altered. However, the probability of size transition and fecundity rate can be set to the Gaussian, Poisson, or negative binomial distributions. In general, if size or fecundity rate is a continuous variable (i.e., not an integer or count variable), then it should be set to the Gaussian distribution. In contrast, if size or fecundity rate is a count, then it should be set to the Poisson distribution. The negative binomial distribution is also provided in cases where the Poisson distribution’s assumptions, such as the mean equaling the variance, are clearly broken. We do not encourage the use of the negative binomial except in such extreme cases, as the extra parameters estimated for the negative binomial distribution reduce the power of the modeling exercises conducted.

Neither the Poisson nor the negative binomial distribution handle 0s. We do not include zero-inflated versions of these distributions because such distributions are technically incorrect, or theoretically sloppy. If 0s occur within the dataset, then the best option is to utilize the probability of observation or the probability of reproduction as vital rates that will also be modeled. This will allow the 0s to be analyzed within these binomial models, and leave the non-zero values for rate estimation under the Poisson or negative binomial distributions. If 0s still exist after doing this, then it is likely that some individuals assigned to visible classes are likely not observable, size class might be determined logarithmically, or reproductive classes are not really reproductive. Further exploration of the dataset should yield the correct answer.

#### Step 3c. Model building and selection

In lefko3, the `modelsearch` function is the main workhorse function that conducts vital rate estimation. Here, we will create a full suite of vital rate models for the Cypripedium candidum dataset. Before proceeding, we need to decide on the linear model building strategy, the correct vital rates to model, the proper statistical distributions for estimated vital rates, the proper parameterizations for each vital rate, and the strategy for determination of the best-fit models.

First, the building strategy. In most cases, the best procedure will be through linear mixed models in which monitoring occasion and individual identity are random terms. We set monitoring occasion as random because we wish to make inferences for the population as a whole and not restrict ourselves to inference only for the years monitored (i.e. our distribution of years sampled is itself a sample of the population in time). We set individual identity to random because we do not wish to pseudoreplicate, since not doing so will result in the assumption that each row in our vertical dataset is statistically independent of all other rows. Thus, we will use `approach = "lme4"` in the parameterization for `modelsearch`, and keep the defaults for `year.as.random`, `indiv`, and `year`, which are set to the default output for whether monitoring occasion is a random or fixed term (random by default), which variable corresponds to individual identity (`individ` by default), and which variable corresponds to time t(`year2` by default).

The mixed modeling approach is usually better, particularly because of its elimination of pseudoreplication. However, a mixed modeling strategy results in lower statistical power and a greater time used in estimating models. `lefko3` users wishing to use a standard generalized linear modeling strategy can set `approach = "glm"`. In this case, individual identity is not used and all observed transitions are treated as independent.

Second, the correct vital rates to model. The function `modelsearch()` estimates up to 5 sets of linear models:

1. survival probability from time t to time t+1 (juvenile and mature),

2. observation probability in time t+1 assuming survival until that time (juvenile and mature),

3. size in time t+1 assuming survival and observation in that time (juvenile and mature),

4. reproduction status in time t+1 assuming survival and observation until that time (individuals leaving the juvenile stages, and those already mature), and

5. fecundity rate assuming survival in time t and observation and reproduction in time t+1 (mature only).

The default settings for `modelsearch` involve the estimation of 1), survival probability, 3) size distribution, and 5) fecundity, which are the minimum required for a full projection matrix. Observation probability (option `obs` in `vitalrates`) should only be included when a life history stage or size exists that cannot be observed. For example, in the case of a plant with vegetative dormancy, the observation probability can be thought of as the sprouting probability, which is a biologically meaningful vital rate (Shefferson et al. 2001). Further, reproduction status (option `repst` in `vitalrates`) should only be modeled if size classification should be stratified by the ability to reproduce, as when 0 fecundity occurs within the dataset. In this latter case, we can imagine that reproductive and non-reproductive individuals of each size class might theoretically exist, and we wish to parameterize transitions allowing individuals to be reproductive or non-reproductive. Since Cypripedium candidum is capable of longs bouts of vegetative dormancy, and since we wish to stratify the population into reproductive and non-reproductive adults, we will set `vitalrates = c("surv", "obs", "size", "repst", "fec")`.

Third, we need to set the proper statistical distributions for each parameter. Survival probability, observation probability, and reproductive status (a probability) are all modeled as binomial variables, and this cannot be changed. In the case of Cypripedium, size was measured as the number of stems and so is a count variable. Likewise, fecundity is actually estimated as the number of fruits produced per plant, and so is also a count variable. Thus, we will use the Poisson distribution for both variables.

Users of `lefko3` who wish to use a Poisson or negative binomial distribution should note that 0s are not allowed, because zero-inflated models are essentially conditional models with different statistical properties at 0 than elsewhere. The `modelsearch` function will subset the dataset used to parameterize size and fecundity to only rows without zeroes in time t+1. If size is potentially 0, then consider using observation probability to catch the probability of becoming size 0 vs. observable sizes.

Fourth, we need the proper model parameterizations for each vital rate, using the `suite` option. The default, `suite = "main"`, under the mixed model setting (`approach = "lme4"`) is for `modelsearch` to estimate a global model that includes sizes in times t and t-1, and reproductive status in times t and t-1, as fixed factors, with individual identity and time t set as random terms in a mixed model framework using R package `lme4` (Bates et al. 2015). However, setting `suite = "full"` will yield a global model that also includes all two-way interactions. We will set to the latter. The default under the GLM setting (`approach = "glm"`) makes time t a fixed term and drops individual identity from consideration. The global model under `suite = "full"` then includes all fixed factors noted before, plus time t and all two-way interactions with it. To eliminate all interactions from the model and only analyze main effects, use `suite = "main"`. If the population is not stratified by reproductive status, then `suite = "size"` will eliminate reproductive status terms and use all others in the global model. If size is not important, then `suite = "rep"` will eiminate size but keep reproductive status and all other terms. Finally, `suite = "cons"` will result in a global model in which both reproductive status and size are not considered.

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

In the model building exercise below, we will use the `suite = "full"` option to run all main effects and their two-way interactions.

    cypmodels3 <- modelsearch(vertdata, historical = TRUE, approach = "lme4",
                              vitalrates = c("surv", "obs", "size", "repst", "fec"),
                              sizedist = "poisson", fecdist = "poisson", suite = "full",
                              size = c("size3added", "size2added", "size1added"),
                              quiet = TRUE)
    #> Warning in modelsearch(vertdata, historical = TRUE, approach = "lme4", vitalrates = c("surv", : WARNING: Fecundity in time t cannot be Poisson-distributed and include 0s. Will develop fecundity models excluding all 0s. Consider adding a reproductive status variable to absorb 0 values.
    #>
    #> Developing global model of survival probability...
    #> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : unable to evaluate scaled gradient
    #> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge: degenerate Hessian with 2
    #> negative eigenvalues
    #>
    #> Developing global model of observation probability...
    #> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0217509 (tol =
    #> 0.002, component 1)
    #>
    #> Developing global model of size (Poisson)...
    #> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0181958 (tol =
    #> 0.002, component 1)
    #> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
    #>  - Rescale variables?
    #>
    #> Developing global model of the probability of reproduction...
    #> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00594956 (tol =
    #> 0.002, component 1)
    #>
    #> Developing global model of fecundity (Poisson)...
    #> fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients
    #> boundary (singular) fit: see ?isSingular
    #>
    #> All global models developed.
    #>
    #>
    #> Commencing dredge of survival probability...
    #>
    #>
    #> Commencing dredge of observation probability...
    #>
    #>
    #> Commencing dredge of size...
    #>
    #>
    #> Commencing dredge of reproduction probability...
    #>
    #>
    #> Commencing dredge of fecundity...
    #>
    #>
    #> Finished dredging all vital rates.

As `modelsearch` works, it produces a good amount of text to allow the user to understand what is going on. It is entirely possible, and actually quite likely, that many warning messages will appear, and these may be of use to users in understanding their data and how well it conforms to their analytical assumptions. We have silenced this with the `quiet = TRUE` option, but we encourage users to allow the function to run unsilenced, in case of large problems in modeling.

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

    summary(cypmodels3)
    #> This LefkoMod object includes 5 linear models.
    #> Best-fit model criterion used: AICc&k
    #>
    #> ────────────────────────────────────────────────────────────────────────────────
    #> Survival model:
    #> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
    #>  Family: binomial  ( logit )
    #> Formula: alive3 ~ size2added + (1 | individ) + (1 | year2)
    #>    Data: surv.data
    #>      AIC      BIC   logLik deviance df.resid
    #> 128.1324 143.2057 -60.0662 120.1324      316
    #> Random effects:
    #>  Groups  Name        Std.Dev.
    #>  individ (Intercept) 1.198378
    #>  year2   (Intercept) 0.008826
    #> Number of obs: 320, groups:  individ, 74; year2, 5
    #> Fixed Effects:
    #> (Intercept)   size2added
    #>      2.0352       0.6344
    #> convergence code 0; 0 optimizer warnings; 1 lme4 warnings
    #>
    #> ────────────────────────────────────────
    #>
    #> Observation model:
    #> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
    #>  Family: binomial  ( logit )
    #> Formula: obsstatus3 ~ size2added + (1 | individ) + (1 | year2)
    #>    Data: obs.data
    #>      AIC      BIC   logLik deviance df.resid
    #> 118.2567 133.1117 -55.1284 110.2567      299
    #> Random effects:
    #>  Groups  Name        Std.Dev.
    #>  individ (Intercept) 1.078e-05
    #>  year2   (Intercept) 8.776e-01
    #> Number of obs: 303, groups:  individ, 70; year2, 5
    #> Fixed Effects:
    #> (Intercept)   size2added
    #>      2.4904       0.3134
    #> convergence code 0; 0 optimizer warnings; 1 lme4 warnings
    #>
    #> ────────────────────────────────────────
    #>
    #> Size model:
    #> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
    #>  Family: poisson  ( log )
    #> Formula: size3added ~ repstatus1 + repstatus2 + size1added + size2added +
    #>     (1 | individ) + (1 | year2) + repstatus1:repstatus2 + size1added:size2added
    #>    Data: size.data
    #>       AIC       BIC    logLik  deviance  df.resid
    #> 1081.6019 1114.5685 -531.8009 1063.6019       279
    #> Random effects:
    #>  Groups  Name        Std.Dev.
    #>  individ (Intercept) 0.06104
    #>  year2   (Intercept) 0.19627
    #> Number of obs: 288, groups:  individ, 70; year2, 5
    #> Fixed Effects:
    #>           (Intercept)             repstatus1             repstatus2             size1added             size2added  repstatus1:repstatus2
    #>              0.213173               0.278913               0.255861               0.131367               0.153437              -0.368510
    #> size1added:size2added
    #>             -0.008696
    #> convergence code 0; 0 optimizer warnings; 2 lme4 warnings
    #>
    #> ────────────────────────────────────────
    #>
    #> Reproductive status model:
    #> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
    #>  Family: binomial  ( logit )
    #> Formula: repstatus3 ~ repstatus2 + size2added + (1 | individ) + (1 | year2)
    #>    Data: repst.data
    #>       AIC       BIC    logLik  deviance  df.resid
    #>  333.6176  351.9324 -161.8088  323.6176       283
    #> Random effects:
    #>  Groups  Name        Std.Dev.
    #>  individ (Intercept) 0.1829
    #>  year2   (Intercept) 0.6250
    #> Number of obs: 288, groups:  individ, 70; year2, 5
    #> Fixed Effects:
    #> (Intercept)   repstatus2   size2added
    #>     -1.4630       1.6457       0.1715
    #>
    #> ────────────────────────────────────────
    #>
    #> Fecundity model:
    #> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
    #>  Family: poisson  ( log )
    #> Formula: feca2 ~ size2added + (1 | individ) + (1 | year2)
    #>    Data: fec.data
    #>      AIC      BIC   logLik deviance df.resid
    #> 156.8531 164.5012 -74.4266 148.8531       46
    #> Random effects:
    #>  Groups  Name        Std.Dev.
    #>  individ (Intercept) 0.000
    #>  year2   (Intercept) 0.172
    #> Number of obs: 50, groups:  individ, 26; year2, 5
    #> Fixed Effects:
    #> (Intercept)   size2added
    #>     0.21931      0.04429
    #> convergence code 0; 0 optimizer warnings; 1 lme4 warnings
    #>
    #> ────────────────────────────────────────────────────────────────────────────────
    #> Juvenile survival model:
    #> [1] 1
    #>
    #> ────────────────────────────────────────
    #>
    #> Juvenile observation model:
    #> [1] 1
    #>
    #> ────────────────────────────────────────
    #>
    #> Juvenile size model:
    #> [1] 1
    #>
    #> ────────────────────────────────────────
    #>
    #> Juvenile reproduction model:
    #> [1] 1
    #>
    #>
    #> ────────────────────────────────────────────────────────────────────────────────
    #>
    #> Number of models in survival table:113
    #>
    #> Number of models in observation table:113
    #>
    #> Number of models in size table:113
    #>
    #> Number of models in reproduction status table:113
    #>
    #> Number of models in fecundity table:113
    #>
    #> 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 reproduction table: 1
    #>
    #>
    #> ────────────────────────────────────────────────────────────────────────────────
    #>
    #> General model parameter names (column 1), and specific names used in these models (column 2):
    #>    mainparams modelparams
    #> 1       year2       year2
    #> 2     individ     individ
    #> 3       patch        <NA>
    #> 4       surv3      alive3
    #> 5        obs3  obsstatus3
    #> 6       size3  size3added
    #> 7      repst3  repstatus3
    #> 8        fec3       feca3
    #> 9        fec2       feca2
    #> 10      size2  size2added
    #> 11      size1  size1added
    #> 12     repst2  repstatus2
    #> 13     repst1  repstatus1
    #> 14        age        <NA>
    #>
    #>
    #> ────────────────────────────────────────────────────────────────────────────────
    #>
    #> Quality control:
    #>
    #> Survival estimated with 74 individuals and 320 individual transitions.
    #> Observation estimated with 70 individuals and 303 individual transitions.
    #> Size estimated with 70 individuals and 288 individual transitions.
    #> Reproductive status estimated with 70 individuals and 288 individual transitions.
    #> Fecundity estimated with 26 individuals and 50 individual transitions.
    #> Juvenile survival not estimated.
    #> Juvenile observation probability not estimated.
    #> Juvenile size transition not estimated.
    #> Juvenile reproduction probability not estimated.
    #> NULL

The `summary` function gives us a great deal of information about all of the models, but also hides a number of more technical details. For example, the complete model selection tables are actually included in the `modelsearch` output, and users can see these by calling them directly, as below.

    cypmodels3$survival_table
    #> Global model call: lme4::glmer(formula = alive3 ~ size1added + size2added + repstatus1 +
    #>     repstatus2 + size1added:size2added + repstatus1:repstatus2 +
    #>     size1added:repstatus1 + size2added:repstatus2 + size1added:repstatus2 +
    #>     size2added:repstatus1 + (1 | individ) + (1 | year2), data = surv.data,
    #>     family = "binomial")
    #> ---
    #> Model selection table
    #>      (Int)      rp1      rp2     sz1    sz2 rp1:rp2 rp1:sz1  rp1:sz2 rp2:sz1    rp2:sz2    sz1:sz2 df  logLik used.criterion delta weight
    #> 9    2.035                           0.6344                                                         4 -60.066          128.3  0.00  0.134
    #> 11   2.006           0.54530         0.5899                                                         5 -59.773          129.7  1.48  0.064
    #> 13   1.538                   0.18620 0.5185                                                         5 -59.781          129.8  1.49  0.063
    #> 10   1.984  0.31570                  0.6240                                                         5 -59.967          130.1  1.87  0.053
    #> 28   2.063 -0.42760 -0.13110         0.5712   46.86                                                 7 -58.039          130.4  2.18  0.045
    #> 74   2.149 -0.77080                  0.5004                  0.81440                                6 -59.407          131.1  2.82  0.033
    #> 15   1.505           0.41530 0.18280 0.4824                                                         6 -59.578          131.4  3.17  0.027
    #> 46   2.042 -1.70300          0.03281 0.5542          1.0750                                         7 -58.546          131.5  3.19  0.027
    #> 12   1.972  0.25460  0.51650         0.5866                                                         6 -59.712          131.7  3.43  0.024
    #> 525  1.453                   0.23930 0.5678                                             -0.0242500  6 -59.745          131.8  3.50  0.023
    #> 267  1.994           0.78800         0.6309                                  -1.287e-01             6 -59.748          131.8  3.50  0.023
    #> 14   1.535  0.04577          0.18200 0.5183                                                         6 -59.779          131.8  3.57  0.022
    #> 64   2.193 -2.84700 -0.13040 0.03241 0.5207  177.70  1.2890                                         9 -56.623          131.8  3.57  0.022
    #> 92   2.161 -0.97690 -0.09133         0.5073   66.09          0.43230                                8 -57.864          132.2  3.93  0.019
    #> 32   1.642 -0.62310 -0.13560 0.16720 0.4722   32.09                                                 8 -57.910          132.3  4.02  0.018
    #> 284  1.999 -0.42880  0.16980         0.6194   78.38                          -1.665e-01             8 -57.991          132.4  4.18  0.017
    #> 160  1.446 -0.80970  0.49640 0.35650 0.4784   29.16                  -0.5133                        9 -57.022          132.6  4.37  0.015
    #> 76   2.126 -0.74500  0.43350         0.4772                  0.75490                                7 -59.228          132.8  4.56  0.014
    #> 143  1.381           0.81620 0.27590 0.4952                          -0.2702                        7 -59.333          133.0  4.76  0.012
    #> 78   1.696 -0.78060          0.15700 0.4225                  0.66390                                7 -59.333          133.0  4.77  0.012
    #> 48   2.022 -1.76700  0.47620 0.03453 0.5164          1.0790                                         8 -58.336          133.1  4.88  0.012
    #> 110  2.112 -2.41400          0.04005 0.5037          1.0360  0.64910                                8 -58.346          133.2  4.89  0.012
    #> 192  1.900 -2.52000  0.32410 0.18860 0.5129   50.38  1.0710          -0.3320                       10 -56.289          133.3  5.03  0.011
    #> 527  1.425           0.41160 0.23320 0.5288                                             -0.0227300  7 -59.546          133.5  5.19  0.010
    #> 271  1.467           0.58300 0.18440 0.5106                                  -9.563e-02             7 -59.561          133.5  5.22  0.010
    #> 16   1.507 -0.05451  0.42620 0.18760 0.4815                                                         7 -59.575          133.5  5.25  0.010
    #> 558  1.937 -1.68200          0.08943 0.6026          1.0560                             -0.0230700  8 -58.524          133.5  5.25  0.010
    #> 268  1.950  0.27090  0.78600         0.6328                                  -1.436e-01             7 -59.680          133.7  5.46  0.009
    #> 526  1.451  0.04333          0.23500 0.5673                                             -0.0240900  7 -59.743          133.8  5.59  0.008
    #> 96   1.988 -1.04200 -0.09154 0.13020 0.4754   67.44          0.37210                                9 -57.649          133.9  5.62  0.008
    #> 576  2.101 -2.79200 -0.13630 0.07267 0.5484   26.46  1.2560                             -0.0149300 10 -56.614          133.9  5.68  0.008
    #> 128  2.161 -2.75000 -0.13640 0.03317 0.5258   42.13  1.3190 -0.12630                               10 -56.618          133.9  5.69  0.008
    #> 320  2.156 -2.81200 -0.07952 0.03480 0.5263   16.11  1.2750                  -2.692e-02            10 -56.623          134.0  5.70  0.008
    #> 544  1.524 -0.63140 -0.14680 0.24430 0.5409   32.12                                     -0.0345200  9 -57.852          134.3  6.02  0.007
    #> 348  2.123 -0.93650  0.04906         0.5338   48.18          0.40450         -7.679e-02             9 -57.855          134.3  6.03  0.007
    #> 288  1.596 -0.62030  0.08667 0.16430 0.5111  102.40                          -1.275e-01             9 -57.877          134.3  6.08  0.006
    #> 80   1.951 -0.83530  0.42180 0.13510 0.4415                  0.70600                                8 -59.007          134.5  6.22  0.006
    #> 5    2.203                   0.42320                                                                4 -63.249          134.6  6.36  0.006
    #> 224  1.483 -0.98140  0.49340 0.34520 0.4569   37.83          0.16060 -0.4970                       10 -56.997          134.7  6.45  0.005
    #> 672  1.444 -0.80970  0.49510 0.35750 0.4797   33.19                  -0.5125            -0.0005682 10 -57.022          134.8  6.50  0.005
    #> 416  1.445 -0.80970  0.49780 0.35650 0.4788   30.05                  -0.5131 -9.738e-04            10 -57.022          134.8  6.50  0.005
    #> 332  2.060 -0.73950  0.75970         0.5324                  0.75250         -1.818e-01             8 -59.169          134.8  6.54  0.005
    #> 176  1.523 -1.75600  0.97530 0.20890 0.4554          1.0470          -0.4025                        9 -58.111          134.8  6.54  0.005
    #> 590  1.546 -0.84910          0.26430 0.5043                  0.69230                    -0.0454200  8 -59.216          134.9  6.64  0.005
    #> 56   2.765 -3.14600  0.26410 0.22520          20.36  1.3760                                         8 -59.221          134.9  6.65  0.005
    #> 112  2.082 -2.34100  0.40920 0.04015 0.4795          1.0320  0.55050                                9 -58.191          135.0  6.70  0.005
    #> 7    2.034           0.80290 0.40000                                                                5 -62.393          135.0  6.72  0.005
    #> 655  1.362           0.80260 0.28700 0.5080                          -0.2618            -0.0063460  8 -59.330          135.1  6.86  0.004
    #> 144  1.382 -0.02074  0.82000 0.27750 0.4949                          -0.2698                        8 -59.332          135.1  6.87  0.004
    #> 399  1.376           0.83470 0.27570 0.4990                          -0.2679 -1.284e-02             8 -59.332          135.1  6.87  0.004
    #> 304  1.959 -1.75900  0.76360 0.03397 0.5650          1.0790                  -1.612e-01             9 -58.292          135.2  6.90  0.004
    #> 560  1.938 -1.75500  0.46630 0.07946 0.5562          1.0650                             -0.0183300  9 -58.319          135.2  6.96  0.004
    #> 622  2.018 -2.36400          0.09076 0.5452          1.0040  0.63730                    -0.0198700  9 -58.327          135.2  6.97  0.004
    #> 448  1.937 -2.55200  0.20510 0.18990 0.4885   21.69  1.0830          -0.3451  7.605e-02            11 -56.280          135.4  7.16  0.004
    #> 256  1.890 -2.44700  0.31900 0.18750 0.5216   30.99  1.1140 -0.13040 -0.3329                       11 -56.282          135.4  7.16  0.004
    #> 704  1.921 -2.52600  0.33430 0.17940 0.5018   74.15  1.0740          -0.3394             0.0051310 11 -56.288          135.4  7.17  0.004
    #> 783  1.401           0.55540 0.22940 0.5477                                  -8.203e-02 -0.0200800  8 -59.533          135.5  7.27  0.004
    #> 528  1.427 -0.05600  0.42280 0.23850 0.5282                                             -0.0228700  8 -59.543          135.5  7.29  0.003
    #> 272  1.470 -0.04976  0.59090 0.18880 0.5094                                  -9.440e-02             8 -59.558          135.6  7.32  0.003
    #> 152  1.913 -0.93410  0.95620 0.65090          20.14                  -0.6027                        8 -59.560          135.6  7.32  0.003
    #> 24   2.154 -0.69120  0.25110 0.39910          18.80                                                 7 -60.686          135.7  7.47  0.003
    #> 608  1.834 -1.11300 -0.09708 0.24060 0.5457   23.46          0.40290                    -0.0405500 10 -57.568          135.8  7.59  0.003
    #>  [ reached getOption("max.print") -- omitted 51 rows ]
    #> Models ranked by used.criterion(x)
    #> Random terms (all models):
    #> '1 | individ', '1 | year2'

Before moving on, we note that the models created above are actually only usable for the construction of historical models. For comparison, we may wish to estimate ahistorical models. In that case, we also need linear models in which the global models tested do not include state at time t-1. Here, we produce these models.

    cypmodels2 <- modelsearch(vertdata, historical = FALSE, approach = "lme4",
                              vitalrates = c("surv", "obs", "size", "repst", "fec"),
                              sizedist = "poisson", fecdist = "poisson", suite = "full",
                              size = c("size3added", "size2added"))
    #> Warning in modelsearch(vertdata, historical = FALSE, approach = "lme4", : WARNING: Fecundity in time t cannot be Poisson-distributed and include 0s. Will develop fecundity models excluding all 0s. Consider adding a reproductive status variable to absorb 0 values.
    #>
    #> Developing global model of survival probability...
    #> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0330463 (tol =
    #> 0.002, component 1)
    #>
    #> Developing global model of observation probability...
    #>
    #>
    #> Developing global model of size (Poisson)...
    #>
    #>
    #> Developing global model of the probability of reproduction...
    #>
    #>
    #> Developing global model of fecundity (Poisson)...
    #> fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
    #> boundary (singular) fit: see ?isSingular
    #>
    #> All global models developed.
    #>
    #>
    #> Commencing dredge of survival probability...
    #> Fixed term is "(Intercept)"
    #> boundary (singular) fit: see ?isSingular
    #> boundary (singular) fit: see ?isSingular
    #> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0268727 (tol =
    #> 0.002, component 1)
    #> boundary (singular) fit: see ?isSingular
    #> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : unable to evaluate scaled gradient
    #> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge: degenerate Hessian with 1
    #> negative eigenvalues
    #> Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
    #> not positive definite or contains NA values: falling back to var-cov estimated from RX
    #> Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
    #> not positive definite or contains NA values: falling back to var-cov estimated from RX
    #>
    #> Commencing dredge of observation probability...
    #> Fixed term is "(Intercept)"
    #> boundary (singular) fit: see ?isSingular
    #>
    #> Commencing dredge of size...
    #> Fixed term is "(Intercept)"
    #>
    #> Commencing dredge of reproduction probability...
    #> Fixed term is "(Intercept)"
    #>
    #> Commencing dredge of fecundity...
    #> Fixed term is "(Intercept)"
    #> boundary (singular) fit: see ?isSingular
    #> fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
    #> boundary (singular) fit: see ?isSingular
    #> boundary (singular) fit: see ?isSingular
    #> fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
    #> boundary (singular) fit: see ?isSingular
    #> fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
    #> boundary (singular) fit: see ?isSingular
    #>
    #> Finished dredging all vital rates.
    #> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0268727 (tol =
    #> 0.002, component 1)
    #> boundary (singular) fit: see ?isSingular
    #> boundary (singular) fit: see ?isSingular
    summary(cypmodels2)
    #> This LefkoMod object includes 5 linear models.
    #> Best-fit model criterion used: AICc&k
    #>
    #> ────────────────────────────────────────────────────────────────────────────────
    #> Survival model:
    #> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
    #>  Family: binomial  ( logit )
    #> Formula: alive3 ~ size2added + (1 | individ) + (1 | year2)
    #>    Data: surv.data
    #>      AIC      BIC   logLik deviance df.resid
    #> 128.1324 143.2057 -60.0662 120.1324      316
    #> Random effects:
    #>  Groups  Name        Std.Dev.
    #>  individ (Intercept) 1.198378
    #>  year2   (Intercept) 0.008826
    #> Number of obs: 320, groups:  individ, 74; year2, 5
    #> Fixed Effects:
    #> (Intercept)   size2added
    #>      2.0352       0.6344
    #> convergence code 0; 0 optimizer warnings; 1 lme4 warnings
    #>
    #> ────────────────────────────────────────
    #>
    #> Observation model:
    #> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
    #>  Family: binomial  ( logit )
    #> Formula: obsstatus3 ~ size2added + (1 | individ) + (1 | year2)
    #>    Data: obs.data
    #>      AIC      BIC   logLik deviance df.resid
    #> 118.2567 133.1117 -55.1284 110.2567      299
    #> Random effects:
    #>  Groups  Name        Std.Dev.
    #>  individ (Intercept) 1.078e-05
    #>  year2   (Intercept) 8.776e-01
    #> Number of obs: 303, groups:  individ, 70; year2, 5
    #> Fixed Effects:
    #> (Intercept)   size2added
    #>      2.4904       0.3134
    #> convergence code 0; 0 optimizer warnings; 1 lme4 warnings
    #>
    #> ────────────────────────────────────────
    #>
    #> Size model:
    #> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
    #>  Family: poisson  ( log )
    #> Formula: size3added ~ size2added + (1 | individ) + (1 | year2)
    #>    Data: size.data
    #>       AIC       BIC    logLik  deviance  df.resid
    #> 1115.3617 1130.0135 -553.6808 1107.3617       284
    #> Random effects:
    #>  Groups  Name        Std.Dev.
    #>  individ (Intercept) 0.55311
    #>  year2   (Intercept) 0.09895
    #> Number of obs: 288, groups:  individ, 70; year2, 5
    #> Fixed Effects:
    #> (Intercept)   size2added
    #>     0.83893      0.04244
    #>
    #> ────────────────────────────────────────
    #>
    #> Reproductive status model:
    #> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
    #>  Family: binomial  ( logit )
    #> Formula: repstatus3 ~ repstatus2 + size2added + (1 | individ) + (1 | year2)
    #>    Data: repst.data
    #>       AIC       BIC    logLik  deviance  df.resid
    #>  333.6176  351.9324 -161.8088  323.6176       283
    #> Random effects:
    #>  Groups  Name        Std.Dev.
    #>  individ (Intercept) 0.1829
    #>  year2   (Intercept) 0.6250
    #> Number of obs: 288, groups:  individ, 70; year2, 5
    #> Fixed Effects:
    #> (Intercept)   repstatus2   size2added
    #>     -1.4630       1.6457       0.1715
    #>
    #> ────────────────────────────────────────
    #>
    #> Fecundity model:
    #> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
    #>  Family: poisson  ( log )
    #> Formula: feca2 ~ size2added + (1 | individ) + (1 | year2)
    #>    Data: fec.data
    #>      AIC      BIC   logLik deviance df.resid
    #> 156.8531 164.5012 -74.4266 148.8531       46
    #> Random effects:
    #>  Groups  Name        Std.Dev.
    #>  individ (Intercept) 0.000
    #>  year2   (Intercept) 0.172
    #> Number of obs: 50, groups:  individ, 26; year2, 5
    #> Fixed Effects:
    #> (Intercept)   size2added
    #>     0.21931      0.04429
    #> convergence code 0; 0 optimizer warnings; 1 lme4 warnings
    #>
    #> ────────────────────────────────────────────────────────────────────────────────
    #> Juvenile survival model:
    #> [1] 1
    #>
    #> ────────────────────────────────────────
    #>
    #> Juvenile observation model:
    #> [1] 1
    #>
    #> ────────────────────────────────────────
    #>
    #> Juvenile size model:
    #> [1] 1
    #>
    #> ────────────────────────────────────────
    #>
    #> Juvenile reproduction model:
    #> [1] 1
    #>
    #>
    #> ────────────────────────────────────────────────────────────────────────────────
    #>
    #> Number of models in survival table:5
    #>
    #> Number of models in observation table:5
    #>
    #> Number of models in size table:5
    #>
    #> Number of models in reproduction status table:5
    #>
    #> Number of models in fecundity table:5
    #>
    #> 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 reproduction table: 1
    #>
    #>
    #> ────────────────────────────────────────────────────────────────────────────────
    #>
    #> General model parameter names (column 1), and specific names used in these models (column 2):
    #>    mainparams modelparams
    #> 1       year2       year2
    #> 2     individ     individ
    #> 3       patch        <NA>
    #> 4       surv3      alive3
    #> 5        obs3  obsstatus3
    #> 6       size3  size3added
    #> 7      repst3  repstatus3
    #> 8        fec3       feca3
    #> 9        fec2       feca2
    #> 10      size2  size2added
    #> 11      size1        <NA>
    #> 12     repst2  repstatus2
    #> 13     repst1        <NA>
    #> 14        age        <NA>
    #>
    #>
    #> ────────────────────────────────────────────────────────────────────────────────
    #>
    #> Quality control:
    #>
    #> Survival estimated with 74 individuals and 320 individual transitions.
    #> Observation estimated with 70 individuals and 303 individual transitions.
    #> Size estimated with 70 individuals and 288 individual transitions.
    #> Reproductive status estimated with 70 individuals and 288 individual transitions.
    #> Fecundity estimated with 26 individuals and 50 individual transitions.
    #> Juvenile survival not estimated.
    #> Juvenile observation probability not estimated.
    #> Juvenile size transition not estimated.
    #> Juvenile reproduction probability not estimated.
    #> NULL

Fewer models were estimated per dredge, since fewer parameters were tested in the global models (size and reproductive status in time t-1 were not included). So, the best-fit models should look a little bit different. However, a more thorough comparison will show that many of the best-fit models are similar between historical and ahistorical analysis - indeed, size appears to be the only major vital rate to differ. This is not guaranteed - in this case, it may be that the relatively small number of years and small overall sample size leaves too little power to find an impact of historical status on most vital rates.

### Step 4. Projection matrix creation

Now we will create population function-based matrices for each year of the study. Function-based matrices have quickly been taking over in population ecology, probably because of their ability to parse out interesting trends and influential factors picked up by the linear modeling of vital rates. Let’s first create a set of ahistorical matrices.

    cypmatrix2 <- flefko2(stageframe = cypframe, repmatrix = rep.assumptions,
                          modelsuite = cypmodels2, overwrite = cypover,
                          data = vertdata, year.as.random = TRUE)

    summary(cypmatrix2)
    #>
    #> This lefkoMat object contains 5 matrices.
    #>
    #> Each matrix is a square matrix with 54 rows and columns, and a total of 2916 elements.
    #> A total of 12035 survival-transitions were estimated, with 2407 per matrix.
    #> A total of 240 fecundity transitions were estimated, with 48 per matrix.
    #>
    #> Vital rate modeling quality control:
    #>
    #> Survival estimated with 74 individuals and 320 individual transitions.
    #> Observation estimated with 70 individuals and 303 individual transitions.
    #> Size estimated with 70 individuals and 288 individual transitions.
    #> Reproductive status estimated with 70 individuals and 288 individual transitions.
    #> Fecundity estimated with 26 individuals and 50 individual transitions.
    #> Juvenile survival estimated with 0 individuals and 0 individual transitions.
    #> Juvenile observation estimated with 0 individuals and 0 individual transitions.
    #> Juvenile size estimated with 0 individuals and 0 individual transitions.
    #> Juvenile reproductive status estimated with 0 individuals and 0 individual transitions.
    #> NULL

A quick glance at the summary output will highlight that many more elements are estimated for function-based matrices. In raw matrices, elements associated with transitions from specific stages are only estimated when individuals actually exist within those particular stages. In function-based matrices, in contrast, the linear models estimated allow the estimation of all elements that are theoretically possible (i.e. only structural 0s are not estimated). Let’s take a look at an example matrix, but only on the top corner to deal with its size.

    cypmatrix2$A[[1]][1:25, 1:8]
    #>       [,1] [,2] [,3] [,4] [,5]         [,6]         [,7]         [,8]
    #>  [1,]  0.1  0.0  0.0 0.00  0.0 0.000000e+00 0.000000e+00 0.000000e+00
    #>  [2,]  0.2  0.0  0.0 0.00  0.0 0.000000e+00 0.000000e+00 0.000000e+00
    #>  [3,]  0.0  0.2  0.0 0.00  0.0 0.000000e+00 0.000000e+00 0.000000e+00
    #>  [4,]  0.0  0.0  0.2 0.00  0.0 0.000000e+00 0.000000e+00 0.000000e+00
    #>  [5,]  0.0  0.0  0.0 0.25  0.4 0.000000e+00 0.000000e+00 0.000000e+00
    #>  [6,]  0.0  0.0  0.0 0.00  0.0 4.701977e-02 3.686912e-02 2.809353e-02
    #>  [7,]  0.0  0.0  0.0 0.00  0.0 1.005860e-01 9.518443e-02 8.657188e-02
    #>  [8,]  0.0  0.0  0.0 0.00  0.0 1.491216e-01 1.472310e-01 1.397142e-01
    #>  [9,]  0.0  0.0  0.0 0.00  0.0 1.473847e-01 1.518244e-01 1.503187e-01
    #> [10,]  0.0  0.0  0.0 0.00  0.0 1.092510e-01 1.174208e-01 1.212961e-01
    #> [11,]  0.0  0.0  0.0 0.00  0.0 6.478709e-02 7.265046e-02 7.830161e-02
    #> [12,]  0.0  0.0  0.0 0.00  0.0 3.201624e-02 3.745852e-02 4.212241e-02
    #> [13,]  0.0  0.0  0.0 0.00  0.0 1.356143e-02 1.655450e-02 1.942267e-02
    #> [14,]  0.0  0.0  0.0 0.00  0.0 5.026301e-03 6.401613e-03 7.836331e-03
    #> [15,]  0.0  0.0  0.0 0.00  0.0 1.655919e-03 2.200444e-03 2.810373e-03
    #> [16,]  0.0  0.0  0.0 0.00  0.0 4.909894e-04 6.807282e-04 9.071055e-04
    #> [17,]  0.0  0.0  0.0 0.00  0.0 1.323465e-04 1.914452e-04 2.661699e-04
    #> [18,]  0.0  0.0  0.0 0.00  0.0 3.270125e-05 4.935449e-05 7.159314e-05
    #> [19,]  0.0  0.0  0.0 0.00  0.0 7.458545e-06 1.174483e-05 1.777550e-05
    #> [20,]  0.0  0.0  0.0 0.00  0.0 1.579644e-06 2.595268e-06 4.098148e-06
    #> [21,]  0.0  0.0  0.0 0.00  0.0 3.122490e-07 5.352471e-07 8.818407e-07
    #> [22,]  0.0  0.0  0.0 0.00  0.0 5.786475e-08 1.034899e-07 1.778951e-07
    #> [23,]  0.0  0.0  0.0 0.00  0.0 1.009249e-08 1.883269e-08 3.377604e-08
    #> [24,]  0.0  0.0  0.0 0.00  0.0 1.662489e-09 3.236706e-09 6.056616e-09
    #> [25,]  0.0  0.0  0.0 0.00  0.0 2.594408e-10 5.270030e-10 1.028893e-09

The matrix is overwhelmingly composed of non-zero elements, unlike in the raw matrix case.

Next, we will create a set of historical Lefkovitch matrices.

    cypmatrix3 <- flefko3(stageframe = cypframe, repmatrix = rep.assumptions,
                          modelsuite = cypmodels3, overwrite = cypover,
                          data = vertdata, yearcol = "year2", year.as.random = TRUE)

    summary(cypmatrix3)
    #>
    #> This lefkoMat object contains 5 matrices.
    #>
    #> Each matrix is a square matrix with 2916 rows and columns, and a total of 8503056 elements.
    #> A total of 588520 survival-transitions were estimated, with 117704 per matrix.
    #> A total of 12960 fecundity transitions were estimated, with 2592 per matrix.
    #>
    #> Vital rate modeling quality control:
    #>
    #> Survival estimated with 74 individuals and 320 individual transitions.
    #> Observation estimated with 70 individuals and 303 individual transitions.
    #> Size estimated with 70 individuals and 288 individual transitions.
    #> Reproductive status estimated with 70 individuals and 288 individual transitions.
    #> Fecundity estimated with 26 individuals and 50 individual transitions.
    #> Juvenile survival estimated with 0 individuals and 0 individual transitions.
    #> Juvenile observation estimated with 0 individuals and 0 individual transitions.
    #> Juvenile size estimated with 0 individuals and 0 individual transitions.
    #> Juvenile reproductive status estimated with 0 individuals and 0 individual transitions.
    #> NULL

Once again, we see many more elements estimated, and many more rows and columns (54 rows and columns in the ahistorical case, and 542 = 2916 rows and columns in the historical case). However, the dominance of structural 0s in historical matrices still yields matrices that are mostly composed of 0s. A quick glance at one matrix will show that. We will focus on only one small section of that matrix.

    cypmatrix3$A[[1]][2001:2050,2036:2045]
    #>               [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
    #>  [1,] 0.000000e+00    0    0    0    0    0    0    0    0     0
    #>  [2,] 0.000000e+00    0    0    0    0    0    0    0    0     0
    #>  [3,] 0.000000e+00    0    0    0    0    0    0    0    0     0
    #>  [4,] 4.551147e-03    0    0    0    0    0    0    0    0     0
    #>  [5,] 3.218013e-05    0    0    0    0    0    0    0    0     0
    #>  [6,] 1.743937e-04    0    0    0    0    0    0    0    0     0
    #>  [7,] 6.300612e-04    0    0    0    0    0    0    0    0     0
    #>  [8,] 1.707245e-03    0    0    0    0    0    0    0    0     0
    #>  [9,] 3.700829e-03    0    0    0    0    0    0    0    0     0
    #> [10,] 6.685300e-03    0    0    0    0    0    0    0    0     0
    #> [11,] 1.035133e-02    0    0    0    0    0    0    0    0     0
    #> [12,] 1.402423e-02    0    0    0    0    0    0    0    0     0
    #> [13,] 1.688923e-02    0    0    0    0    0    0    0    0     0
    #> [14,] 1.830555e-02    0    0    0    0    0    0    0    0     0
    #> [15,] 1.803696e-02    0    0    0    0    0    0    0    0     0
    #> [16,] 1.629128e-02    0    0    0    0    0    0    0    0     0
    #> [17,] 1.358267e-02    0    0    0    0    0    0    0    0     0
    #> [18,] 1.051550e-02    0    0    0    0    0    0    0    0     0
    #> [19,] 7.598222e-03    0    0    0    0    0    0    0    0     0
    #> [20,] 5.147130e-03    0    0    0    0    0    0    0    0     0
    #> [21,] 3.281628e-03    0    0    0    0    0    0    0    0     0
    #> [22,] 1.976014e-03    0    0    0    0    0    0    0    0     0
    #> [23,] 1.127222e-03    0    0    0    0    0    0    0    0     0
    #> [24,] 6.108754e-04    0    0    0    0    0    0    0    0     0
    #> [25,] 3.152873e-04    0    0    0    0    0    0    0    0     0
    #> [26,] 1.553306e-04    0    0    0    0    0    0    0    0     0
    #> [27,] 7.319850e-05    0    0    0    0    0    0    0    0     0
    #> [28,] 3.305705e-05    0    0    0    0    0    0    0    0     0
    #> [29,] 1.794085e-04    0    0    0    0    0    0    0    0     0
    #> [30,] 9.722686e-04    0    0    0    0    0    0    0    0     0
    #> [31,] 3.512676e-03    0    0    0    0    0    0    0    0     0
    #> [32,] 9.518122e-03    0    0    0    0    0    0    0    0     0
    #> [33,] 2.063262e-02    0    0    0    0    0    0    0    0     0
    #> [34,] 3.727145e-02    0    0    0    0    0    0    0    0     0
    #> [35,] 5.771004e-02    0    0    0    0    0    0    0    0     0
    #> [36,] 7.818699e-02    0    0    0    0    0    0    0    0     0
    #> [37,] 9.415971e-02    0    0    0    0    0    0    0    0     0
    #> [38,] 1.020559e-01    0    0    0    0    0    0    0    0     0
    #> [39,] 1.005585e-01    0    0    0    0    0    0    0    0     0
    #> [40,] 9.082608e-02    0    0    0    0    0    0    0    0     0
    #> [41,] 7.572519e-02    0    0    0    0    0    0    0    0     0
    #> [42,] 5.862535e-02    0    0    0    0    0    0    0    0     0
    #> [43,] 4.236111e-02    0    0    0    0    0    0    0    0     0
    #> [44,] 2.869594e-02    0    0    0    0    0    0    0    0     0
    #> [45,] 1.829552e-02    0    0    0    0    0    0    0    0     0
    #> [46,] 1.101654e-02    0    0    0    0    0    0    0    0     0
    #> [47,] 6.284415e-03    0    0    0    0    0    0    0    0     0
    #> [48,] 3.405712e-03    0    0    0    0    0    0    0    0     0
    #> [49,] 1.757769e-03    0    0    0    0    0    0    0    0     0
    #> [50,] 8.659888e-04    0    0    0    0    0    0    0    0     0

Now let’s estimate some mean matrices. In this case, with all relevant matrix elements holding non-zero values, we may use the geometric mean across time, as follows for the ahistorical matrices.

    tmeans2r <- lmean(cypmatrix2, time = "geometric")

And now the historical mean. We will use this and the ahistorical mean in analyses in the next step.

    tmeans3r <- lmean(cypmatrix3, time = "geometric")

### Step 5. Analyses

Let’s now conduct some quick population analyses. First, we can estimate the deterministic growth rate for these means. First the ahistorical case.

    lambda3(tmeans2r)
    #>   pop patch    lambda
    #> 1   1     1 0.8868898
    #> 2   1   All 0.8868898

Now the historical case.

    lambda3(tmeans3r)
    #>   pop patch    lambda
    #> 1   1     1 0.8159135
    #> 2   1   All 0.8159135

We find that the historical lambda is lower than the ahistorical lambda, suggesting that trade-offs are working across time to offset survival and/or reproduction.

Now let’s take a look at the stable stage distributions. First the ahistorical case, with summary.

    tm2ss <- stablestage3(tmeans2r)
    summary(tm2ss)
    #>      matrix     new_stage_id  orig_stage_id      original_size      ss_prop
    #>  Min.   :1.0   Min.   : 1.0   Length:108         Min.   : 0.00   Min.   :0.000e+00
    #>  1st Qu.:1.0   1st Qu.:14.0   Class :character   1st Qu.: 4.00   1st Qu.:0.000e+00
    #>  Median :1.5   Median :27.5   Mode  :character   Median :11.00   Median :9.451e-05
    #>  Mean   :1.5   Mean   :27.5                      Mean   :11.11   Mean   :1.852e-02
    #>  3rd Qu.:2.0   3rd Qu.:41.0                      3rd Qu.:18.00   3rd Qu.:1.934e-02
    #>  Max.   :2.0   Max.   :54.0                      Max.   :24.00   Max.   :2.283e-01

Now the first 6 rows.

    head(tm2ss)
    #>   matrix new_stage_id orig_stage_id original_size     ss_prop
    #> 1      1            1            SD             0 0.205168634
    #> 2      1            2            P1             0 0.228302122
    #> 3      1            3            P2             0 0.051483779
    #> 4      1            4            P3             0 0.011609947
    #> 5      1            5            SL             0 0.005961283
    #> 6      1            6             D             0 0.019337881

Next, the historical case. We will skip the historical stage-pair distributions and look directly at the predicted distributions of ahistorical stages.

    tm3ss <- stablestage3(tmeans3r)
    summary(tm3ss$ahist)
    #>      matrix     new_stage_id     ss_prop
    #>  Min.   :1.0   Min.   : 1.0   Min.   :0.0000755
    #>  1st Qu.:1.0   1st Qu.:14.0   1st Qu.:0.0014256
    #>  Median :1.5   Median :27.5   Median :0.0065553
    #>  Mean   :1.5   Mean   :27.5   Mean   :0.0185185
    #>  3rd Qu.:2.0   3rd Qu.:41.0   3rd Qu.:0.0127742
    #>  Max.   :2.0   Max.   :54.0   Max.   :0.3245035

And the first 6 rows.

    head(tm3ss$ahist)
    #>   matrix new_stage_id     ss_prop
    #> 1      1            1 0.289074019
    #> 2      1            2 0.324503508
    #> 3      1            3 0.017369247
    #> 4      1            4 0.004257621
    #> 5      1            5 0.002559200
    #> 6      1            6 0.009451046

Overall, these look like fairly large shifts. For comparison, we can try graphing these together.

    plot(tm2ss$ss_prop, ylab = "Stable Stage Proportion", xlab = "Stage Index", ylim = c(0, 0.35), type = "l", lty = 1, lwd = 2)
    lines(tm3ss$ahist[,3], col = "red", lty = 2, lwd = 2)
    legend("topright", c("ahistorical", "historical"), lty = c(1, 2), col = c("black", "red"), lwd = 2, bty = "n")

plot of chunk Ch3an2st5ln710

Although the two have much general agreement, there is nonetheless evidence of some disagreement. We leave it to the user to assess these differences.

Finally, let’s take a peek at the reproductive values associated with both ahistorical and historical approaches. First, the ahistorical set.

    repvalue3(tmeans2r)
    #>     matrix new_stage_id orig_stage_id original_size left_vector rep_value
    #> 1        1            1            SD             0   0.0000000  0.000000
    #> 2        1            2            P1             0   0.0000000  0.000000
    #> 3        1            3            P2             0   0.0000000  0.000000
    #> 4        1            4            P3             0   0.0000000  0.000000
    #> 5        1            5            SL             0   0.0000000  0.000000
    #> 6        1            6             D             0  -0.1149892  1.000000
    #> 7        1            7            V1             1  -0.1230093  1.069747
    #> 8        1            8            V2             2  -0.1281946  1.114841
    #> 9        1            9            V3             3  -0.1315767  1.144253
    #> 10       1           10            V4             4  -0.1338763  1.164251
    #> 11       1           11            V5             5  -0.1355452  1.178765
    #> 12       1           12            V6             6  -0.1368542  1.190149
    #> 13       1           13            V7             7  -0.1379633  1.199794
    #> 14       1           14            V8             8  -0.1389671  1.208524
    #> 15       1           15            V9             9  -0.1399212  1.216821
    #> 16       1           16           V10            10  -0.1408574  1.224963
    #> 17       1           17           V11            11  -0.1417924  1.233094
    #> 18       1           18           V12            12  -0.1427327  1.241271
    #> 19       1           19           V13            13  -0.1436783  1.249495
    #> 20       1           20           V14            14  -0.1446251  1.257728
    #> 21       1           21           V15            15  -0.1455662  1.265913
    #> 22       1           22           V16            16  -0.1464943  1.273983
    #> 23       1           23           V17            17  -0.1474018  1.281876
    #> 24       1           24           V18            18  -0.1482822  1.289532
    #> 25       1           25           V19            19  -0.1491304  1.296908
    #> 26       1           26           V20            20  -0.1499427  1.303972
    #> 27       1           27           V21            21  -0.1507169  1.310705
    #> 28       1           28           V22            22  -0.1514522  1.317100
    #> 29       1           29           V23            23  -0.1521491  1.323160
    #> 30       1           30           V24            24  -0.1528086  1.328896
    #> 31       1           31            F1             1  -0.1218207  1.059411
    #> 32       1           32            F2             2  -0.1274049  1.107974
    #> 33       1           33            F3             3  -0.1312803  1.141676
    #> 34       1           34            F4             4  -0.1341254  1.166417
    #> 35       1           35            F5             5  -0.1363550  1.185807
    #> 36       1           36            F6             6  -0.1382082  1.201924
    #> 37       1           37            F7             7  -0.1398184  1.215927
    #> 38       1           38            F8             8  -0.1412592  1.228456
    #> 39       1           39            F9             9  -0.1425710  1.239864
    #> 40       1           40           F10            10  -0.1437769  1.250352
    #> 41       1           41           F11            11  -0.1448909  1.260040
    #> 42       1           42           F12            12  -0.1459225  1.269011
    #> 43       1           43           F13            13  -0.1468789  1.277328
    #> 44       1           44           F14            14  -0.1477666  1.285048
    #> 45       1           45           F15            15  -0.1485914  1.292221
    #> 46       1           46           F16            16  -0.1493591  1.298898
    #> 47       1           47           F17            17  -0.1500753  1.305125
    #> 48       1           48           F18            18  -0.1507451  1.310951
    #> 49       1           49           F19            19  -0.1513738  1.316418
    #> 50       1           50           F20            20  -0.1519658  1.321567
    #> 51       1           51           F21            21  -0.1525257  1.326435
    #> 52       1           52           F22            22  -0.1530571  1.331057
    #> 53       1           53           F23            23  -0.1535637  1.335462
    #> 54       1           54           F24            24  -0.1540484  1.339677
    #> 55       2            1            SD             0   0.0000000  0.000000
    #> 56       2            2            P1             0   0.0000000  0.000000
    #> 57       2            3            P2             0   0.0000000  0.000000
    #> 58       2            4            P3             0   0.0000000  0.000000
    #> 59       2            5            SL             0   0.0000000  0.000000
    #> 60       2            6             D             0  -0.1149892  1.000000
    #> 61       2            7            V1             1  -0.1230093  1.069747
    #> 62       2            8            V2             2  -0.1281946  1.114841
    #> 63       2            9            V3             3  -0.1315767  1.144253
    #> 64       2           10            V4             4  -0.1338763  1.164251
    #> 65       2           11            V5             5  -0.1355452  1.178765
    #> 66       2           12            V6             6  -0.1368542  1.190149
    #> 67       2           13            V7             7  -0.1379633  1.199794
    #> 68       2           14            V8             8  -0.1389671  1.208524
    #> 69       2           15            V9             9  -0.1399212  1.216821
    #> 70       2           16           V10            10  -0.1408574  1.224963
    #> 71       2           17           V11            11  -0.1417924  1.233094
    #> 72       2           18           V12            12  -0.1427327  1.241271
    #> 73       2           19           V13            13  -0.1436783  1.249495
    #> 74       2           20           V14            14  -0.1446251  1.257728
    #> 75       2           21           V15            15  -0.1455662  1.265913
    #> 76       2           22           V16            16  -0.1464943  1.273983
    #> 77       2           23           V17            17  -0.1474018  1.281876
    #> 78       2           24           V18            18  -0.1482822  1.289532
    #> 79       2           25           V19            19  -0.1491304  1.296908
    #> 80       2           26           V20            20  -0.1499427  1.303972
    #> 81       2           27           V21            21  -0.1507169  1.310705
    #> 82       2           28           V22            22  -0.1514522  1.317100
    #> 83       2           29           V23            23  -0.1521491  1.323160
    #> 84       2           30           V24            24  -0.1528086  1.328896
    #> 85       2           31            F1             1  -0.1218207  1.059411
    #> 86       2           32            F2             2  -0.1274049  1.107974
    #> 87       2           33            F3             3  -0.1312803  1.141676
    #> 88       2           34            F4             4  -0.1341254  1.166417
    #> 89       2           35            F5             5  -0.1363550  1.185807
    #> 90       2           36            F6             6  -0.1382082  1.201924
    #> 91       2           37            F7             7  -0.1398184  1.215927
    #> 92       2           38            F8             8  -0.1412592  1.228456
    #> 93       2           39            F9             9  -0.1425710  1.239864
    #> 94       2           40           F10            10  -0.1437769  1.250352
    #> 95       2           41           F11            11  -0.1448909  1.260040
    #> 96       2           42           F12            12  -0.1459225  1.269011
    #> 97       2           43           F13            13  -0.1468789  1.277328
    #> 98       2           44           F14            14  -0.1477666  1.285048
    #> 99       2           45           F15            15  -0.1485914  1.292221
    #> 100      2           46           F16            16  -0.1493591  1.298898
    #> 101      2           47           F17            17  -0.1500753  1.305125
    #> 102      2           48           F18            18  -0.1507451  1.310951
    #> 103      2           49           F19            19  -0.1513738  1.316418
    #> 104      2           50           F20            20  -0.1519658  1.321567
    #> 105      2           51           F21            21  -0.1525257  1.326435
    #> 106      2           52           F22            22  -0.1530571  1.331057
    #> 107      2           53           F23            23  -0.1535637  1.335462
    #> 108      2           54           F24            24  -0.1540484  1.339677

And let’s compare these to the historical case.

    repvalue3(tmeans3r)$ahist
    #>     matrix new_stage_id  rep_value
    #> 1        1            1 0.00000000
    #> 2        1            2 0.00000000
    #> 3        1            3 0.00000000
    #> 4        1            4 0.00000000
    #> 5        1            5 0.00000000
    #> 6        1            6 1.00000000
    #> 7        1            7 1.11337163
    #> 8        1            8 1.25989440
    #> 9        1            9 1.36130980
    #> 10       1           10 1.40413383
    #> 11       1           11 1.37804439
    #> 12       1           12 1.28111047
    #> 13       1           13 1.12415468
    #> 14       1           14 0.93071691
    #> 15       1           15 0.73091305
    #> 16       1           16 0.55172715
    #> 17       1           17 0.40917911
    #> 18       1           18 0.30650513
    #> 19       1           19 0.23796212
    #> 20       1           20 0.19442578
    #> 21       1           21 0.16729948
    #> 22       1           22 0.15005974
    #> 23       1           23 0.13827432
    #> 24       1           24 0.12913590
    #> 25       1           25 0.12102582
    #> 26       1           26 0.11319538
    #> 27       1           27 0.10544336
    #> 28       1           28 0.09785310
    #> 29       1           29 0.09055702
    #> 30       1           30 0.08363655
    #> 31       1           31 1.24482371
    #> 32       1           32 1.35373424
    #> 33       1           33 1.40137667
    #> 34       1           34 1.37966649
    #> 35       1           35 1.28654989
    #> 36       1           36 1.13098199
    #> 37       1           37 0.93537863
    #> 38       1           38 0.73167196
    #> 39       1           39 0.55062241
    #> 40       1           40 0.41041837
    #> 41       1           41 0.31324037
    #> 42       1           42 0.25087072
    #> 43       1           43 0.21246065
    #> 44       1           44 0.18901643
    #> 45       1           45 0.17439549
    #> 46       1           46 0.16473622
    #> 47       1           47 0.15766164
    #> 48       1           48 0.15170294
    #> 49       1           49 0.14596621
    #> 50       1           50 0.13995298
    #> 51       1           51 0.13345689
    #> 52       1           52 0.12647446
    #> 53       1           53 0.11910600
    #> 54       1           54 0.11147341
    #> 55       2            1 0.00000000
    #> 56       2            2 0.00000000
    #> 57       2            3 0.00000000
    #> 58       2            4 0.00000000
    #> 59       2            5 0.00000000
    #> 60       2            6 1.00000000
    #> 61       2            7 1.11337163
    #> 62       2            8 1.25989440
    #> 63       2            9 1.36130980
    #> 64       2           10 1.40413383
    #> 65       2           11 1.37804439
    #> 66       2           12 1.28111047
    #> 67       2           13 1.12415468
    #> 68       2           14 0.93071691
    #> 69       2           15 0.73091305
    #> 70       2           16 0.55172715
    #> 71       2           17 0.40917911
    #> 72       2           18 0.30650513
    #> 73       2           19 0.23796212
    #> 74       2           20 0.19442578
    #> 75       2           21 0.16729948
    #> 76       2           22 0.15005974
    #> 77       2           23 0.13827432
    #> 78       2           24 0.12913590
    #> 79       2           25 0.12102582
    #> 80       2           26 0.11319538
    #> 81       2           27 0.10544336
    #> 82       2           28 0.09785310
    #> 83       2           29 0.09055702
    #> 84       2           30 0.08363655
    #> 85       2           31 1.24482371
    #> 86       2           32 1.35373424
    #> 87       2           33 1.40137667
    #> 88       2           34 1.37966649
    #> 89       2           35 1.28654989
    #> 90       2           36 1.13098199
    #> 91       2           37 0.93537863
    #> 92       2           38 0.73167196
    #> 93       2           39 0.55062241
    #> 94       2           40 0.41041837
    #> 95       2           41 0.31324037
    #> 96       2           42 0.25087072
    #> 97       2           43 0.21246065
    #> 98       2           44 0.18901643
    #> 99       2           45 0.17439549
    #> 100      2           46 0.16473622
    #> 101      2           47 0.15766164
    #> 102      2           48 0.15170294
    #> 103      2           49 0.14596621
    #> 104      2           50 0.13995298
    #> 105      2           51 0.13345689
    #> 106      2           52 0.12647446
    #> 107      2           53 0.11910600
    #> 108      2           54 0.11147341

Note that the historical case predicts a greater importance of small adults and lesser importance of large adults than the ahistorical matrix does.

## Acknowledgements

The project resulting in this package and this tutorial was funded by Grant-In-Aid 19H03298 from the Japan Society for the Promotion of Science.

## Literature cited

Bartoń, Kamil A. 2014. “MuMIn: Multi-Model Inference.” https://CRAN.R-project.org/package=MuMIn.

Bates, Douglas, Martin Maechler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using Lme4.” Journal of Statistical Software 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.

Burnham, Kenneth P., and David R. Anderson. 2002. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. New York, New York, USA: Springer-Verlag New York, Inc.

Ehrlén, Johan. 2000. “The Dynamics of Plant Populations: Does the History of Individuals Matter?” Ecology 81 (6): 1675–84. https://doi.org/10.1890/0012-9658(2000)081[1675:TDOPPD]2.0.CO;2.

Ellner, Stephen P., and Mark Rees. 2006. “Integral Projection Models for Species with Complex Demography.” The American Naturalist 167 (3): 410–28. https://doi.org/10.1086/499438.

Salguero-Gómez, Roberto, and J. B. Plotkin. 2010. “Matrix Dimensions Bias Demographic Inferences: Implications for Comparative Plant Demography.” The American Naturalist 176 (6): 710–22. https://doi.org/10.1086/657044.

Shefferson, Richard P., Ryo Mizuta, and Michael J. Hutchings. 2017. “Predicting Evolution in Response to Climate Change: The Example of Sprouting Probability in Three Dormancy-Prone Orchid Species.” Royal Society Open Science 4 (1): 160647. https://doi.org/10.1098/rsos.160647.

Shefferson, Richard P., and Deborah A. Roach. 2010. “Longitudinal Analysis of Plantago: Adaptive Benefits of Iteroparity in a Short-Lived, Herbaceous Perennial.” Ecology 91 (2): 441–47. https://doi.org/10.1890/09-0423.1.

Shefferson, Richard P., Brett K. Sandercock, Joyce Proper, and Steven R. Beissinger. 2001. “Estimating Dormancy and Survival of a Rare Herbaceous Perennial Using Mark-Recapture Models.” Ecology 82 (1): 145–56. https://doi.org/10.1890/0012-9658(2001)082[0145:EDASOA]2.0.CO;2.

Shefferson, Richard P., Robert J. Warren II, and H. Ronald Pulliam. 2014. “Life History Costs Make Perfect Sprouting Maladaptive in Two Herbaceous Perennials.” Journal of Ecology 102 (5): 1318–28. https://doi.org/10.1111/1365-2745.12281.