---
title: "Lathyrus vernus case studies"
author: Richard P. Shefferson
output: rmarkdown::html_vignette
bibliography: Lefko3Tutorial.bib
vignette: >
  %\VignetteIndexEntry{Lathyrus vernus case studies}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---



This document was built in Markdown in R 4.0.2.

## CASE STUDIES OF SWEDISH *Lathyrus vernus* POPULATION

### ORGANISM AND POPULATION

  <p style="text-indent: 20px">*Lathyrus vernus* (family Fabaceae) is a long-lived forest herb, native to Europe and large parts of northern Asia. Individuals increase slowly in size and usually flower only after 10-15 years of vegetative growth. Flowering individuals have an average conditional life span of 44.3 years [@ehrlen_how_2002]. *Lathyrus vernus* lacks organs for vegetative spread and individuals are well delimited [@ehrlen_assessing_2002]. One or several erect shoots of up to 40 cm height emerge from a subterranean rhizome in March-April. Flowering occurs about four weeks after shoot emergence. Shoot growth is determinate, and the number of flowers is determined in the previous year [@ehrlen_storage_2001]. Individuals may not produce above-ground structures every year but can remain dormant in one season. *Lathyrus vernus* is self-compatible but requires visits from bumble-bees to produce seeds. Individuals produce few, large seeds and establishment from seeds is relatively frequent [@ehrlen_seedling_1996,]. The pre-dispersal seed predator *Bruchus atomarius* often consumes a large fraction of developing seeds, and roe deers (*Capreolus capreolus*) sometimes consume the shoots [@ehrlen_timing_2009].</p>
  
  <p style="text-indent: 20px">Data for this study was collected from individuals in six permanent plots in a population of *L. vernus* in a deciduous forest in the Tullgarn area, SE Sweden (58.9496 N, 17.6097 E), during 1988–1991 [@ehrlen_demography_1995]. The six plots were relatively similar with regard to soil type, elevation, slope and canopy cover. Within each plot, all individuals were marked with numbered tags that remained over the study period, and their locations were carefully mapped. New individuals were included in the study in each year. Individuals were recorded at least three times every growth season. At the time of shoot emergence, we recorded whether individuals were alive and produced above-ground shoots, and if shoots had been grazed. During flowering, we recorded flower number and the height and diameter of all shoots. At fruit maturation, we counted the number of intact and damaged seeds. To derive a measure of above-ground size for each individual, we calculated the volume of each shoot as π × (diameter/2)<sup>2</sup> × height, and summed the volumes of all shoots. This measure is closely correlated with the dry mass of aboveground tissues (r<sup>2</sup> = 0.924, P < 0.001, n = 50, log-transformed values; Ehrlen 1995). Size of individuals that had been grazed was estimated based on measures of shoot diameter in grazed shoots, and the relationship between shoot diameter and shoot height in non-grazed individuals. Only individuals with an aboveground volume of more than 230 mm<sup>3</sup> flowered and produced fruits during this study. Individuals that lacked above-ground structures in one season but reappeared in the following year were considered dormant. Individuals that lacked above-ground structures in two subsequent seasons were considered dead from the year in which they first lacked above-ground structures. Probabilities of seeds surviving to the next year, and of being present as seedlings or seeds in the soil seed bank, were derived from separate yearly sowing experiments in separate plots adjacent to each subplot [@ehrlen_seedling_1996].</p>

### ANALYSES WITH *LATHYRUS* DATA

  <p style="text-indent: 20px">We will analyze these data in three different ways to illustrate the power and flexiblity of package `lefko3`:</p> 
  
  1) through the estimation of **raw matrices**, with the intention of producing matrices similar to those published in Ehrlén [-@ehrlen_dynamics_2000]; 
  
  2) through the estimation of function-derived matrices using a stage classification different from Ehrlén [-@ehrlen_dynamics_2000], developed using the natural logarithm of the size measure used in that study; and
  
  3) through the construction of an integral projection model.
  
### Analysis 1: Raw matrix estimation

  <p style="text-indent: 20px">Here, we will attempt to estimate matrices similar to and based on the dataset used in Ehrlén [-@ehrlen_dynamics_2000]. These matrices will not be the same, as the dataset currently includes more individuals for those years as well as an extra year of data, and also includes differences in assumptions regarding transitions to unobservable life history stages. However, they will be very similar.</p>
  
#### Step 1. Data characterization and reorganization

  <p style="text-indent: 20px">The dataset that we have provided is organized in horizontal format, meaning that rows correspond to unique individuals and columns correspond to individual condition in particular years. Looking at the original Excel spreadsheet, you will note a repeating pattern in the names of the columns. Package `lefko3` includes functions to handle data in horizontal format based on these patterns, 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*).</p>

<img src = "lathscreen1.png" width = "700"></img>

  <p style="text-indent: 20px">This dataset includes information on 1119 individuals, so there are 1119 rows with data (not counting the header). There are 38 columns. The first two columns are variables giving identifying information about each individual, with each individual's data entirely restricted to one row. This is followed by four sets of nine columns, each named `VolumeXX`, `lnVolXX`, `FCODEXX`, `FlowXX`, `IntactseedXX`, `Dead19XX`, `DormantXX`, `Missing19XX`, and `SeedlingXX`, with the `XX` in each case corresponding to the year of observation and with years organized consecutively. Thus, columns 3-11 refer to year 1988, columns 12-20 refer to year 1989, etc. To properly conduct this exercise, we need to know the exact number of years used, which is 4 years here (includes all years from and including 1988 to 1991), we need the columns to be repeated in the same order in each year, and we need the years in consecutive order with no extra columns between them.</p>

  <p style="text-indent: 20px">First, we load the dataset, and look at its dimensions.</p>

```r
rm(list=ls(all=TRUE))

library(lefko3)

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

  <p style="text-indent: 20px">After looking over the dataset, we need to create a **stageframe** describing the life history of the species and linking it to the data. A stageframe is a data frame that describes all stages in the life history of the organism, in a way usable by the functions in this package. It needs to include complete descriptions of all stages that occur in the dataset, with each stage defined uniquely. Since this object can be used for automated classification of individuals, all sizes, reproductive states, and other characteristics defining each stage in the dataset need to be accounted for. This can be difficult if a few data points exist outside the range of sizes specified in the stageframe, so great care must be taken to include all size values and values of other descriptor variables occurring within the dataset. The final description of each stage occurring in the dataset may not overlap completely with any other stage also found in the dataset, although partial overlap is expected.</p>
  
  <p style="text-indent: 20px">Here, we create a stageframe named `lathframe` based on Ehrlén [-@ehrlen_dynamics_2000]. We build this by creating vectors of the values describing each stage, always in the same order. Of particular note are the vectors input as `sizes` and `binhalfwidth` in the `sf_create` function. In the case where sizes are binned, the values input in the former are the central values of each bin, while the latter represents one-half of the width of the bin. If size values are not to be binned, then narrow binwidths can be used. For example, in this dataset, vegetatively dormant individuals necessarily have a size of 0, and so we can set the `halfbinwidth` for this stage to 0.5.</p>

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

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

To be useful to others reading your work, or to yourself in the future, it helps to add text descriptions of the stages. Here, we also add some descriptions to this stageframe as comments.

```r
lathframe$comments <- c("Dormant seed", "Seedling", "Very small vegetative", 
                        "Small vegetative", "Very large vegetative", 
                        "Flowering", "Dormant")
lathframe
#>   stagenames size repstatus obsstatus propstatus immstatus matstatus indataset binhalfwidth_raw sizebin_min sizebin_max sizebin_center
#> 1         Sd    0         0         0          1         1         0         0              0.0         0.0         0.0              0
#> 2        Sdl  100         0         1          0         1         0         1            100.0         0.0       200.0            100
#> 3        VSm   13         0         1          0         0         1         1             11.0         2.0        24.0             13
#> 4         Sm  127         0         1          0         0         1         1            103.0        24.0       230.0            127
#> 5        VLa 3730         0         1          0         0         1         1           3500.0       230.0      7230.0           3730
#> 6        Flo 3800         1         1          0         0         1         1           3800.0         0.0      7600.0           3800
#> 7       Dorm    0         0         0          0         0         1         1              0.5        -0.5         0.5              0
#>   sizebin_width              comments
#> 1             0          Dormant seed
#> 2           200              Seedling
#> 3            22 Very small vegetative
#> 4           206      Small vegetative
#> 5          7000 Very large vegetative
#> 6          7600             Flowering
#> 7             1               Dormant
```
  
  <p style="text-indent: 20px">A further important point has to do with the meaning of 0 sizes. In most cases, a size of 0 will mean that the individual is alive but unobservable. However, a size of 0 may have different meanings in other cases, such as when the size metric used is a logarithm and so values of 0 and lower are possible in observable individuals in the dataset. The dataset should be explored carefully for these situations, particularly in the case of the creation of a function-based matrix, such as an IPM, in order to make sure that the stageframe accurately describes and matches the values actually occurring in the dataset.</p>
  
  <p style="text-indent: 20px">Once the stageframe is created, we can reorganize the dataset into vertical format. Here, 'vertical' format is a way of organizing demogaphic data in which each row corresponds to the state of a single individual in two (if ahistorical) or three (if historical) consecutive time intervals. To handle this, we use the `verticalize3()` function, which creates historically-formatted vertical datasets, as below.</p>

```r
lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988, patchidcol = "SUBPLOT", 
                         individcol = "GENET", blocksize = 9, juvcol = "Seedling1988", 
                         size1col = "Volume88", repstr1col = "FCODE88", 
                         fec1col = "Intactseed88", dead1col = "Dead1988", 
                         nonobs1col = "Dormant1988", stageassign = lathframe, 
                         stagesize = "sizea", censorcol = "Missing1988", 
                         censorkeep = NA, censor = TRUE)
summary(lathvert)
#>      rowid         popid         patchid    individ           year2          sizea1           sizea2           sizea3      
#>  Min.   :   1.0   Mode:logical   1:670   Min.   :  1.00   Min.   :1988   Min.   :   3.4   Min.   :   1.8   Min.   :   2.1  
#>  1st Qu.: 237.0   NA's:2515      2:201   1st Qu.: 39.00   1st Qu.:1988   1st Qu.:  45.0   1st Qu.:  15.6   1st Qu.:  21.0  
#>  Median : 523.0                  3:552   Median : 78.00   Median :1989   Median : 439.5   Median : 127.5   Median : 131.2  
#>  Mean   : 537.5                  4:469   Mean   : 92.82   Mean   :1989   Mean   : 699.5   Mean   : 508.6   Mean   : 466.9  
#>  3rd Qu.: 819.5                  5:305   3rd Qu.:139.50   3rd Qu.:1990   3rd Qu.:1025.5   3rd Qu.: 732.5   3rd Qu.: 696.0  
#>  Max.   :1118.0                  6:318   Max.   :276.00   Max.   :1990   Max.   :7032.0   Max.   :7032.0   Max.   :6645.8  
#>                                                                          NA's   :1125     NA's   :126      NA's   :406     
#>     repstra1         repstra2         repstra3          feca1            feca2            feca3          size1added       size2added    
#>  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   : 0.000   Min.   : 0.000   Min.   : 0.000   Min.   :   0.0   Min.   :   0.0  
#>  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.: 0.000   1st Qu.: 0.000   1st Qu.: 0.000   1st Qu.:   0.0   1st Qu.:  12.6  
#>  Median :0.0000   Median :0.0000   Median :0.0000   Median : 0.000   Median : 0.000   Median : 2.000   Median :   9.0   Median : 106.5  
#>  Mean   :0.3259   Mean   :0.2508   Mean   :0.2614   Mean   : 5.472   Mean   : 4.802   Mean   : 5.969   Mean   : 386.6   Mean   : 483.1  
#>  3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.: 7.000   3rd Qu.: 6.000   3rd Qu.: 9.000   3rd Qu.: 621.7   3rd Qu.: 732.5  
#>  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :66.000   Max.   :66.000   Max.   :66.000   Max.   :7032.0   Max.   :7032.0  
#>  NA's   :1125     NA's   :127      NA's   :407      NA's   :2061     NA's   :1915     NA's   :1969                                      
#>    size3added       fec1added         fec2added        fec3added        obsstatus1       obsstatus2       obsstatus3       repstatus1    
#>  Min.   :   0.0   Min.   : 0.0000   Min.   : 0.000   Min.   : 0.000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
#>  1st Qu.:  10.0   1st Qu.: 0.0000   1st Qu.: 0.000   1st Qu.: 0.000   1st Qu.:0.0000   1st Qu.:1.0000   1st Qu.:1.0000   1st Qu.:0.0000  
#>  Median :  73.2   Median : 0.0000   Median : 0.000   Median : 0.000   Median :1.0000   Median :1.0000   Median :1.0000   Median :0.0000  
#>  Mean   : 391.5   Mean   : 0.9877   Mean   : 1.146   Mean   : 1.296   Mean   :0.5527   Mean   :0.9499   Mean   :0.8386   Mean   :0.1801  
#>  3rd Qu.: 545.0   3rd Qu.: 0.0000   3rd Qu.: 0.000   3rd Qu.: 0.000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:0.0000  
#>  Max.   :6645.8   Max.   :66.0000   Max.   :66.000   Max.   :66.000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
#>                                                                                                                                          
#>    repstatus2       repstatus3       fecstatus1        fecstatus2       fecstatus3       firstseen       lastseen        alive1      
#>  Min.   :0.0000   Min.   :0.0000   Min.   :0.00000   Min.   :0.0000   Min.   :0.0000   Min.   :1988   Min.   :1988   Min.   :0.0000  
#>  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:1988   1st Qu.:1991   1st Qu.:0.0000  
#>  Median :0.0000   Median :0.0000   Median :0.00000   Median :0.0000   Median :0.0000   Median :1988   Median :1991   Median :1.0000  
#>  Mean   :0.2382   Mean   :0.2191   Mean   :0.08946   Mean   :0.1058   Mean   :0.1121   Mean   :1988   Mean   :1991   Mean   :0.5813  
#>  3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:1988   3rd Qu.:1991   3rd Qu.:1.0000  
#>  Max.   :1.0000   Max.   :1.0000   Max.   :1.00000   Max.   :1.0000   Max.   :1.0000   Max.   :1990   Max.   :1991   Max.   :1.0000  
#>                                                                                                                                      
#>      alive2      alive3           obsage        obslifespan       stage1             stage2             stage3            matstatus1    
#>  Min.   :1   Min.   :0.0000   Min.   :0.0000   Min.   :0.000   Length:2515        Length:2515        Length:2515        Min.   :0.0000  
#>  1st Qu.:1   1st Qu.:1.0000   1st Qu.:0.0000   1st Qu.:2.000   Class :character   Class :character   Class :character   1st Qu.:1.0000  
#>  Median :1   Median :1.0000   Median :1.0000   Median :3.000   Mode  :character   Mode  :character   Mode  :character   Median :1.0000  
#>  Mean   :1   Mean   :0.8887   Mean   :0.8258   Mean   :2.449                                                            Mean   :0.9368  
#>  3rd Qu.:1   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:3.000                                                            3rd Qu.:1.0000  
#>  Max.   :1   Max.   :1.0000   Max.   :2.0000   Max.   :3.000                                                            Max.   :1.0000  
#>                                                                                                                                         
#>    matstatus2       matstatus3
#>  Min.   :0.0000   Min.   :1   
#>  1st Qu.:1.0000   1st Qu.:1   
#>  Median :1.0000   Median :1   
#>  Mean   :0.8883   Mean   :1   
#>  3rd Qu.:1.0000   3rd Qu.:1   
#>  Max.   :1.0000   Max.   :1   
#> 
dim(lathvert)
#> [1] 2515   42
```
  
  <p style="text-indent: 20px">It generally pays to explore a dataset thoroughly prior to full analysis. One way to do this is to look at summaries of the dataset after it has been subset to specific cases of interest. For example, after subsetting to only cases in which individuals are vegetatively dormant in time *t*, we can see that this particular stage is actually designated not by a size of 0 but with NA entries. We can also see, looking at the stage designations and at the dimensions of the aforementioned subset, that vegetative dormancy is reasonably plentiful in this dataset.</p>

```r
summary(lathvert[which(lathvert$stage2 == "Dorm"),])
#>      rowid         popid         patchid    individ           year2          sizea1           sizea2        sizea3          repstra1     
#>  Min.   :   1.0   Mode:logical   1:26    Min.   :  1.00   Min.   :1989   Min.   :   3.4   Min.   : NA   Min.   :   4.0   Min.   :0.0000  
#>  1st Qu.: 309.8   NA's:126       2:13    1st Qu.: 28.25   1st Qu.:1989   1st Qu.: 110.5   1st Qu.: NA   1st Qu.:  54.4   1st Qu.:0.0000  
#>  Median : 502.0                  3:38    Median : 66.50   Median :1989   Median : 580.0   Median : NA   Median : 123.0   Median :0.0000  
#>  Mean   : 529.7                  4:13    Mean   : 74.79   Mean   :1989   Mean   : 646.9   Mean   :NaN   Mean   : 285.8   Mean   :0.3248  
#>  3rd Qu.: 863.2                  5:28    3rd Qu.:113.75   3rd Qu.:1990   3rd Qu.: 795.8   3rd Qu.: NA   3rd Qu.: 385.2   3rd Qu.:1.0000  
#>  Max.   :1079.0                  6: 8    Max.   :234.00   Max.   :1990   Max.   :2941.0   Max.   : NA   Max.   :2637.6   Max.   :1.0000  
#>                                                                          NA's   :9        NA's   :126   NA's   :9        NA's   :9       
#>     repstra2      repstra3          feca1            feca2         feca3       size1added        size2added   size3added     
#>  Min.   : NA   Min.   :0.0000   Min.   : 0.000   Min.   : NA   Min.   : 0    Min.   :   0.00   Min.   :0    Min.   :   0.00  
#>  1st Qu.: NA   1st Qu.:0.0000   1st Qu.: 0.000   1st Qu.: NA   1st Qu.: 0    1st Qu.:  66.45   1st Qu.:0    1st Qu.:  32.02  
#>  Median : NA   Median :0.0000   Median : 0.000   Median : NA   Median : 0    Median : 481.60   Median :0    Median : 109.80  
#>  Mean   :NaN   Mean   :0.1111   Mean   : 4.526   Mean   :NaN   Mean   : 3    Mean   : 600.68   Mean   :0    Mean   : 265.41  
#>  3rd Qu.: NA   3rd Qu.:0.0000   3rd Qu.: 5.750   3rd Qu.: NA   3rd Qu.: 4    3rd Qu.: 732.50   3rd Qu.:0    3rd Qu.: 360.05  
#>  Max.   : NA   Max.   :1.0000   Max.   :30.000   Max.   : NA   Max.   :17    Max.   :2941.00   Max.   :0    Max.   :2637.60  
#>  NA's   :126   NA's   :9        NA's   :88       NA's   :126   NA's   :113                                                   
#>    fec1added        fec2added   fec3added         obsstatus1       obsstatus2   obsstatus3       repstatus1       repstatus2
#>  Min.   : 0.000   Min.   :0   Min.   : 0.0000   Min.   :0.0000   Min.   :0    Min.   :0.0000   Min.   :0.0000   Min.   :0   
#>  1st Qu.: 0.000   1st Qu.:0   1st Qu.: 0.0000   1st Qu.:1.0000   1st Qu.:0    1st Qu.:1.0000   1st Qu.:0.0000   1st Qu.:0   
#>  Median : 0.000   Median :0   Median : 0.0000   Median :1.0000   Median :0    Median :1.0000   Median :0.0000   Median :0   
#>  Mean   : 1.365   Mean   :0   Mean   : 0.3095   Mean   :0.9286   Mean   :0    Mean   :0.9286   Mean   :0.3016   Mean   :0   
#>  3rd Qu.: 0.000   3rd Qu.:0   3rd Qu.: 0.0000   3rd Qu.:1.0000   3rd Qu.:0    3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:0   
#>  Max.   :30.000   Max.   :0   Max.   :17.0000   Max.   :1.0000   Max.   :0    Max.   :1.0000   Max.   :1.0000   Max.   :0   
#>                                                                                                                             
#>    repstatus3       fecstatus1      fecstatus2   fecstatus3        firstseen       lastseen        alive1      alive2      alive3 
#>  Min.   :0.0000   Min.   :0.000   Min.   :0    Min.   :0.00000   Min.   :1988   Min.   :1990   Min.   :1   Min.   :1   Min.   :1  
#>  1st Qu.:0.0000   1st Qu.:0.000   1st Qu.:0    1st Qu.:0.00000   1st Qu.:1988   1st Qu.:1991   1st Qu.:1   1st Qu.:1   1st Qu.:1  
#>  Median :0.0000   Median :0.000   Median :0    Median :0.00000   Median :1988   Median :1991   Median :1   Median :1   Median :1  
#>  Mean   :0.1032   Mean   :0.119   Mean   :0    Mean   :0.03175   Mean   :1988   Mean   :1991   Mean   :1   Mean   :1   Mean   :1  
#>  3rd Qu.:0.0000   3rd Qu.:0.000   3rd Qu.:0    3rd Qu.:0.00000   3rd Qu.:1988   3rd Qu.:1991   3rd Qu.:1   3rd Qu.:1   3rd Qu.:1  
#>  Max.   :1.0000   Max.   :1.000   Max.   :0    Max.   :1.00000   Max.   :1989   Max.   :1991   Max.   :1   Max.   :1   Max.   :1  
#>                                                                                                                                   
#>      obsage       obslifespan       stage1             stage2             stage3            matstatus1       matstatus2   matstatus3
#>  Min.   :1.000   Min.   :2.000   Length:126         Length:126         Length:126         Min.   :0.0000   Min.   :1    Min.   :1   
#>  1st Qu.:1.000   1st Qu.:3.000   Class :character   Class :character   Class :character   1st Qu.:1.0000   1st Qu.:1    1st Qu.:1   
#>  Median :1.000   Median :3.000   Mode  :character   Mode  :character   Mode  :character   Median :1.0000   Median :1    Median :1   
#>  Mean   :1.357   Mean   :2.794                                                            Mean   :0.8968   Mean   :1    Mean   :1   
#>  3rd Qu.:2.000   3rd Qu.:3.000                                                            3rd Qu.:1.0000   3rd Qu.:1    3rd Qu.:1   
#>  Max.   :2.000   Max.   :3.000                                                            Max.   :1.0000   Max.   :1    Max.   :1   
#> 
dim(lathvert[which(lathvert$stage2 == "Dorm"),])
#> [1] 126  42
```

The fact that dormant individuals have been assigned NA for size does not affect us in this exercise, because the construction of raw matrices depends on the stage designations rather than on the size classifications. However, if we wished to build function-derived matrices, as in Analyses 2 and 3, then we would need to change these NAs to 0 or another appropriate number.

  <p style="text-indent: 20px">Voilá! Now we can move on to creating a few extra objects that will help us estimate full population projection matrices.</p>
  
#### Step 2. Provide supplemental information for matrix estimation

  <p style="text-indent: 20px">For our next step, we need to create a **reproductive matrix**. This matrix shows which stages are reproductive, which stages they lead to the reproduction of, and at what level reproduction occurs. 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 `lefko3`. This matrix has rows and columns equal to 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, this matrix looks like a nearly empty population matrix, but notes the per-individual mean modifiers on fecundity for each stage that actually reproduces.</p>
  
  <p style="text-indent: 20px">Here, we first create a 0 matrix with dimensionality equal to the number of rows in `lathframe`. Then we modify elements corresponding to fecundity by the expected mean seed dormancy probability (row 1), and by the germination rate (row 2).</p>

```r
lathrepm <- matrix(0, 7, 7)
lathrepm[1, 6] <- 0.345
lathrepm[2, 6] <- 0.054

lathrepm
#>      [,1] [,2] [,3] [,4] [,5]  [,6] [,7]
#> [1,]    0    0    0    0    0 0.345    0
#> [2,]    0    0    0    0    0 0.054    0
#> [3,]    0    0    0    0    0 0.000    0
#> [4,]    0    0    0    0    0 0.000    0
#> [5,]    0    0    0    0    0 0.000    0
#> [6,]    0    0    0    0    0 0.000    0
#> [7,]    0    0    0    0    0 0.000    0
```

  <p style="text-indent: 20px">Next we will provide some given transitions, via an **overwrite table**. 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 a germinated seedling, respectively. Let's start with the ahistorical case.</p>

```r
lathover2 <- overwrite(stage3 = c("Sd", "Sdl"), stage2 = c("Sd", "Sd"), 
                       givenrate = c(0.345, 0.054))
lathover2
#>   stage3 stage2 stage1 eststage3 eststage2 eststage1 givenrate convtype
#> 1     Sd     Sd     NA        NA        NA        NA     0.345        1
#> 2    Sdl     Sd     NA        NA        NA        NA     0.054        1
```

And now the historical case. Here we need to show the stages in time step *t*-1 for this to work properly. Note the use of the `"rep"` designation for `Stage1` - this is shorthand telling R to use all reproductive stages for the time interval. Here, there is only one reproductive stage, but in other cases, such as in an IPM, this shorthand notation can represent many stages in a single statement.

```r
lathover3 <- overwrite(stage3 = c("Sd", "Sd", "Sdl"), stage2 = c("Sd", "Sd", "Sd"), 
                       stage1 = c("Sd", "rep", "rep"), givenrate = c(0.345, 0.345, 0.054))
lathover3
#>   stage3 stage2 stage1 eststage3 eststage2 eststage1 givenrate convtype
#> 1     Sd     Sd     Sd        NA        NA        NA     0.345        1
#> 2     Sd     Sd    rep        NA        NA        NA     0.345        1
#> 3    Sdl     Sd    rep        NA        NA        NA     0.054        1
```

  <p style="text-indent: 20px">These overwrite tables show us that we have survival-transition probabilities (`convtype = 1`), that the given transitions originate from the dormany seed stage (Sd) in time *t* (and reproductive stages in time *t*-1 in the historical case), and the specific values to be used in overwriting: 0.345 and 0.054. If we wished, we could have used the values of transitions to be estimated within this matrix as proxies for these values, in which case the `eststageX` columns would have entries and the `givenrate` column would be blank.</p>

  <p style="text-indent: 20px">Now we are read to create some matrices!</p>
  
#### Step 3. Estimate matrices

  <p style="text-indent: 20px">Now let's try to create the raw Lefkovitch matrices starting with the ahistorical case, based on those shown in Ehrlén [-@ehrlen_dynamics_2000]. That study shows a mean matrix covering years 1989 and 1990 as time *t*. We will utilize the entire dataset instead, covering 1988 to 1991, as follows.</p>

```r
ehrlen2 <- rlefko2(data = lathvert, stageframe = lathframe, year = "all", 
                   stages = c("stage3", "stage2"), repmatrix = lathrepm, 
                   overwrite = lathover2, yearcol = "year2", 
                   indivcol = "individ")
ehrlen2
#> $A
#> $A[[1]]
#>       [,1]       [,2]    [,3]       [,4]      [,5]       [,6] [,7]
#> [1,] 0.345 0.00000000 0.00000 0.00000000 0.0000000 1.02941848    0
#> [2,] 0.054 0.00000000 0.00000 0.00000000 0.0000000 0.16112637    0
#> [3,] 0.000 0.61855670 0.78125 0.07692308 0.0000000 0.00000000    0
#> [4,] 0.000 0.03092784 0.12500 0.67307692 0.2044199 0.07024793    0
#> [5,] 0.000 0.00000000 0.00000 0.03205128 0.3535912 0.21487603    0
#> [6,] 0.000 0.00000000 0.00000 0.03846154 0.2928177 0.61983471    0
#> [7,] 0.000 0.10309278 0.03125 0.12179487 0.1160221 0.09090909    0
#> 
#> $A[[2]]
#>       [,1]       [,2]        [,3]       [,4]       [,5]       [,6]       [,7]
#> [1,] 0.345 0.00000000 0.000000000 0.00000000 0.00000000 2.83876712 0.00000000
#> [2,] 0.054 0.00000000 0.000000000 0.00000000 0.00000000 0.44432877 0.00000000
#> [3,] 0.000 0.70085470 0.784615385 0.17592593 0.02205882 0.00000000 0.15068493
#> [4,] 0.000 0.02564103 0.069230769 0.61574074 0.32352941 0.18721461 0.50684932
#> [5,] 0.000 0.00000000 0.000000000 0.02314815 0.29411765 0.28767123 0.09589041
#> [6,] 0.000 0.00000000 0.000000000 0.02314815 0.21323529 0.42922374 0.12328767
#> [7,] 0.000 0.02564103 0.007692308 0.06481481 0.07352941 0.07305936 0.12328767
#> 
#> $A[[3]]
#>       [,1]      [,2]        [,3]       [,4]       [,5]       [,6]      [,7]
#> [1,] 0.345 0.0000000 0.000000000 0.00000000 0.00000000 0.86750000 0.0000000
#> [2,] 0.054 0.0000000 0.000000000 0.00000000 0.00000000 0.13578261 0.0000000
#> [3,] 0.000 0.6716418 0.606299213 0.03533569 0.00000000 0.00000000 0.1132075
#> [4,] 0.000 0.0000000 0.173228346 0.45229682 0.05785124 0.04347826 0.3584906
#> [5,] 0.000 0.0000000 0.003937008 0.17667845 0.28925620 0.19565217 0.4528302
#> [6,] 0.000 0.0000000 0.000000000 0.15901060 0.52892562 0.66666667 0.0754717
#> [7,] 0.000 0.0000000 0.000000000 0.00000000 0.00000000 0.00000000 0.0000000
#> 
#> 
#> $U
#> $U[[1]]
#>       [,1]       [,2]    [,3]       [,4]      [,5]       [,6] [,7]
#> [1,] 0.345 0.00000000 0.00000 0.00000000 0.0000000 0.00000000    0
#> [2,] 0.054 0.00000000 0.00000 0.00000000 0.0000000 0.00000000    0
#> [3,] 0.000 0.61855670 0.78125 0.07692308 0.0000000 0.00000000    0
#> [4,] 0.000 0.03092784 0.12500 0.67307692 0.2044199 0.07024793    0
#> [5,] 0.000 0.00000000 0.00000 0.03205128 0.3535912 0.21487603    0
#> [6,] 0.000 0.00000000 0.00000 0.03846154 0.2928177 0.61983471    0
#> [7,] 0.000 0.10309278 0.03125 0.12179487 0.1160221 0.09090909    0
#> 
#> $U[[2]]
#>       [,1]       [,2]        [,3]       [,4]       [,5]       [,6]       [,7]
#> [1,] 0.345 0.00000000 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [2,] 0.054 0.00000000 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [3,] 0.000 0.70085470 0.784615385 0.17592593 0.02205882 0.00000000 0.15068493
#> [4,] 0.000 0.02564103 0.069230769 0.61574074 0.32352941 0.18721461 0.50684932
#> [5,] 0.000 0.00000000 0.000000000 0.02314815 0.29411765 0.28767123 0.09589041
#> [6,] 0.000 0.00000000 0.000000000 0.02314815 0.21323529 0.42922374 0.12328767
#> [7,] 0.000 0.02564103 0.007692308 0.06481481 0.07352941 0.07305936 0.12328767
#> 
#> $U[[3]]
#>       [,1]      [,2]        [,3]       [,4]       [,5]       [,6]      [,7]
#> [1,] 0.345 0.0000000 0.000000000 0.00000000 0.00000000 0.00000000 0.0000000
#> [2,] 0.054 0.0000000 0.000000000 0.00000000 0.00000000 0.00000000 0.0000000
#> [3,] 0.000 0.6716418 0.606299213 0.03533569 0.00000000 0.00000000 0.1132075
#> [4,] 0.000 0.0000000 0.173228346 0.45229682 0.05785124 0.04347826 0.3584906
#> [5,] 0.000 0.0000000 0.003937008 0.17667845 0.28925620 0.19565217 0.4528302
#> [6,] 0.000 0.0000000 0.000000000 0.15901060 0.52892562 0.66666667 0.0754717
#> [7,] 0.000 0.0000000 0.000000000 0.00000000 0.00000000 0.00000000 0.0000000
#> 
#> 
#> $F
#> $F[[1]]
#>      [,1] [,2] [,3] [,4] [,5]      [,6] [,7]
#> [1,]    0    0    0    0    0 1.0294185    0
#> [2,]    0    0    0    0    0 0.1611264    0
#> [3,]    0    0    0    0    0 0.0000000    0
#> [4,]    0    0    0    0    0 0.0000000    0
#> [5,]    0    0    0    0    0 0.0000000    0
#> [6,]    0    0    0    0    0 0.0000000    0
#> [7,]    0    0    0    0    0 0.0000000    0
#> 
#> $F[[2]]
#>      [,1] [,2] [,3] [,4] [,5]      [,6] [,7]
#> [1,]    0    0    0    0    0 2.8387671    0
#> [2,]    0    0    0    0    0 0.4443288    0
#> [3,]    0    0    0    0    0 0.0000000    0
#> [4,]    0    0    0    0    0 0.0000000    0
#> [5,]    0    0    0    0    0 0.0000000    0
#> [6,]    0    0    0    0    0 0.0000000    0
#> [7,]    0    0    0    0    0 0.0000000    0
#> 
#> $F[[3]]
#>      [,1] [,2] [,3] [,4] [,5]      [,6] [,7]
#> [1,]    0    0    0    0    0 0.8675000    0
#> [2,]    0    0    0    0    0 0.1357826    0
#> [3,]    0    0    0    0    0 0.0000000    0
#> [4,]    0    0    0    0    0 0.0000000    0
#> [5,]    0    0    0    0    0 0.0000000    0
#> [6,]    0    0    0    0    0 0.0000000    0
#> [7,]    0    0    0    0    0 0.0000000    0
#> 
#> 
#> $hstages
#> NULL
#> 
#> $ahstages
#>   new_stage_id orig_stage_id original_size bin_size_ctr bin_size_min bin_size_max repstatus obsstatus propstatus immstatus matstatus
#> 1            1            Sd             0            0          0.0          0.0         0         0          1         1         0
#> 2            2           Sdl           100          100          0.0        200.0         0         1          0         1         0
#> 3            3           VSm            13           13          2.0         24.0         0         1          0         0         1
#> 4            4            Sm           127          127         24.0        230.0         0         1          0         0         1
#> 5            5           VLa          3730         3730        230.0       7230.0         0         1          0         0         1
#> 6            6           Flo          3800         3800          0.0       7600.0         1         1          0         0         1
#> 7            7          Dorm             0            0         -0.5          0.5         0         0          0         0         1
#>   indataset bin_size_width bin_raw_halfwidth alive              comments stageno
#> 1         0              0               0.0     1          Dormant seed       1
#> 2         1            200             100.0     1              Seedling       2
#> 3         1             22              11.0     1 Very small vegetative       3
#> 4         1            206             103.0     1      Small vegetative       4
#> 5         1           7000            3500.0     1 Very large vegetative       5
#> 6         1           7600            3800.0     1             Flowering       6
#> 7         1              1               0.5     1               Dormant       7
#> 
#> $labels
#>   pop patch year2
#> 1  NA    NA  1988
#> 2  NA    NA  1989
#> 3  NA    NA  1990
#> 
#> $matrixqc
#> [1] 68  6  3
#> 
#> $dataqc
#> [1]  276 2515
#> 
#> attr(,"class")
#> [1] "lefkoMat"
```

The input for the `rlefko2()` function includes `year = "all"`, which can be changed to `year = c(1989, 1990)` to focus just on years 1989 and 1990, as in the paper, or `year = 1989` to focus exclusively on the transition from 1989 to 1990 (the year entered is the year in time *t*). 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.

  <p style="text-indent: 20px">The output from this analysis includes is a `lefkoMat` object, which is a list object with the following elements:</p>
  
  **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
  
  
  <p style="text-indent: 20px">We can understand `lefkoMat` objects through the summary `summary` function.</p>

```r
summary(ehrlen2)
#> 
#> This lefkoMat object contains 3 matrices.
#> 
#> Each matrix is a square matrix with 7 rows and columns, and a total of 49 elements.
#> A total of 68 survival-transitions were estimated, with 22.6666666666667 per matrix.
#> A total of 6 fecundity transitions were estimated, with 2 per matrix.
#> 
#> The dataset contains a total of 276 unique individuals and 2515 unique transitions.
#> NULL
```

  <p style="text-indent: 20px">We start off learning that 3 matrices were estimated, and 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 size of the dataset. It is typical for population ecologists to consider the total number of transitions in a dataset, 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, among other problems. So, this output can be very helpful to understand the degree to which estimated matrices might be overparameterized.</p>

  <p style="text-indent: 20px">Now let's look at historical matrices. Because of the size of historical matrices, we will only show the summary.</p>

```r
ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = c(1989, 1990), 
                   stages = c("stage3", "stage2", "stage1"), repmatrix = lathrepm, 
                   overwrite = lathover3, yearcol = "year2", 
                   indivcol = "individ")

summary(ehrlen3)
#> 
#> This lefkoMat object contains 2 matrices.
#> 
#> Each matrix is a square matrix with 49 rows and columns, and a total of 2401 elements.
#> A total of 131 survival-transitions were estimated, with 65.5 per matrix.
#> A total of 14 fecundity transitions were estimated, with 7 per matrix.
#> 
#> The dataset contains a total of 276 unique individuals and 2515 unique transitions.
#> NULL
```

Since we are estimating raw matrices and three consecutive time steps of data are needed to estimate each transition, there is one less matrix estimated in the historical case than in the ahistorical case. These matrices are also quite a bit bigger than ahistorical matrices, and so require more parameters to be estimated. However, they are also composed mostly of structural 0s, which works against overparameterization and helps us to make more realistic matrices.

  <p style="text-indent: 20px">We can see the impact of structural zeroes by eliminating some of them in the process of matrix estimation. The easy way to do that is to set `reduce = TRUE` within the `rlefko3()` call, which eliminates stage pairs in which both column and row are zero vectors. This will end up giving us matrices with 19 fewer rows and columns.</p>

```r
ehrlen3red <- rlefko3(data = lathvert, stageframe = lathframe, year = c(1989, 1990), 
                   stages = c("stage3", "stage2", "stage1"), repmatrix = lathrepm, 
                   overwrite = lathover3, yearcol = "year2", 
                   indivcol = "individ", reduce = TRUE)

summary(ehrlen3red)
#> 
#> This lefkoMat object contains 2 matrices.
#> 
#> Each matrix is a square matrix with 30 rows and columns, and a total of 900 elements.
#> A total of 131 survival-transitions were estimated, with 65.5 per matrix.
#> A total of 14 fecundity transitions were estimated, with 7 per matrix.
#> 
#> The dataset contains a total of 276 unique individuals and 2515 unique transitions.
#> NULL
```

  <p style="text-indent: 20px">Next we will create a mean ahistorical matrix. 

```r
ehrlen2mean <- lmean(ehrlen2)
ehrlen2mean
#> $A
#> $A[[1]]
#>       [,1]       [,2]        [,3]       [,4]        [,5]       [,6]       [,7]
#> [1,] 0.345 0.00000000 0.000000000 0.00000000 0.000000000 1.57856187 0.00000000
#> [2,] 0.054 0.00000000 0.000000000 0.00000000 0.000000000 0.24707925 0.00000000
#> [3,] 0.000 0.66368440 0.724054866 0.09606156 0.007352941 0.00000000 0.08796416
#> [4,] 0.000 0.01885629 0.122486372 0.58037149 0.195266847 0.10031360 0.28844663
#> [5,] 0.000 0.00000000 0.001312336 0.07729263 0.312321669 0.23273315 0.18290687
#> [6,] 0.000 0.00000000 0.000000000 0.07354010 0.344992865 0.57190837 0.06625312
#> [7,] 0.000 0.04291127 0.012980769 0.06220323 0.063183837 0.05465615 0.04109589
#> 
#> $A[[2]]
#>       [,1]       [,2]        [,3]       [,4]        [,5]       [,6]       [,7]
#> [1,] 0.345 0.00000000 0.000000000 0.00000000 0.000000000 1.57856187 0.00000000
#> [2,] 0.054 0.00000000 0.000000000 0.00000000 0.000000000 0.24707925 0.00000000
#> [3,] 0.000 0.66368440 0.724054866 0.09606156 0.007352941 0.00000000 0.08796416
#> [4,] 0.000 0.01885629 0.122486372 0.58037149 0.195266847 0.10031360 0.28844663
#> [5,] 0.000 0.00000000 0.001312336 0.07729263 0.312321669 0.23273315 0.18290687
#> [6,] 0.000 0.00000000 0.000000000 0.07354010 0.344992865 0.57190837 0.06625312
#> [7,] 0.000 0.04291127 0.012980769 0.06220323 0.063183837 0.05465615 0.04109589
#> 
#> 
#> $U
#> $U[[1]]
#>       [,1]       [,2]        [,3]       [,4]        [,5]       [,6]       [,7]
#> [1,] 0.345 0.00000000 0.000000000 0.00000000 0.000000000 0.00000000 0.00000000
#> [2,] 0.054 0.00000000 0.000000000 0.00000000 0.000000000 0.00000000 0.00000000
#> [3,] 0.000 0.66368440 0.724054866 0.09606156 0.007352941 0.00000000 0.08796416
#> [4,] 0.000 0.01885629 0.122486372 0.58037149 0.195266847 0.10031360 0.28844663
#> [5,] 0.000 0.00000000 0.001312336 0.07729263 0.312321669 0.23273315 0.18290687
#> [6,] 0.000 0.00000000 0.000000000 0.07354010 0.344992865 0.57190837 0.06625312
#> [7,] 0.000 0.04291127 0.012980769 0.06220323 0.063183837 0.05465615 0.04109589
#> 
#> $U[[2]]
#>       [,1]       [,2]        [,3]       [,4]        [,5]       [,6]       [,7]
#> [1,] 0.345 0.00000000 0.000000000 0.00000000 0.000000000 0.00000000 0.00000000
#> [2,] 0.054 0.00000000 0.000000000 0.00000000 0.000000000 0.00000000 0.00000000
#> [3,] 0.000 0.66368440 0.724054866 0.09606156 0.007352941 0.00000000 0.08796416
#> [4,] 0.000 0.01885629 0.122486372 0.58037149 0.195266847 0.10031360 0.28844663
#> [5,] 0.000 0.00000000 0.001312336 0.07729263 0.312321669 0.23273315 0.18290687
#> [6,] 0.000 0.00000000 0.000000000 0.07354010 0.344992865 0.57190837 0.06625312
#> [7,] 0.000 0.04291127 0.012980769 0.06220323 0.063183837 0.05465615 0.04109589
#> 
#> 
#> $F
#> $F[[1]]
#>      [,1] [,2] [,3] [,4] [,5]      [,6] [,7]
#> [1,]    0    0    0    0    0 1.5785619    0
#> [2,]    0    0    0    0    0 0.2470792    0
#> [3,]    0    0    0    0    0 0.0000000    0
#> [4,]    0    0    0    0    0 0.0000000    0
#> [5,]    0    0    0    0    0 0.0000000    0
#> [6,]    0    0    0    0    0 0.0000000    0
#> [7,]    0    0    0    0    0 0.0000000    0
#> 
#> $F[[2]]
#>      [,1] [,2] [,3] [,4] [,5]      [,6] [,7]
#> [1,]    0    0    0    0    0 1.5785619    0
#> [2,]    0    0    0    0    0 0.2470792    0
#> [3,]    0    0    0    0    0 0.0000000    0
#> [4,]    0    0    0    0    0 0.0000000    0
#> [5,]    0    0    0    0    0 0.0000000    0
#> [6,]    0    0    0    0    0 0.0000000    0
#> [7,]    0    0    0    0    0 0.0000000    0
#> 
#> 
#> $hstages
#> NULL
#> 
#> $ahstages
#>   new_stage_id orig_stage_id original_size bin_size_ctr bin_size_min bin_size_max repstatus obsstatus propstatus immstatus matstatus
#> 1            1            Sd             0            0          0.0          0.0         0         0          1         1         0
#> 2            2           Sdl           100          100          0.0        200.0         0         1          0         1         0
#> 3            3           VSm            13           13          2.0         24.0         0         1          0         0         1
#> 4            4            Sm           127          127         24.0        230.0         0         1          0         0         1
#> 5            5           VLa          3730         3730        230.0       7230.0         0         1          0         0         1
#> 6            6           Flo          3800         3800          0.0       7600.0         1         1          0         0         1
#> 7            7          Dorm             0            0         -0.5          0.5         0         0          0         0         1
#>   indataset bin_size_width bin_raw_halfwidth alive              comments stageno
#> 1         0              0               0.0     1          Dormant seed       1
#> 2         1            200             100.0     1              Seedling       2
#> 3         1             22              11.0     1 Very small vegetative       3
#> 4         1            206             103.0     1      Small vegetative       4
#> 5         1           7000            3500.0     1 Very large vegetative       5
#> 6         1           7600            3800.0     1             Flowering       6
#> 7         1              1               0.5     1               Dormant       7
#> 
#> $labels
#>   pop patch
#> 1   1     1
#> 2   1   All
#> 
#> $matrixqc
#> [1] 56  4  2
#> 
#> $dataqc
#> [1]  276 2515
#> 
#> attr(,"class")
#> [1] "lefkoMat"
```
  
  <p style="text-indent: 20px">First of all, let's take a look at the structure of the above output. Function `lmean()` creates a `lefkoMat` object, just as `rlefko()` does, and so we have the main composite mean matrix (shown in element `A`), as well as the mean survival-transition matrix (`U`) and the mean fecundity matrix (`F`), followed by a section outlining the definitions and order of historical paired stages (`hstages`, shown here as `NA` because the matrices are ahistorical), a section outlining the actual stages as outlined in the `stageframe` object used to create these matrices (`ahstages`), a section outlining the definitions and order of the matrices (`labels`), and then two quality control sections used in output for the `summary()` function (`matrixqc` and `dataqc`). So, all necessary information is retained for ease of use, and to focus in on just the main mean matrix itself, we can call the first matrix of object `A`.</p>

```r
ehrlen2mean$A[[1]]
#>       [,1]       [,2]        [,3]       [,4]        [,5]       [,6]       [,7]
#> [1,] 0.345 0.00000000 0.000000000 0.00000000 0.000000000 1.57856187 0.00000000
#> [2,] 0.054 0.00000000 0.000000000 0.00000000 0.000000000 0.24707925 0.00000000
#> [3,] 0.000 0.66368440 0.724054866 0.09606156 0.007352941 0.00000000 0.08796416
#> [4,] 0.000 0.01885629 0.122486372 0.58037149 0.195266847 0.10031360 0.28844663
#> [5,] 0.000 0.00000000 0.001312336 0.07729263 0.312321669 0.23273315 0.18290687
#> [6,] 0.000 0.00000000 0.000000000 0.07354010 0.344992865 0.57190837 0.06625312
#> [7,] 0.000 0.04291127 0.012980769 0.06220323 0.063183837 0.05465615 0.04109589
```

  <p style="text-indent: 20px">Matrix designations make an impact on what we see in the output to `lmean`. Here, we see two matrices, but they are actually completely equal, because the second matrix is a grand mean of patch-level means, and there is only one patch-level mean matrix. We would see additional matrices if we had split our data by population and/or patch and asked for means for each. To see the exact specifications for each matrix, we can look at the `labels` component of our mean `lefkoMat` object.</p>

```r
ehrlen2mean$labels
#>   pop patch
#> 1   1     1
#> 2   1   All
```

The designations A and 1 are the default symbols used for population and patch, respectively, when no designations are supplied. This scenario is generally the case when no such division is made in the dataset.

  <p style="text-indent: 20px">We may wish to error check by assessing the survival of ahistorical stages. We can try the following, which shows us the survival as the column sum for the main survival-transition matrix.</p>

```r
colSums(ehrlen2mean$U[[1]])
#> [1] 0.3990000 0.7254520 0.8608343 0.8894690 0.9231182 0.9596113 0.6666667
```

All of the values are within the realm of possibility for probabilities, so this matrix appears to be fine.

  <p style="text-indent: 20px">And now the historical mean matrix. We will use the reduced matrix to simplify later analyses, and show only the top-left corner of the rather large matrix.</p>

```r
ehrlen3mean <- lmean(ehrlen3red)
ehrlen3mean$A[[1]][1:20,1:8]
#>        [,1] [,2]       [,3]      [,4] [,5]       [,6]      [,7] [,8]
#>  [1,] 0.345    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
#>  [2,] 0.000    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
#>  [3,] 0.000    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
#>  [4,] 0.000    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
#>  [5,] 0.000    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
#>  [6,] 0.000    0 0.68796296 0.0000000    0 0.79373737 0.0000000    0
#>  [7,] 0.000    0 0.03302469 0.0000000    0 0.07555556 0.0000000    0
#>  [8,] 0.000    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
#>  [9,] 0.000    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
#> [10,] 0.000    0 0.00000000 0.3333333    0 0.00000000 0.5416667    0
#> [11,] 0.000    0 0.00000000 0.4166667    0 0.00000000 0.4583333    0
#> [12,] 0.000    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
#> [13,] 0.000    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
#> [14,] 0.000    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
#> [15,] 0.000    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
#> [16,] 0.000    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
#> [17,] 0.000    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
#> [18,] 0.000    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
#> [19,] 0.000    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
#> [20,] 0.000    0 0.00000000 0.0000000    0 0.00000000 0.0000000    0
```

Do not fear the prevalence of 0's in this matrix - this is normal, both because most 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.

  <p style="text-indent: 20px">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.</p>

```r
ehrlen3mean$hstages
#>    stcod3 stcod2 stage3 stage2
#> 1      Sd     Sd      1      1
#> 2     Sdl     Sd      2      1
#> 11    VSm    Sdl      3      2
#> 12     Sm    Sdl      4      2
#> 15   Dorm    Sdl      7      2
#> 19    VSm    VSm      3      3
#> 20     Sm    VSm      4      3
#> 21    VLa    VSm      5      3
#> 23   Dorm    VSm      7      3
#> 27    VSm     Sm      3      4
#> 28     Sm     Sm      4      4
#> 29    VLa     Sm      5      4
#> 30    Flo     Sm      6      4
#> 31   Dorm     Sm      7      4
#> 35    VSm    VLa      3      5
#> 36     Sm    VLa      4      5
#> 37    VLa    VLa      5      5
#> 38    Flo    VLa      6      5
#> 39   Dorm    VLa      7      5
#> 41     Sd    Flo      1      6
#> 42    Sdl    Flo      2      6
#> 44     Sm    Flo      4      6
#> 45    VLa    Flo      5      6
#> 46    Flo    Flo      6      6
#> 47   Dorm    Flo      7      6
#> 51    VSm   Dorm      3      7
#> 52     Sm   Dorm      4      7
#> 53    VLa   Dorm      5      7
#> 54    Flo   Dorm      6      7
#> 55   Dorm   Dorm      7      7
```

  <p style="text-indent: 20px">There are 30 pairs of ahistorical stages. 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 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 dormant seed at 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 seedling in time *t*-1 and very small adult 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*!</p>

  <p style="text-indent: 20px">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.</p>

```r
colSums(ehrlen3mean$U[[1]])
#>  [1] 0.3450000 0.0000000 0.7209877 0.7500000 1.0000000 0.8692929 1.0000000 0.0000000 1.0000000 0.8793860 0.8652597 1.0000000 0.9000000
#> [14] 1.0000000 0.5000000 0.9381313 0.9052419 0.9027326 1.0000000 0.3990000 0.0000000 0.9045911 0.9229692 0.9560545 1.0000000 0.4500000
#> [27] 0.3243243 0.3571429 0.4444444 0.5000000
```

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. [-@shefferson_life_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.

```r
summary(colSums(ehrlen3mean$U[[1]]))
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.0000  0.4458  0.8743  0.6945  0.9516  1.0000
```

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

  <p style="text-indent: 20px">We may wish to conduct all sorts of analyses with these matrices, and some are possible with `lefko3` itself. Let's estimate the deterministic population growth rate in each case via eigen analysis. 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. Note that each lambda estimate also includes a data frame describing the matrices in order (this is the `labels` object within the output list). Here is the set of ahistorical annual lambda estimates.</p>

```r
lambda3(ehrlen2)
#>   pop patch year2    lambda
#> 1  NA    NA  1988 0.8952585
#> 2  NA    NA  1989 0.9069721
#> 3  NA    NA  1990 0.9297366
```

Here is the lambda associated with the mean matrix.

```r
lambda3(ehrlen2mean)
#>   pop patch    lambda
#> 1   1     1 0.9369727
#> 2   1   All 0.9369727
```

Readers may be surprised to see that lambda for the mean matrix is on the high side relative to the annual matrices. This is most likely a result of the fact that some elements are 0's due simply to a lack of individuals transitioning through them. You may consider this a likely artefact of low sample size.

  <p style="text-indent: 20px">And we will now look at the same numbers for the historical analyses. First the annual matrices.</p>

```r
lambda3(ehrlen3red)
#>   pop patch year2    lambda
#> 1  NA    NA  1989 0.8859389
#> 2  NA    NA  1990 0.9045258
```

Now the mean matrix.

```r
lambda3(ehrlen3mean)
#>   pop patch    lambda
#> 1   1     1 0.8847749
#> 2   1   All 0.8847749
```

  <p style="text-indent: 20px">Readers will likely observe both that there are fewer lambda estimates in the historical case, and that the mean lambda is lower. Both are expected. First, because there are 4 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_longitudinal_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_life_2014]. While we do not argue that the drop in lambda is due to this specific trade-off, and it is certainly possible for historical matrix analysis to lead to higher lambda estimates, we do argue that this lambda is likely to be more realistic than the higher lambda estimated in the ahistorical case.</p>

  <p style="text-indent: 20px">We can also take a peek at the stable stage distributions, as follows for the ahistorical case. Our numbers here match those produced by package `popbio`'s `stable.stage()` function.</p>

```r
stablestage3(ehrlen2mean)
#>    matrix new_stage_id orig_stage_id original_size    ss_prop
#> 1       1            1            Sd             0 0.30247427
#> 2       1            2           Sdl           100 0.04734382
#> 3       1            3           VSm            13 0.24631667
#> 4       1            4            Sm           127 0.18488072
#> 5       1            5           VLa          3730 0.07469276
#> 6       1            6           Flo          3800 0.11343017
#> 7       1            7          Dorm             0 0.03086158
#> 8       2            1            Sd             0 0.30247427
#> 9       2            2           Sdl           100 0.04734382
#> 10      2            3           VSm            13 0.24631667
#> 11      2            4            Sm           127 0.18488072
#> 12      2            5           VLa          3730 0.07469276
#> 13      2            6           Flo          3800 0.11343017
#> 14      2            7          Dorm             0 0.03086158
```

  <p style="text-indent: 20px">The data frame output shows us the stages themselves, which matrix they refer to, and the stable stage distribution (in the `ss_prop` column). We can do this for the historical case as well. Because historical output for the `stablestage3()` function is a list with two data frames, let's take a look at each of these data frames. The first will be the historical output.</p>

```r
ehrlen3mss <- stablestage3(ehrlen3mean)
ehrlen3mss$hist[,c(1,2,3,5)]
#>    matrix orig_stage_id_2 orig_stage_id_1 new_stage_id_1
#> 1       1              Sd              Sd              1
#> 2       1             Sdl              Sd              1
#> 3       1             VSm             Sdl              2
#> 4       1              Sm             Sdl              2
#> 5       1            Dorm             Sdl              2
#> 6       1             VSm             VSm              3
#> 7       1              Sm             VSm              3
#> 8       1             VLa             VSm              3
#> 9       1            Dorm             VSm              3
#> 10      1             VSm              Sm              4
#> 11      1              Sm              Sm              4
#> 12      1             VLa              Sm              4
#> 13      1             Flo              Sm              4
#> 14      1            Dorm              Sm              4
#> 15      1             VSm             VLa              5
#> 16      1              Sm             VLa              5
#> 17      1             VLa             VLa              5
#> 18      1             Flo             VLa              5
#> 19      1            Dorm             VLa              5
#> 20      1              Sd             Flo              6
#> 21      1             Sdl             Flo              6
#> 22      1              Sm             Flo              6
#> 23      1             VLa             Flo              6
#> 24      1             Flo             Flo              6
#> 25      1            Dorm             Flo              6
#> 26      1             VSm            Dorm              7
#> 27      1              Sm            Dorm              7
#> 28      1             VLa            Dorm              7
#> 29      1             Flo            Dorm              7
#> 30      1            Dorm            Dorm              7
#> 31      2              Sd              Sd              1
#> 32      2             Sdl              Sd              1
#> 33      2             VSm             Sdl              2
#> 34      2              Sm             Sdl              2
#> 35      2            Dorm             Sdl              2
#> 36      2             VSm             VSm              3
#> 37      2              Sm             VSm              3
#> 38      2             VLa             VSm              3
#> 39      2            Dorm             VSm              3
#> 40      2             VSm              Sm              4
#> 41      2              Sm              Sm              4
#> 42      2             VLa              Sm              4
#> 43      2             Flo              Sm              4
#> 44      2            Dorm              Sm              4
#> 45      2             VSm             VLa              5
#> 46      2              Sm             VLa              5
#> 47      2             VLa             VLa              5
#> 48      2             Flo             VLa              5
#> 49      2            Dorm             VLa              5
#> 50      2              Sd             Flo              6
#> 51      2             Sdl             Flo              6
#> 52      2              Sm             Flo              6
#> 53      2             VLa             Flo              6
#> 54      2             Flo             Flo              6
#> 55      2            Dorm             Flo              6
#> 56      2             VSm            Dorm              7
#> 57      2              Sm            Dorm              7
#> 58      2             VLa            Dorm              7
#> 59      2             Flo            Dorm              7
#> 60      2            Dorm            Dorm              7
```

This is hard to read because it is by stage pair. It may make more sense if we look at the `$ahist` portion, which shows us the stable stage distribtuion according to the original stages given in our stageframe.

```r
ehrlen3mss$ahist
#>    matrix new_stage_id    ss_prop
#> 1       1            1 0.33108399
#> 2       1            2 0.04394258
#> 3       1            3 0.23880523
#> 4       1            4 0.19396392
#> 5       1            5 0.07119439
#> 6       1            6 0.10503696
#> 7       1            7 0.01597293
#> 8       2            1 0.33108399
#> 9       2            2 0.04394258
#> 10      2            3 0.23880523
#> 11      2            4 0.19396392
#> 12      2            5 0.07119439
#> 13      2            6 0.10503696
#> 14      2            7 0.01597293
```

  <p style="text-indent: 20px">Notice that this stable stage distribution is different from the ahistorical case, suggesting a strong influence of individual history.</p>
  
  <p style="text-indent: 20px">Let's take a look at the reproductive values now, in similar order to the stable stage distribution case.</p>

```r
repvalue3(ehrlen2mean)
#>    matrix new_stage_id orig_stage_id original_size left_vector rep_value
#> 1       1            1            Sd             0  -0.0165976   1.00000
#> 2       1            2           Sdl           100  -0.1819506  10.96246
#> 3       1            3           VSm            13  -0.2278100  13.72548
#> 4       1            4            Sm           127  -0.3596400  21.66819
#> 5       1            5           VLa          3730  -0.5107549  30.77282
#> 6       1            6           Flo          3800  -0.6629890  39.94487
#> 7       1            7          Dorm             0  -0.2914706  17.56101
#> 8       2            1            Sd             0  -0.0165976   1.00000
#> 9       2            2           Sdl           100  -0.1819506  10.96246
#> 10      2            3           VSm            13  -0.2278100  13.72548
#> 11      2            4            Sm           127  -0.3596400  21.66819
#> 12      2            5           VLa          3730  -0.5107549  30.77282
#> 13      2            6           Flo          3800  -0.6629890  39.94487
#> 14      2            7          Dorm             0  -0.2914706  17.56101
```

These values are also the same as those produced by `popbio`'s `reproductive.value()` function. But there are differences when we look at the historical case. Note that `popbio` typically fails in the historical case, as that package is not generally made to handle extremely large matrices.

```r
ehrlen3mrv <- repvalue3(ehrlen3mean)
ehrlen3mrv$hist[,c(1,2,3,5,6)]
#>    matrix orig_stage_id_2 orig_stage_id_1 new_stage_id_1 left_vector
#> 1       1              Sd              Sd              1   0.0000000
#> 2       1             Sdl              Sd              1   0.0000000
#> 3       1             VSm             Sdl              2   0.1585647
#> 4       1              Sm             Sdl              2   0.1738053
#> 5       1            Dorm             Sdl              2   0.1070537
#> 6       1             VSm             VSm              3   0.1927766
#> 7       1              Sm             VSm              3   0.2322782
#> 8       1             VLa             VSm              3   0.0000000
#> 9       1            Dorm             VSm              3   0.1022135
#> 10      1             VSm              Sm              4   0.2077572
#> 11      1              Sm              Sm              4   0.2028627
#> 12      1             VLa              Sm              4   0.2672135
#> 13      1             Flo              Sm              4   0.3022418
#> 14      1            Dorm              Sm              4   0.1074098
#> 15      1             VSm             VLa              5   0.0875093
#> 16      1              Sm             VLa              5   0.2461279
#> 17      1             VLa             VLa              5   0.2701883
#> 18      1             Flo             VLa              5   0.2931642
#> 19      1            Dorm             VLa              5   0.1144807
#> 20      1              Sd             Flo              6   0.0000000
#> 21      1             Sdl             Flo              6   0.0000000
#> 22      1              Sm             Flo              6   0.2498206
#> 23      1             VLa             Flo              6   0.2852715
#> 24      1             Flo             Flo              6   0.3365592
#> 25      1            Dorm             Flo              6   0.1260929
#> 26      1             VSm            Dorm              7   0.1025115
#> 27      1              Sm            Dorm              7   0.0904359
#> 28      1             VLa            Dorm              7   0.1108299
#> 29      1             Flo            Dorm              7   0.1571749
#> 30      1            Dorm            Dorm              7   0.0639748
#> 31      2              Sd              Sd              1   0.0000000
#> 32      2             Sdl              Sd              1   0.0000000
#> 33      2             VSm             Sdl              2   0.1585647
#> 34      2              Sm             Sdl              2   0.1738053
#> 35      2            Dorm             Sdl              2   0.1070537
#> 36      2             VSm             VSm              3   0.1927766
#> 37      2              Sm             VSm              3   0.2322782
#> 38      2             VLa             VSm              3   0.0000000
#> 39      2            Dorm             VSm              3   0.1022135
#> 40      2             VSm              Sm              4   0.2077572
#> 41      2              Sm              Sm              4   0.2028627
#> 42      2             VLa              Sm              4   0.2672135
#> 43      2             Flo              Sm              4   0.3022418
#> 44      2            Dorm              Sm              4   0.1074098
#> 45      2             VSm             VLa              5   0.0875093
#> 46      2              Sm             VLa              5   0.2461279
#> 47      2             VLa             VLa              5   0.2701883
#> 48      2             Flo             VLa              5   0.2931642
#> 49      2            Dorm             VLa              5   0.1144807
#> 50      2              Sd             Flo              6   0.0000000
#> 51      2             Sdl             Flo              6   0.0000000
#> 52      2              Sm             Flo              6   0.2498206
#> 53      2             VLa             Flo              6   0.2852715
#> 54      2             Flo             Flo              6   0.3365592
#> 55      2            Dorm             Flo              6   0.1260929
#> 56      2             VSm            Dorm              7   0.1025115
#> 57      2              Sm            Dorm              7   0.0904359
#> 58      2             VLa            Dorm              7   0.1108299
#> 59      2             Flo            Dorm              7   0.1571749
#> 60      2            Dorm            Dorm              7   0.0639748
```

  <p style="text-indent: 20px">The final column is the reproductive value of the stage pair. Once again, this is hard to interpret, so we can take a look at the ahistorical summary.</p>

```r
ehrlen3mrv$ahist
#>    matrix new_stage_id rep_value
#> 1       1            1 0.0000000
#> 2       1            2 0.0000000
#> 3       1            3 1.0000000
#> 4       1            4 1.0787534
#> 5       1            5 1.3549138
#> 6       1            6 1.6366109
#> 7       1            7 0.5780928
#> 8       2            1 0.0000000
#> 9       2            2 0.0000000
#> 10      2            3 1.0000000
#> 11      2            4 1.0787534
#> 12      2            5 1.3549138
#> 13      2            6 1.6366109
#> 14      2            7 0.5780928
```

Notice that the reproductive values are quite different here, once again suggesting a strong importance to individual history.

  <p style="text-indent: 20px">Further analytical tools are being planned for `lefko3`, but all packages that handle projection matrices can handle the individual matrices produced and saved in `lefkoMat` objects in this package.
  
#### Other issues

  <p style="text-indent: 20px">The `lmean()` function gives the option of choosing the geometric or arithmetic mean, and also of treating 0's slightly differently depending on the situation. In principle, we advocate taking the geometric mean across time because survival and reproduction are geometric processes. Indeed, projecting an arithmetic mean matrix forward will lead to a higher projected population size over time than projecting the matrices themselves would, while the geometric mean does not suffer from this upward bias. However, a number of issues do arise when geometric mean matrices are estimated, and these need to be considered prior to a decision as to which method to take.</p>
  
  <p style="text-indent: 20px">In datasets sparse enough to yield many constituent matrices with 0's caused by individuals being entirely but temporarily absent from a particular stage, the use of the geometric mean may result in a downward bias since a single 0 will yield a geometric mean of 0, regardless of the other values being averaged. Here, the dataset is sufficiently large and the stage definitions are sufficiently broad that we generally avoid this problem. Conversely, it is possible that some mean probabilities estimated through arithmetic rather geometric approaches will end up greater than 1.0, leading to problems with analysis. Though the literature does not exist on this, one of us (RPS) suspects that relationships with time might account for the tendency of matrix projection models to overestimate lambda and future population size [@crone_ability_2013].</p>
  
  <p style="text-indent: 20px">Here, we create a geometric mean matrix using the `sparse` option.</p>

```r
ehrlen2geomean <- lmean(ehrlen2, time = "geometric", sparse = TRUE)
ehrlen2geomean$A[[1]]
#>       [,1]       [,2]        [,3]       [,4]       [,5]       [,6]       [,7]
#> [1,] 0.345 0.00000000 0.000000000 0.00000000 0.00000000 1.36352698 0.00000000
#> [2,] 0.054 0.00000000 0.000000000 0.00000000 0.00000000 0.21342162 0.00000000
#> [3,] 0.000 0.66279884 0.718970853 0.07819879 0.02205882 0.00000000 0.13060885
#> [4,] 0.000 0.02816064 0.114448308 0.57230654 0.15640472 0.08300071 0.42626365
#> [5,] 0.000 0.00000000 0.003937008 0.05079821 0.31100573 0.22953894 0.20837964
#> [6,] 0.000 0.00000000 0.000000000 0.05211823 0.32083642 0.56185311 0.09646103
#> [7,] 0.000 0.05141405 0.015504342 0.08884882 0.09236361 0.08149699 0.12328767
```
  
  <p style="text-indent: 20px">The `sparse = TRUE` option is vital to these matrices. This option tells R how to interpret 0s occurring within matrices. It is quite typical for raw matrices to include 0s in some elements simply because individuals did not exist in some stages within a particular year. In such cases, a 0 is a misrepresentation of that transition, because it reflects a lack of data (including low sample size) rather than a genuine transition impossibility. In the extreme case, treating a 0 just like another transition may make the mean element 0 as well when `time = "geometric"`. Thus, setting `sparse = TRUE` tells R to ignore 0's except when they are structural, and structural 0s are those 0s that are 0 in every matrix. This will be particularly useful in large matrices such as historical matrices, but we can still see the effects by comparing the last mean matrix above to one in which the `sparse` option is set to `FALSE` (the default), as below.</p>

```r
ehrlen2geomeanns <- lmean(ehrlen2, time = "geometric")
ehrlen2geomeanns$A[[1]]
#>       [,1]      [,2]      [,3]       [,4]      [,5]       [,6] [,7]
#> [1,] 0.345 0.0000000 0.0000000 0.00000000 0.0000000 1.36352698    0
#> [2,] 0.054 0.0000000 0.0000000 0.00000000 0.0000000 0.21342162    0
#> [3,] 0.000 0.6627988 0.7189709 0.07819879 0.0000000 0.00000000    0
#> [4,] 0.000 0.0000000 0.1144483 0.57230654 0.1564047 0.08300071    0
#> [5,] 0.000 0.0000000 0.0000000 0.05079821 0.3110057 0.22953894    0
#> [6,] 0.000 0.0000000 0.0000000 0.05211823 0.3208364 0.56185311    0
#> [7,] 0.000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000    0
```

  <p style="text-indent: 20px">With `sparse = FALSE`, which is the default, 13 transitions dropped to 0 in the mean matrix. This includes all of the transitions to and from vegetative dormancy, the transition from seedling to small vegetative, the transition from very small vegetative to large vegetative, and a number of others with typically lower probabilities due to fewer individuals experiencing them. We will use the former mean matrix, with `sparse = TRUE`. Note, however, that the sparseness of data used in most situations leading to the construction of raw matrices will likely lead to the use of arithmetic means.

  <p style="text-indent: 20px">As a comparison, let's now check the survival probabilities estimated using the geometric approach, and compare with the original arithmetic mean estimates as well as with the `sparse` option off.

```r
colSums(ehrlen2geomean$U[[1]])
#> [1] 0.3990000 0.7423735 0.8528605 0.8422706 0.9026693 0.9558898 0.9850008
colSums(ehrlen2mean$U[[1]])
#> [1] 0.3990000 0.7254520 0.8608343 0.8894690 0.9231182 0.9596113 0.6666667
colSums(ehrlen2geomeanns$U[[1]])
#> [1] 0.3990000 0.6627988 0.8334192 0.7534218 0.7882469 0.8743928 0.0000000
```

  <p style="text-indent: 20px">This is an interesting comparison. The lowest estimates come from the true geometric mean, which turns all elements to 0 if a single 0 exists. The highest is the arithmetic approach, where all survival estimates are higher except for the dormant seed stage (equal everywhere), and vegetative dormancy (highest in the geometric sparse case).</p>
  
  <p style="text-indent: 20px">Let's compare lambda in each case. First, the geometric sparse case</p>

```r
lambda3(ehrlen2geomean)
#>   pop patch    lambda
#> 1   1     1 0.9250929
#> 2   1   All 0.9250929
```

Now the arithmetic case.

```r
lambda3(ehrlen2mean)
#>   pop patch    lambda
#> 1   1     1 0.9369727
#> 2   1   All 0.9369727
```

Finally, the true geometric case.

```r
lambda3(ehrlen2geomeanns)
#>   pop patch    lambda
#> 1   1     1 0.8542391
#> 2   1   All 0.8542391
```

The true geometric mean has the lowest lambda, while the highest is the true arithmetic case.


### Analysis 2: Function-derived matrix estimation

  <p style="text-indent: 20px">In this analysis, we will build a set of function-derived matrices using the *Lathyrus* dataset. To spice things up, we will use a slightly different approach to size classification, using the natural log sizes instead of the normal sizes. This has to do with the fact that volume is used as the size metric here, and so should have an allometric relationship to some vital rates (note that all size metrics have allometric relationships by default, but this is clearer when size is based on something more strongly related to mass, as is volume). We will also categorize individuals as reproductive vs. non-reproductive.</p>

#### Step 1. Data characterization and reorganization

  <p style="text-indent: 20px">First, we will create a stageframe for the dataset. For this purpose, let's look back at a summary of the main dataset, focusing at the distribution of log sizes.</p>

```r
summary(c(lathyrus$Volume88, lathyrus$Volume89, lathyrus$Volume90, lathyrus$Volume91))
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>     1.8    14.7   123.0   484.2   732.5  7032.0    1248
summary(c(lathyrus$lnVol88, lathyrus$lnVol89, lathyrus$lnVol90, lathyrus$lnVol91))
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>   0.600   2.700   4.800   4.777   6.600   8.900    1248
```
  
  <p style="text-indent: 20px">It is important to note the size minima and maxima, because we have been using 0 as the size of vegetatively dormant individuals. What we note here is that apart from the many NAs that occur in the dataset, the smallest log sizes are still above 0, meaning that we are still able to use the value 0 as an indicator of vegetative dormancy. This would not be the case if the smallest size were negative, as might happen when the volume (not log volume) is between 0 and 1.</p>
  
  <p style="text-indent: 20px">It can also help to take a look at plots of these distributions. First we will look at the raw volume data.</p>

```r
plot(density(c(lathyrus$Volume88, lathyrus$Volume89, lathyrus$Volume90, lathyrus$Volume91), 
             na.rm = TRUE), main = "Volume")
```

![plot of chunk Ch2an2st1ln383](Ch2an2st1ln383-1.png)

Next we will look at the log volume data.

```r
plot(density(c(lathyrus$lnVol88, lathyrus$lnVol89, lathyrus$lnVol90, lathyrus$lnVol91), 
             na.rm = TRUE), main = "Log volume")
```

![plot of chunk Ch2an2st1ln389](Ch2an2st1ln389-1.png)

  <p style="text-indent: 20px">We see two very differently shaped distributions here. The log volume distribution looks 'better' than the volume distribution, in the sense that it is closer to some semblance of normality. This is helpful since our size metrics are decimals, and so cannot be treated as integers (in the latter case, we could try modeling them as either Poisson- or negative binomial-distributed). We will work with log volume in this example, and treat it as Gaussian-distributed.</p>
  
  <p style="text-indent: 20px">We need to cover all log volumes actually occurring in the dataset, because our approach will include developing estimates of vital rates given inputs including all possible sizes, reproductive status, and so forth. We build this by creating vectors of the values describing each stage, always in the same order.</p>

```r
sizevector <- c(0, 4.6, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9)
stagevector <- c("Sd", "Sdl", "Dorm", "Sz1nr", "Sz2nr", "Sz3nr", "Sz4nr", "Sz5nr",
                 "Sz6nr", "Sz7nr", "Sz8nr", "Sz9nr", "Sz1r", "Sz2r", "Sz3r", "Sz4r",
                 "Sz5r", "Sz6r", "Sz7r", "Sz8r", "Sz9r")
repvector <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1)
obsvector <- c(0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
matvector <- c(0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
immvector <- c(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
indataset <- c(0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
binvec <- c(0, 4.6, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 
            0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5)

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

  <p style="text-indent: 20px">To be useful to others reading your work, or to yourself in the future, it helps to add text descriptions of the stages. Here, we also add some descriptions to this stageframe as comments.</p>

```r
lathframeln$comments <- c("Dormant seed", "Seedling", "Dormant", "Size 1 Non-reprod", "Size 2 Non-reprod",
                        "Size 3 Non-reprod", "Size 4 Non-reprod", "Size 5 Non-reprod", "Size 6 Non-reprod",
                        "Size 7 Non-reprod", "Size 8 Non-reprod", "Size 9 Non-reprod", "Size 1 Reprod",
                        "Size 2 Reprod", "Size 3 Reprod", "Size 4 Reprod", "Size 5 Reprod", "Size 6 Reprod",
                        "Size 7 Reprod", "Size 8 Reprod", "Size 9 Reprod")
lathframeln
#>    stagenames size repstatus obsstatus propstatus immstatus matstatus indataset binhalfwidth_raw sizebin_min sizebin_max sizebin_center
#> 1          Sd  0.0         0         0          1         1         0         0              0.0         0.0         0.0            0.0
#> 2         Sdl  4.6         0         1          0         1         0         1              4.6         0.0         9.2            4.6
#> 3        Dorm  0.0         0         0          0         0         1         1              0.5        -0.5         0.5            0.0
#> 4       Sz1nr  1.0         0         1          0         0         1         1              0.5         0.5         1.5            1.0
#> 5       Sz2nr  2.0         0         1          0         0         1         1              0.5         1.5         2.5            2.0
#> 6       Sz3nr  3.0         0         1          0         0         1         1              0.5         2.5         3.5            3.0
#> 7       Sz4nr  4.0         0         1          0         0         1         1              0.5         3.5         4.5            4.0
#> 8       Sz5nr  5.0         0         1          0         0         1         1              0.5         4.5         5.5            5.0
#> 9       Sz6nr  6.0         0         1          0         0         1         1              0.5         5.5         6.5            6.0
#> 10      Sz7nr  7.0         0         1          0         0         1         1              0.5         6.5         7.5            7.0
#> 11      Sz8nr  8.0         0         1          0         0         1         1              0.5         7.5         8.5            8.0
#> 12      Sz9nr  9.0         0         1          0         0         1         1              0.5         8.5         9.5            9.0
#> 13       Sz1r  1.0         1         1          0         0         1         1              0.5         0.5         1.5            1.0
#> 14       Sz2r  2.0         1         1          0         0         1         1              0.5         1.5         2.5            2.0
#> 15       Sz3r  3.0         1         1          0         0         1         1              0.5         2.5         3.5            3.0
#> 16       Sz4r  4.0         1         1          0         0         1         1              0.5         3.5         4.5            4.0
#> 17       Sz5r  5.0         1         1          0         0         1         1              0.5         4.5         5.5            5.0
#> 18       Sz6r  6.0         1         1          0         0         1         1              0.5         5.5         6.5            6.0
#> 19       Sz7r  7.0         1         1          0         0         1         1              0.5         6.5         7.5            7.0
#> 20       Sz8r  8.0         1         1          0         0         1         1              0.5         7.5         8.5            8.0
#> 21       Sz9r  9.0         1         1          0         0         1         1              0.5         8.5         9.5            9.0
#>    sizebin_width          comments
#> 1            0.0      Dormant seed
#> 2            9.2          Seedling
#> 3            1.0           Dormant
#> 4            1.0 Size 1 Non-reprod
#> 5            1.0 Size 2 Non-reprod
#> 6            1.0 Size 3 Non-reprod
#> 7            1.0 Size 4 Non-reprod
#> 8            1.0 Size 5 Non-reprod
#> 9            1.0 Size 6 Non-reprod
#> 10           1.0 Size 7 Non-reprod
#> 11           1.0 Size 8 Non-reprod
#> 12           1.0 Size 9 Non-reprod
#> 13           1.0     Size 1 Reprod
#> 14           1.0     Size 2 Reprod
#> 15           1.0     Size 3 Reprod
#> 16           1.0     Size 4 Reprod
#> 17           1.0     Size 5 Reprod
#> 18           1.0     Size 6 Reprod
#> 19           1.0     Size 7 Reprod
#> 20           1.0     Size 8 Reprod
#> 21           1.0     Size 9 Reprod
```

  <p style="text-indent: 20px">Once the stageframe is created, we can reorganize the dataset into vertical format. Here, 'vertical' format is a way of organizing demographic data in which each row corresponds to the state of a single individual in two (if ahistorical) or three (if historical) consecutive time steps. To handle this, we use the `verticalize3` function, which creates historically-formatted vertical datasets, as below. Please note that we need to get rid of NAs for `modelsearch` to work properly when we actually build models of vital rates, so we will now use the `NAas0 = TRUE` option.</p>

```r
lathvertln <- verticalize3(lathyrus, noyears = 4, firstyear = 1988, patchidcol = "SUBPLOT", 
                         individcol = "GENET", blocksize = 9, juvcol = "Seedling1988", 
                         size1col = "lnVol88", repstr1col = "FCODE88", 
                         fec1col = "Intactseed88", dead1col = "Dead1988", 
                         nonobs1col = "Dormant1988", stageassign = lathframeln, 
                         stagesize = "sizea", censorcol = "Missing1988", 
                         censorkeep = NA, NAas0 = TRUE, censor = TRUE)
summary(lathvertln)
#>      rowid         popid         patchid    individ           year2          sizea1          sizea2          sizea3         repstra1     
#>  Min.   :   1.0   Mode:logical   1:670   Min.   :  1.00   Min.   :1988   Min.   :0.000   Min.   :0.000   Min.   :0.000   Min.   :0.0000  
#>  1st Qu.: 237.0   NA's:2515      2:201   1st Qu.: 39.00   1st Qu.:1988   1st Qu.:0.000   1st Qu.:2.500   1st Qu.:2.300   1st Qu.:0.0000  
#>  Median : 523.0                  3:552   Median : 78.00   Median :1989   Median :2.200   Median :4.700   Median :4.300   Median :0.0000  
#>  Mean   : 537.5                  4:469   Mean   : 92.82   Mean   :1989   Mean   :2.944   Mean   :4.587   Mean   :4.062   Mean   :0.1801  
#>  3rd Qu.: 819.5                  5:305   3rd Qu.:139.50   3rd Qu.:1990   3rd Qu.:6.400   3rd Qu.:6.600   3rd Qu.:6.300   3rd Qu.:0.0000  
#>  Max.   :1118.0                  6:318   Max.   :276.00   Max.   :1990   Max.   :8.900   Max.   :8.900   Max.   :8.800   Max.   :1.0000  
#>     repstra2         repstra3          feca1             feca2            feca3          obsstatus1       obsstatus2       obsstatus3    
#>  Min.   :0.0000   Min.   :0.0000   Min.   : 0.0000   Min.   : 0.000   Min.   : 0.000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
#>  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.: 0.0000   1st Qu.: 0.000   1st Qu.: 0.000   1st Qu.:0.0000   1st Qu.:1.0000   1st Qu.:1.0000  
#>  Median :0.0000   Median :0.0000   Median : 0.0000   Median : 0.000   Median : 0.000   Median :1.0000   Median :1.0000   Median :1.0000  
#>  Mean   :0.2382   Mean   :0.2191   Mean   : 0.9877   Mean   : 1.146   Mean   : 1.296   Mean   :0.5527   Mean   :0.9499   Mean   :0.8386  
#>  3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.: 0.0000   3rd Qu.: 0.000   3rd Qu.: 0.000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1.0000  
#>  Max.   :1.0000   Max.   :1.0000   Max.   :66.0000   Max.   :66.000   Max.   :66.000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000  
#>    repstatus1       repstatus2       repstatus3       fecstatus1        fecstatus2       fecstatus3       firstseen       lastseen   
#>  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00000   Min.   :0.0000   Min.   :0.0000   Min.   :1988   Min.   :1988  
#>  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:1988   1st Qu.:1991  
#>  Median :0.0000   Median :0.0000   Median :0.0000   Median :0.00000   Median :0.0000   Median :0.0000   Median :1988   Median :1991  
#>  Mean   :0.1801   Mean   :0.2382   Mean   :0.2191   Mean   :0.08946   Mean   :0.1058   Mean   :0.1121   Mean   :1988   Mean   :1991  
#>  3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:1988   3rd Qu.:1991  
#>  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00000   Max.   :1.0000   Max.   :1.0000   Max.   :1990   Max.   :1991  
#>      alive1           alive2      alive3           obsage        obslifespan       stage1             stage2             stage3         
#>  Min.   :0.0000   Min.   :1   Min.   :0.0000   Min.   :0.0000   Min.   :0.000   Length:2515        Length:2515        Length:2515       
#>  1st Qu.:0.0000   1st Qu.:1   1st Qu.:1.0000   1st Qu.:0.0000   1st Qu.:2.000   Class :character   Class :character   Class :character  
#>  Median :1.0000   Median :1   Median :1.0000   Median :1.0000   Median :3.000   Mode  :character   Mode  :character   Mode  :character  
#>  Mean   :0.5813   Mean   :1   Mean   :0.8887   Mean   :0.8258   Mean   :2.449                                                           
#>  3rd Qu.:1.0000   3rd Qu.:1   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:3.000                                                           
#>  Max.   :1.0000   Max.   :1   Max.   :1.0000   Max.   :2.0000   Max.   :3.000                                                           
#>    matstatus1       matstatus2       matstatus3
#>  Min.   :0.0000   Min.   :0.0000   Min.   :1   
#>  1st Qu.:1.0000   1st Qu.:1.0000   1st Qu.:1   
#>  Median :1.0000   Median :1.0000   Median :1   
#>  Mean   :0.9368   Mean   :0.8883   Mean   :1   
#>  3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1   
#>  Max.   :1.0000   Max.   :1.0000   Max.   :1
dim(lathvertln)
#> [1] 2515   36
```
  
  <p style="text-indent: 20px">Something to note here is that fecundity is typically but not always an integer. We can see this below, where we show each non-integer in the dataset. These are not errors, but due to the fact that fecundity was estimated in two ways in this dataset, one way leading to integers, and one leading to decimals.</p>

```r
lathvertln$feca2[which(lathvertln$feca2 != round(lathvertln$feca2))]
#> [1] 12.500000  3.818182 26.600000 23.333333  5.333333 24.500000
length(lathvertln$feca2)
#> [1] 2515
```
  
In this case, we can either treat fecundity as a continuous variable, or round the values and treat them as count variables. The choice of which approach to take will have major repercussions on the analysis, and may severely alter projected population dynamics. In this dataset, fecundity is mostly counts of intact seeds, and only differs in a few cases where the seed output was estimated based on other models. Indeed, only 6 of 2515 entries are not integers! So, we will round fecundity so that we can assume count variables in the analysis, as below.

```r
lathvertln$feca2 <- round(lathvertln$feca2)
lathvertln$feca1 <- round(lathvertln$feca1)
lathvertln$feca3 <- round(lathvertln$feca3)
```

  <p style="text-indent: 20px">All set! Now to move on to supplemental descriptive information.</p>
  
#### Step 2. Provide supplemental information for matrix estimation

  <p style="text-indent: 20px">Next we need to create a **reproductive matrix**, which desribes which stages are reproductive, which stages they lead to the reproduction of, and at what level they reproduce. 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 `lefko3`. This matrix has rows and columns equal to the number of stages described in the stageframe for this dataset, and the rows and columns refer to these stages in the exact 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 `lathframeln`. Then we modify elements corresponding to fecundity by the expected mean seed dormancy probability (row 1), and by the germination rate (row 2).</p>

```r
lathrepmln <- matrix(0, 21, 21)
lathrepmln[1, c(13:21)] <- 0.345
lathrepmln[2, c(13:21)] <- 0.054
```

  <p style="text-indent: 20px">Next we need to provide some given transitions. However, we will use the same overwrite tables as in the last example (object `lathover`), because these relate only to stages that are the same in both analyses. Here they are again.</p>

```r
lathover2 <- overwrite(stage3 = c("Sd", "Sdl"), stage2 = c("Sd", "Sd"), 
                      givenrate = c(0.345, 0.054))

lathover3 <- overwrite(stage3 = c("Sd", "Sd", "Sdl"), stage2 = c("Sd", "Sd", "Sd"), 
                      stage1 = c("Sd", "rep", "rep"), givenrate = c(0.345, 0.345, 0.054))
```
  
  <p style="text-indent: 20px">Next, we move to a new step: modeling vital rates.</p>

#### Step 3. Develop models of demographic parameters

  <p style="text-indent: 20px">Here, we will run the `modelsearch` function with the new vertical dataset, first for the ahistorical model set and then for the historical model set. This function will develop our models for us, provided we tell the function how this should all be handled. This function looks simple, but it is actually quite complicated in that it automates a crucial portion of demographic analysis. Specifically, it automates 1) the building of global models for each vital rate requested, 2) the exhaustive construction of all reduced models, and 3) the evaluation of models leading to the determination of the best-fit model. Let's look at some of the options that we will utilize for this purpose (please note that this list is only of the options shown in the block below - further options are shown in the documentation for `modelsearch`, and further theoretical details are shown in the Chapter 1 tutorial).</p>
  
  **historical**: Setting this to TRUE or FALSE tells R whether to include state in time *t*-1 in the global models.
  
  **approach**: This option tells R whether to use generalized linear models (`glm`), in which all factors are treated as fixed, or mixed effects models (`lme4`), in which factors are treated as either fixed or random, most notably time and individual identity. We encourage users to use the latter option, as it accounts for pseudoreplication, but the former approach is more common.
  
  **suite**: This tells R which factors to include in the global model. The options include `size`, which includes size only; `rep`, which includes reproductive status only; `main`, which includes both size and reproductive status as main effects only; `full`, which includes both size and reproductive status and all two-way interactions; and `const`, which includes only an intercept. These factors are in addition to individual identity and time.
  
  **vitalrates**: This option tells R which vital rates to model. The default is to model survival, size, and fecundity, but users can also model observation status and reproduction status.
  
  **juvestimate**: This optional term tells R the name of the juvenile stage transitioning to maturity, in cases where the dataset includes data on juveniles.
  
  **juvsize**: This optional term denotes whether size should be used as an independent term in models involving transition from the juvenile stage.
  
  **bestfit**: This describes whether the best-fit model was chosen as the model with the lowest AICc (`AICc`) or as the most parsimonious model (`AICc&k`), where the latter would be the model that has the fewest estimated parameters and is within 2 AICc units of the model with the lowest AICc. We encourage users to choose the latter, although many, perhaps most situations will lead to the same model.
  
  **sizedist**: This designates the distribution used for size. The options include `gaussian`, `poisson`, and `negbin`, the last of which refers to the negative binomial distribution.
  
  **fecdist**: This designates the distribution used for fecundity. The options include `gaussian`, `poisson`, and `negbin`, the last of which refers to the negative binomial distribution.
  
  **fectime**: This designates whether fecundity is estimated within time *t* (the default) or time *t*+1.
  
  **indiv**: This should include the name of the variable corresponding to individual identity in the vertical dataset.
  
  **patch**: This should include the name of the variable corresponding to patch identity in the vertical dataset.
  
  **year**: This should include the name of the variable corresponding to time *t* in the vertical dataset.
  
  **year.as.random**: This tells R whether to treat year as a random factor, and is only used when `approach = "lme4"`.
  
  **patch.as.random**: This tells R whether to treat patch identity as a random factor, and is only used when `approach = "lme4"`.
  
  **show.model.tables**: This tells R to include the full dredge model tables showing all models developed and their associated AICc values.
  
In addition, there are other options that provide flexibility in handling datasets with different designations for key variables, and allowing manual designation of stages. Here we create the ahistorical model set.

```r
lathmodelsln2 <- modelsearch(lathvertln, historical = FALSE, approach = "lme4", suite = "full",
                             vitalrates = c("surv", "obs", "size", "repst", "fec"), juvestimate = "Sdl",
                             bestfit = "AICc&k", sizedist = "gaussian", fecdist = "poisson",
                             indiv = "individ", patch = "patchid", year = "year2",
                             year.as.random = TRUE, patch.as.random = TRUE, quiet = TRUE)
#> Warning in modelsearch(lathvertln, 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...
#> 
#> 
#> Developing global model of observation probability...
#> 
#> 
#> Developing global model of size (Gaussian)...
#> 
#> 
#> 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
#> 
#> Developing global model of juvenile survival probability...
#> boundary (singular) fit: see ?isSingular
#> 
#> Developing global model of juvenile observation probability...
#> 
#> 
#> Developing global model of juvenile size (Gaussian)...
#> boundary (singular) fit: see ?isSingular
#> 
#> Juvenile reproduction response is constant so will not model it.
#> 
#> 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...
#> 
#> 
#> Commencing dredge of juvenile survival probability...
#> 
#> 
#> Commencing dredge of juvenile observation probability...
#> 
#> 
#> Commencing dredge of juvenile size...
#> 
#> 
#> Finished dredging all vital rates.
```

  <p style="text-indent: 20px">The output can be rather verbose, and so we have limited it with the \code{quiet = TRUE} option. The function was developed to provide text marker posts of what the function is doing and what it has accomplished, as well as to show all warnings from all workhorse functions used. Because we used the `lme4` approach here, this includes warnings originating from modeling mixed linear models with package `lme4` [@bates_fitting_2015]. It also shows warnings originating from the `dredge` function of package `MuMIn` [@barton_mumin_2014], which is the core function used in model building and AICc estimation. We encourage users to get familiar with interpreting these warnings and assessing the degree to which they impact their own analyses.</p>

  <p style="text-indent: 20px">In developing these linear models, the `modelsearch()` function set up a series of nested datasets for use in estimating vital rates as conditional rates and probabilities. The exact subsets created depend on which parameters are called for. It is rather important to consider which are required, and particular attention needs to be paid if there are size classes with sizes of 0. In these situations, it is best to consider a size of 0 as unobservable, and so to introduce observation status as a vital rate. This sets up datasets for size estimation that do not include 0 in the response term, because all 0 responses are already absorbed by observation status.</p>
  
  <p style="text-indent: 20px">A further point of interest is in the style of model selection allowed. Package `lefko3` differs from many other approaches in that it allows a more intelligent approach to selection. Indeed, instead of simply selecting the model with the lowest AICc, the reader can set `modelsearch()` to choose the model with fewest parameters within 2.0 AICc units of the model with the lowest AICc. This is the default setting, and is more in-line with model selection protocols preferred by information theorists [@burnham_model_2002].</p>
  
  <p style="text-indent: 20px">Let's take a peek at the summary for the model selection output. This condenses the output to the necessities.</p>

```r
summary(lathmodelsln2)
#> This LefkoMod object includes 8 linear models.
#> Best-fit model criterion used: AICc&k
#> 
#> ────────────────────────────────────────────────────────────────────────────────
#> Survival model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: binomial  ( logit )
#> Formula: alive3 ~ repstatus2 + (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: surv.data
#>       AIC       BIC    logLik  deviance  df.resid 
#> 1278.4711 1307.0289 -634.2356 1268.4711      2229 
#> Random effects:
#>  Groups  Name        Std.Dev.
#>  individ (Intercept) 0.3685  
#>  patchid (Intercept) 0.3160  
#>  year2   (Intercept) 0.6642  
#> Number of obs: 2234, groups:  individ, 257; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)   repstatus2  
#>       2.391        1.045  
#> 
#> ────────────────────────────────────────
#> 
#> Observation model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: binomial  ( logit )
#> Formula: obsstatus3 ~ repstatus2 + sizea2 + (1 | individ) + (1 | year2) +      (1 | patchid) + repstatus2:sizea2
#>    Data: obs.data
#>       AIC       BIC    logLik  deviance  df.resid 
#>  773.3746  812.6817 -379.6873  759.3746      2022 
#> Random effects:
#>  Groups  Name        Std.Dev.
#>  individ (Intercept) 0.3439  
#>  patchid (Intercept) 0.6078  
#>  year2   (Intercept) 2.9018  
#> Number of obs: 2029, groups:  individ, 254; patchid, 6; year2, 3
#> Fixed Effects:
#>       (Intercept)         repstatus2             sizea2  repstatus2:sizea2  
#>           4.90971           -6.27869           -0.09714            0.91857  
#> 
#> ────────────────────────────────────────
#> 
#> Size model:
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: sizea3 ~ repstatus2 + sizea2 + (1 | individ) + (1 | year2) +      (1 | patchid) + repstatus2:sizea2
#>    Data: size.data
#> REML criterion at convergence: 5986.416
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  individ  (Intercept) 0.6379  
#>  patchid  (Intercept) 0.1652  
#>  year2    (Intercept) 0.3511  
#>  Residual             1.0576  
#> Number of obs: 1916, groups:  individ, 254; patchid, 6; year2, 3
#> Fixed Effects:
#>       (Intercept)         repstatus2             sizea2  repstatus2:sizea2  
#>            2.6531            -1.5411             0.4095             0.2967  
#> 
#> ────────────────────────────────────────
#> 
#> Reproductive status model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: binomial  ( logit )
#> Formula: repstatus3 ~ repstatus2 + sizea2 + (1 | individ) + (1 | year2) +      (1 | patchid)
#>    Data: repst.data
#>       AIC       BIC    logLik  deviance  df.resid 
#> 1613.2218 1646.5698 -800.6109 1601.2218      1910 
#> Random effects:
#>  Groups  Name        Std.Dev. 
#>  individ (Intercept) 4.414e-05
#>  patchid (Intercept) 3.292e-01
#>  year2   (Intercept) 6.736e-01
#> Number of obs: 1916, groups:  individ, 254; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)   repstatus2       sizea2  
#>     -5.9117       0.6675       0.8264  
#> convergence code 0; 0 optimizer warnings; 1 lme4 warnings 
#> 
#> ────────────────────────────────────────
#> 
#> Fecundity model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: poisson  ( log )
#> Formula: feca2 ~ sizea2 + (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: fec.data
#>       AIC       BIC    logLik  deviance  df.resid 
#>  2036.563  2054.462 -1013.282  2026.563       260 
#> Random effects:
#>  Groups  Name        Std.Dev.
#>  individ (Intercept) 0.5422  
#>  patchid (Intercept) 0.1372  
#>  year2   (Intercept) 0.1573  
#> Number of obs: 265, groups:  individ, 101; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)       sizea2  
#>     -2.4599       0.6374  
#> 
#> ────────────────────────────────────────────────────────────────────────────────
#> Juvenile survival model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: binomial  ( logit )
#> Formula: alive3 ~ (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: juvsurv.data
#>       AIC       BIC    logLik  deviance  df.resid 
#>  333.7298  348.2833 -162.8649  325.7298       277 
#> Random effects:
#>  Groups  Name        Std.Dev.
#>  individ (Intercept) 0.0000  
#>  patchid (Intercept) 0.2104  
#>  year2   (Intercept) 0.0000  
#> Number of obs: 281, groups:  individ, 187; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)  
#>       1.001  
#> convergence code 0; 0 optimizer warnings; 1 lme4 warnings 
#> 
#> ────────────────────────────────────────
#> 
#> Juvenile observation model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: binomial  ( logit )
#> Formula: obsstatus3 ~ (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: juvobs.data
#>      AIC      BIC   logLik deviance df.resid 
#>  98.1375 111.4490 -45.0688  90.1375      202 
#> Random effects:
#>  Groups  Name        Std.Dev.
#>  individ (Intercept) 1.1474  
#>  patchid (Intercept) 0.8774  
#>  year2   (Intercept) 1.1929  
#> Number of obs: 206, groups:  individ, 151; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)  
#>       4.005  
#> 
#> ────────────────────────────────────────
#> 
#> Juvenile size model:
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: sizea3 ~ (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: juvsize.data
#> REML criterion at convergence: 243.3232
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  individ  (Intercept) 0.08956 
#>  patchid  (Intercept) 0.00000 
#>  year2    (Intercept) 0.09066 
#>  Residual             0.43789 
#> Number of obs: 193, groups:  individ, 144; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)  
#>        2.29  
#> convergence code 0; 0 optimizer warnings; 1 lme4 warnings 
#> 
#> ────────────────────────────────────────
#> 
#> Juvenile reproduction model:
#> [1] 0
#> 
#> 
#> ────────────────────────────────────────────────────────────────────────────────
#> 
#> 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     patchid
#> 4       surv3      alive3
#> 5        obs3  obsstatus3
#> 6       size3      sizea3
#> 7      repst3  repstatus3
#> 8        fec3       feca3
#> 9        fec2       feca2
#> 10      size2      sizea2
#> 11      size1        <NA>
#> 12     repst2  repstatus2
#> 13     repst1        <NA>
#> 14        age        <NA>
#> 
#> 
#> ────────────────────────────────────────────────────────────────────────────────
#> 
#> Quality control:
#> 
#> Survival estimated with 257 individuals and 2234 individual transitions.
#> Observation estimated with 254 individuals and 2029 individual transitions.
#> Size estimated with 254 individuals and 1916 individual transitions.
#> Reproductive status estimated with 254 individuals and 1916 individual transitions.
#> Fecundity estimated with 101 individuals and 265 individual transitions.
#> Juvenile survival estimated with 187 individuals and 281 individual transitions.
#> Juvenile observation estimated with 151 individuals and 206 individual transitions.
#> Juvenile size estimated with 144 individuals and 193 individual transitions.
#> Juvenile reproduction probability not estimated.
#> NULL
```

  <p style="text-indent: 20px">The best-fit models vary a bit in complexity. For example, survival is influenced by reproductive status in the current year, as well as by patch, individual identity, and year, while size and observation status are influenced by all of the as well as size in the current year. We can see these models explicitly, as well as the model tables developed, by calling them directly from the lefkoMod object. Here is the best-fit model.</p>

```r
lathmodelsln2$survival_model
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: binomial  ( logit )
#> Formula: alive3 ~ repstatus2 + (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: surv.data
#>       AIC       BIC    logLik  deviance  df.resid 
#> 1278.4711 1307.0289 -634.2356 1268.4711      2229 
#> Random effects:
#>  Groups  Name        Std.Dev.
#>  individ (Intercept) 0.3685  
#>  patchid (Intercept) 0.3160  
#>  year2   (Intercept) 0.6642  
#> Number of obs: 2234, groups:  individ, 257; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)   repstatus2  
#>       2.391        1.045
```

  <p style="text-indent: 20px">And now we move on to the historical model set. Note that this search will take more time, because the inclusion of status in time *t*-1 (and all two-way interactions when `suite = "full"`) increases the number of models that the `dredge` function will build and compare. The inclusion of these terms may also create more warnings as the function proceeds. Please also note that we include `quiet = TRUE` here, which prevents the model selection protocol from issuing warnings and diagnostic messages. We do this only because model selection generates a particularly verbose list of warnings in this case. Normally it is best to stick with the default, `quiet = FALSE`.</p>

```r
lathmodelsln3 <- modelsearch(lathvertln, historical = TRUE, approach = "lme4", suite = "full", 
                          vitalrates = c("surv", "obs", "size", "repst", "fec"), juvestimate = "Sdl",
                          bestfit = "AICc&k", sizedist = "gaussian", fecdist = "poisson",
                          indiv = "individ", patch = "patchid", year = "year2",
                          year.as.random = TRUE, patch.as.random = TRUE,
                          show.model.tables = TRUE, quiet = TRUE)
#> Warning in modelsearch(lathvertln, historical = TRUE, 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.00620144 (tol =
#> 0.002, component 1)
#> 
#> 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.0622003 (tol =
#> 0.002, component 1)
#> 
#> Developing global model of size (Gaussian)...
#> 
#> 
#> 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.0363529 (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
#> 
#> Developing global model of juvenile survival probability...
#> boundary (singular) fit: see ?isSingular
#> 
#> Developing global model of juvenile observation probability...
#> 
#> 
#> Developing global model of juvenile size (Gaussian)...
#> boundary (singular) fit: see ?isSingular
#> 
#> Juvenile reproduction response is constant so will not model it.
#> 
#> 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...
#> 
#> 
#> Commencing dredge of juvenile survival probability...
#> 
#> 
#> Commencing dredge of juvenile observation probability...
#> 
#> 
#> Commencing dredge of juvenile size...
#> 
#> 
#> Finished dredging all vital rates.
summary(lathmodelsln3)
#> This LefkoMod object includes 8 linear models.
#> Best-fit model criterion used: AICc&k
#> 
#> ────────────────────────────────────────────────────────────────────────────────
#> Survival model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: binomial  ( logit )
#> Formula: alive3 ~ repstatus2 + sizea1 + sizea2 + (1 | individ) + (1 |      year2) + (1 | patchid) + repstatus2:sizea1
#>    Data: surv.data
#>       AIC       BIC    logLik  deviance  df.resid 
#> 1248.3337 1294.0261 -616.1668 1232.3337      2226 
#> Random effects:
#>  Groups  Name        Std.Dev.
#>  individ (Intercept) 0.2729  
#>  patchid (Intercept) 0.2949  
#>  year2   (Intercept) 1.0027  
#> Number of obs: 2234, groups:  individ, 257; patchid, 6; year2, 3
#> Fixed Effects:
#>       (Intercept)         repstatus2             sizea1             sizea2  repstatus2:sizea1  
#>           2.27341            2.18473            0.21316           -0.09658           -0.24630  
#> 
#> ────────────────────────────────────────
#> 
#> Observation model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: binomial  ( logit )
#> Formula: obsstatus3 ~ repstatus1 + sizea1 + (1 | individ) + (1 | year2) +      (1 | patchid)
#>    Data: obs.data
#>       AIC       BIC    logLik  deviance  df.resid 
#>  772.7908  806.4825 -380.3954  760.7908      2023 
#> Random effects:
#>  Groups  Name        Std.Dev.
#>  individ (Intercept) 0.2163  
#>  patchid (Intercept) 0.5783  
#>  year2   (Intercept) 3.1318  
#> Number of obs: 2029, groups:  individ, 254; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)   repstatus1       sizea1  
#>      4.9068       0.7758      -0.1760  
#> 
#> ────────────────────────────────────────
#> 
#> Size model:
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: sizea3 ~ repstatus1 + repstatus2 + sizea1 + sizea2 + (1 | individ) +  
#>     (1 | year2) + (1 | patchid) + repstatus1:repstatus2 + repstatus1:sizea1 +      sizea1:sizea2
#>    Data: size.data
#> REML criterion at convergence: 5210.081
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  individ  (Intercept) 0.1645  
#>  patchid  (Intercept) 0.1144  
#>  year2    (Intercept) 0.4532  
#>  Residual             0.9168  
#> Number of obs: 1916, groups:  individ, 254; patchid, 6; year2, 3
#> Fixed Effects:
#>           (Intercept)             repstatus1             repstatus2                 sizea1                 sizea2  repstatus1:repstatus2  
#>               0.53719               -3.17885                0.46813                0.60500                0.77221               -0.37342  
#>     repstatus1:sizea1          sizea1:sizea2  
#>               0.56141               -0.09019  
#> 
#> ────────────────────────────────────────
#> 
#> Reproductive status model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: binomial  ( logit )
#> Formula: repstatus3 ~ repstatus1 + repstatus2 + sizea1 + sizea2 + (1 |      individ) + (1 | year2) + (1 | patchid) + repstatus1:sizea1 +  
#>     repstatus1:sizea2 + sizea1:sizea2
#>    Data: repst.data
#>       AIC       BIC    logLik  deviance  df.resid 
#> 1531.1030 1592.2409 -754.5515 1509.1030      1905 
#> Random effects:
#>  Groups  Name        Std.Dev. 
#>  individ (Intercept) 0.0004744
#>  patchid (Intercept) 0.3566851
#>  year2   (Intercept) 0.5721380
#> Number of obs: 1916, groups:  individ, 254; patchid, 6; year2, 3
#> Fixed Effects:
#>       (Intercept)         repstatus1         repstatus2             sizea1             sizea2  repstatus1:sizea1  repstatus1:sizea2  
#>          -8.21370           -0.12889            0.61756            0.39851            1.15713            0.57036           -0.55935  
#>     sizea1:sizea2  
#>          -0.05787  
#> convergence code 0; 0 optimizer warnings; 1 lme4 warnings 
#> 
#> ────────────────────────────────────────
#> 
#> Fecundity model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: poisson  ( log )
#> Formula: feca2 ~ repstatus1 + sizea1 + sizea2 + (1 | individ) + (1 | year2) +      (1 | patchid) + repstatus1:sizea2 + sizea1:sizea2
#>    Data: fec.data
#>       AIC       BIC    logLik  deviance  df.resid 
#>  2019.045  2051.262 -1000.522  2001.045       256 
#> Random effects:
#>  Groups  Name        Std.Dev. 
#>  individ (Intercept) 0.5170899
#>  patchid (Intercept) 0.1316614
#>  year2   (Intercept) 0.0001251
#> Number of obs: 265, groups:  individ, 101; patchid, 6; year2, 3
#> Fixed Effects:
#>       (Intercept)         repstatus1             sizea1             sizea2  repstatus1:sizea2      sizea1:sizea2  
#>           -3.0775            -4.9835             0.7627             0.7020             0.7301            -0.1051  
#> convergence code 0; 0 optimizer warnings; 3 lme4 warnings 
#> 
#> ────────────────────────────────────────────────────────────────────────────────
#> Juvenile survival model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: binomial  ( logit )
#> Formula: alive3 ~ (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: juvsurv.data
#>       AIC       BIC    logLik  deviance  df.resid 
#>  333.7298  348.2833 -162.8649  325.7298       277 
#> Random effects:
#>  Groups  Name        Std.Dev.
#>  individ (Intercept) 0.0000  
#>  patchid (Intercept) 0.2104  
#>  year2   (Intercept) 0.0000  
#> Number of obs: 281, groups:  individ, 187; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)  
#>       1.001  
#> convergence code 0; 0 optimizer warnings; 1 lme4 warnings 
#> 
#> ────────────────────────────────────────
#> 
#> Juvenile observation model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: binomial  ( logit )
#> Formula: obsstatus3 ~ (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: juvobs.data
#>      AIC      BIC   logLik deviance df.resid 
#>  98.1375 111.4490 -45.0688  90.1375      202 
#> Random effects:
#>  Groups  Name        Std.Dev.
#>  individ (Intercept) 1.1474  
#>  patchid (Intercept) 0.8774  
#>  year2   (Intercept) 1.1929  
#> Number of obs: 206, groups:  individ, 151; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)  
#>       4.005  
#> 
#> ────────────────────────────────────────
#> 
#> Juvenile size model:
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: sizea3 ~ (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: juvsize.data
#> REML criterion at convergence: 243.3232
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  individ  (Intercept) 0.08956 
#>  patchid  (Intercept) 0.00000 
#>  year2    (Intercept) 0.09066 
#>  Residual             0.43789 
#> Number of obs: 193, groups:  individ, 144; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)  
#>        2.29  
#> convergence code 0; 0 optimizer warnings; 1 lme4 warnings 
#> 
#> ────────────────────────────────────────
#> 
#> Juvenile reproduction model:
#> [1] 0
#> 
#> 
#> ────────────────────────────────────────────────────────────────────────────────
#> 
#> 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     patchid
#> 4       surv3      alive3
#> 5        obs3  obsstatus3
#> 6       size3      sizea3
#> 7      repst3  repstatus3
#> 8        fec3       feca3
#> 9        fec2       feca2
#> 10      size2      sizea2
#> 11      size1      sizea1
#> 12     repst2  repstatus2
#> 13     repst1  repstatus1
#> 14        age        <NA>
#> 
#> 
#> ────────────────────────────────────────────────────────────────────────────────
#> 
#> Quality control:
#> 
#> Survival estimated with 257 individuals and 2234 individual transitions.
#> Observation estimated with 254 individuals and 2029 individual transitions.
#> Size estimated with 254 individuals and 1916 individual transitions.
#> Reproductive status estimated with 254 individuals and 1916 individual transitions.
#> Fecundity estimated with 101 individuals and 265 individual transitions.
#> Juvenile survival estimated with 187 individuals and 281 individual transitions.
#> Juvenile observation estimated with 151 individuals and 206 individual transitions.
#> Juvenile size estimated with 144 individuals and 193 individual transitions.
#> Juvenile reproduction probability not estimated.
#> NULL
```
The summary will no doubt reveal some similarities, but also many differences. Historical models incorporate a further time step back in the global model, so some of the best-fit models will have these historical terms in them. Note that model sets that include historical terms should not be used to create ahistorical matrices, since the best-fit models can only be used to predict rates and probabilities using all terms in those final models.

  <p style="text-indent: 20px">Now we are ready to estimate matrices.</p>
  
#### Step 4. Estimate matrices
  
  <p style="text-indent: 20px">We will start by estimating the ahistorical set of matrices. Notice that we are matching the ahistorical matrix creation function, `flefko2`, with the appropriate ahistorical input, including the ahistorical lefkoMod object `lathmodelsln2`.</p>

```r
lathmat2ln <- flefko2(stageframe = lathframeln, modelsuite = lathmodelsln2,
                      data = lathvertln, repmatrix = lathrepmln, overwrite = lathover2,
                      year.as.random = FALSE, patch.as.random = FALSE)
summary(lathmat2ln)
#> 
#> This lefkoMat object contains 18 matrices.
#> 
#> Each matrix is a square matrix with 21 rows and columns, and a total of 441 elements.
#> A total of 6714 survival-transitions were estimated, with 373 per matrix.
#> A total of 324 fecundity transitions were estimated, with 18 per matrix.
#> 
#> Vital rate modeling quality control:
#> 
#> Survival estimated with 257 individuals and 2234 individual transitions.
#> Observation estimated with 254 individuals and 2029 individual transitions.
#> Size estimated with 254 individuals and 1916 individual transitions.
#> Reproductive status estimated with 254 individuals and 1916 individual transitions.
#> Fecundity estimated with 101 individuals and 265 individual transitions.
#> Juvenile survival estimated with 187 individuals and 281 individual transitions.
#> Juvenile observation estimated with 151 individuals and 206 individual transitions.
#> Juvenile size estimated with 144 individuals and 193 individual transitions.
#> Juvenile reproductive status estimated with 0 individuals and 0 individual transitions.
#> NULL
```

  <p style="text-indent: 20px">Let's take a look at the first of these matrices.</p>

```r
lathmat2ln$A[[1]]
#>        [,1]         [,2]         [,3]         [,4]         [,5]         [,6]         [,7]         [,8]         [,9]        [,10]
#>  [1,] 0.345 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#>  [2,] 0.054 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
#>  [3,] 0.000 1.308258e-02 5.431474e-02 5.951592e-02 6.517972e-02 7.134032e-02 7.803300e-02 8.529400e-02 9.316011e-02 1.016683e-01
#>  [4,] 0.000 8.527702e-03 7.994576e-02 3.791326e-02 1.540099e-02 5.329694e-03 1.553032e-03 3.722605e-04 7.052148e-05 1.004242e-05
#>  [5,] 0.000 5.253179e-01 2.574400e-01 1.760678e-01 1.031444e-01 5.147629e-02 2.163183e-02 7.477701e-03 2.042917e-03 4.195425e-04
#>  [6,] 0.000 1.758073e-01 3.390616e-01 3.344187e-01 2.825298e-01 2.033453e-01 1.232333e-01 6.143433e-02 2.420482e-02 7.168614e-03
#>  [7,] 0.000 3.196508e-04 1.826430e-01 2.597903e-01 3.165227e-01 3.285359e-01 2.871343e-01 2.064315e-01 1.172937e-01 5.009753e-02
#>  [8,] 0.000 3.157470e-09 4.023919e-02 8.254245e-02 1.450331e-01 2.170968e-01 2.736300e-01 2.837023e-01 2.324715e-01 1.431922e-01
#>  [9,] 0.000 1.694443e-16 3.625913e-03 1.072639e-02 2.718013e-02 5.867401e-02 1.066508e-01 1.594672e-01 1.884458e-01 1.673956e-01
#> [10,] 0.000 4.940139e-26 1.336310e-04 5.701008e-04 2.083329e-03 6.485748e-03 1.700148e-02 3.666080e-02 6.247770e-02 8.003700e-02
#> [11,] 0.000 7.824846e-38 2.014277e-06 1.239286e-05 6.531092e-05 2.932218e-04 1.108488e-03 3.447103e-03 8.471991e-03 1.565162e-02
#> [12,] 0.000 6.733437e-52 1.241804e-08 1.101827e-07 8.374054e-07 5.421934e-06 2.955950e-05 1.325649e-04 4.698593e-04 1.251843e-03
#> [13,] 0.000 0.000000e+00 2.131152e-04 2.309559e-04 2.143904e-04 1.695422e-04 1.128950e-04 6.183867e-05 2.677032e-05 8.711431e-06
#> [14,] 0.000 0.000000e+00 6.862702e-04 1.072551e-03 1.435827e-03 1.637506e-03 1.572489e-03 1.242171e-03 7.755017e-04 3.639376e-04
#> [15,] 0.000 0.000000e+00 9.038527e-04 2.037175e-03 3.932973e-03 6.468591e-03 8.958239e-03 1.020527e-02 9.188273e-03 6.218507e-03
#> [16,] 0.000 0.000000e+00 4.868801e-04 1.582563e-03 4.406173e-03 1.045101e-02 2.087274e-02 3.429172e-02 4.452530e-02 4.345775e-02
#> [17,] 0.000 0.000000e+00 1.072675e-04 5.028231e-04 2.018943e-03 6.906037e-03 1.989107e-02 4.712768e-02 8.824737e-02 1.242139e-01
#> [18,] 0.000 0.000000e+00 9.665769e-06 6.534185e-05 3.783627e-04 1.866471e-03 7.752801e-03 2.649016e-02 7.153499e-02 1.452095e-01
#> [19,] 0.000 0.000000e+00 3.562266e-07 3.472878e-06 2.900110e-05 2.063173e-04 1.235893e-03 6.089970e-03 2.371686e-02 6.942913e-02
#> [20,] 0.000 0.000000e+00 5.369554e-09 7.549349e-08 9.091647e-07 9.327640e-06 8.057964e-05 5.726213e-04 3.216011e-03 1.357720e-02
#> [21,] 0.000 0.000000e+00 3.310335e-11 6.711990e-10 1.165715e-08 1.724765e-07 2.148777e-06 2.202124e-05 1.783610e-04 1.085927e-03
#>              [,11]        [,12]        [,13]        [,14]        [,15]        [,16]        [,17]        [,18]        [,19]        [,20]
#>  [1,] 0.000000e+00 0.000000e+00 4.685856e-02 8.863230e-02 1.676468e-01 3.171015e-01 5.997931e-01 1.134500e+00 2.145891e+00 4.058923e+00
#>  [2,] 0.000000e+00 0.000000e+00 7.334383e-03 1.387288e-02 2.624036e-02 4.963328e-02 9.388065e-02 1.775739e-01 3.358786e-01 6.353097e-01
#>  [3,] 1.108554e-01 1.207574e-01 9.230245e-01 8.503261e-01 7.211756e-01 5.360521e-01 3.384881e-01 1.841610e-01 9.042234e-02 4.191361e-02
#>  [4,] 1.042441e-06 8.014322e-08 1.621162e-02 1.447885e-02 7.482686e-03 2.086890e-03 2.934903e-04 2.008612e-05 6.750008e-07 1.173045e-08
#>  [5,] 6.280543e-05 6.963384e-06 2.474932e-02 4.156147e-02 4.038635e-02 2.117858e-02 5.600300e-03 7.206655e-04 4.553676e-05 1.487965e-06
#>  [6,] 1.547621e-03 2.474547e-04 1.545332e-02 4.879437e-02 8.915253e-02 8.790551e-02 4.370696e-02 1.057531e-02 1.256440e-03 7.719558e-05
#>  [7,] 1.559745e-02 3.596608e-03 3.946412e-03 2.342990e-02 8.049234e-02 1.492303e-01 1.395121e-01 6.347089e-02 1.417892e-02 1.638001e-03
#>  [8,] 6.429316e-02 2.138023e-02 4.121970e-04 4.601432e-03 2.972332e-02 1.036144e-01 1.821356e-01 1.558037e-01 6.544353e-02 1.421535e-02
#>  [9,] 1.083921e-01 5.198204e-02 1.760877e-05 3.696047e-04 4.489127e-03 2.942421e-02 9.725229e-02 1.564238e-01 1.235412e-01 5.045720e-02
#> [10,] 7.473990e-02 5.169115e-02 3.076630e-07 1.214238e-05 2.772990e-04 3.417523e-03 2.123864e-02 6.423178e-02 9.538478e-02 7.325047e-02
#> [11,] 2.107799e-02 2.102331e-02 2.198585e-09 1.631518e-07 7.005783e-06 1.623454e-04 1.897035e-03 1.078745e-02 3.012091e-02 4.349302e-02
#> [12,] 2.431240e-03 3.497097e-03 6.425887e-12 8.966067e-10 7.239138e-08 3.154208e-06 6.930201e-05 7.409851e-04 3.890263e-03 1.056211e-02
#> [13,] 2.066434e-06 3.630407e-07 1.925200e-04 3.929183e-04 4.640286e-04 2.957370e-04 9.504269e-05 1.486414e-05 1.141476e-06 4.533107e-08
#> [14,] 1.244993e-04 3.154343e-05 2.939089e-04 1.127871e-03 2.504505e-03 3.001255e-03 1.813578e-03 5.333072e-04 7.700601e-05 5.750081e-06
#> [15,] 3.067851e-03 1.120945e-03 1.835149e-04 1.324153e-03 5.528672e-03 1.245725e-02 1.415388e-02 7.825948e-03 2.124733e-03 2.983141e-04
#> [16,] 3.091886e-02 1.629227e-02 4.686536e-05 6.358268e-04 4.991622e-03 2.114770e-02 4.517901e-02 4.696976e-02 2.397761e-02 6.329880e-03
#> [17,] 1.274484e-01 9.685028e-02 4.895018e-06 1.248709e-04 1.843251e-03 1.468338e-02 5.898204e-02 1.152979e-01 1.106698e-01 5.493373e-02
#> [18,] 2.148659e-01 2.354734e-01 2.091118e-07 1.003011e-05 2.783871e-04 4.169758e-03 3.149377e-02 1.157568e-01 2.089174e-01 1.949865e-01
#> [19,] 1.481570e-01 2.341558e-01 3.653631e-09 3.295127e-07 1.719632e-05 4.843034e-04 6.877831e-03 4.753283e-02 1.613027e-01 2.830687e-01
#> [20,] 4.178294e-02 9.523347e-02 2.610915e-11 4.427517e-09 4.344541e-07 2.300626e-05 6.143278e-04 7.982933e-03 5.093668e-02 1.680742e-01
#> [21,] 4.819450e-03 1.584150e-02 7.631020e-14 2.433158e-11 4.489253e-09 4.469885e-07 2.244247e-05 5.483441e-04 6.578722e-03 4.081615e-02
#>              [,21]
#>  [1,] 7.677397e+00
#>  [2,] 1.201679e+00
#>  [3,] 1.888172e-02
#>  [4,] 1.134509e-10
#>  [5,] 2.705864e-08
#>  [6,] 2.639526e-06
#>  [7,] 1.053096e-04
#>  [8,] 1.718429e-03
#>  [9,] 1.146878e-02
#> [10,] 3.130579e-02
#> [11,] 3.495055e-02
#> [12,] 1.595898e-02
#> [13,] 1.001862e-09
#> [14,] 2.389495e-07
#> [15,] 2.330914e-05
#> [16,] 9.299681e-04
#> [17,] 1.517511e-02
#> [18,] 1.012785e-01
#> [19,] 2.764553e-01
#> [20,] 3.086414e-01
#> [21,] 1.409306e-01
```

The number of elements estimated is vastly greater now than in the previous case, because each matrix element must be estimated based on the vital rate linear models supplied. And the matrix is overwhelmingly composed of elements that are to be estimated, whereas in the raw matrix case, elements would only be estimated if individuals actually passed through that transition. These are interesting and important differences between raw and function-derived matrices.

  <p style="text-indent: 20px">Now we will estimate the historical matrices.</p>

```r
lathmat3ln <- flefko3(stageframe = lathframeln, modelsuite = lathmodelsln3,
                      data = lathvertln, repmatrix = lathrepmln, overwrite = lathover3,
                      year.as.random = FALSE, patch.as.random = FALSE)
summary(lathmat3ln)
#> 
#> This lefkoMat object contains 18 matrices.
#> 
#> Each matrix is a square matrix with 441 rows and columns, and a total of 194481 elements.
#> A total of 133902 survival-transitions were estimated, with 7439 per matrix.
#> A total of 6804 fecundity transitions were estimated, with 378 per matrix.
#> 
#> Vital rate modeling quality control:
#> 
#> Survival estimated with 257 individuals and 2234 individual transitions.
#> Observation estimated with 254 individuals and 2029 individual transitions.
#> Size estimated with 254 individuals and 1916 individual transitions.
#> Reproductive status estimated with 254 individuals and 1916 individual transitions.
#> Fecundity estimated with 101 individuals and 265 individual transitions.
#> Juvenile survival estimated with 187 individuals and 281 individual transitions.
#> Juvenile observation estimated with 151 individuals and 206 individual transitions.
#> Juvenile size estimated with 144 individuals and 193 individual transitions.
#> Juvenile reproductive status estimated with 0 individuals and 0 individual transitions.
#> NULL
```

Once again, the matrices are vastly bigger. Indeed, the ahistorical matrices here are 21 x 21, so contain 441 elements each. The historical matrices are 21<sup>2</sup> x 21<sup>2</sup> (441 x 441), or 194,481 elements large. However, because of the prevalence of structural 0's in the latter, only 7,439 elements are actually estimated. Another point of interest here is that the number of **A** matrices estimated is 18, which is the same as the number of ahistorical matrices. In the raw matrix case, there would actually only be 12 historical matrices estimated, because of the lack of information for the time step previous to the start of the dataset in each of the 6 patches. Here, linear modeling allows us to make assumptions leading to estimated values for even those first transitions.

  <p style="text-indent: 20px">Because of the size of the matrices and the proliferation of structural zeroes, we may wish to analyze these matrices after removing some unneeded rows and columns. Let's re-do the above, but overwrite the models with smaller matrices using the `reduce = TRUE` option, which eliminates stages or stage-pairs with corresponding row and column vectors both equal to zero vectors.</p>

```r
lathmat3ln <- flefko3(stageframe = lathframeln, modelsuite = lathmodelsln3, 
                    data = lathvertln, repmatrix = lathrepmln, overwrite = lathover3, 
                    year.as.random = FALSE, patch.as.random = FALSE, reduce = TRUE)

summary(lathmat3ln)
#> 
#> This lefkoMat object contains 18 matrices.
#> 
#> Each matrix is a square matrix with 420 rows and columns, and a total of 176400 elements.
#> A total of 133902 survival-transitions were estimated, with 7439 per matrix.
#> A total of 6804 fecundity transitions were estimated, with 378 per matrix.
#> 
#> Vital rate modeling quality control:
#> 
#> Survival estimated with 257 individuals and 2234 individual transitions.
#> Observation estimated with 254 individuals and 2029 individual transitions.
#> Size estimated with 254 individuals and 1916 individual transitions.
#> Reproductive status estimated with 254 individuals and 1916 individual transitions.
#> Fecundity estimated with 101 individuals and 265 individual transitions.
#> Juvenile survival estimated with 187 individuals and 281 individual transitions.
#> Juvenile observation estimated with 151 individuals and 206 individual transitions.
#> Juvenile size estimated with 144 individuals and 193 individual transitions.
#> Juvenile reproductive status estimated with 0 individuals and 0 individual transitions.
#> NULL
```
This exercise eliminated 21 rows and columns, so the matrices definitely got smaller.

  <p style="text-indent: 20px">Now the mean matrices. Note that these include 6 mean matrices for the patches, followed by a grand mean matrix.</p>

```r
lathmat2lnmean <- lmean(lathmat2ln)
lathmat3lnmean <- lmean(lathmat3ln)

summary(lathmat2lnmean)
#> 
#> This lefkoMat object contains 7 matrices.
#> 
#> Each matrix is a square matrix with 21 rows and columns, and a total of 441 elements.
#> A total of 2611 survival-transitions were estimated, with 373 per matrix.
#> A total of 126 fecundity transitions were estimated, with 18 per matrix.
#> 
#> Vital rate modeling quality control:
#> 
#> Survival estimated with 257 individuals and 2234 individual transitions.
#> Observation estimated with 254 individuals and 2029 individual transitions.
#> Size estimated with 254 individuals and 1916 individual transitions.
#> Reproductive status estimated with 254 individuals and 1916 individual transitions.
#> Fecundity estimated with 101 individuals and 265 individual transitions.
#> Juvenile survival estimated with 187 individuals and 281 individual transitions.
#> Juvenile observation estimated with 151 individuals and 206 individual transitions.
#> Juvenile size estimated with 144 individuals and 193 individual transitions.
#> Juvenile reproductive status estimated with 0 individuals and 0 individual transitions.
#> NULL
summary(lathmat3lnmean)
#> 
#> This lefkoMat object contains 7 matrices.
#> 
#> Each matrix is a square matrix with 420 rows and columns, and a total of 176400 elements.
#> A total of 52073 survival-transitions were estimated, with 7439 per matrix.
#> A total of 2646 fecundity transitions were estimated, with 378 per matrix.
#> 
#> Vital rate modeling quality control:
#> 
#> Survival estimated with 257 individuals and 2234 individual transitions.
#> Observation estimated with 254 individuals and 2029 individual transitions.
#> Size estimated with 254 individuals and 1916 individual transitions.
#> Reproductive status estimated with 254 individuals and 1916 individual transitions.
#> Fecundity estimated with 101 individuals and 265 individual transitions.
#> Juvenile survival estimated with 187 individuals and 281 individual transitions.
#> Juvenile observation estimated with 151 individuals and 206 individual transitions.
#> Juvenile size estimated with 144 individuals and 193 individual transitions.
#> Juvenile reproductive status estimated with 0 individuals and 0 individual transitions.
#> NULL
```

And here we do some quality control. We will focus on the grand mean U matrix in each case.

```r
writeLines("Stage survival in ahistorical grand mean matrix")
#> Stage survival in ahistorical grand mean matrix
summary(colSums(lathmat2lnmean$U[[7]]))
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.3990  0.8949  0.8961  0.8870  0.9551  0.9598

writeLines("\nStage survival in historical grand mean matrix")
#> 
#> Stage survival in historical grand mean matrix
summary(colSums(lathmat3lnmean$U[[7]]))
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.0000  0.8726  0.9429  0.8657  0.9635  0.9776
```
Quick scans of these numbers show no real problems.

#### Step 5. Analyses

  <p style="text-indent: 20px">Finally, let's estimate the deterministic population growth rate in each case. First the ahistorical case.</p>

```r
lambda3(lathmat2lnmean)
#>   pop patch    lambda
#> 1   1     1 0.9623922
#> 2   1     2 0.9166536
#> 3   1     3 0.9382975
#> 4   1     4 0.8957659
#> 5   1     5 0.8834887
#> 6   1     6 0.9533979
#> 7   1   All 0.9250544
```

Now the historical case.

```r
lambda3(lathmat3lnmean)
#>   pop patch    lambda
#> 1   1     1 0.9645367
#> 2   1     2 0.9397586
#> 3   1     3 0.9459287
#> 4   1     4 0.9221450
#> 5   1     5 0.8917885
#> 6   1     6 0.9862835
#> 7   1   All 0.9422742
```
Note that the population growth rate values differ, and that in this case they are higher in the historical case. These are also both higher than in the raw matrix case, which probably reflects the case that raw matrices can have many 0's simply because of low sample size. These lambda values are potentially more realistic, though only if the assumptions under which the matrices were constructed are realistic.

  <p style="text-indent: 20px">Now let's look at the stable stage distributions, as follows for the ahistorical case. Our numbers here match those produced by package `popbio`'s `stable.stage()` function.</p>

```r
stablestage3(lathmat2lnmean)
#>     matrix new_stage_id orig_stage_id original_size      ss_prop
#> 1        1            1            Sd           0.0 1.926355e-01
#> 2        1            2           Sdl           4.6 3.015164e-02
#> 3        1            3          Dorm           0.0 5.517375e-02
#> 4        1            4         Sz1nr           1.0 9.465746e-03
#> 5        1            5         Sz2nr           2.0 5.788582e-02
#> 6        1            6         Sz3nr           3.0 1.170173e-01
#> 7        1            7         Sz4nr           4.0 1.778001e-01
#> 8        1            8         Sz5nr           5.0 1.571412e-01
#> 9        1            9         Sz6nr           6.0 7.310088e-02
#> 10       1           10         Sz7nr           7.0 1.759963e-02
#> 11       1           11         Sz8nr           8.0 2.235860e-03
#> 12       1           12         Sz9nr           9.0 1.598921e-04
#> 13       1           13          Sz1r           1.0 1.490597e-04
#> 14       1           14          Sz2r           2.0 1.500553e-03
#> 15       1           15          Sz3r           3.0 8.268151e-03
#> 16       1           16          Sz4r           4.0 2.386070e-02
#> 17       1           17          Sz5r           5.0 3.533135e-02
#> 18       1           18          Sz6r           6.0 2.689423e-02
#> 19       1           19          Sz7r           7.0 1.083238e-02
#> 20       1           20          Sz8r           8.0 2.457276e-03
#> 21       1           21          Sz9r           9.0 3.390302e-04
#> 22       2            1            Sd           0.0 1.634621e-01
#> 23       2            2           Sdl           4.6 2.558535e-02
#> 24       2            3          Dorm           0.0 8.326414e-02
#> 25       2            4         Sz1nr           1.0 1.433495e-02
#> 26       2            5         Sz2nr           2.0 7.096926e-02
#> 27       2            6         Sz3nr           3.0 1.380675e-01
#> 28       2            7         Sz4nr           4.0 1.883963e-01
#> 29       2            8         Sz5nr           5.0 1.496962e-01
#> 30       2            9         Sz6nr           6.0 6.299523e-02
#> 31       2           10         Sz7nr           7.0 1.372097e-02
#> 32       2           11         Sz8nr           8.0 1.570285e-03
#> 33       2           12         Sz9nr           9.0 1.004937e-04
#> 34       2           13          Sz1r           1.0 1.743944e-04
#> 35       2           14          Sz2r           2.0 1.587675e-03
#> 36       2           15          Sz3r           3.0 7.978266e-03
#> 37       2           16          Sz4r           4.0 2.104865e-02
#> 38       2           17          Sz5r           5.0 2.846314e-02
#> 39       2           18          Sz6r           6.0 1.972697e-02
#> 40       2           19          Sz7r           7.0 7.201399e-03
#> 41       2           20          Sz8r           8.0 1.473500e-03
#> 42       2           21          Sz9r           9.0 1.832471e-04
#> 43       3            1            Sd           0.0 1.998832e-01
#> 44       3            2           Sdl           4.6 3.128605e-02
#> 45       3            3          Dorm           0.0 7.355740e-02
#> 46       3            4         Sz1nr           1.0 8.281027e-03
#> 47       3            5         Sz2nr           2.0 5.327696e-02
#> 48       3            6         Sz3nr           3.0 1.014798e-01
#> 49       3            7         Sz4nr           4.0 1.601106e-01
#> 50       3            8         Sz5nr           5.0 1.574062e-01
#> 51       3            9         Sz6nr           6.0 8.416290e-02
#> 52       3           10         Sz7nr           7.0 2.369910e-02
#> 53       3           11         Sz8nr           8.0 3.568005e-03
#> 54       3           12         Sz9nr           9.0 3.066495e-04
#> 55       3           13          Sz1r           1.0 7.939909e-05
#> 56       3           14          Sz2r           2.0 8.714847e-04
#> 57       3           15          Sz3r           3.0 5.454738e-03
#> 58       3           16          Sz4r           4.0 1.826626e-02
#> 59       3           17          Sz5r           5.0 3.173116e-02
#> 60       3           18          Sz6r           6.0 2.854976e-02
#> 61       3           19          Sz7r           7.0 1.368827e-02
#> 62       3           20          Sz8r           8.0 3.722924e-03
#> 63       3           21          Sz9r           9.0 6.181143e-04
#> 64       4            1            Sd           0.0 1.457089e-01
#> 65       4            2           Sdl           4.6 2.280660e-02
#> 66       4            3          Dorm           0.0 3.762152e-02
#> 67       4            4         Sz1nr           1.0 1.346386e-02
#> 68       4            5         Sz2nr           2.0 7.373736e-02
#> 69       4            6         Sz3nr           3.0 1.615466e-01
#> 70       4            7         Sz4nr           4.0 2.242193e-01
#> 71       4            8         Sz5nr           5.0 1.677929e-01
#> 72       4            9         Sz6nr           6.0 6.408351e-02
#> 73       4           10         Sz7nr           7.0 1.240002e-02
#> 74       4           11         Sz8nr           8.0 1.235381e-03
#> 75       4           12         Sz9nr           9.0 6.693246e-05
#> 76       4           13          Sz1r           1.0 2.082509e-04
#> 77       4           14          Sz2r           2.0 1.832518e-03
#> 78       4           15          Sz3r           3.0 8.502509e-03
#> 79       4           16          Sz4r           4.0 2.021375e-02
#> 80       4           17          Sz5r           5.0 2.427979e-02
#> 81       4           18          Sz6r           6.0 1.475324e-02
#> 82       4           19          Sz7r           7.0 4.643048e-03
#> 83       4           20          Sz8r           8.0 8.014379e-04
#> 84       4           21          Sz9r           9.0 8.262208e-05
#> 85       5            1            Sd           0.0 7.710097e-02
#> 86       5            2           Sdl           4.6 1.206800e-02
#> 87       5            3          Dorm           0.0 1.044259e-01
#> 88       5            4         Sz1nr           1.0 2.567431e-02
#> 89       5            5         Sz2nr           2.0 9.775685e-02
#> 90       5            6         Sz3nr           3.0 1.900308e-01
#> 91       5            7         Sz4nr           4.0 2.256583e-01
#> 92       5            8         Sz5nr           5.0 1.485491e-01
#> 93       5            9         Sz6nr           6.0 5.091582e-02
#> 94       5           10         Sz7nr           7.0 8.895344e-03
#> 95       5           11         Sz8nr           8.0 7.988081e-04
#> 96       5           12         Sz9nr           9.0 3.873100e-05
#> 97       5           13          Sz1r           1.0 2.371339e-04
#> 98       5           14          Sz2r           2.0 1.833367e-03
#> 99       5           15          Sz3r           3.0 7.672476e-03
#> 100      5           16          Sz4r           4.0 1.662393e-02
#> 101      5           17          Sz5r           5.0 1.823781e-02
#> 102      5           18          Sz6r           6.0 1.010464e-02
#> 103      5           19          Sz7r           7.0 2.886581e-03
#> 104      5           20          Sz8r           8.0 4.494441e-04
#> 105      5           21          Sz9r           9.0 4.168442e-05
#> 106      6            1            Sd           0.0 2.150166e-01
#> 107      6            2           Sdl           4.6 3.365477e-02
#> 108      6            3          Dorm           0.0 4.326665e-02
#> 109      6            4         Sz1nr           1.0 8.724746e-03
#> 110      6            5         Sz2nr           2.0 5.761931e-02
#> 111      6            6         Sz3nr           3.0 1.127029e-01
#> 112      6            7         Sz4nr           4.0 1.671843e-01
#> 113      6            8         Sz5nr           5.0 1.428120e-01
#> 114      6            9         Sz6nr           6.0 6.379700e-02
#> 115      6           10         Sz7nr           7.0 1.476912e-02
#> 116      6           11         Sz8nr           8.0 1.818030e-03
#> 117      6           12         Sz9nr           9.0 1.276185e-04
#> 118      6           13          Sz1r           1.0 2.304757e-04
#> 119      6           14          Sz2r           2.0 2.215994e-03
#> 120      6           15          Sz3r           3.0 1.156495e-02
#> 121      6           16          Sz4r           4.0 3.160724e-02
#> 122      6           17          Sz5r           5.0 4.451157e-02
#> 123      6           18          Sz6r           6.0 3.247525e-02
#> 124      6           19          Sz7r           7.0 1.268547e-02
#> 125      6           20          Sz8r           8.0 2.829220e-03
#> 126      6           21          Sz9r           9.0 3.868074e-04
#> 127      7            1            Sd           0.0 1.671435e-01
#> 128      7            2           Sdl           4.6 2.616158e-02
#> 129      7            3          Dorm           0.0 6.614619e-02
#> 130      7            4         Sz1nr           1.0 1.273137e-02
#> 131      7            5         Sz2nr           2.0 6.687129e-02
#> 132      7            6         Sz3nr           3.0 1.342981e-01
#> 133      7            7         Sz4nr           4.0 1.903432e-01
#> 134      7            8         Sz5nr           5.0 1.559146e-01
#> 135      7            9         Sz6nr           6.0 6.766094e-02
#> 136      7           10         Sz7nr           7.0 1.529048e-02
#> 137      7           11         Sz8nr           8.0 1.829913e-03
#> 138      7           12         Sz9nr           9.0 1.231824e-04
#> 139      7           13          Sz1r           1.0 1.766866e-04
#> 140      7           14          Sz2r           2.0 1.637860e-03
#> 141      7           15          Sz3r           3.0 8.335812e-03
#> 142      7           16          Sz4r           4.0 2.231679e-02
#> 143      7           17          Sz5r           5.0 3.080017e-02
#> 144      7           18          Sz6r           6.0 2.193995e-02
#> 145      7           19          Sz7r           7.0 8.287115e-03
#> 146      7           20          Sz8r           8.0 1.762796e-03
#> 147      7           21          Sz9r           9.0 2.283664e-04
```

The output is a dataframe with stage information and stable stage distributions for each matrix, with matrix info appended. Here, we have 7 sets of stable stage distributions because there are 7 matrices. Note that the last matrix, the 7th, is the grand mean matrix.

  <p style="text-indent: 20px">Now the historical case. In this case, the historical output includes 2 data frames, so we will first create an object to hold the output and then look at the first data frame. However, because of the size of the output, we will first only inspect the summary.</p>

```r
ehrlen3ss <- stablestage3(lathmat3lnmean)
summary(ehrlen3ss$hist)
#>      matrix  orig_stage_id_2 orig_stage_id_1 new_stage_id_2  new_stage_id_1     ss_prop         
#>  Min.   :1   Sdl    : 147    Sz1r   : 147    Min.   : 1.00   Min.   : 1.00   Min.   :0.000e+00  
#>  1st Qu.:2   Sz1r   : 147    Sz2r   : 147    1st Qu.: 6.00   1st Qu.: 6.00   1st Qu.:1.170e-06  
#>  Median :4   Sz2r   : 147    Sz3r   : 147    Median :11.00   Median :11.00   Median :5.047e-05  
#>  Mean   :4   Sz3r   : 147    Sz4r   : 147    Mean   :11.35   Mean   :11.34   Mean   :2.381e-03  
#>  3rd Qu.:6   Sz4r   : 147    Sz5r   : 147    3rd Qu.:16.25   3rd Qu.:16.25   3rd Qu.:9.793e-04  
#>  Max.   :7   Sz5r   : 147    Sz6r   : 147    Max.   :21.00   Max.   :21.00   Max.   :1.428e-01  
#>              (Other):2058    (Other):2058
```

Now the first 6 entries.

```r
head(ehrlen3ss$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.10390398
#> 2      1             Sdl              Sd              2              1 0.01044612
#> 3      1            Sz1r              Sd             13              1 0.00000000
#> 4      1            Sz2r              Sd             14              1 0.00000000
#> 5      1            Sz3r              Sd             15              1 0.00000000
#> 6      1            Sz4r              Sd             16              1 0.00000000
```
  
Here we see the historical portion of the stable stage distribution, where the distribution of stage pairs is given. If we focus on the `$ahist` elements, then we can see the distribution in terms of the original stages shown in our stageframe.

```r
ehrlen3ss$ahist
#>     matrix new_stage_id      ss_prop
#> 1        1            1 2.904905e-01
#> 2        1            2 3.965097e-02
#> 3        1            3 5.086115e-02
#> 4        1            4 1.378194e-02
#> 5        1            5 5.192587e-02
#> 6        1            6 8.452315e-02
#> 7        1            7 1.315953e-01
#> 8        1            8 1.389902e-01
#> 9        1            9 8.011630e-02
#> 10       1           10 2.286952e-02
#> 11       1           11 3.114642e-03
#> 12       1           12 2.080753e-04
#> 13       1           13 6.216238e-05
#> 14       1           14 4.901771e-04
#> 15       1           15 3.176462e-03
#> 16       1           16 1.323398e-02
#> 17       1           17 2.881058e-02
#> 18       1           18 2.938736e-02
#> 19       1           19 1.354000e-02
#> 20       1           20 2.873936e-03
#> 21       1           21 2.978350e-04
#> 22       2            1 3.293560e-01
#> 23       2            2 4.460362e-02
#> 24       2            3 7.503005e-02
#> 25       2            4 1.377543e-02
#> 26       2            5 4.929159e-02
#> 27       2            6 6.958748e-02
#> 28       2            7 1.047371e-01
#> 29       2            8 1.192559e-01
#> 30       2            9 7.994358e-02
#> 31       2           10 2.798095e-02
#> 32       2           11 4.852653e-03
#> 33       2           12 4.274685e-04
#> 34       2           13 3.581766e-05
#> 35       2           14 2.641169e-04
#> 36       2           15 1.733290e-03
#> 37       2           16 8.185246e-03
#> 38       2           17 2.169553e-02
#> 39       2           18 2.779340e-02
#> 40       2           19 1.633745e-02
#> 41       2           20 4.497577e-03
#> 42       2           21 6.157321e-04
#> 43       3            1 3.030531e-01
#> 44       3            2 4.112462e-02
#> 45       3            3 7.068055e-02
#> 46       3            4 1.282831e-02
#> 47       3            5 4.698760e-02
#> 48       3            6 6.996636e-02
#> 49       3            7 1.101013e-01
#> 50       3            8 1.293136e-01
#> 51       3            9 8.839602e-02
#> 52       3           10 3.117940e-02
#> 53       3           11 5.404644e-03
#> 54       3           12 4.730567e-04
#> 55       3           13 3.586078e-05
#> 56       3           14 2.750070e-04
#> 57       3           15 1.872245e-03
#> 58       3           16 9.062924e-03
#> 59       3           17 2.426829e-02
#> 60       3           18 3.111277e-02
#> 61       3           19 1.821202e-02
#> 62       3           20 4.977748e-03
#> 63       3           21 6.745912e-04
#> 64       4            1 3.013603e-01
#> 65       4            2 4.056705e-02
#> 66       4            3 3.639390e-02
#> 67       4            4 1.028321e-02
#> 68       4            5 4.709258e-02
#> 69       4            6 8.124459e-02
#> 70       4            7 1.384844e-01
#> 71       4            8 1.543841e-01
#> 72       4            9 9.033601e-02
#> 73       4           10 2.525692e-02
#> 74       4           11 3.256448e-03
#> 75       4           12 2.005654e-04
#> 76       4           13 3.552841e-05
#> 77       4           14 3.249514e-04
#> 78       4           15 2.362749e-03
#> 79       4           16 1.042560e-02
#> 80       4           17 2.296365e-02
#> 81       4           18 2.288661e-02
#> 82       4           19 9.995029e-03
#> 83       4           20 1.960375e-03
#> 84       4           21 1.854356e-04
#> 85       5            1 1.906478e-01
#> 86       5            2 2.537448e-02
#> 87       5            3 1.129607e-01
#> 88       5            4 2.881927e-02
#> 89       5            5 7.100550e-02
#> 90       5            6 1.096652e-01
#> 91       5            7 1.496255e-01
#> 92       5            8 1.468541e-01
#> 93       5            9 8.275200e-02
#> 94       5           10 2.325891e-02
#> 95       5           11 3.081804e-03
#> 96       5           12 1.989711e-04
#> 97       5           13 5.164441e-05
#> 98       5           14 3.351805e-04
#> 99       5           15 1.931907e-03
#> 100      5           16 7.898268e-03
#> 101      5           17 1.749632e-02
#> 102      5           18 1.801820e-02
#> 103      5           19 8.182190e-03
#> 104      5           20 1.674532e-03
#> 105      5           21 1.674055e-04
#> 106      6            1 4.081751e-01
#> 107      6            2 5.607100e-02
#> 108      6            3 2.732529e-02
#> 109      6            4 4.710593e-03
#> 110      6            5 3.613981e-02
#> 111      6            6 4.651838e-02
#> 112      6            7 7.993770e-02
#> 113      6            8 1.018675e-01
#> 114      6            9 7.330415e-02
#> 115      6           10 2.765458e-02
#> 116      6           11 5.263497e-03
#> 117      6           12 5.111613e-04
#> 118      6           13 3.269480e-05
#> 119      6           14 3.127559e-04
#> 120      6           15 2.426999e-03
#> 121      6           16 1.222737e-02
#> 122      6           17 3.352641e-02
#> 123      6           18 4.515433e-02
#> 124      6           19 2.875833e-02
#> 125      6           20 8.764769e-03
#> 126      6           21 1.317572e-03
#> 127      7            1 3.101625e-01
#> 128      7            2 4.203916e-02
#> 129      7            3 5.975082e-02
#> 130      7            4 1.277985e-02
#> 131      7            5 4.894671e-02
#> 132      7            6 7.480581e-02
#> 133      7            7 1.171315e-01
#> 134      7            8 1.315670e-01
#> 135      7            9 8.358260e-02
#> 136      7           10 2.696603e-02
#> 137      7           11 4.235194e-03
#> 138      7           12 3.329049e-04
#> 139      7           13 4.378201e-05
#> 140      7           14 3.418966e-04
#> 141      7           15 2.289362e-03
#> 142      7           16 1.037110e-02
#> 143      7           17 2.532016e-02
#> 144      7           18 2.941429e-02
#> 145      7           19 1.559678e-02
#> 146      7           20 3.852606e-03
#> 147      7           21 4.699761e-04
```

These are quite different values, showing the impact of history once again.

  <p style="text-indent: 20px">Now let's take a peek at the reproductive values. Let's start with the ahistorical case.</p>

```r
repvalue3(lathmat2lnmean)
#>     matrix new_stage_id orig_stage_id original_size left_vector rep_value
#> 1        1            1            Sd           0.0  -0.0090146   1.00000
#> 2        1            2           Sdl           4.6  -0.1030660  11.43323
#> 3        1            3          Dorm           0.0  -0.1318171  14.62262
#> 4        1            4         Sz1nr           1.0  -0.1345282  14.92337
#> 5        1            5         Sz2nr           2.0  -0.1367046  15.16480
#> 6        1            6         Sz3nr           3.0  -0.1391253  15.43333
#> 7        1            7         Sz4nr           4.0  -0.1426341  15.82257
#> 8        1            8         Sz5nr           5.0  -0.1485511  16.47895
#> 9        1            9         Sz6nr           6.0  -0.1588004  17.61591
#> 10       1           10         Sz7nr           7.0  -0.1752023  19.43539
#> 11       1           11         Sz8nr           8.0  -0.1982553  21.99269
#> 12       1           12         Sz9nr           9.0  -0.2266261  25.13990
#> 13       1           13          Sz1r           1.0  -0.1340193  14.86692
#> 14       1           14          Sz2r           2.0  -0.1386382  15.37930
#> 15       1           15          Sz3r           3.0  -0.1443168  16.00923
#> 16       1           16          Sz4r           4.0  -0.1544172  17.12968
#> 17       1           17          Sz5r           5.0  -0.1741441  19.31801
#> 18       1           18          Sz6r           6.0  -0.2121657  23.53579
#> 19       1           19          Sz7r           7.0  -0.2814259  31.21890
#> 20       1           20          Sz8r           8.0  -0.3978079  44.12929
#> 21       1           21          Sz9r           9.0  -0.5763184  63.93167
#> 22       2            1            Sd           0.0   0.0093332   1.00000
#> 23       2            2           Sdl           4.6   0.0988024  10.58612
#> 24       2            3          Dorm           0.0   0.1201485  12.87324
#> 25       2            4         Sz1nr           1.0   0.1227827  13.15548
#> 26       2            5         Sz2nr           2.0   0.1248224  13.37402
#> 27       2            6         Sz3nr           3.0   0.1270459  13.61226
#> 28       2            7         Sz4nr           4.0   0.1302921  13.96007
#> 29       2            8         Sz5nr           5.0   0.1358690  14.55760
#> 30       2            9         Sz6nr           6.0   0.1457258  15.61370
#> 31       2           10         Sz7nr           7.0   0.1617930  17.33521
#> 32       2           11         Sz8nr           8.0   0.1847227  19.79200
#> 33       2           12         Sz9nr           9.0   0.2132903  22.85286
#> 34       2           13          Sz1r           1.0   0.1262342  13.52529
#> 35       2           14          Sz2r           2.0   0.1306857  14.00224
#> 36       2           15          Sz3r           3.0   0.1361939  14.59241
#> 37       2           16          Sz4r           4.0   0.1462195  15.66660
#> 38       2           17          Sz5r           5.0   0.1664015  17.82899
#> 39       2           18          Sz6r           6.0   0.2064649  22.12156
#> 40       2           19          Sz7r           7.0   0.2813907  30.14943
#> 41       2           20          Sz8r           8.0   0.4103551  43.96725
#> 42       2           21          Sz9r           9.0   0.6128799  65.66664
#> 43       3            1            Sd           0.0   0.0095307   1.00000
#> 44       3            2           Sdl           4.6   0.1047140  10.98702
#> 45       3            3          Dorm           0.0   0.1309757  13.74251
#> 46       3            4         Sz1nr           1.0   0.1333239  13.98889
#> 47       3            5         Sz2nr           2.0   0.1353931  14.20600
#> 48       3            6         Sz3nr           3.0   0.1378376  14.46248
#> 49       3            7         Sz4nr           4.0   0.1414246  14.83885
#> 50       3            8         Sz5nr           5.0   0.1474539  15.47147
#> 51       3            9         Sz6nr           6.0   0.1579366  16.57135
#> 52       3           10         Sz7nr           7.0   0.1748870  18.34986
#> 53       3           11         Sz8nr           8.0   0.1988899  20.86834
#> 54       3           12         Sz9nr           9.0   0.2282753  23.95158
#> 55       3           13          Sz1r           1.0   0.1357881  14.24744
#> 56       3           14          Sz2r           2.0   0.1396581  14.65350
#> 57       3           15          Sz3r           3.0   0.1447880  15.19175
#> 58       3           16          Sz4r           4.0   0.1544211  16.20249
#> 59       3           17          Sz5r           5.0   0.1739107  18.24742
#> 60       3           18          Sz6r           6.0   0.2123649  22.28219
#> 61       3           19          Sz7r           7.0   0.2830692  29.70078
#> 62       3           20          Sz8r           8.0   0.4004786  42.01985
#> 63       3           21          Sz9r           9.0   0.5737043  60.19540
#> 64       4            1            Sd           0.0   0.0091859   1.00000
#> 65       4            2           Sdl           4.6   0.0936908  10.19941
#> 66       4            3          Dorm           0.0   0.1105711  12.03705
#> 67       4            4         Sz1nr           1.0   0.1135397  12.36022
#> 68       4            5         Sz2nr           2.0   0.1156691  12.59203
#> 69       4            6         Sz3nr           3.0   0.1177915  12.82308
#> 70       4            7         Sz4nr           4.0   0.1207450  13.14460
#> 71       4            8         Sz5nr           5.0   0.1257888  13.69368
#> 72       4            9         Sz6nr           6.0   0.1349549  14.69153
#> 73       4           10         Sz7nr           7.0   0.1506832  16.40375
#> 74       4           11         Sz8nr           8.0   0.1744973  18.99621
#> 75       4           12         Sz9nr           9.0   0.2058241  22.40653
#> 76       4           13          Sz1r           1.0   0.1168015  12.71530
#> 77       4           14          Sz2r           2.0   0.1226285  13.34964
#> 78       4           15          Sz3r           3.0   0.1295179  14.09964
#> 79       4           16          Sz4r           4.0   0.1408020  15.32806
#> 80       4           17          Sz5r           5.0   0.1620299  17.63898
#> 81       4           18          Sz6r           6.0   0.2031427  22.11462
#> 82       4           19          Sz7r           7.0   0.2804364  30.52901
#> 83       4           20          Sz8r           8.0   0.4169110  45.38597
#> 84       4           21          Sz9r           9.0   0.6398222  69.65264
#> 85       5            1            Sd           0.0  -0.0111981   1.00000
#> 86       5            2           Sdl           4.6  -0.1116679   9.97204
#> 87       5            3          Dorm           0.0  -0.1309525  11.69417
#> 88       5            4         Sz1nr           1.0  -0.1341790  11.98230
#> 89       5            5         Sz2nr           2.0  -0.1361570  12.15894
#> 90       5            6         Sz3nr           3.0  -0.1378181  12.30728
#> 91       5            7         Sz4nr           4.0  -0.1399680  12.49926
#> 92       5            8         Sz5nr           5.0  -0.1436470  12.82780
#> 93       5            9         Sz6nr           6.0  -0.1504446  13.43483
#> 94       5           10         Sz7nr           7.0  -0.1622292  14.48721
#> 95       5           11         Sz8nr           8.0  -0.1801128  16.08423
#> 96       5           12         Sz9nr           9.0  -0.2035902  18.18078
#> 97       5           13          Sz1r           1.0  -0.1402118  12.52103
#> 98       5           14          Sz2r           2.0  -0.1453164  12.97688
#> 99       5           15          Sz3r           3.0  -0.1506177  13.45029
#> 100      5           16          Sz4r           4.0  -0.1594040  14.23491
#> 101      5           17          Sz5r           5.0  -0.1766314  15.77334
#> 102      5           18          Sz6r           6.0  -0.2110002  18.84250
#> 103      5           19          Sz7r           7.0  -0.2771176  24.74684
#> 104      5           20          Sz8r           8.0  -0.3962593  35.38630
#> 105      5           21          Sz9r           9.0  -0.5943450  53.07552
#> 106      6            1            Sd           0.0   0.0092260   1.00000
#> 107      6            2           Sdl           4.6   0.1039461  11.26665
#> 108      6            3          Dorm           0.0   0.1307526  14.17219
#> 109      6            4         Sz1nr           1.0   0.1338417  14.50701
#> 110      6            5         Sz2nr           2.0   0.1364577  14.79056
#> 111      6            6         Sz3nr           3.0   0.1394728  15.11736
#> 112      6            7         Sz4nr           4.0   0.1438862  15.59573
#> 113      6            8         Sz5nr           5.0   0.1511659  16.38477
#> 114      6            9         Sz6nr           6.0   0.1631110  17.67949
#> 115      6           10         Sz7nr           7.0   0.1809620  19.61435
#> 116      6           11         Sz8nr           8.0   0.2044965  22.16524
#> 117      6           12         Sz9nr           9.0   0.2321974  25.16772
#> 118      6           13          Sz1r           1.0   0.1330824  14.42471
#> 119      6           14          Sz2r           2.0   0.1384249  15.00378
#> 120      6           15          Sz3r           3.0   0.1450823  15.72537
#> 121      6           16          Sz4r           4.0   0.1565517  16.96853
#> 122      6           17          Sz5r           5.0   0.1780228  19.29577
#> 123      6           18          Sz6r           6.0   0.2174579  23.57012
#> 124      6           19          Sz7r           7.0   0.2858893  30.98735
#> 125      6           20          Sz8r           8.0   0.3966273  42.99017
#> 126      6           21          Sz9r           9.0   0.5629653  61.01943
#> 127      7            1            Sd           0.0  -0.0095329   1.00000
#> 128      7            2           Sdl           4.6  -0.1023998  10.74173
#> 129      7            3          Dorm           0.0  -0.1255361  13.16872
#> 130      7            4         Sz1nr           1.0  -0.1283692  13.46591
#> 131      7            5         Sz2nr           2.0  -0.1305512  13.69480
#> 132      7            6         Sz3nr           3.0  -0.1328951  13.94068
#> 133      7            7         Sz4nr           4.0  -0.1362611  14.29377
#> 134      7            8         Sz5nr           5.0  -0.1419483  14.89036
#> 135      7            9         Sz6nr           6.0  -0.1518504  15.92909
#> 136      7           10         Sz7nr           7.0  -0.1678411  17.60651
#> 137      7           11         Sz8nr           8.0  -0.1905847  19.99231
#> 138      7           12         Sz9nr           9.0  -0.2189275  22.96547
#> 139      7           13          Sz1r           1.0  -0.1306521  13.70539
#> 140      7           14          Sz2r           2.0  -0.1355285  14.21692
#> 141      7           15          Sz3r           3.0  -0.1414061  14.83348
#> 142      7           16          Sz4r           4.0  -0.1516681  15.90996
#> 143      7           17          Sz5r           5.0  -0.1716750  18.00869
#> 144      7           18          Sz6r           6.0  -0.2104921  22.08059
#> 145      7           19          Sz7r           7.0  -0.2819936  29.58109
#> 146      7           20          Sz8r           8.0  -0.4039580  42.37514
#> 147      7           21          Sz9r           9.0  -0.5944450  62.35721
```

As before, these are similar to those estimated by `popbio`'s `reproductive.value()` function. Let's now see the historical case.

```r
repvalue3(lathmat3lnmean)$ahist
#>     matrix new_stage_id rep_value
#> 1        1            1 0.0000000
#> 2        1            2 1.0000000
#> 3        1            3 1.3975164
#> 4        1            4 1.2280362
#> 5        1            5 1.6324483
#> 6        1            6 1.7964686
#> 7        1            7 1.9803999
#> 8        1            8 2.1754203
#> 9        1            9 2.3944457
#> 10       1           10 2.6305321
#> 11       1           11 2.8495758
#> 12       1           12 3.0285688
#> 13       1           13 0.7514644
#> 14       1           14 1.5057453
#> 15       1           15 1.9744964
#> 16       1           16 2.3254100
#> 17       1           17 2.6827236
#> 18       1           18 3.1367012
#> 19       1           19 3.8207087
#> 20       1           20 4.9436458
#> 21       1           21 6.7599245
#> 22       2            1 0.0000000
#> 23       2            2 1.0000000
#> 24       2            3 1.3593648
#> 25       2            4 1.1823853
#> 26       2            5 1.5876436
#> 27       2            6 1.7326227
#> 28       2            7 1.9338790
#> 29       2            8 2.1676129
#> 30       2            9 2.4290436
#> 31       2           10 2.7026341
#> 32       2           11 2.9480588
#> 33       2           12 3.1375191
#> 34       2           13 0.7929432
#> 35       2           14 1.5063771
#> 36       2           15 1.9954792
#> 37       2           16 2.4158040
#> 38       2           17 2.8508869
#> 39       2           18 3.3996557
#> 40       2           19 4.2045608
#> 41       2           20 5.4739947
#> 42       2           21 7.4070470
#> 43       3            1 0.0000000
#> 44       3            2 1.0000000
#> 45       3            3 1.3997350
#> 46       3            4 1.2207032
#> 47       3            5 1.6089864
#> 48       3            6 1.7477457
#> 49       3            7 1.9315053
#> 50       3            8 2.1366035
#> 51       3            9 2.3606762
#> 52       3           10 2.5934273
#> 53       3           11 2.8021534
#> 54       3           12 2.9630963
#> 55       3           13 0.8069945
#> 56       3           14 1.5154643
#> 57       3           15 1.9675418
#> 58       3           16 2.3351232
#> 59       3           17 2.7069953
#> 60       3           18 3.1710132
#> 61       3           19 3.8462391
#> 62       3           20 4.9022245
#> 63       3           21 6.5003847
#> 64       4            1 0.0000000
#> 65       4            2 1.0000000
#> 66       4            3 1.3232467
#> 67       4            4 1.1109167
#> 68       4            5 1.5523237
#> 69       4            6 1.7325941
#> 70       4            7 1.9558045
#> 71       4            8 2.1877438
#> 72       4            9 2.4442189
#> 73       4           10 2.7298066
#> 74       4           11 3.0083144
#> 75       4           12 3.2403892
#> 76       4           13 0.7515776
#> 77       4           14 1.5355680
#> 78       4           15 2.0649666
#> 79       4           16 2.4791483
#> 80       4           17 2.9037197
#> 81       4           18 3.4406105
#> 82       4           19 4.2516446
#> 83       4           20 5.6101633
#> 84       4           21 7.8839651
#> 85       5            1 0.0000000
#> 86       5            2 1.0000000
#> 87       5            3 1.2126869
#> 88       5            4 1.1165393
#> 89       5            5 1.4812835
#> 90       5            6 1.6630908
#> 91       5            7 1.8424863
#> 92       5            8 2.0222644
#> 93       5            9 2.2126141
#> 94       5           10 2.4197105
#> 95       5           11 2.6196204
#> 96       5           12 2.7824580
#> 97       5           13 0.7070892
#> 98       5           14 1.4426179
#> 99       5           15 1.9367040
#> 100      5           16 2.3164916
#> 101      5           17 2.6671235
#> 102      5           18 3.0778297
#> 103      5           19 3.6846194
#> 104      5           20 4.7133177
#> 105      5           21 6.4460108
#> 106      6            1 0.0000000
#> 107      6            2 1.0000000
#> 108      6            3 1.5342309
#> 109      6            4 1.2560986
#> 110      6            5 1.7124331
#> 111      6            6 1.8297668
#> 112      6            7 2.0439000
#> 113      6            8 2.3030656
#> 114      6            9 2.5836872
#> 115      6           10 2.8541081
#> 116      6           11 3.0797718
#> 117      6           12 3.2469241
#> 118      6           13 0.8448615
#> 119      6           14 1.5376210
#> 120      6           15 1.9815261
#> 121      6           16 2.3676317
#> 122      6           17 2.7966272
#> 123      6           18 3.3396137
#> 124      6           19 4.0912334
#> 125      6           20 5.1742827
#> 126      6           21 6.7107023
#> 127      7            1 0.0000000
#> 128      7            2 1.0000000
#> 129      7            3 1.3671875
#> 130      7            4 1.1767219
#> 131      7            5 1.5898560
#> 132      7            6 1.7477774
#> 133      7            7 1.9508850
#> 134      7            8 2.1727877
#> 135      7            9 2.4152731
#> 136      7           10 2.6688684
#> 137      7           11 2.8994148
#> 138      7           12 3.0819562
#> 139      7           13 0.7678584
#> 140      7           14 1.5020675
#> 141      7           15 1.9849468
#> 142      7           16 2.3741501
#> 143      7           17 2.7732833
#> 144      7           18 3.2729429
#> 145      7           19 4.0048979
#> 146      7           20 5.1691564
#> 147      7           21 6.9866105
```

Once again, we see strong differences here. We leave it to the reader to consider the importance of such differences.


### Analysis 3: Integral projection models (IPMs)

  <p style="text-indent: 20px">In this analysis, we will build an integral projection model. This is essentially a form of function-derived matrix, but one in which size must be modeled as Gaussian, and in which the size distribution is broken up into many fine-scale classes or size bins. Although the number of these bins varies, package `lefko3` uses a default of 100 (this can be changed as an option). Note also that, because of the high dimensionality of IPMs, we will not create a historical version. The ahistorical version under the default setting will have at least 10,000 elements per matrix (dimensions of 100 by 100 at minimum), while a historical version would have at least 100,000,000 - clearly such a matrix would be overparameterized!</p>
  
#### Step 1. Data characterization and reorganization

  <p style="text-indent: 20px">To begin, we need to create a stageframe for this dataset, in which we identify the key stages used in analysis. We will base all of this on Ehrlén [-@ehrlen_dynamics_2000], but modify the size bits to allow IPM construction and make all mature stages other than vegetative dormancy reproductive.</p>
  
  <p style="text-indent: 20px">In the stageframe code below, we show that we want an IPM by choosing two stages that serve as the size limits for IPM size classification. These two size classes should have exactly the same characteristics in the stageframe table other than size. We mark these in the vector that we load into the `stagenames` option using the string `ipm`. Package `lefko3` will then rename all IPM size classes according to its own conventions. Although we leave the number of size classes to the default setting here (100 bins), we may alter that using the `ipmbins` option, for example by setting `ipmbins = 25`, which would create 25 IPM size classes rather than 100.</p>

```r
sizevector <- c(0, 100, 0, 1, 7100)
stagevector <- c("Sd", "Sdl", "Dorm", "ipm", "ipm")
repvector <- c(0, 0, 0, 1, 1)
obsvector <- c(0, 1, 0, 1, 1)
matvector <- c(0, 0, 1, 1, 1)
immvector <- c(1, 1, 0, 0, 0)
propvector <- c(1, 0, 0, 0, 0)
indataset <- c(0, 1, 1, 1, 1)
binvec <- c(0, 100, 0.5, 1, 1)

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

  <p style="text-indent: 20px">We will also add some descriptive comments to this stageframe so that we know what each of these stages is.</p>

```r
lathframeipm$comments <- c("Dormant seed", "Seedling", "Dormant", rep("ipm stage", 100))
lathframeipm
#>                  stagenames       size repstatus obsstatus propstatus immstatus matstatus indataset binhalfwidth_raw sizebin_min
#> 1                        Sd    0.00000         0         0          1         1         0         0          0.00000     0.00000
#> 2                       Sdl  100.00000         0         1          0         1         0         1        100.00000     0.00000
#> 3                      Dorm    0.00000         0         0          0         0         1         1          0.50000    -0.50000
#> 4    sz35.85354 rp1 mt1 ob1   35.85354         1         1          0         0         1         1         35.67751     0.17603
#> 5   sz107.20855 rp1 mt1 ob1  107.20855         1         1          0         0         1         1         35.67751    71.53104
#> 6   sz178.56356 rp1 mt1 ob1  178.56356         1         1          0         0         1         1         35.67751   142.88606
#> 7   sz249.91858 rp1 mt1 ob1  249.91858         1         1          0         0         1         1         35.67751   214.24107
#> 8   sz321.27359 rp1 mt1 ob1  321.27359         1         1          0         0         1         1         35.67751   285.59609
#> 9   sz392.62861 rp1 mt1 ob1  392.62861         1         1          0         0         1         1         35.67751   356.95110
#> 10  sz463.98362 rp1 mt1 ob1  463.98362         1         1          0         0         1         1         35.67751   428.30612
#> 11  sz535.33864 rp1 mt1 ob1  535.33864         1         1          0         0         1         1         35.67751   499.66113
#> 12  sz606.69365 rp1 mt1 ob1  606.69365         1         1          0         0         1         1         35.67751   571.01615
#> 13  sz678.04867 rp1 mt1 ob1  678.04867         1         1          0         0         1         1         35.67751   642.37116
#> 14  sz749.40368 rp1 mt1 ob1  749.40368         1         1          0         0         1         1         35.67751   713.72618
#> 15   sz820.7587 rp1 mt1 ob1  820.75870         1         1          0         0         1         1         35.67751   785.08119
#> 16  sz892.11371 rp1 mt1 ob1  892.11371         1         1          0         0         1         1         35.67751   856.43621
#> 17  sz963.46873 rp1 mt1 ob1  963.46873         1         1          0         0         1         1         35.67751   927.79122
#> 18 sz1034.82374 rp1 mt1 ob1 1034.82374         1         1          0         0         1         1         35.67751   999.14624
#> 19 sz1106.17876 rp1 mt1 ob1 1106.17876         1         1          0         0         1         1         35.67751  1070.50125
#> 20 sz1177.53377 rp1 mt1 ob1 1177.53377         1         1          0         0         1         1         35.67751  1141.85626
#> 21 sz1248.88879 rp1 mt1 ob1 1248.88879         1         1          0         0         1         1         35.67751  1213.21128
#> 22  sz1320.2438 rp1 mt1 ob1 1320.24380         1         1          0         0         1         1         35.67751  1284.56629
#> 23 sz1391.59882 rp1 mt1 ob1 1391.59882         1         1          0         0         1         1         35.67751  1355.92131
#> 24 sz1462.95383 rp1 mt1 ob1 1462.95383         1         1          0         0         1         1         35.67751  1427.27632
#> 25 sz1534.30885 rp1 mt1 ob1 1534.30885         1         1          0         0         1         1         35.67751  1498.63134
#> 26 sz1605.66386 rp1 mt1 ob1 1605.66386         1         1          0         0         1         1         35.67751  1569.98635
#> 27 sz1677.01888 rp1 mt1 ob1 1677.01888         1         1          0         0         1         1         35.67751  1641.34137
#> 28 sz1748.37389 rp1 mt1 ob1 1748.37389         1         1          0         0         1         1         35.67751  1712.69638
#> 29 sz1819.72891 rp1 mt1 ob1 1819.72891         1         1          0         0         1         1         35.67751  1784.05140
#> 30 sz1891.08392 rp1 mt1 ob1 1891.08392         1         1          0         0         1         1         35.67751  1855.40641
#> 31 sz1962.43893 rp1 mt1 ob1 1962.43893         1         1          0         0         1         1         35.67751  1926.76143
#> 32 sz2033.79395 rp1 mt1 ob1 2033.79395         1         1          0         0         1         1         35.67751  1998.11644
#> 33 sz2105.14896 rp1 mt1 ob1 2105.14896         1         1          0         0         1         1         35.67751  2069.47146
#> 34 sz2176.50398 rp1 mt1 ob1 2176.50398         1         1          0         0         1         1         35.67751  2140.82647
#> 35 sz2247.85899 rp1 mt1 ob1 2247.85899         1         1          0         0         1         1         35.67751  2212.18149
#> 36 sz2319.21401 rp1 mt1 ob1 2319.21401         1         1          0         0         1         1         35.67751  2283.53650
#> 37 sz2390.56902 rp1 mt1 ob1 2390.56902         1         1          0         0         1         1         35.67751  2354.89152
#> 38 sz2461.92404 rp1 mt1 ob1 2461.92404         1         1          0         0         1         1         35.67751  2426.24653
#> 39 sz2533.27905 rp1 mt1 ob1 2533.27905         1         1          0         0         1         1         35.67751  2497.60155
#> 40 sz2604.63407 rp1 mt1 ob1 2604.63407         1         1          0         0         1         1         35.67751  2568.95656
#> 41 sz2675.98908 rp1 mt1 ob1 2675.98908         1         1          0         0         1         1         35.67751  2640.31158
#> 42  sz2747.3441 rp1 mt1 ob1 2747.34410         1         1          0         0         1         1         35.67751  2711.66659
#> 43 sz2818.69911 rp1 mt1 ob1 2818.69911         1         1          0         0         1         1         35.67751  2783.02160
#> 44 sz2890.05413 rp1 mt1 ob1 2890.05413         1         1          0         0         1         1         35.67751  2854.37662
#> 45 sz2961.40914 rp1 mt1 ob1 2961.40914         1         1          0         0         1         1         35.67751  2925.73163
#> 46 sz3032.76416 rp1 mt1 ob1 3032.76416         1         1          0         0         1         1         35.67751  2997.08665
#> 47 sz3104.11917 rp1 mt1 ob1 3104.11917         1         1          0         0         1         1         35.67751  3068.44166
#> 48 sz3175.47419 rp1 mt1 ob1 3175.47419         1         1          0         0         1         1         35.67751  3139.79668
#> 49  sz3246.8292 rp1 mt1 ob1 3246.82920         1         1          0         0         1         1         35.67751  3211.15169
#> 50 sz3318.18422 rp1 mt1 ob1 3318.18422         1         1          0         0         1         1         35.67751  3282.50671
#> 51 sz3389.53923 rp1 mt1 ob1 3389.53923         1         1          0         0         1         1         35.67751  3353.86172
#> 52 sz3460.89425 rp1 mt1 ob1 3460.89425         1         1          0         0         1         1         35.67751  3425.21674
#> 53 sz3532.24926 rp1 mt1 ob1 3532.24926         1         1          0         0         1         1         35.67751  3496.57175
#> 54 sz3603.60428 rp1 mt1 ob1 3603.60428         1         1          0         0         1         1         35.67751  3567.92677
#> 55 sz3674.95929 rp1 mt1 ob1 3674.95929         1         1          0         0         1         1         35.67751  3639.28178
#> 56  sz3746.3143 rp1 mt1 ob1 3746.31430         1         1          0         0         1         1         35.67751  3710.63680
#> 57 sz3817.66932 rp1 mt1 ob1 3817.66932         1         1          0         0         1         1         35.67751  3781.99181
#> 58 sz3889.02433 rp1 mt1 ob1 3889.02433         1         1          0         0         1         1         35.67751  3853.34683
#> 59 sz3960.37935 rp1 mt1 ob1 3960.37935         1         1          0         0         1         1         35.67751  3924.70184
#> 60 sz4031.73436 rp1 mt1 ob1 4031.73436         1         1          0         0         1         1         35.67751  3996.05686
#> 61 sz4103.08938 rp1 mt1 ob1 4103.08938         1         1          0         0         1         1         35.67751  4067.41187
#> 62 sz4174.44439 rp1 mt1 ob1 4174.44439         1         1          0         0         1         1         35.67751  4138.76689
#> 63 sz4245.79941 rp1 mt1 ob1 4245.79941         1         1          0         0         1         1         35.67751  4210.12190
#> 64 sz4317.15442 rp1 mt1 ob1 4317.15442         1         1          0         0         1         1         35.67751  4281.47692
#> 65 sz4388.50944 rp1 mt1 ob1 4388.50944         1         1          0         0         1         1         35.67751  4352.83193
#> 66 sz4459.86445 rp1 mt1 ob1 4459.86445         1         1          0         0         1         1         35.67751  4424.18695
#> 67 sz4531.21947 rp1 mt1 ob1 4531.21947         1         1          0         0         1         1         35.67751  4495.54196
#> 68 sz4602.57448 rp1 mt1 ob1 4602.57448         1         1          0         0         1         1         35.67751  4566.89697
#> 69  sz4673.9295 rp1 mt1 ob1 4673.92950         1         1          0         0         1         1         35.67751  4638.25199
#> 70 sz4745.28451 rp1 mt1 ob1 4745.28451         1         1          0         0         1         1         35.67751  4709.60700
#> 71 sz4816.63953 rp1 mt1 ob1 4816.63953         1         1          0         0         1         1         35.67751  4780.96202
#>    sizebin_max sizebin_center sizebin_width     comments
#> 1      0.00000        0.00000       0.00000 Dormant seed
#> 2    200.00000      100.00000     200.00000     Seedling
#> 3      0.50000        0.00000       1.00000      Dormant
#> 4     71.53104       35.85354      71.35501    ipm stage
#> 5    142.88606      107.20855      71.35502    ipm stage
#> 6    214.24107      178.56356      71.35501    ipm stage
#> 7    285.59609      249.91858      71.35502    ipm stage
#> 8    356.95110      321.27359      71.35501    ipm stage
#> 9    428.30612      392.62861      71.35502    ipm stage
#> 10   499.66113      463.98362      71.35501    ipm stage
#> 11   571.01615      535.33864      71.35502    ipm stage
#> 12   642.37116      606.69365      71.35501    ipm stage
#> 13   713.72618      678.04867      71.35502    ipm stage
#> 14   785.08119      749.40368      71.35501    ipm stage
#> 15   856.43621      820.75870      71.35502    ipm stage
#> 16   927.79122      892.11371      71.35501    ipm stage
#> 17   999.14624      963.46873      71.35502    ipm stage
#> 18  1070.50125     1034.82374      71.35501    ipm stage
#> 19  1141.85626     1106.17876      71.35501    ipm stage
#> 20  1213.21128     1177.53377      71.35502    ipm stage
#> 21  1284.56629     1248.88879      71.35501    ipm stage
#> 22  1355.92131     1320.24380      71.35502    ipm stage
#> 23  1427.27632     1391.59882      71.35501    ipm stage
#> 24  1498.63134     1462.95383      71.35502    ipm stage
#> 25  1569.98635     1534.30885      71.35501    ipm stage
#> 26  1641.34137     1605.66386      71.35502    ipm stage
#> 27  1712.69638     1677.01888      71.35501    ipm stage
#> 28  1784.05140     1748.37389      71.35502    ipm stage
#> 29  1855.40641     1819.72891      71.35501    ipm stage
#> 30  1926.76143     1891.08392      71.35502    ipm stage
#> 31  1998.11644     1962.43893      71.35501    ipm stage
#> 32  2069.47146     2033.79395      71.35502    ipm stage
#> 33  2140.82647     2105.14896      71.35501    ipm stage
#> 34  2212.18149     2176.50398      71.35502    ipm stage
#> 35  2283.53650     2247.85899      71.35501    ipm stage
#> 36  2354.89152     2319.21401      71.35502    ipm stage
#> 37  2426.24653     2390.56902      71.35501    ipm stage
#> 38  2497.60155     2461.92404      71.35502    ipm stage
#> 39  2568.95656     2533.27905      71.35501    ipm stage
#> 40  2640.31158     2604.63407      71.35502    ipm stage
#> 41  2711.66659     2675.98908      71.35501    ipm stage
#> 42  2783.02160     2747.34410      71.35501    ipm stage
#> 43  2854.37662     2818.69911      71.35502    ipm stage
#> 44  2925.73163     2890.05413      71.35501    ipm stage
#> 45  2997.08665     2961.40914      71.35502    ipm stage
#> 46  3068.44166     3032.76416      71.35501    ipm stage
#> 47  3139.79668     3104.11917      71.35502    ipm stage
#> 48  3211.15169     3175.47419      71.35501    ipm stage
#> 49  3282.50671     3246.82920      71.35502    ipm stage
#> 50  3353.86172     3318.18422      71.35501    ipm stage
#> 51  3425.21674     3389.53923      71.35502    ipm stage
#> 52  3496.57175     3460.89425      71.35501    ipm stage
#> 53  3567.92677     3532.24926      71.35502    ipm stage
#> 54  3639.28178     3603.60428      71.35501    ipm stage
#> 55  3710.63680     3674.95929      71.35502    ipm stage
#> 56  3781.99181     3746.31430      71.35501    ipm stage
#> 57  3853.34683     3817.66932      71.35502    ipm stage
#> 58  3924.70184     3889.02433      71.35501    ipm stage
#> 59  3996.05686     3960.37935      71.35502    ipm stage
#> 60  4067.41187     4031.73436      71.35501    ipm stage
#> 61  4138.76689     4103.08938      71.35502    ipm stage
#> 62  4210.12190     4174.44439      71.35501    ipm stage
#> 63  4281.47692     4245.79941      71.35502    ipm stage
#> 64  4352.83193     4317.15442      71.35501    ipm stage
#> 65  4424.18695     4388.50944      71.35502    ipm stage
#> 66  4495.54196     4459.86445      71.35501    ipm stage
#> 67  4566.89697     4531.21947      71.35501    ipm stage
#> 68  4638.25199     4602.57448      71.35502    ipm stage
#> 69  4709.60700     4673.92950      71.35501    ipm stage
#> 70  4780.96202     4745.28451      71.35502    ipm stage
#> 71  4852.31703     4816.63953      71.35501    ipm stage
#>  [ reached 'max' / getOption("max.print") -- omitted 32 rows ]
```

A look at this stageframe shows that the IPM portion technically starts with the fourth stage an keeps going until the 103rd stage. Stage names within this range are concatenations of the central size (designated with "sz"), reproductive status (designated with "rp"), maturity status (designated with "mt"), and observation status (designated with "ob"). The first three stages, which fall outside of the IPM classification, are left unaltered.

  <p style="text-indent: 20px">To work with this dataset, we first need to format the data into 'vertical' format, in which each row corresponds to the state of a single individual in two (if ahistorical) or three (if historical) consecutive time intervals. Because this is an IPM, we will need to estimate linear models of vital rates, and that will require that NAs are avoided in key terms used in estimation. For this purpose, we will set `NAas0 = TRUE`. We will also set `NRasRep = TRUE` because there are mature individuals that do not reproduce that we wish to include in reproductive stages (setting this option to TRUE makes sure that the reproductive status of non-reproductive individuals is set to 1, although the actual fecundity is not altered).</p>

```r
lathvertipm <- verticalize3(lathyrus, noyears = 4, firstyear = 1988, patchidcol = "SUBPLOT", 
                            individcol = "GENET", blocksize = 9, juvcol = "Seedling1988", 
                            size1col = "Volume88", repstr1col = "FCODE88",
                            fec1col = "Intactseed88", dead1col = "Dead1988", 
                            nonobs1col = "Dormant1988", stageassign = lathframeipm, 
                            stagesize = "sizea", censorcol = "Missing1988",censorkeep = NA, 
                            censor = TRUE, NAas0 = TRUE, NRasRep = TRUE)
summary(lathvertipm)
#>      rowid         popid         patchid    individ           year2          sizea1           sizea2           sizea3      
#>  Min.   :   1.0   Mode:logical   1:670   Min.   :  1.00   Min.   :1988   Min.   :   0.0   Min.   :   0.0   Min.   :   0.0  
#>  1st Qu.: 237.0   NA's:2515      2:201   1st Qu.: 39.00   1st Qu.:1988   1st Qu.:   0.0   1st Qu.:  12.6   1st Qu.:  10.0  
#>  Median : 523.0                  3:552   Median : 78.00   Median :1989   Median :   9.0   Median : 106.5   Median :  73.2  
#>  Mean   : 537.5                  4:469   Mean   : 92.82   Mean   :1989   Mean   : 386.6   Mean   : 483.1   Mean   : 391.5  
#>  3rd Qu.: 819.5                  5:305   3rd Qu.:139.50   3rd Qu.:1990   3rd Qu.: 621.7   3rd Qu.: 732.5   3rd Qu.: 545.0  
#>  Max.   :1118.0                  6:318   Max.   :276.00   Max.   :1990   Max.   :7032.0   Max.   :7032.0   Max.   :6645.8  
#>     repstra1         repstra2         repstra3          feca1             feca2            feca3          obsstatus1       obsstatus2    
#>  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   : 0.0000   Min.   : 0.000   Min.   : 0.000   Min.   :0.0000   Min.   :0.0000  
#>  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.: 0.0000   1st Qu.: 0.000   1st Qu.: 0.000   1st Qu.:0.0000   1st Qu.:1.0000  
#>  Median :0.0000   Median :0.0000   Median :0.0000   Median : 0.0000   Median : 0.000   Median : 0.000   Median :1.0000   Median :1.0000  
#>  Mean   :0.1801   Mean   :0.2382   Mean   :0.2191   Mean   : 0.9877   Mean   : 1.146   Mean   : 1.296   Mean   :0.5527   Mean   :0.9499  
#>  3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.: 0.0000   3rd Qu.: 0.000   3rd Qu.: 0.000   3rd Qu.:1.0000   3rd Qu.:1.0000  
#>  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :66.0000   Max.   :66.000   Max.   :66.000   Max.   :1.0000   Max.   :1.0000  
#>    obsstatus3       repstatus1       repstatus2       repstatus3       fecstatus1        fecstatus2       fecstatus3       firstseen   
#>  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.00000   Min.   :0.0000   Min.   :0.0000   Min.   :1988  
#>  1st Qu.:1.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.00000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:1988  
#>  Median :1.0000   Median :0.0000   Median :0.0000   Median :0.0000   Median :0.00000   Median :0.0000   Median :0.0000   Median :1988  
#>  Mean   :0.8386   Mean   :0.1801   Mean   :0.2382   Mean   :0.2191   Mean   :0.08946   Mean   :0.1058   Mean   :0.1121   Mean   :1988  
#>  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:0.00000   3rd Qu.:0.0000   3rd Qu.:0.0000   3rd Qu.:1988  
#>  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   :1.00000   Max.   :1.0000   Max.   :1.0000   Max.   :1990  
#>     lastseen        alive1           alive2      alive3           obsage        obslifespan       stage1             stage2         
#>  Min.   :1988   Min.   :0.0000   Min.   :1   Min.   :0.0000   Min.   :0.0000   Min.   :0.000   Length:2515        Length:2515       
#>  1st Qu.:1991   1st Qu.:0.0000   1st Qu.:1   1st Qu.:1.0000   1st Qu.:0.0000   1st Qu.:2.000   Class :character   Class :character  
#>  Median :1991   Median :1.0000   Median :1   Median :1.0000   Median :1.0000   Median :3.000   Mode  :character   Mode  :character  
#>  Mean   :1991   Mean   :0.5813   Mean   :1   Mean   :0.8887   Mean   :0.8258   Mean   :2.449                                        
#>  3rd Qu.:1991   3rd Qu.:1.0000   3rd Qu.:1   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:3.000                                        
#>  Max.   :1991   Max.   :1.0000   Max.   :1   Max.   :1.0000   Max.   :2.0000   Max.   :3.000                                        
#>     stage3            matstatus1       matstatus2       matstatus3
#>  Length:2515        Min.   :0.0000   Min.   :0.0000   Min.   :1   
#>  Class :character   1st Qu.:1.0000   1st Qu.:1.0000   1st Qu.:1   
#>  Mode  :character   Median :1.0000   Median :1.0000   Median :1   
#>                     Mean   :0.9368   Mean   :0.8883   Mean   :1   
#>                     3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.:1   
#>                     Max.   :1.0000   Max.   :1.0000   Max.   :1
```

  <p style="text-indent: 20px">Voila! Now we move on to create the extra bits of information needed for matrix estimation.</p>
  
#### Step 2: Develop supplemental information for matrix estimation

  <p style="text-indent: 20px">Now we will provide some given transitions. This is the same overwrite object as used before in ahistorical matrix estimation.</p>

```r
lathover2 <- overwrite(stage3 = c("Sd", "Sdl"), stage2 = c("Sd", "Sd"), givenrate = c(0.345, 0.054))
```

  <p style="text-indent: 20px">We will also create a reproductive matrix, which describes not only which stages are reproductive, but which stages they lead to the reproduction of, and at what level. This matrix is composed mostly of 0s, but fecundity is noted as non-zero entries equal to a modifier to the full fecundity estimated by R. This matrix has rows and columns equal to 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 `lathframeipm`. Then we modify elements corresponding to fecundity by the expected mean seed dormancy probability (row 1), and by the germination rate (row 2).</p>

```r
lathrepmipm <- matrix(0, 103, 103)
lathrepmipm[1, 4:103] <- 0.345
lathrepmipm[2, 4:103] <- 0.054
```

  <p style="text-indent: 20px">Let's now use the same rounding procedure as in Analysis 2 for fecundity in our vertical dataset, since we know that fecundity is all integers except for 6 data points that were estimated differently from the rest. Then we will take a peek at the degree to which 0's exist in fecundity.</p>

```r
lathvertipm$feca2 <- round(lathvertipm$feca2)
lathvertipm$feca1 <- round(lathvertipm$feca1)
lathvertipm$feca3 <- round(lathvertipm$feca3)

summary(subset(lathvertipm, repstatus2 == 1)$feca2)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   0.000   0.000   0.000   4.791   6.000  66.000
length(subset(lathvertipm, repstatus2 == 1)$feca2)
#> [1] 599
length(which(subset(lathvertipm, repstatus2 == 1)$feca2 > 0))
#> [1] 265
```
Nearly half of reproductive individuals do not actually reproduce in this example. Under the circumstance, we will not be able to use the Poisson or negative binomial distribution to model fecundity, and so must rely on the Gaussian.

#### Step 3: Develop models of demographic parameters

  <p style="text-indent: 20px">Integral projection models (IPMs) require functions of vital rates to populate them. Here, we will create these models using `modelsearch`. This looks similar to the `modelsearch` call in the last Analysis, but differs in that we are not including models of reproductive status. Note that we will only create ahistorical models, since IPMs already have high dimensionality.</p>

```r
lathmodels2ipm <- modelsearch(lathvertipm, historical = FALSE, approach = "lme4", suite = "size",
                              vitalrates = c("surv", "obs", "size", "fec"), juvestimate = "Sdl",
                              bestfit = "AICc&k", sizedist = "gaussian", fecdist = "gaussian",
                              indiv = "individ", patch = "patchid", year = "year2",
                              year.as.random = TRUE, patch.as.random = TRUE,
                              show.model.tables = TRUE)
#> 
#> 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.0278817 (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?;Model is nearly unidentifiable: large eigenvalue ratio
#>  - Rescale variables?
#> 
#> 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.0233345 (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?;Model is nearly unidentifiable: large eigenvalue ratio
#>  - Rescale variables?
#> 
#> Developing global model of size (Gaussian)...
#> boundary (singular) fit: see ?isSingular
#> 
#> Developing global model of fecundity (Gaussian)...
#> 
#> 
#> Developing global model of juvenile survival probability...
#> boundary (singular) fit: see ?isSingular
#> 
#> Developing global model of juvenile observation probability...
#> 
#> 
#> Developing global model of juvenile size (Gaussian)...
#> 
#> 
#> All global models developed.
#> 
#> 
#> Commencing dredge of survival probability...
#> Fixed term is "(Intercept)"
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0278817 (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?;Model is nearly unidentifiable: large eigenvalue ratio
#>  - Rescale variables?
#> 
#> Commencing dredge of observation probability...
#> Fixed term is "(Intercept)"
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0233345 (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?;Model is nearly unidentifiable: large eigenvalue ratio
#>  - Rescale variables?
#> 
#> Commencing dredge of size...
#> Warning in dredge(size.global.model, rank = used.criterion): comparing models fitted by REML
#> Fixed term is "(Intercept)"
#> boundary (singular) fit: see ?isSingular
#> 
#> Commencing dredge of fecundity...
#> Warning in dredge(fec.global.model, rank = used.criterion): comparing models fitted by REML
#> Fixed term is "(Intercept)"
#> 
#> Commencing dredge of juvenile survival probability...
#> Fixed term is "(Intercept)"
#> boundary (singular) fit: see ?isSingular
#> 
#> Commencing dredge of juvenile observation probability...
#> Fixed term is "(Intercept)"
#> 
#> Commencing dredge of juvenile size...
#> Warning in dredge(juv.size.global.model, rank = used.criterion): comparing models fitted by REML
#> Fixed term is "(Intercept)"
#> 
#> Finished dredging all vital rates.
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0278817 (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?;Model is nearly unidentifiable: large eigenvalue ratio
#>  - Rescale variables?
#> boundary (singular) fit: see ?isSingular
#> boundary (singular) fit: see ?isSingular
```

  <p style="text-indent: 20px">Now we can take a look at the summary of these models.</p>

```r
summary(lathmodels2ipm)
#> This LefkoMod object includes 7 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 ~ sizea2 + (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: surv.data
#>       AIC       BIC    logLik  deviance  df.resid 
#> 1283.1696 1311.7274 -636.5848 1273.1696      2229 
#> Random effects:
#>  Groups  Name        Std.Dev.
#>  individ (Intercept) 0.4130  
#>  patchid (Intercept) 0.3205  
#>  year2   (Intercept) 0.6041  
#> Number of obs: 2234, groups:  individ, 257; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)       sizea2  
#>   2.3072081    0.0006649  
#> convergence code 0; 0 optimizer warnings; 3 lme4 warnings 
#> 
#> ────────────────────────────────────────
#> 
#> Observation model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: binomial  ( logit )
#> Formula: obsstatus3 ~ (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: obs.data
#>       AIC       BIC    logLik  deviance  df.resid 
#>  775.5882  798.0494 -383.7941  767.5882      2025 
#> Random effects:
#>  Groups  Name        Std.Dev.
#>  individ (Intercept) 0.2226  
#>  patchid (Intercept) 0.5822  
#>  year2   (Intercept) 2.8777  
#> Number of obs: 2029, groups:  individ, 254; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)  
#>       4.409  
#> 
#> ────────────────────────────────────────
#> 
#> Size model:
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: sizea3 ~ sizea2 + (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: size.data
#> REML criterion at convergence: 29285.58
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  individ  (Intercept)   0.00  
#>  patchid  (Intercept)  57.13  
#>  year2    (Intercept) 210.37  
#>  Residual             502.47  
#> Number of obs: 1916, groups:  individ, 254; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)       sizea2  
#>    176.0197       0.6113  
#> convergence code 0; 0 optimizer warnings; 1 lme4 warnings 
#> 
#> ────────────────────────────────────────
#> 
#> Reproductive status model:
#> [1] 1
#> 
#> ────────────────────────────────────────
#> 
#> Fecundity model:
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: feca2 ~ sizea2 + (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: fec.data
#> REML criterion at convergence: 12729.05
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  individ  (Intercept) 0.7057  
#>  patchid  (Intercept) 0.1606  
#>  year2    (Intercept) 0.9824  
#>  Residual             4.1016  
#> Number of obs: 2234, groups:  individ, 257; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)       sizea2  
#>    -0.50396      0.00319  
#> 
#> ────────────────────────────────────────────────────────────────────────────────
#> Juvenile survival model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: binomial  ( logit )
#> Formula: alive3 ~ (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: juvsurv.data
#>       AIC       BIC    logLik  deviance  df.resid 
#>  333.7298  348.2833 -162.8649  325.7298       277 
#> Random effects:
#>  Groups  Name        Std.Dev.
#>  individ (Intercept) 0.0000  
#>  patchid (Intercept) 0.2104  
#>  year2   (Intercept) 0.0000  
#> Number of obs: 281, groups:  individ, 187; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)  
#>       1.001  
#> convergence code 0; 0 optimizer warnings; 1 lme4 warnings 
#> 
#> ────────────────────────────────────────
#> 
#> Juvenile observation model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#>  Family: binomial  ( logit )
#> Formula: obsstatus3 ~ (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: juvobs.data
#>      AIC      BIC   logLik deviance df.resid 
#>  98.1375 111.4490 -45.0688  90.1375      202 
#> Random effects:
#>  Groups  Name        Std.Dev.
#>  individ (Intercept) 1.1474  
#>  patchid (Intercept) 0.8774  
#>  year2   (Intercept) 1.1929  
#> Number of obs: 206, groups:  individ, 151; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)  
#>       4.005  
#> 
#> ────────────────────────────────────────
#> 
#> Juvenile size model:
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: sizea3 ~ (1 | individ) + (1 | year2) + (1 | patchid)
#>    Data: juvsize.data
#> REML criterion at convergence: 1303.003
#> Random effects:
#>  Groups   Name        Std.Dev.
#>  individ  (Intercept) 0.988829
#>  patchid  (Intercept) 0.001148
#>  year2    (Intercept) 1.182503
#>  Residual             6.997522
#> Number of obs: 193, groups:  individ, 144; patchid, 6; year2, 3
#> Fixed Effects:
#> (Intercept)  
#>       11.02  
#> 
#> ────────────────────────────────────────
#> 
#> Juvenile reproduction model:
#> [1] 1
#> 
#> 
#> ────────────────────────────────────────────────────────────────────────────────
#> 
#> Number of models in survival table:2
#> 
#> Number of models in observation table:2
#> 
#> Number of models in size table:2
#> 
#> Number of models in reproduction status table: 1
#> 
#> Number of models in fecundity table:2
#> 
#> 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     patchid
#> 4       surv3      alive3
#> 5        obs3  obsstatus3
#> 6       size3      sizea3
#> 7      repst3  repstatus3
#> 8        fec3       feca3
#> 9        fec2       feca2
#> 10      size2      sizea2
#> 11      size1        <NA>
#> 12     repst2  repstatus2
#> 13     repst1        <NA>
#> 14        age        <NA>
#> 
#> 
#> ────────────────────────────────────────────────────────────────────────────────
#> 
#> Quality control:
#> 
#> Survival estimated with 257 individuals and 2234 individual transitions.
#> Observation estimated with 254 individuals and 2029 individual transitions.
#> Size estimated with 254 individuals and 1916 individual transitions.
#> Reproduction probability not estimated.
#> Fecundity estimated with 257 individuals and 2234 individual transitions.
#> Juvenile survival estimated with 187 individuals and 281 individual transitions.
#> Juvenile observation estimated with 151 individuals and 206 individual transitions.
#> Juvenile size estimated with 144 individuals and 193 individual transitions.
#> Juvenile reproduction probability not estimated.
#> NULL
```

  <p style="text-indent: 20px">Let's move on now to the matrices themselves.</p>

#### Step 4. Estimate matrices

  <p style="text-indent: 20px">We will now create the suite of matrices covering the patches and years of study. Note that we will use the `negfec = TRUE` option. This option prevents fecundity estimates added to the matrix to be negative - a situation that can happen when fecundity is treated as Gaussian.</p>

```r
lathmat2ipm <- flefko2(stageframe = lathframeipm, repmatrix = lathrepmipm, 
                       modelsuite = lathmodels2ipm, overwrite = lathover2,
                       data = lathvertipm, year.as.random = FALSE, 
                       patch.as.random = FALSE, reduce = FALSE, negfec = TRUE)
summary(lathmat2ipm)
#> 
#> This lefkoMat object contains 18 matrices.
#> 
#> Each matrix is a square matrix with 103 rows and columns, and a total of 10609 elements.
#> A total of 183744 survival-transitions were estimated, with 10208 per matrix.
#> A total of 3600 fecundity transitions were estimated, with 200 per matrix.
#> 
#> Vital rate modeling quality control:
#> 
#> Survival estimated with 257 individuals and 2234 individual transitions.
#> Observation estimated with 254 individuals and 2029 individual transitions.
#> Size estimated with 254 individuals and 1916 individual transitions.
#> Reproductive status estimated with 0 individuals and 0 individual transitions.
#> Fecundity estimated with 257 individuals and 2234 individual transitions.
#> Juvenile survival estimated with 187 individuals and 281 individual transitions.
#> Juvenile observation estimated with 151 individuals and 206 individual transitions.
#> Juvenile size estimated with 144 individuals and 193 individual transitions.
#> Juvenile reproductive status estimated with 0 individuals and 0 individual transitions.
#> NULL
```

  <p style="text-indent: 20px">One thing to note is that the IPM is a large matrix with most elements estimated. This makes it likely to become overparameterized if the total sample size in the dataset and the number of stages needed are not properly considered. For example, here survival is estimated based on data for 257 individuals, yielding 2234 individual transitions actually used in modeling, while each matrix includes over 10,000 estimated elements, most of which are survival-transitions. Caution is advised.</p>
  
  <p style="text-indent: 20px">Let's take a look at the top-left corner of the first matrix (the matrix is too huge to inspect in full here).</p>

```r
lathmat2ipm$A[[1]][1:25,1:5]
#>        [,1]          [,2]        [,3]        [,4]        [,5]
#>  [1,] 0.345  0.000000e+00 0.000000000 0.083054724 0.104285820
#>  [2,] 0.054  0.000000e+00 0.000000000 0.012999870 0.016322998
#>  [3,] 0.000  1.308258e-02 0.087052521 0.087123573 0.087260363
#>  [4,] 0.000  5.390777e-03 0.044957756 0.044076854 0.042133183
#>  [5,] 0.000  2.743738e-41 0.047447401 0.046806766 0.045297724
#>  [6,] 0.000 9.682695e-125 0.049075172 0.048713382 0.047727656
#>  [7,] 0.000 2.369143e-253 0.049745417 0.049685513 0.049283968
#>  [8,] 0.000  0.000000e+00 0.049418087 0.049665278 0.049874992
#>  [9,] 0.000  0.000000e+00 0.048112799 0.048653917 0.049465437
#> [10,] 0.000  0.000000e+00 0.045906790 0.046711558 0.048079779
#> [11,] 0.000  0.000000e+00 0.042927448 0.043951403 0.045799941
#> [12,] 0.000  0.000000e+00 0.039340042 0.040528706 0.042757173
#> [13,] 0.000  0.000000e+00 0.035332666 0.036626428 0.039119644
#> [14,] 0.000  0.000000e+00 0.031099944 0.032439041 0.035076997
#> [15,] 0.000  0.000000e+00 0.026827774 0.028156800 0.030824196
#> [16,] 0.000  0.000000e+00 0.022680430 0.023951912 0.026546222
#> [17,] 0.000  0.000000e+00 0.018791427 0.019968199 0.022405545
#> [18,] 0.000  0.000000e+00 0.015258430 0.016314706 0.018533179
#> [19,] 0.000  0.000000e+00 0.012142321 0.013063554 0.015024018
#> [20,] 0.000  0.000000e+00 0.009469682 0.010251449 0.011936145
#> [21,] 0.000  0.000000e+00 0.007237868 0.007884075 0.009293593
#> [22,] 0.000  0.000000e+00 0.005421604 0.005942348 0.007091614
#> [23,] 0.000  0.000000e+00 0.003980031 0.004389419 0.005303323
#> [24,] 0.000  0.000000e+00 0.002863433 0.003177590 0.003886807
#> [25,] 0.000  0.000000e+00 0.002018967 0.002254396 0.002791769
```

This matrix is very large, of course, so is difficult to read properly. We can get another handle on quality control by checking the column sums of the first U matrix, to make sure that all column sums look like survival probabilities.

```r
summary(colSums(lathmat2ipm$U[[1]]))
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> 0.01847 0.98280 0.99642 0.95570 0.99893 0.99968
```

  <p style="text-indent: 20px">Let's estimate the mean IPM matrices. This code will estimate mean matrices for each of the 6 patches, followed by a grand mean.</p>

```r
lathipmmean <- lmean(lathmat2ipm)
summary(lathipmmean)
#> 
#> This lefkoMat object contains 7 matrices.
#> 
#> Each matrix is a square matrix with 103 rows and columns, and a total of 10609 elements.
#> A total of 71456 survival-transitions were estimated, with 10208 per matrix.
#> A total of 1400 fecundity transitions were estimated, with 200 per matrix.
#> 
#> Vital rate modeling quality control:
#> 
#> Survival estimated with 257 individuals and 2234 individual transitions.
#> Observation estimated with 254 individuals and 2029 individual transitions.
#> Size estimated with 254 individuals and 1916 individual transitions.
#> Reproductive status estimated with 0 individuals and 0 individual transitions.
#> Fecundity estimated with 257 individuals and 2234 individual transitions.
#> Juvenile survival estimated with 187 individuals and 281 individual transitions.
#> Juvenile observation estimated with 151 individuals and 206 individual transitions.
#> Juvenile size estimated with 144 individuals and 193 individual transitions.
#> Juvenile reproductive status estimated with 0 individuals and 0 individual transitions.
#> NULL
```

  <p style="text-indent: 20px">As a final check, let's take a look at the column sums of the grand mean survival-transition matrix.</p>

```r
colSums(lathipmmean$U[[7]])
#>   [1] 0.39900000 0.01847336 0.58222244 0.59650891 0.62444329 0.65158362 0.67778757 0.70293110 0.72691001 0.74964093 0.77106158 0.79113064
#>  [13] 0.80982697 0.82714842 0.84311028 0.85774340 0.87109210 0.88321195 0.89416756 0.90403028 0.91287611 0.92078370 0.92783256 0.93410146
#>  [25] 0.93966710 0.94460301 0.94897869 0.95285900 0.95630373 0.95936739 0.96209915 0.96454292 0.96673756 0.96871712 0.97051124 0.97214548
#>  [37] 0.97364172 0.97501860 0.97629188 0.97747483 0.97857856 0.97961236 0.98058399 0.98149991 0.98236550 0.98318531 0.98396311 0.98470215
#>  [49] 0.98540517 0.98607454 0.98671234 0.98732040 0.98790034 0.98845364 0.98898163 0.98948556 0.98996657 0.99042573 0.99086405 0.99128249
#>  [61] 0.99168193 0.99206325 0.99242725 0.99277470 0.99310636 0.99342292 0.99372507 0.99401345 0.99428868 0.99455135 0.99480202 0.99504124
#>  [73] 0.99526952 0.99548736 0.99569522 0.99589357 0.99608282 0.99626340 0.99643569 0.99660007 0.99675691 0.99690654 0.99704930 0.99718549
#>  [85] 0.99731543 0.99743938 0.99755763 0.99767044 0.99777805 0.99788070 0.99797863 0.99807204 0.99816115 0.99824614 0.99832722 0.99840455
#>  [97] 0.99847832 0.99854867 0.99861578 0.99867977 0.99874080 0.99879899 0.99885446
summary(colSums(lathipmmean$U[[7]]))
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> 0.01847 0.94679 0.98732 0.93055 0.99617 0.99885
```

All looks fine!

#### Step 5. Analyses

  <p style="text-indent: 20px">Now let's estimate the deterministic population growth rate for the IPM analysis.</p>

```r
lambda3(lathipmmean)
#>   pop patch    lambda
#> 1   1     1 0.8720046
#> 2   1     2 0.8867581
#> 3   1     3 0.8452272
#> 4   1     4 0.7979244
#> 5   1     5 0.7709210
#> 6   1     6 0.9503947
#> 7   1   All 0.8974453
```

These lambda estimates are in line with those already estimated, though perhaps most so with the original lambda estimates from Analysis 1. So, all seems well. Let's also compare against lambda estimates evaluated with the geometric mean matrix, since there should not be any singleton 0's causing mean matrix elements to equal 0.

```r
lambda3(lmean(lathmat2ipm, time = "geometric"))
#>   pop patch    lambda
#> 1   1     1 0.8088286
#> 2   1     2 0.7596964
#> 3   1     3 0.7539980
#> 4   1     4 0.8196549
#> 5   1     5 0.7115295
#> 6   1     6 0.8512272
#> 7   1   All 0.7658357
```

The geometric mean matrix lambdas are substantially lower, suggesting a strong influence of our assumptions about time. Though we are aware that the typical population matrix study would not use the geometric mean approach, we would encourage its consideration in this case, because of the strong chance that inference of population viability would be too 'rosy' using the arithmetic approach.

  <p style="text-indent: 20px">Now let's look at the stable stage distribution. We will only look at the summary and first 6 entries, given the size of the output.</p>

```r
summary(stablestage3(lathipmmean))
#>      matrix   new_stage_id orig_stage_id      original_size     ss_prop         
#>  Min.   :1   Min.   :  1   Length:721         Min.   :   0   Min.   :0.0000000  
#>  1st Qu.:2   1st Qu.: 26   Class :character   1st Qu.:1606   1st Qu.:0.0000000  
#>  Median :4   Median : 52   Mode  :character   Median :3461   Median :0.0000001  
#>  Mean   :4   Mean   : 52                      Mean   :3465   Mean   :0.0097087  
#>  3rd Qu.:6   3rd Qu.: 78                      3rd Qu.:5316   3rd Qu.:0.0003023  
#>  Max.   :7   Max.   :103                      Max.   :7100   Max.   :0.8501904
head(stablestage3(lathipmmean))
#>   matrix new_stage_id           orig_stage_id original_size      ss_prop
#> 1      1            1                      Sd       0.00000 0.8415146109
#> 2      1            2                     Sdl     100.00000 0.1317153400
#> 3      1            3                    Dorm       0.00000 0.0034542301
#> 4      1            4  sz35.85354 rp1 mt1 ob1      35.85354 0.0017273280
#> 5      1            5 sz107.20855 rp1 mt1 ob1     107.20855 0.0009965845
#> 6      1            6 sz178.56356 rp1 mt1 ob1     178.56356 0.0010726486
```

The stable stage distribution appears to be dominated by seeds, followed by seedlings, and adults get progressively less common as they get larger.

  <p style="text-indent: 20px">Finally, the reproductive values.</p>

```r
summary(repvalue3(lathipmmean))
#>      matrix   new_stage_id orig_stage_id      original_size   left_vector           rep_value      
#>  Min.   :1   Min.   :  1   Length:721         Min.   :   0   Min.   :-0.6046683   Min.   :      0  
#>  1st Qu.:2   1st Qu.: 26   Class :character   1st Qu.:1606   1st Qu.:-0.0000238   1st Qu.:      2  
#>  Median :4   Median : 52   Mode  :character   Median :3461   Median : 0.0000000   Median :     42  
#>  Mean   :4   Mean   : 52                      Mean   :3465   Mean   :-0.0041197   Mean   : 212950  
#>  3rd Qu.:6   3rd Qu.: 78                      3rd Qu.:5316   3rd Qu.: 0.0000015   3rd Qu.:  13061  
#>  Max.   :7   Max.   :103                      Max.   :7100   Max.   : 0.6046321   Max.   :6046683
head(repvalue3(lathipmmean))
#>   matrix new_stage_id           orig_stage_id original_size left_vector rep_value
#> 1      1            1                      Sd       0.00000       0e+00         0
#> 2      1            2                     Sdl     100.00000       0e+00         0
#> 3      1            3                    Dorm       0.00000      -2e-07         1
#> 4      1            4  sz35.85354 rp1 mt1 ob1      35.85354      -2e-07         1
#> 5      1            5 sz107.20855 rp1 mt1 ob1     107.20855      -2e-07         1
#> 6      1            6 sz178.56356 rp1 mt1 ob1     178.56356      -2e-07         1
```

A quick scan through these values shows that the highest reproductive values are for the largest adults, all the way at the bottom of the data frame output.

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

<div id="refs"></div>
