CASE STUDIES OF SWEDISH Lathyrus vernus POPULATION
ORGANISM AND POPULATION
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 (Ehrlén and Lehtilä 2002). Lathyrus vernus lacks organs for vegetative spread and individuals are well delimited (Ehrlén 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 (Ehrlén and Van Groenendael 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 (Ehrlén and Eriksson 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 (Ehrlén and Münzbergová 2009).
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 (Ehrlén 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)2 × height, and summed the volumes of all shoots. This measure is closely correlated with the dry mass of aboveground tissues (r2 = 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 mm3 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 (Ehrlén and Eriksson 1996).
ANALYSES WITH LATHYRUS DATA
We will analyze these data in three different ways to illustrate the power and flexiblity of package lefko3:
through the estimation of raw matrices, with the intention of producing matrices similar to those published in Ehrlén (2000);
through the estimation of function-derived matrices using a stage classification different from Ehrlén (2000), developed using the natural logarithm of the size measure used in that study; and
through the construction of an integral projection model.
Analysis 1: Raw matrix estimation
Here, we will attempt to estimate matrices similar to and based on the dataset used in Ehrlén (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.
Step 1. Data characterization and reorganization
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).
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.
First, we load the dataset, and look at its dimensions.
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.
Here, we create a stageframe named lathframe based on Ehrlén (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.
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.
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 DormantA 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.
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.
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 42It 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.
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 42The 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.
Voilá! Now we can move on to creating a few extra objects that will help us estimate full population projection matrices.
Step 2. Provide supplemental information for matrix estimation
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.
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).
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 0Next 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.
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 1And 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.
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
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.
Now we are read to create some matrices!
Step 3. Estimate matrices
Now let’s try to create the raw Lefkovitch matrices starting with the ahistorical case, based on those shown in Ehrlén (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.
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.
The output from this analysis includes is a lefkoMat object, which is a list object with the following elements:
A: a list of full population projection matrices, in order of population, patch, and year
U: a list of matrices showing only survival-transition elements, in the same order as A
F: a list of matrices showing only fecundity elements, in the same order as A
hstages: a data frame showing the order of paired stages (given if matrices are historical, otherwise NA)
ahstages: this is the stageframe used in analysis, with stages reordered and edited as they occur in the matrix
labels: a table showing the order of matrices, according to population, patch, and year
matrixqc: a short vector used in summary statements to describe the overall quality of each matrix
dataqc: a short vector used in summary statements to describe key sampling aspects of the dataset
We can understand lefkoMat objects through the summary summary function.
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.
#> NULLWe 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.
Now let’s look at historical matrices. Because of the size of historical matrices, we will only show the summary.
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.
#> NULLSince 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.
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.
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
Next we will create a mean ahistorical matrix.
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"
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.
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
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.
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.
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.
colSums(ehrlen2mean$U[[1]])
#> [1] 0.3990000 0.7254520 0.8608343 0.8894690 0.9231182 0.9596113 0.6666667All of the values are within the realm of possibility for probabilities, so this matrix appears to be fine.
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.
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 0Do 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.
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.
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
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!
Error-checking is more difficult with historical matrices because they are typically one or two orders of magnitude bigger than their ahistorical counterparts. Take a look at the column sums here to see the difficulty.
colSums(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.5000000While not too bad here, other historical matrices often have more than 100 columns, sometimes many, many more (some historical matrices used in Shefferson et al. (2014) had dimensions of over 2500 x 2500!). In these cases we can assess the distribution of survival estimates for historical stages, which is given as the set of column sums in the survival-transition matrix, as below.
summary(colSums(ehrlen3mean$U[[1]]))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.0000 0.4458 0.8743 0.6945 0.9516 1.0000As 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
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.
lambda3(ehrlen2)
#> pop patch year2 lambda
#> 1 NA NA 1988 0.8952585
#> 2 NA NA 1989 0.9069721
#> 3 NA NA 1990 0.9297366Here is the lambda associated with the mean matrix.
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.
And we will now look at the same numbers for the historical analyses. First the annual matrices.
Now the mean matrix.
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 and Roach 2010). One particularly common such trade-off is the cost of growth: an individual that grows a great deal in one time step due to great environmental conditions in that year might pay a large cost of survival, growth, or reproduction in the next if those environmental conditions deteriorate (Shefferson, Warren II, and Pulliam 2014). While we do not argue that the drop in lambda is due to this specific trade-off, 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.
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.
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
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.
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 7This 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.
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.01597293Notice that this stable stage distribution is different from the ahistorical case, suggesting a strong influence of individual history.
Let’s take a look at the reproductive values now, in similar order to the stable stage distribution case.
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.56101These 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.
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.0639748The 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.
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.5780928Notice that the reproductive values are quite different here, once again suggesting a strong importance to individual history.
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
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.
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 et al. 2013).
Here, we create a geometric mean matrix using the sparse option.
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
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.
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
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.
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.
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.0000000This 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).
Let’s compare lambda in each case. First, the geometric sparse case
Now the arithmetic case.
Finally, the true geometric case.
The true geometric mean has the lowest lambda, while the highest is the true arithmetic case.
Analysis 2: Function-derived matrix estimation
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.
Step 1. Data characterization and reorganization
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.
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 1248It 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.
It can also help to take a look at plots of these distributions. First we will look at the raw volume data.
plot(density(c(lathyrus$Volume88, lathyrus$Volume89, lathyrus$Volume90, lathyrus$Volume91),
na.rm = TRUE), main = "Volume")
plot of chunk Ch2an2st1ln383
Next we will look at the log volume data.
plot(density(c(lathyrus$lnVol88, lathyrus$lnVol89, lathyrus$lnVol90, lathyrus$lnVol91),
na.rm = TRUE), main = "Log volume")
plot of chunk Ch2an2st1ln389
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.
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.
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)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.
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
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.
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 36Something 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.
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] 2515In 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.
lathvertln$feca2 <- round(lathvertln$feca2)
lathvertln$feca1 <- round(lathvertln$feca1)
lathvertln$feca3 <- round(lathvertln$feca3)All set! Now to move on to supplemental descriptive information.
Step 2. Provide supplemental information for matrix estimation
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).
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.
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))Next, we move to a new step: modeling vital rates.
Step 3. Develop models of demographic parameters
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).
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.
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.
The output can be rather verbose, and so we have limited it with the 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 et al. 2015). It also shows warnings originating from the dredge function of package MuMIn (Bartoń 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.
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.
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 and Anderson 2002).
Let’s take a peek at the summary for the model selection output. This condenses the output to the necessities.
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.
#> NULLThe 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.
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
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.
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.
#> NULLThe 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.
Now we are ready to estimate matrices.
Step 4. Estimate matrices
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.
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.
#> NULLLet’s take a look at the first of these matrices.
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-01The 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.
Now we will estimate the historical matrices.
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.
#> NULLOnce again, the matrices are vastly bigger. Indeed, the ahistorical matrices here are 21 x 21, so contain 441 elements each. The historical matrices are 212 x 212 (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.
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.
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.
#> NULLThis exercise eliminated 21 rows and columns, so the matrices definitely got smaller.
Now the mean matrices. Note that these include 6 mean matrices for the patches, followed by a grand mean matrix.
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.
#> NULLAnd here we do some quality control. We will focus on the grand mean U matrix in each case.
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.9776Quick scans of these numbers show no real problems.
Step 5. Analyses
Finally, let’s estimate the deterministic population growth rate in each case. First the ahistorical case.
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.9250544Now the historical case.
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.9422742Note 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.
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.
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-04The 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.
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.
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):2058Now the first 6 entries.
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.00000000Here 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.
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-04These are quite different values, showing the impact of history once again.
Now let’s take a peek at the reproductive values. Let’s start with the ahistorical case.
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.35721As before, these are similar to those estimated by popbio’s reproductive.value() function. Let’s now see the historical case.
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.9866105Once 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)
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!
Step 1. Data characterization and reorganization
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 (2000), but modify the size bits to allow IPM construction and make all mature stages other than vegetative dormancy reproductive.
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.
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)We will also add some descriptive comments to this stageframe so that we know what each of these stages is.
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.
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).
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. :1Voila! Now we move on to create the extra bits of information needed for matrix estimation.
Step 2: Develop supplemental information for matrix estimation
Now we will provide some given transitions. This is the same overwrite object as used before in ahistorical matrix estimation.
lathover2 <- overwrite(stage3 = c("Sd", "Sdl"), stage2 = c("Sd", "Sd"), givenrate = c(0.345, 0.054))
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).
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.
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] 265Nearly 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
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.
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 ?isSingularNow we can take a look at the summary of these models.
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.
#> NULLLet’s move on now to the matrices themselves.
Step 4. Estimate matrices
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.
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.
#> NULLOne 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.
Let’s take a look at the top-left corner of the first matrix (the matrix is too huge to inspect in full here).
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.002791769This 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.
summary(colSums(lathmat2ipm$U[[1]]))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.01847 0.98280 0.99642 0.95570 0.99893 0.99968Let’s estimate the mean IPM matrices. This code will estimate mean matrices for each of the 6 patches, followed by a grand mean.
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.
#> NULLAs a final check, let’s take a look at the column sums of the grand mean survival-transition matrix.
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.99885All looks fine!
Step 5. Analyses
Now let’s estimate the deterministic population growth rate for the IPM analysis.
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.8974453These 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.
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.7658357The 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.
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.
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.0010726486The stable stage distribution appears to be dominated by seeds, followed by seedlings, and adults get progressively less common as they get larger.
Finally, the reproductive values.
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 1A 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.