CASE STUDIES OF AMERICAN Cypripedium candidum POPULATION
ORGANISM AND POPULATION
In this tutorial, we will focus on dataset cypdata, which is a demographic dataset for an American population of the white lady’s slipper, Cypripedium candidum. This population was monitored annually from 2004 to 2009. This species is of conservation concern, and the population is located within a state nature preserve located in northwestern Illinois, USA. More about this population and its characteristics is given in Shefferson et al. (2001) and Shefferson et al. (2017).
Population matrix projection modeling requires an appropriate life history model showing how stages are related. The figure below shows a very general life history model detailing how stages are related to one another in Cypripedium candidum. The first stage of life is a dormant seed stage, although an individual may germinate in the year following seed production. The first germinated stage is a protocorm, which is an underground, mycoheterotrophic stage unique to the families Orchidaceae and Pyrolaceae. There are three years of protocorm stages, followed by a seedling stage, and finally a set of stages that comprise the size-classified adult portion of life. The figure shows 49 such stages, each for a different number of stems and one of two reproductive statuses. These stages may be compressed for different circumstances (more on this later).
We can see a variety of transitions within this figure. The juvenile stages have fairly simple transitions. New recruits may enter the population directly from germination of a seed produced the previous year, in which case they start in the protocorm 1 stage, or they may begin as dormant seed. Dormant seed may remain dormant, die, or germinate into the protocorm 1 stage. Protocorms exist for up to 3 years, yielding the protocorm 1, 2, and 3 stages, without any possibility of staying within each of these stages for more than a single year. This leads to a seedling stage, in which the plant may persist for many years, before becoming mature. Here, maturity does not really refer to reproduction per se, but rather to a morphology indistinguishable from a reproductive plant except for the lack of a flower. The first mature stage is usually either vegetative dormancy (dorm), during which time the plant does not sprout, or a small, non-flowering adult (1V). Once in this portion of the life history, the plant may move among 49 mature stages, including vegetative dormancy, 1-24 shoots without flowers, or 1-24 shoots with at least one flower.
The dataset cypdata, and the equivalent dataset cypvert, which is simply structured differently, both include only data for the adult stages, and so later we will need to set juvenile transitions to constants.
ANALYSES WITH CYPRIPEDIUM DATA
We will analyze these data in two different ways to illustrate the power and flexibility of package lefko3:
through the estimation of raw matrices using a simplified life history; and
through the estimation of function-derived matrices using a count-based size metric and the general life history model shown above.
Analysis 1. Raw matrix estimation
In this portion of the tutorial, we will create raw matrices with these data. Here, we use the term ‘raw’ to refer to the fact that we will estimate matrix elements as exact proportions of individuals surviving and transitioning to different stages. This requires us to develop a life history model that is both biologically realistic, and numerically parsimonious. The former requirement means that stages need to be defined in biologically meaningful ways, while the latter requirement means that we need to design our life stages in such a way that most years include some individuals in each stage. We also need to consider the fact that very low numbers of stages result in biased matrix analyses, so we want to make sure that we have at least 7 stages in the final model (Salguero-Gómez and Plotkin 2010).
Step 1a. Data characterization and reorganization (horizontal file)
First let’s wipe the memory, load lefko3, and then load the data.
rm(list=ls(all=TRUE))
library(lefko3)
data(cypdata)
summary(cypdata)
#> plantid patch censor Inf2.04 Inf.04 Veg.04 Pod.04 Inf2.05 Inf.05
#> Min. : 164.0 A:23 Min. :1 Min. :0 Min. :0.0000 Min. : 0.000 Min. :0.0 Min. :0.00000 Min. : 0.000
#> 1st Qu.: 265.0 B:35 1st Qu.:1 1st Qu.:0 1st Qu.:0.0000 1st Qu.: 1.000 1st Qu.:0.0 1st Qu.:0.00000 1st Qu.: 0.000
#> Median : 455.0 C:19 Median :1 Median :0 Median :0.0000 Median : 2.000 Median :0.0 Median :0.00000 Median : 0.000
#> Mean : 669.1 Mean :1 Mean :0 Mean :0.6923 Mean : 2.923 Mean :0.2 Mean :0.04478 Mean : 1.537
#> 3rd Qu.: 829.0 3rd Qu.:1 3rd Qu.:0 3rd Qu.:1.0000 3rd Qu.: 4.000 3rd Qu.:0.0 3rd Qu.:0.00000 3rd Qu.: 2.000
#> Max. :1560.0 Max. :1 Max. :0 Max. :8.0000 Max. :12.000 Max. :3.0 Max. :1.00000 Max. :18.000
#> NA's :12 NA's :12 NA's :12 NA's :12 NA's :10 NA's :10
#> Veg.05 Pod.05 Inf2.06 Inf.06 Veg.06 Pod.06 Inf2.07 Inf.07
#> Min. :0.000 Min. :0.0000 Min. :0 Min. : 0.0000 Min. : 0.000 Min. :0.0000 Min. :0 Min. :0.0000
#> 1st Qu.:1.000 1st Qu.:0.0000 1st Qu.:0 1st Qu.: 0.0000 1st Qu.: 1.000 1st Qu.:0.0000 1st Qu.:0 1st Qu.:0.0000
#> Median :1.000 Median :0.0000 Median :0 Median : 0.0000 Median : 2.000 Median :0.0000 Median :0 Median :0.0000
#> Mean :2.134 Mean :0.6567 Mean :0 Mean : 0.9016 Mean : 2.213 Mean :0.3934 Mean :0 Mean :0.6271
#> 3rd Qu.:3.000 3rd Qu.:1.0000 3rd Qu.:0 3rd Qu.: 1.0000 3rd Qu.: 3.000 3rd Qu.:0.0000 3rd Qu.:0 3rd Qu.:1.0000
#> Max. :9.000 Max. :7.0000 Max. :0 Max. :18.0000 Max. :13.000 Max. :4.0000 Max. :0 Max. :7.0000
#> NA's :10 NA's :10 NA's :16 NA's :16 NA's :16 NA's :16 NA's :18
#> Veg.07 Pod.07 Inf2.08 Inf.08 Veg.08 Pod.08 Inf2.09 Inf.09
#> Min. : 0.000 Min. :0.0000 Min. :0 Min. : 0.0000 Min. : 0.00 Min. :0.0000 Min. :0 Min. : 0.000
#> 1st Qu.: 1.000 1st Qu.:0.0000 1st Qu.:0 1st Qu.: 0.0000 1st Qu.: 1.00 1st Qu.:0.0000 1st Qu.:0 1st Qu.: 0.000
#> Median : 2.000 Median :0.0000 Median :0 Median : 0.0000 Median : 2.00 Median :0.0000 Median :0 Median : 1.000
#> Mean : 2.627 Mean :0.0678 Mean :0 Mean : 0.8868 Mean : 2.83 Mean :0.1509 Mean :0 Mean : 1.833
#> 3rd Qu.: 4.000 3rd Qu.:0.0000 3rd Qu.:0 3rd Qu.: 1.0000 3rd Qu.: 4.00 3rd Qu.:0.0000 3rd Qu.:0 3rd Qu.: 2.000
#> Max. :13.000 Max. :1.0000 Max. :0 Max. :11.0000 Max. :13.00 Max. :2.0000 Max. :0 Max. :11.000
#> NA's :18 NA's :18 NA's :24 NA's :24 NA's :24 NA's :24 NA's :17 NA's :17
#> Veg.09 Pod.09
#> Min. : 0.000 Min. :0.000
#> 1st Qu.: 1.000 1st Qu.:0.000
#> Median : 1.000 Median :1.000
#> Mean : 2.233 Mean :1.133
#> 3rd Qu.: 3.000 3rd Qu.:1.000
#> Max. :10.000 Max. :8.000
#> NA's :17 NA's :17
The dataset that we have provided is organized in horizontal format, meaning that rows correspond to unique individuals and columns correspond to state in particular years. Looking at the original Excel spreadsheet (below), you will note a repeating pattern in the names of the columns. Package lefko3 includes functions to handle data in horizontal format, as well as functions to handle vertically formatted data (i.e. data for individuals is broken up across rows, where each row is a unique combination of individual and year in time t).
In this dataset, there are 77 individuals, so there are 77 rows with data (not counting the header). There are 27 columns. Note the the first 3 columns are variables giving specific information about each individual, with each individual’s data entirely restricted to one row. This is followed by a number of sets of 4 columns, each named Inf2.XX, Inf.XX, Veg.XX, and Pod.XX. The XX in each case corresponds to a specific year, which are organized consecutively. Thus, columns 4-7 refer to year 04 (short for 2004), columns 8-11 refer to year 05, columns 12-15 refer to year 06, etc. To properly conduct this exercise, we need to know the exact number of years used, which is six years here (includes all years from 2004 to 2009). Note that each year MUST utilize exactly the same pattern of columns. Also note that this package also includes dataset cypvert, which is the same dataset but set in ahistorical vertical format.
Now on to the assessment of size. The full sizes are actually the sums of columns (representing sprouts) governing different years. We will take these sums, and then assess the distribution of sizes across years. We will look at all years and look for general patterns and abnormalities.
size.04 <- cypdata$Inf2.04 + cypdata$Inf.04 + cypdata$Veg.04
size.05 <- cypdata$Inf2.05 + cypdata$Inf.05 + cypdata$Veg.05
size.06 <- cypdata$Inf2.06 + cypdata$Inf.06 + cypdata$Veg.06
size.07 <- cypdata$Inf2.07 + cypdata$Inf.07 + cypdata$Veg.07
size.08 <- cypdata$Inf2.08 + cypdata$Inf.08 + cypdata$Veg.08
size.09 <- cypdata$Inf2.09 + cypdata$Inf.09 + cypdata$Veg.09
summary(c(size.04, size.05, size.06, size.07, size.08, size.09))
#> Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
#> 1.000 1.000 2.000 3.581 5.000 24.000 97Now let’s quickly plot this.
plot of chunk Ch3an1st1aln79pl
This exercise gives us a reasonable idea of size classes to use for adults. We will have a dormant class (size = 0 shoots), extra small class (1 shoot), small class (2-3 shoots), medium class (4-5 shoots), large class (6-10 shoots), and extra large class (>10 shoots). Let’s define a stageframe that shows this.
sizevector <- c(0, 0, 0, 0, 0, 0, 1, 2.5, 4.5, 8, 17.5)
stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "XSm", "Sm", "Md", "Lg", "XLg")
repvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
obsvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
matvector <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)
immvector <- c(0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0)
propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
indataset <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)
binvec <- c(0, 0, 0, 0, 0, 0.5, 0.5, 1, 1, 2.5, 7)
cypframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector,
repstatus = repvector, obsstatus = obsvector,
matstatus = matvector, propstatus = propvector,
immstatus = immvector, indataset = indataset,
binhalfwidth = binvec)Now we will add some comments to the stageframe for our later use in interpretation.
cypframe_raw$comments[(cypframe_raw$stagenames == "SD")] <- "Dormant seed"
cypframe_raw$comments[(cypframe_raw$stagenames == "P1")] <- "1st yr protocorm"
cypframe_raw$comments[(cypframe_raw$stagenames == "P2")] <- "2nd yr protocorm"
cypframe_raw$comments[(cypframe_raw$stagenames == "P3")] <- "3rd yr protocorm"
cypframe_raw$comments[(cypframe_raw$stagenames == "SL")] <- "Seedling"
cypframe_raw$comments[(cypframe_raw$stagenames == "D")] <- "Dormant adult"
cypframe_raw$comments[(cypframe_raw$stagenames == "XSm")] <- "Extra small adult (1 shoot)"
cypframe_raw$comments[(cypframe_raw$stagenames == "Sm")] <- "Small adult (2-3 shoots)"
cypframe_raw$comments[(cypframe_raw$stagenames == "Md")] <- "Medium adult (4-5 shoots)"
cypframe_raw$comments[(cypframe_raw$stagenames == "Lg")] <- "Large adult (6-10 shoots)"
cypframe_raw$comments[(cypframe_raw$stagenames == "XLg")] <- "Extra large adult (>10 shoots)"
cypframe_raw
#> stagenames size repstatus obsstatus propstatus immstatus matstatus indataset binhalfwidth_raw sizebin_min sizebin_max sizebin_center
#> 1 SD 0.0 0 0 1 0 0 0 0.0 0.0 0.0 0.0
#> 2 P1 0.0 0 0 0 1 0 0 0.0 0.0 0.0 0.0
#> 3 P2 0.0 0 0 0 1 0 0 0.0 0.0 0.0 0.0
#> 4 P3 0.0 0 0 0 1 0 0 0.0 0.0 0.0 0.0
#> 5 SL 0.0 0 0 0 1 0 0 0.0 0.0 0.0 0.0
#> 6 D 0.0 0 0 0 0 1 1 0.5 -0.5 0.5 0.0
#> 7 XSm 1.0 1 1 0 0 1 1 0.5 0.5 1.5 1.0
#> 8 Sm 2.5 1 1 0 0 1 1 1.0 1.5 3.5 2.5
#> 9 Md 4.5 1 1 0 0 1 1 1.0 3.5 5.5 4.5
#> 10 Lg 8.0 1 1 0 0 1 1 2.5 5.5 10.5 8.0
#> 11 XLg 17.5 1 1 0 0 1 1 7.0 10.5 24.5 17.5
#> sizebin_width comments
#> 1 0 Dormant seed
#> 2 0 1st yr protocorm
#> 3 0 2nd yr protocorm
#> 4 0 3rd yr protocorm
#> 5 0 Seedling
#> 6 1 Dormant adult
#> 7 1 Extra small adult (1 shoot)
#> 8 2 Small adult (2-3 shoots)
#> 9 2 Medium adult (4-5 shoots)
#> 10 5 Large adult (6-10 shoots)
#> 11 14 Extra large adult (>10 shoots)
Next we will create the vertical dataset. Note that because we are lumping reproductive and non-reproductive individuals into the non-dormant adult classes, we need to set NRasRP = TRUE. Otherwise, verticalize3() will attempt to use the reproductive status of individuals in classification, and will fail due to the presence of non-reproductive adults.
cypraw_v1 <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004,
patchidcol = "patch", individcol = "plantid",
blocksize = 4, size1col = "Inf2.04", size2col = "Inf.04",
size3col = "Veg.04", repstr1col = "Inf.04",
repstr2col = "Inf2.04", fec1col = "Pod.04",
stageassign = cypframe_raw, stagesize = "sizeadded",
NAas0 = TRUE, NRasRep = TRUE)
summary(cypraw_v1)
#> rowid popid patchid individ year2 sizea1 sizea2 sizea3
#> Min. : 1.00 Mode:logical A: 93 Min. : 164.0 Min. :2004 Min. :0.000000 Min. :0.000000 Min. :0.000000
#> 1st Qu.:21.00 NA's:320 B:154 1st Qu.: 391.0 1st Qu.:2005 1st Qu.:0.000000 1st Qu.:0.000000 1st Qu.:0.000000
#> Median :37.50 C: 73 Median : 453.0 Median :2006 Median :0.000000 Median :0.000000 Median :0.000000
#> Mean :38.45 Mean : 651.5 Mean :2006 Mean :0.009375 Mean :0.009375 Mean :0.009375
#> 3rd Qu.:56.00 3rd Qu.: 476.0 3rd Qu.:2007 3rd Qu.:0.000000 3rd Qu.:0.000000 3rd Qu.:0.000000
#> Max. :77.00 Max. :1560.0 Max. :2008 Max. :1.000000 Max. :1.000000 Max. :1.000000
#> repstra1 repstra2 repstra3 feca1 feca2 feca3 sizeb1
#> Min. : 0.0000 Min. : 0.0000 Min. : 0.000 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. : 0.0000
#> 1st Qu.: 0.0000 1st Qu.: 0.0000 1st Qu.: 0.000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.: 0.0000
#> Median : 0.0000 Median : 0.0000 Median : 0.000 Median :0.0000 Median :0.0000 Median :0.0000 Median : 0.0000
#> Mean : 0.7469 Mean : 0.8969 Mean : 1.069 Mean :0.2656 Mean :0.2906 Mean :0.4562 Mean : 0.7469
#> 3rd Qu.: 1.0000 3rd Qu.: 1.0000 3rd Qu.: 1.000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.: 1.0000
#> Max. :18.0000 Max. :18.0000 Max. :18.000 Max. :7.0000 Max. :7.0000 Max. :8.0000 Max. :18.0000
#> sizeb2 sizeb3 repstrb1 repstrb2 repstrb3 sizec1 sizec2
#> Min. : 0.0000 Min. : 0.000 Min. :0.000000 Min. :0.000000 Min. :0.000000 Min. : 0.0 Min. : 0.000
#> 1st Qu.: 0.0000 1st Qu.: 0.000 1st Qu.:0.000000 1st Qu.:0.000000 1st Qu.:0.000000 1st Qu.: 0.0 1st Qu.: 1.000
#> Median : 0.0000 Median : 0.000 Median :0.000000 Median :0.000000 Median :0.000000 Median : 1.0 Median : 2.000
#> Mean : 0.8969 Mean : 1.069 Mean :0.009375 Mean :0.009375 Mean :0.009375 Mean : 1.9 Mean : 2.416
#> 3rd Qu.: 1.0000 3rd Qu.: 1.000 3rd Qu.:0.000000 3rd Qu.:0.000000 3rd Qu.:0.000000 3rd Qu.: 3.0 3rd Qu.: 3.000
#> Max. :18.0000 Max. :18.000 Max. :1.000000 Max. :1.000000 Max. :1.000000 Max. :13.0 Max. :13.000
#> sizec3 size1added size2added size3added obsstatus1 obsstatus2 obsstatus3 repstatus1
#> Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. :0.0000 Min. :0.0000 Min. :0.0 Min. :0.0000
#> 1st Qu.: 1.000 1st Qu.: 0.000 1st Qu.: 1.000 1st Qu.: 1.000 1st Qu.:0.0000 1st Qu.:1.0000 1st Qu.:1.0 1st Qu.:0.0000
#> Median : 1.000 Median : 2.000 Median : 2.000 Median : 2.000 Median :1.0000 Median :1.0000 Median :1.0 Median :0.0000
#> Mean : 2.209 Mean : 2.656 Mean : 3.322 Mean : 3.288 Mean :0.7469 Mean :0.9531 Mean :0.9 Mean :0.2875
#> 3rd Qu.: 3.000 3rd Qu.: 4.000 3rd Qu.: 4.000 3rd Qu.: 4.000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0 3rd Qu.:1.0000
#> Max. :13.000 Max. :21.000 Max. :24.000 Max. :24.000 Max. :1.0000 Max. :1.0000 Max. :1.0 Max. :1.0000
#> repstatus2 repstatus3 fecstatus1 fecstatus2 fecstatus3 firstseen lastseen alive1
#> Min. :0.0000 Min. :0.0 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :2004 Min. :2004 Min. :0.0000
#> 1st Qu.:0.0000 1st Qu.:0.0 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:2004 1st Qu.:2009 1st Qu.:1.0000
#> Median :0.0000 Median :0.0 Median :0.0000 Median :0.0000 Median :0.0000 Median :2004 Median :2009 Median :1.0000
#> Mean :0.3688 Mean :0.4 Mean :0.1344 Mean :0.1562 Mean :0.2219 Mean :2004 Mean :2009 Mean :0.7688
#> 3rd Qu.:1.0000 3rd Qu.:1.0 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:2004 3rd Qu.:2009 3rd Qu.:1.0000
#> Max. :1.0000 Max. :1.0 Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :2008 Max. :2009 Max. :1.0000
#> alive2 alive3 obsage obslifespan stage1 stage2 stage3 matstatus1
#> Min. :1 Min. :0.0000 Min. :0.000 Min. :0.000 Length:320 Length:320 Length:320 Min. :1
#> 1st Qu.:1 1st Qu.:1.0000 1st Qu.:1.000 1st Qu.:5.000 Class :character Class :character Class :character 1st Qu.:1
#> Median :1 Median :1.0000 Median :2.000 Median :5.000 Mode :character Mode :character Mode :character Median :1
#> Mean :1 Mean :0.9469 Mean :1.853 Mean :4.556 Mean :1
#> 3rd Qu.:1 3rd Qu.:1.0000 3rd Qu.:3.000 3rd Qu.:5.000 3rd Qu.:1
#> Max. :1 Max. :1.0000 Max. :4.000 Max. :5.000 Max. :1
#> matstatus2 matstatus3
#> Min. :1 Min. :1
#> 1st Qu.:1 1st Qu.:1
#> Median :1 Median :1
#> Mean :1 Mean :1
#> 3rd Qu.:1 3rd Qu.:1
#> Max. :1 Max. :1
dim(cypraw_v1)
#> [1] 320 48Step 1b. Data characterization and reorganization (vertical file)
We may also wish to see how to proceed if our original dataset is already in vertical, but ahistorical, format. Here, we see such a situation and use the historicalize3() function. First, let’s load the ahistorical vertical raw data file.
data(cypvert)
summary(cypvert)
#> plantid patch censor year2 Inf2.2 Inf.2 Veg.2 Pod.2
#> Min. : 164.0 Length:331 Min. :1 Min. :2004 Min. :0.000000 Min. : 0.000 Min. : 0.000 Min. :0.0000
#> 1st Qu.: 391.0 Class :character 1st Qu.:1 1st Qu.:2005 1st Qu.:0.000000 1st Qu.: 0.000 1st Qu.: 1.000 1st Qu.:0.0000
#> Median : 454.0 Mode :character Median :1 Median :2006 Median :0.000000 Median : 0.000 Median : 2.000 Median :0.0000
#> Mean : 662.8 Mean :1 Mean :2006 Mean :0.009646 Mean : 0.941 Mean : 2.534 Mean :0.3049
#> 3rd Qu.: 664.0 3rd Qu.:1 3rd Qu.:2007 3rd Qu.:0.000000 3rd Qu.: 1.000 3rd Qu.: 3.000 3rd Qu.:0.0000
#> Max. :1560.0 Max. :1 Max. :2008 Max. :1.000000 Max. :18.000 Max. :13.000 Max. :7.0000
#> NA's :20 NA's :26 NA's :26 NA's :26
#> Inf2.3 Inf.3 Veg.3 Pod.3
#> Min. :0.000000 Min. : 0.000 Min. : 0.00 Min. :0.0000
#> 1st Qu.:0.000000 1st Qu.: 0.000 1st Qu.: 1.00 1st Qu.:0.0000
#> Median :0.000000 Median : 0.000 Median : 1.00 Median :0.0000
#> Mean :0.009804 Mean : 1.173 Mean : 2.39 Mean :0.4933
#> 3rd Qu.:0.000000 3rd Qu.: 1.000 3rd Qu.: 3.00 3rd Qu.:0.0000
#> Max. :1.000000 Max. :18.000 Max. :13.00 Max. :8.0000
#> NA's :25 NA's :31 NA's :31 NA's :31And let’s also look at its dimensions.
This dataset appears to be longer and narrower, with more rows and fewer columns. This is because we now split data for each individual across multiple columns. A single column designates time in year t, given as year2. This dataset then includes columns for only pairs of consecutive years corresponding to times t and t+1. State in time t-1 is not presented because this is an ahistorical dataset. Fortunately, this dataset includes the plantid variable, which is an individual identity term and must be supplied for conversion. The historicalize3() function can use individual identity variables to reorganize datasets into historical vertical format.
cypraw_v2 <- historicalize3(data = cypvert, patchidcol = "patch", individcol = "plantid",
year2col = "year2", sizea2col = "Inf2.2", sizea3col = "Inf2.3",
sizeb2col = "Inf.2", sizeb3col = "Inf.3", sizec2col = "Veg.2",
sizec3col = "Veg.3", repstra2col = "Inf2.2", repstra3col = "Inf2.3",
repstrb2col = "Inf.2", repstrb3col = "Inf.3", feca2col = "Pod.2",
feca3col = "Pod.3", repstrrel = 2, stageassign = cypframe_raw,
stagesize = "sizeadded", censorcol = "censor", censor = FALSE,
NAas0 = TRUE, NRasRep = TRUE, reduce = TRUE)
summary(cypraw_v2)
#> popid patchid rowid sizea1 sizeb1 sizec1 repstra1
#> Mode:logical Length:320 Min. : 1.0 Min. :0.000000 Min. : 0.0000 Min. : 0.000 Min. :0.000000
#> NA's:320 Class :character 1st Qu.: 41.0 1st Qu.:0.000000 1st Qu.: 0.0000 1st Qu.: 0.000 1st Qu.:0.000000
#> Mode :character Median : 95.5 Median :0.000000 Median : 0.0000 Median : 1.000 Median :0.000000
#> Mean :111.3 Mean :0.009375 Mean : 0.7375 Mean : 1.881 Mean :0.009375
#> 3rd Qu.:178.2 3rd Qu.:0.000000 3rd Qu.: 1.0000 3rd Qu.: 3.000 3rd Qu.:0.000000
#> Max. :311.0 Max. :1.000000 Max. :18.0000 Max. :13.000 Max. :1.000000
#> repstrb1 feca1 stage1 year2 sizea2 sizeb2 sizec2
#> Min. : 0.0000 Min. :0.0000 Length:320 Min. :2004 Min. :0.000000 Min. : 0.0000 Min. : 0.000
#> 1st Qu.: 0.0000 1st Qu.:0.0000 Class :character 1st Qu.:2005 1st Qu.:0.000000 1st Qu.: 0.0000 1st Qu.: 1.000
#> Median : 0.0000 Median :0.0000 Mode :character Median :2006 Median :0.000000 Median : 0.0000 Median : 1.500
#> Mean : 0.7375 Mean :0.2625 Mean :2006 Mean :0.009375 Mean : 0.8812 Mean : 2.378
#> 3rd Qu.: 1.0000 3rd Qu.:0.0000 3rd Qu.:2007 3rd Qu.:0.000000 3rd Qu.: 1.0000 3rd Qu.: 3.000
#> Max. :18.0000 Max. :7.0000 Max. :2008 Max. :1.000000 Max. :18.0000 Max. :13.000
#> repstra2 repstrb2 feca2 stage2 sizea3 sizeb3 sizec3
#> Min. :0.000000 Min. : 0.0000 Min. :0.0000 Length:320 Min. :0.000000 Min. : 0.000 Min. : 0.000
#> 1st Qu.:0.000000 1st Qu.: 0.0000 1st Qu.:0.0000 Class :character 1st Qu.:0.000000 1st Qu.: 0.000 1st Qu.: 1.000
#> Median :0.000000 Median : 0.0000 Median :0.0000 Mode :character Median :0.000000 Median : 0.000 Median : 1.000
#> Mean :0.009375 Mean : 0.8812 Mean :0.2875 Mean :0.009375 Mean : 1.059 Mean : 2.209
#> 3rd Qu.:0.000000 3rd Qu.: 1.0000 3rd Qu.:0.0000 3rd Qu.:0.000000 3rd Qu.: 1.000 3rd Qu.: 3.000
#> Max. :1.000000 Max. :18.0000 Max. :7.0000 Max. :1.000000 Max. :18.000 Max. :13.000
#> repstra3 repstrb3 feca3 stage3 firstseen lastseen alive1 alive2
#> Min. :0.000000 Min. : 0.000 Min. :0.0000 Length:320 Min. :2004 Min. :2005 Min. :0.0000 Min. :1
#> 1st Qu.:0.000000 1st Qu.: 0.000 1st Qu.:0.0000 Class :character 1st Qu.:2004 1st Qu.:2009 1st Qu.:1.0000 1st Qu.:1
#> Median :0.000000 Median : 0.000 Median :0.0000 Mode :character Median :2004 Median :2009 Median :1.0000 Median :1
#> Mean :0.009375 Mean : 1.059 Mean :0.4437 Mean :2004 Mean :2009 Mean :0.7688 Mean :1
#> 3rd Qu.:0.000000 3rd Qu.: 1.000 3rd Qu.:0.0000 3rd Qu.:2004 3rd Qu.:2009 3rd Qu.:1.0000 3rd Qu.:1
#> Max. :1.000000 Max. :18.000 Max. :8.0000 Max. :2008 Max. :2009 Max. :1.0000 Max. :1
#> alive3 obsstatus1 obsstatus2 obsstatus3 matstatus1 matstatus2 matstatus3 repstatus1 repstatus2
#> Min. :1 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :1 Min. :1 Min. :1 Min. :0.0000 Min. :0.0000
#> 1st Qu.:1 1st Qu.:0.0000 1st Qu.:1.0000 1st Qu.:1.0000 1st Qu.:1 1st Qu.:1 1st Qu.:1 1st Qu.:0.0000 1st Qu.:0.0000
#> Median :1 Median :1.0000 Median :1.0000 Median :1.0000 Median :1 Median :1 Median :1 Median :0.0000 Median :0.0000
#> Mean :1 Mean :0.7344 Mean :0.9281 Mean :0.9125 Mean :1 Mean :1 Mean :1 Mean :0.2812 Mean :0.3563
#> 3rd Qu.:1 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1 3rd Qu.:1 3rd Qu.:1 3rd Qu.:1.0000 3rd Qu.:1.0000
#> Max. :1 Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1 Max. :1 Max. :1 Max. :1.0000 Max. :1.0000
#> repstatus3 fecstatus1 fecstatus2 fecstatus3 individ obsage obslifespan size1added
#> Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. : 164.0 Min. :0.000 Min. :1.000 Min. : 0.000
#> 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.: 391.0 1st Qu.:1.000 1st Qu.:5.000 1st Qu.: 0.000
#> Median :0.0000 Median :0.0000 Median :0.0000 Median :0.0000 Median : 453.0 Median :2.000 Median :5.000 Median : 1.000
#> Mean :0.4031 Mean :0.1313 Mean :0.1531 Mean :0.2125 Mean : 651.5 Mean :1.853 Mean :4.706 Mean : 2.628
#> 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.: 476.0 3rd Qu.:3.000 3rd Qu.:5.000 3rd Qu.: 4.000
#> Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1560.0 Max. :4.000 Max. :5.000 Max. :21.000
#> size2added size3added
#> Min. : 0.000 Min. : 0.000
#> 1st Qu.: 1.000 1st Qu.: 1.000
#> Median : 2.000 Median : 2.000
#> Mean : 3.269 Mean : 3.278
#> 3rd Qu.: 4.000 3rd Qu.: 4.000
#> Max. :24.000 Max. :24.000We can compare the dimensions of these datasets.
The lengths of the datasets are the same in terms of rows and columns, and the variables and data are the same although the order of the rows might not exactly match.
Step 2. Provide supplemental information for matrix estimation
For our next step, we need to create a reproductive matrix, which essentially tells R not only which stages are reproductive, but which stages they lead to the reproduction of, and at what level. This matrix is mostly composed of 0s, but fecundity is noted as non-zero entries equal to a modifier to the full fecundity estimated by R. This matrix has as many rows and columns as the number of stages described in the stageframe for this dataset, and the rows and columns refer to these stages in the same order as in the stageframe. In many ways, it looks like a nearly empty population matrix, but notes the per-individual mean modifiers on fecundity for each stage that actually reproduces. Here, we first create a 0 matrix with dimensionality equal to the number of rows in cypframe_raw. Then we modify elements corresponding to fecundity by dividing fecundity evenly between dormant seeds (row 1) and germinating seeds (row 2).
rep_cyp_raw <- matrix(0, 11, 11)
rep_cyp_raw[1:2,7:11] <- 0.5
rep_cyp_raw
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
#> [1,] 0 0 0 0 0 0 0.5 0.5 0.5 0.5 0.5
#> [2,] 0 0 0 0 0 0 0.5 0.5 0.5 0.5 0.5
#> [3,] 0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0
#> [4,] 0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0
#> [5,] 0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0
#> [6,] 0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0
#> [7,] 0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0
#> [8,] 0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0
#> [9,] 0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0
#> [10,] 0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0
#> [11,] 0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0Next we will provide some given transitions. In this case, we are providing the seed dormancy probability and germination rate, which in this case are provided as transitions from the dormant seed stage to another year of seed dormancy or to a first-year protocorm, respectively. We are also providing the survival-transition probabilities between different year protocorm stages (P1, P2, and P3), to the seedling stage (SL), and from the seedling stage to some of the adult stages (XSm, Sm, and D). Let’s start with the ahistorical case.
cypover2r <- overwrite(stage3 = c("SD", "P1", "P2", "P3", "D", "XSm", "Sm", "SL", "SL"),
stage2 = c("SD", "SD", "P1", "P2", "P3", "P3", "P3", "P3", "SL"),
eststage3 = c(NA, NA, NA, NA, "D", "XSm", "Sm", NA, NA),
eststage2 = c(NA, NA, NA, NA, "D", "D", "D", NA, NA),
givenrate = c(0.1, 0.2, 0.2, 0.2, NA, NA, NA, 0.25, 0.4),
type = c("S", "S", "S", "S", "S", "S", "S", "S", "S"))
cypover2r
#> stage3 stage2 stage1 eststage3 eststage2 eststage1 givenrate convtype
#> 1 SD SD NA <NA> <NA> NA 0.10 1
#> 2 P1 SD NA <NA> <NA> NA 0.20 1
#> 3 P2 P1 NA <NA> <NA> NA 0.20 1
#> 4 P3 P2 NA <NA> <NA> NA 0.20 1
#> 5 D P3 NA D D NA NA 1
#> 6 XSm P3 NA XSm D NA NA 1
#> 7 Sm P3 NA Sm D NA NA 1
#> 8 SL P3 NA <NA> <NA> NA 0.25 1
#> 9 SL SL NA <NA> <NA> NA 0.40 1This overwrite table shows us that we have survival-transition probabilities (convtype = 1), that the given transitions are ahistorical, and particularly outlines probabilities for transitions that we cannot estimate with our dataset, which in this case refers to the immature stages of life. While 6 of these survival-transitions are given in the givenrate column, we also mark 3 of them as survival-transitions that we wish to use other estimates as proxies for. This is indicated via the eststageX columns, which have entries corresponding the stages to use as proxies (note that the givenrate entries are blank for these cases).
And now the historical case. Here we need to show the states in time step t-1 for this to work properly. Because of the extra time step, we will use the short-hand term "rep" to code for reproductive stages leading to the seeds and first-year protocorms that must survive to the next year.< /p>
cypover3r <- overwrite(stage3 = c("SD", "SD", "P1", "P1", "P2", "P3", "D", "XSm", "Sm",
"SL", "SL", "SL"), stage2 = c("SD", "SD", "SD", "SD", "P1", "P2",
"P3", "P3", "P3", "P3", "SL", "SL"), stage1 = c("SD", "rep", "SD",
"rep", "SD", "P1", "P2", "P2", "P2", "P2", "P3", "SL"),
eststage3 = c(NA, NA, NA, NA, NA, NA, "D", "XSm", "Sm", NA, NA, NA),
eststage2 = c(NA, NA, NA, NA, NA, NA, "D", "D", "D", NA, NA, NA),
eststage1 = c(NA, NA, NA, NA, NA, NA, "D", "D", "D", NA, NA, NA),
givenrate = c(0.1, 0.1, 0.2, 0.2, 0.2, 0.2, NA, NA, NA, 0.25, 0.4, 0.4),
type = c("S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S", "S"))
cypover3r
#> stage3 stage2 stage1 eststage3 eststage2 eststage1 givenrate convtype
#> 1 SD SD SD <NA> <NA> <NA> 0.10 1
#> 2 SD SD rep <NA> <NA> <NA> 0.10 1
#> 3 P1 SD SD <NA> <NA> <NA> 0.20 1
#> 4 P1 SD rep <NA> <NA> <NA> 0.20 1
#> 5 P2 P1 SD <NA> <NA> <NA> 0.20 1
#> 6 P3 P2 P1 <NA> <NA> <NA> 0.20 1
#> 7 D P3 P2 D D D NA 1
#> 8 XSm P3 P2 XSm D D NA 1
#> 9 Sm P3 P2 Sm D D NA 1
#> 10 SL P3 P2 <NA> <NA> <NA> 0.25 1
#> 11 SL SL P3 <NA> <NA> <NA> 0.40 1
#> 12 SL SL SL <NA> <NA> <NA> 0.40 1Now we are read to create some matrices!
Step 3. Estimate matrices
We will begin with the creation of a set of ahistorical matrices for the Cypripedium candidum dataset. This should not take terribly long. The rlefko2 function was created to deal with the construction of standard ahistorical matrices using raw data. Matrices may strongly differ, particularly if the demographic dataset is somewhat sparse.
cypmatrix2r <- rlefko2(data = cypraw_v1, stageframe = cypframe_raw, year = "all",
patch = "all", stages = c("stage3", "stage2", "stage1"),
size = c("size3added", "size2added"),
repmatrix = rep_cyp_raw, overwrite = cypover2r,
yearcol = "year2", patchcol = "patchid",
indivcol = "individ")The input for the rlefko2() function includes year = "all", but can be set to focus on any set of years included within the data. Package lefko3 actually includes a great deal of flexibility here, and can estimate many matrices covering all of the populations, patches, and years occurring in a specific dataset. The function-based matrix approach in the next section will showcase some more of this flexibility.
The output from this analysis includes is a lefkoMat object, which is a list object with the following elements:
A: a list of full population projection matrices, in order of population, patch, and year
U: a list of matrices showing only survival-transition elements, in the same order as A
F: a list of matrices showing only fecundity elements, in the same order as A
hstages: a data frame showing the order of paired stages (given if matrices are historical, otherwise NA)
ahstages: this is the stageframe used in analysis, with stages reordered and edited as they occur in the matrix
labels: a table showing the order of matrices, according to population, patch, and year
matrixqc: a short vector used in summary statements to describe the overall quality of each matrix
dataqc: a short vector used in summary statements to describe key sampling aspects of the dataset
Objects of class lefkoMat have their own summary statements, and we can use this to understand more about this particular object.
summary(cypmatrix2r)
#>
#> This lefkoMat object contains 15 matrices.
#>
#> Each matrix is a square matrix with 11 rows and columns, and a total of 121 elements.
#> A total of 258 survival-transitions were estimated, with 17.2 per matrix.
#> A total of 66 fecundity transitions were estimated, with 4.4 per matrix.
#>
#> The dataset contains a total of 74 unique individuals and 320 unique transitions.
#> NULLWe start off learning that 15 matrices were estimated, and we learn the dimensionality of those matrices. Of note here is the output telling us how many elements were actually estimated, both overall and per matrix, and the number of individuals and transitions the matrices are based on. Matrices are often overparameterized in population ecology, meaning that the number of elements estimated is quite high given the length of the data. It is typical for population ecologists to consider the total number of transitions in a dataset as a measure of the statistical power of a matrix, but the number of individuals used is just as important because each transition that an individual experiences is dependent on the other transitions that it also experiences. Indeed, this is the fundamental point that led to the development of historical matrices and of this package - the assumption that the status of an individual in the next time step is dependent only on its current state is too simplistic and leads to pseudoreplication. So, this output can be very helpful to understand the degree to which estimated matrices might be overparameterized.
Let’s look at a sample matrix here.
cypmatrix2r$A[[1]]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
#> [1,] 0.1 0.0 0.0 0.00 0.0 0 0.0000000 0.1 0.5 0.0 0
#> [2,] 0.2 0.0 0.0 0.00 0.0 0 0.0000000 0.1 0.5 0.0 0
#> [3,] 0.0 0.2 0.0 0.00 0.0 0 0.0000000 0.0 0.0 0.0 0
#> [4,] 0.0 0.0 0.2 0.00 0.0 0 0.0000000 0.0 0.0 0.0 0
#> [5,] 0.0 0.0 0.0 0.25 0.4 0 0.0000000 0.0 0.0 0.0 0
#> [6,] 0.0 0.0 0.0 0.00 0.0 0 0.0000000 0.0 0.0 0.0 0
#> [7,] 0.0 0.0 0.0 0.00 0.0 0 0.6363636 0.2 0.0 0.0 0
#> [8,] 0.0 0.0 0.0 0.00 0.0 0 0.2727273 0.6 0.5 0.0 0
#> [9,] 0.0 0.0 0.0 0.00 0.0 0 0.0000000 0.0 0.5 0.5 0
#> [10,] 0.0 0.0 0.0 0.00 0.0 0 0.0000000 0.2 0.0 0.5 0
#> [11,] 0.0 0.0 0.0 0.00 0.0 0 0.0000000 0.0 0.0 0.0 0The reader will note that although this is an ahistorical matrix, it is predominantly composed of 0 elements. This is a result of the sparseness of the data, and will likely lead to different elements shifting between 0 and positive elements across time.
Now we will create some historical matrices. Historical matrix construction parses the data much more finely across many more stages than ahistorical matrix construction, so historical matrices are even more likely to differ strongly across time, particularly as the number of individuals in a dataset decreases. Let’s see what these matrices look like here.
cypmatrix3r <- rlefko3(data = cypraw_v1, stageframe = cypframe_raw, year = "all",
patch = "all", stages = c("stage3", "stage2", "stage1"),
size = c("size3added", "size2added", "size1added"),
repmatrix = rep_cyp_raw, overwrite = cypover3r,
yearcol = "year2", patchcol = "patchid", indivcol = "individ")
summary(cypmatrix3r)
#>
#> This lefkoMat object contains 12 matrices.
#>
#> Each matrix is a square matrix with 121 rows and columns, and a total of 14641 elements.
#> A total of 388 survival-transitions were estimated, with 32.3333333333333 per matrix.
#> A total of 72 fecundity transitions were estimated, with 6 per matrix.
#>
#> The dataset contains a total of 74 unique individuals and 320 unique transitions.
#> NULLThere are at least two things to note here. First, there are 3 fewer matrices here than in the ahistorical case. There are 3 patches that we are estimating matrices for, and 6 years of data for each patch, leading to 5 possible ahistorical time transitions and 15 possible ahistorical matrices. Since historical matrices require 3 years of transitions, that means that only 4 historical transitions are possible per patch, leading to 12 total matrices. Second, the dimensionality of the matrices is the square of the dimensionality of the ahistorical matrices. This leads to vastly more matrix elements within each matrix, although it turns out that most of these matrix elements are structural 0s because they reflect impossible transitions. Indeed, in this case, although there are 14641 elements in each matrix, on average only 38.333 are actually estimated, with the rest set to 0.
Let’s look at the first matrix, corresponding to the transition from 2004 and 2005 to 2006 in the first patch. Because this is a huge matrix, we will only look at the top corner, and a middle section.
cypmatrix3r$A[[1]][1:25,1:10]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 0.1 0.0 0 0 0 0 0 0 0 0
#> [2,] 0.2 0.0 0 0 0 0 0 0 0 0
#> [3,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [4,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [5,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [6,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [7,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [8,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [9,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [10,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [11,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [12,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [13,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [14,] 0.0 0.2 0 0 0 0 0 0 0 0
#> [15,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [16,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [17,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [18,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [19,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [20,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [21,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [22,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [23,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [24,] 0.0 0.0 0 0 0 0 0 0 0 0
#> [25,] 0.0 0.0 0 0 0 0 0 0 0 0
cypmatrix3r$A[[1]][66:90,73:81]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [2,] 0.07142857 0.0000000 0 0 0 0 0 0 0
#> [3,] 0.07142857 0.0000000 0 0 0 0 0 0 0
#> [4,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [5,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [6,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [7,] 0.14285714 0.0000000 0 0 0 0 0 0 0
#> [8,] 0.71428571 0.0000000 0 0 0 0 0 0 0
#> [9,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [10,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [11,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [12,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [13,] 0.00000000 0.3333333 0 0 0 0 0 0 0
#> [14,] 0.00000000 0.3333333 0 0 0 0 0 0 0
#> [15,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [16,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [17,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [18,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [19,] 0.00000000 0.6666667 0 0 0 0 0 0 0
#> [20,] 0.00000000 0.3333333 0 0 0 0 0 0 0
#> [21,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [22,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [23,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [24,] 0.00000000 0.0000000 0 0 0 0 0 0 0
#> [25,] 0.00000000 0.0000000 0 0 0 0 0 0 0The full matrix is actually not shown here, but we can focus on portions of it if we wish. These matrices may also be exported to Excel or another spreadsheet program to look over in detail.
Next we will create a mean ahistorical matrix. Note that we will treat time arithmetically in this mean, meaning that matrix elements in the mean matrix are arithmetic means across time. With regards to which should be estimated, arithmetic or geometric mean: Our own feeling is that the geometric mean is a better choice in principle when taking a mean matrix elements across time, but the literature is generally full of the former case. The reason that the geometric mean may be better is that survival and reproduction are multiplicative over time, and the overall impact of variability in multiplicative processes is not symmetrical with respect to outcomes. For example, if I take the arithmetic mean matrix of all 4 matrices resulting from a 5 year study, and multiply the matrix four times over by the original vector of individuals across their respective stages, then I will end up projecting more individuals alive in the 5th year than actually exist! In multiplicative (or geometric) processes, the impact of lower values on the outcome is greater than the impact of higher values. Thus, in principle, the geometric mean is preferable. However, in practice, occasional lack of individuals moving through some stages may lead to sparse instances of elements dropping to 0. These situations cause the geometric mean to drop to 0, as well. In raw matrix projection, such situations may end up leading to geometric mean matrices overwhelmingly composed of 0s if the stages are not designed to absorb as many individuals as possible in each year. This can lead to a situation in the geometric mean projects a far lower number of individuals alive after some period.
The lmean() function in lefko3 incorporates a number of interesting options to help deal with these situations. Let’s do a comparison. First the pure arithmetic grand mean matrix, which is the last A matrix in this object.
cypr2mean <- lmean(cypmatrix2r, time = "arithmetic")
cypr2mean$A[[4]]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
#> [1,] 0.1 0.0 0.0 0.00000000 0.0 0.00000000 0.007870370 0.07174603 0.18898148 0.29238095 0.46666667
#> [2,] 0.2 0.0 0.0 0.00000000 0.0 0.00000000 0.007870370 0.07174603 0.18898148 0.29238095 0.46666667
#> [3,] 0.0 0.2 0.0 0.00000000 0.0 0.00000000 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [4,] 0.0 0.0 0.2 0.00000000 0.0 0.00000000 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [5,] 0.0 0.0 0.0 0.25000000 0.4 0.00000000 0.000000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [6,] 0.0 0.0 0.0 0.02222222 0.0 0.02222222 0.067255892 0.02814815 0.01944444 0.00000000 0.00000000
#> [7,] 0.0 0.0 0.0 0.21111111 0.0 0.21111111 0.519745070 0.22907407 0.06944444 0.06666667 0.00000000
#> [8,] 0.0 0.0 0.0 0.15555556 0.0 0.15555556 0.254737855 0.40314815 0.19611111 0.18222222 0.02222222
#> [9,] 0.0 0.0 0.0 0.00000000 0.0 0.05555556 0.007407407 0.21351852 0.48814815 0.16793651 0.00000000
#> [10,] 0.0 0.0 0.0 0.00000000 0.0 0.00000000 0.037037037 0.06500000 0.22685185 0.54349206 0.08888889
#> [11,] 0.0 0.0 0.0 0.00000000 0.0 0.02222222 0.000000000 0.00000000 0.00000000 0.03968254 0.48888889
writeLines("\nColumn sums (survival probabilities) for grand arithmetic mean matrix")
#>
#> Column sums (survival probabilities) for grand arithmetic mean matrix
summary(colSums(cypr2mean$U[[4]]))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.2000 0.3500 0.6000 0.6028 0.9125 1.0000For comparison, let’s look first at a geometric mean matrix.
lmean(cypmatrix2r, time = "geometric")$A[[4]]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
#> [1,] 0.1 0.0 0.0 0.00 0.0 0 0.0000000 0.00000000 0.0000000 0.1123006 0.0000000
#> [2,] 0.2 0.0 0.0 0.00 0.0 0 0.0000000 0.00000000 0.0000000 0.1123006 0.0000000
#> [3,] 0.0 0.2 0.0 0.00 0.0 0 0.0000000 0.00000000 0.0000000 0.0000000 0.0000000
#> [4,] 0.0 0.0 0.2 0.00 0.0 0 0.0000000 0.00000000 0.0000000 0.0000000 0.0000000
#> [5,] 0.0 0.0 0.0 0.25 0.4 0 0.0000000 0.00000000 0.0000000 0.0000000 0.0000000
#> [6,] 0.0 0.0 0.0 0.00 0.0 0 0.0000000 0.00000000 0.0000000 0.0000000 0.0000000
#> [7,] 0.0 0.0 0.0 0.00 0.0 0 0.4952254 0.08359622 0.0000000 0.0000000 0.0000000
#> [8,] 0.0 0.0 0.0 0.00 0.0 0 0.0000000 0.31081295 0.0000000 0.0000000 0.0000000
#> [9,] 0.0 0.0 0.0 0.00 0.0 0 0.0000000 0.00000000 0.1469767 0.0000000 0.0000000
#> [10,] 0.0 0.0 0.0 0.00 0.0 0 0.0000000 0.00000000 0.0000000 0.3689368 0.0000000
#> [11,] 0.0 0.0 0.0 0.00 0.0 0 0.0000000 0.00000000 0.0000000 0.0000000 0.2901835
writeLines("\nColumn sums (survival probabilities) for grand geometric mean matrix")
#>
#> Column sums (survival probabilities) for grand geometric mean matrix
summary(colSums(lmean(cypmatrix2r, time = "geometric")$U[[4]]))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.0000 0.2000 0.2902 0.2769 0.3817 0.4952Note that this matrix suggests extremely low survival probabilities relative to the arithmetic case. This is due to the fact that our data is relatively sparse. So, a single year in which no individuals start in a particular stage will result in a zero entry in transition probabilities starting from that stage, leading the geometric mean matrix to have zero entries in that column. In comparison, we can look at a geometric mean that scraps 0s unless they are structural.
lmean(cypmatrix2r, time = "geometric", sparse = TRUE)$A[[4]]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
#> [1,] 0.1 0.0 0.0 0.0000000 0.0 0.0000000 0.03935185 0.1793293 0.28972162 0.38437657 0.5206737
#> [2,] 0.2 0.0 0.0 0.0000000 0.0 0.0000000 0.03935185 0.1793293 0.28972162 0.38437657 0.5206737
#> [3,] 0.0 0.2 0.0 0.0000000 0.0 0.0000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
#> [4,] 0.0 0.0 0.2 0.0000000 0.0 0.0000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
#> [5,] 0.0 0.0 0.0 0.2500000 0.4 0.0000000 0.00000000 0.0000000 0.00000000 0.00000000 0.0000000
#> [6,] 0.0 0.0 0.0 0.1111111 0.0 0.1111111 0.20004424 0.1407407 0.09722222 0.00000000 0.0000000
#> [7,] 0.0 0.0 0.0 0.5088982 0.0 0.5088982 0.49522542 0.2967068 0.34722222 0.33333333 0.0000000
#> [8,] 0.0 0.0 0.0 0.5257834 0.0 0.5257834 0.36473518 0.4117641 0.43129808 0.51356794 0.1111111
#> [9,] 0.0 0.0 0.0 0.0000000 0.0 0.2777778 0.03703704 0.2888765 0.60017161 0.39258730 0.0000000
#> [10,] 0.0 0.0 0.0 0.0000000 0.0 0.0000000 0.14814815 0.1933207 0.39360704 0.56387133 0.3027494
#> [11,] 0.0 0.0 0.0 0.0000000 0.0 0.1111111 0.00000000 0.0000000 0.00000000 0.08908708 0.5031646
writeLines("\nColumn sums (survival probabilities) for grand geometric mean matrix")
#>
#> Column sums (survival probabilities) for grand geometric mean matrix
summary(colSums(lmean(cypmatrix2r, time = "geometric", sparse = TRUE)$U[[4]]))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.200 0.350 1.245 1.026 1.465 1.892Clearly, no stage can have a probability of survival over 1.0, but eliminating 0 entries from the estimation of mean matrices can sometimes result in this effect under both arithmetic and geometric scenarios. So, here we scrap this approach in favor of the arithmetic mean.
And now the historical grand mean matrix, with a peek at a middle portion with some non-zero values.
cypr3mean <- lmean(cypmatrix3r, time = "arithmetic")
cypr3mean$A[[4]][66:90,73:80]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] 0.000000000 0.00000000 0 0 0 0 0 0
#> [2,] 0.005952381 0.00000000 0 0 0 0 0 0
#> [3,] 0.005952381 0.00000000 0 0 0 0 0 0
#> [4,] 0.000000000 0.00000000 0 0 0 0 0 0
#> [5,] 0.000000000 0.00000000 0 0 0 0 0 0
#> [6,] 0.000000000 0.00000000 0 0 0 0 0 0
#> [7,] 0.136904762 0.00000000 0 0 0 0 0 0
#> [8,] 0.537301587 0.00000000 0 0 0 0 0 0
#> [9,] 0.225000000 0.00000000 0 0 0 0 0 0
#> [10,] 0.000000000 0.00000000 0 0 0 0 0 0
#> [11,] 0.000000000 0.00000000 0 0 0 0 0 0
#> [12,] 0.000000000 0.00000000 0 0 0 0 0 0
#> [13,] 0.000000000 0.02777778 0 0 0 0 0 0
#> [14,] 0.000000000 0.02777778 0 0 0 0 0 0
#> [15,] 0.000000000 0.00000000 0 0 0 0 0 0
#> [16,] 0.000000000 0.00000000 0 0 0 0 0 0
#> [17,] 0.000000000 0.00000000 0 0 0 0 0 0
#> [18,] 0.000000000 0.00000000 0 0 0 0 0 0
#> [19,] 0.000000000 0.30138889 0 0 0 0 0 0
#> [20,] 0.000000000 0.09861111 0 0 0 0 0 0
#> [21,] 0.000000000 0.10000000 0 0 0 0 0 0
#> [22,] 0.000000000 0.02777778 0 0 0 0 0 0
#> [23,] 0.000000000 0.00000000 0 0 0 0 0 0
#> [24,] 0.000000000 0.00000000 0 0 0 0 0 0
#> [25,] 0.000000000 0.00000000 0 0 0 0 0 0Do not fear the prevalence of 0s in this matrix - this is normal, both because many elements are structural 0s and so cannot equal anything else, and because this is a raw matrix, meaning that transitions that do not actually occur in the dataset cannot equal anything other than 0.
To understand the dominance of structural 0s in the historical case, let’s take a look at the hstages object associated with this mean matrix.
cypr3mean$hstages
#> stcod3 stcod2 stage3 stage2
#> 1 SD SD 1 1
#> 2 P1 SD 2 1
#> 3 P2 SD 3 1
#> 4 P3 SD 4 1
#> 5 SL SD 5 1
#> 6 D SD 6 1
#> 7 XSm SD 7 1
#> 8 Sm SD 8 1
#> 9 Md SD 9 1
#> 10 Lg SD 10 1
#> 11 XLg SD 11 1
#> 13 SD P1 1 2
#> 14 P1 P1 2 2
#> 15 P2 P1 3 2
#> 16 P3 P1 4 2
#> 17 SL P1 5 2
#> 18 D P1 6 2
#> 19 XSm P1 7 2
#> 20 Sm P1 8 2
#> 21 Md P1 9 2
#> 22 Lg P1 10 2
#> 23 XLg P1 11 2
#> 25 SD P2 1 3
#> 26 P1 P2 2 3
#> 27 P2 P2 3 3
#> 28 P3 P2 4 3
#> 29 SL P2 5 3
#> 30 D P2 6 3
#> 31 XSm P2 7 3
#> 32 Sm P2 8 3
#> 33 Md P2 9 3
#> 34 Lg P2 10 3
#> 35 XLg P2 11 3
#> 37 SD P3 1 4
#> 38 P1 P3 2 4
#> 39 P2 P3 3 4
#> 40 P3 P3 4 4
#> 41 SL P3 5 4
#> 42 D P3 6 4
#> 43 XSm P3 7 4
#> 44 Sm P3 8 4
#> 45 Md P3 9 4
#> 46 Lg P3 10 4
#> 47 XLg P3 11 4
#> 49 SD SL 1 5
#> 50 P1 SL 2 5
#> 51 P2 SL 3 5
#> 52 P3 SL 4 5
#> 53 SL SL 5 5
#> 54 D SL 6 5
#> 55 XSm SL 7 5
#> 56 Sm SL 8 5
#> 57 Md SL 9 5
#> 58 Lg SL 10 5
#> 59 XLg SL 11 5
#> 61 SD D 1 6
#> 62 P1 D 2 6
#> 63 P2 D 3 6
#> 64 P3 D 4 6
#> 65 SL D 5 6
#> 66 D D 6 6
#> 67 XSm D 7 6
#> 68 Sm D 8 6
#> 69 Md D 9 6
#> 70 Lg D 10 6
#> 71 XLg D 11 6
#> 73 SD XSm 1 7
#> 74 P1 XSm 2 7
#> 75 P2 XSm 3 7
#> 76 P3 XSm 4 7
#> 77 SL XSm 5 7
#> 78 D XSm 6 7
#> 79 XSm XSm 7 7
#> 80 Sm XSm 8 7
#> 81 Md XSm 9 7
#> 82 Lg XSm 10 7
#> 83 XLg XSm 11 7
#> 85 SD Sm 1 8
#> 86 P1 Sm 2 8
#> 87 P2 Sm 3 8
#> 88 P3 Sm 4 8
#> 89 SL Sm 5 8
#> 90 D Sm 6 8
#> 91 XSm Sm 7 8
#> 92 Sm Sm 8 8
#> 93 Md Sm 9 8
#> 94 Lg Sm 10 8
#> 95 XLg Sm 11 8
#> 97 SD Md 1 9
#> 98 P1 Md 2 9
#> 99 P2 Md 3 9
#> 100 P3 Md 4 9
#> 101 SL Md 5 9
#> 102 D Md 6 9
#> 103 XSm Md 7 9
#> 104 Sm Md 8 9
#> 105 Md Md 9 9
#> 106 Lg Md 10 9
#> 107 XLg Md 11 9
#> 109 SD Lg 1 10
#> 110 P1 Lg 2 10
#> 111 P2 Lg 3 10
#> 112 P3 Lg 4 10
#> 113 SL Lg 5 10
#> 114 D Lg 6 10
#> 115 XSm Lg 7 10
#> 116 Sm Lg 8 10
#> 117 Md Lg 9 10
#> 118 Lg Lg 10 10
#> 119 XLg Lg 11 10
#> 121 SD XLg 1 11
#> 122 P1 XLg 2 11
#> 123 P2 XLg 3 11
#> 124 P3 XLg 4 11
#> 125 SL XLg 5 11
#> 126 D XLg 6 11
#> 127 XSm XLg 7 11
#> 128 Sm XLg 8 11
#> 129 Md XLg 9 11
#> 130 Lg XLg 10 11
#> 131 XLg XLg 11 11There are 121 pairs of ahistorical stages, and these pairs correspond to the rows and columns of the historical matrices output by rlefko3 in this case. The pairs are interpreted so that matrix columns represent the states of the individual in times t-1 and t, and the rows represent states in times t and t+1. For an element in the matrix to contain a number other than 0, it must represent the same stage at time t in both the column stage pairs and the row stage pairs. The element [1, 1], for example, represents the transition probability from dormant seed at times t-1 and t (column pair) to times t and t+1 (row pair) - the time t stages match, and so this element is possible. However, element [1, 2] represents the transition probability from seed in time t-1 and protocorm 1 in time t (column pair), to dormant seed in time t and in time t+1 (row pair). Clearly [1, 2] is a structural 0 because it is impossible for an individual to be in two stages at once in time t!
Error-checking is more difficult with historical matrices because they are typically one or two orders of magnitude bigger than their ahistorical counterparts. Take a look at the column sums here to see the difficulty.
colSums(cypr3mean$U[[4]])
#> [1] 0.30000000 0.20000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [13] 0.00000000 0.20000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [25] 0.00000000 0.33333333 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [37] 0.00000000 0.40000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [49] 0.40000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [61] 0.08333333 0.16666667 0.16666667 0.00000000 0.00000000 0.08333333 0.30000000 0.00000000 0.00000000 0.00000000 0.00000000 0.50000000
#> [73] 0.89920635 0.52777778 0.08333333 0.25000000 0.00000000 0.30000000 0.00000000 0.00000000 0.00000000 0.00000000 0.16666667 0.77500000
#> [85] 0.88888889 0.50000000 0.41666667 0.00000000 0.30000000 0.00000000 0.00000000 0.00000000 0.00000000 0.16666667 0.25000000 0.50000000
#> [97] 0.75000000 0.50000000 0.00000000 0.30000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.08333333 0.41666667 0.50000000
#> [109] 0.83333333 0.08333333 0.30000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.08333333 0.00000000 0.16666667
#> [121] 0.58333333While not too bad here, other historical matrices often have more than 100 columns, sometimes many, many more (some historical matrices used in Shefferson et al. (2014) had dimensions of over 2500 x 2500!). In these cases we can assess the distribution of survival estimates for historical stages, which is given as the set of column sums in the survival-transition matrix, as below.
summary(colSums(cypr3mean$U[[4]]))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.0000 0.0000 0.0000 0.1137 0.1667 0.8992As long as all of the numbers above are between 0 and 1, then all is probably well. Fine-scale error-checking would require outputting the matrix into a spreadsheet and assessing it using the hstages output as a guide to what the elements refer to.
Step 4. Analyses
Now let’s estimate the deterministic population growth rate in each case. We will start by looking at the annual population growth rate estimated from the ahistorical analyses, followed by the population growth rate associated with the mean matrix from that analysis. First the annual ahistorical matrices.
lambda3(cypmatrix2r)
#> pop patch year2 lambda
#> 1 NA A 2004 0.9782070
#> 2 NA A 2005 1.0000000
#> 3 NA A 2006 1.0000000
#> 4 NA A 2007 1.0000000
#> 5 NA A 2008 1.0000000
#> 6 NA B 2004 0.8964607
#> 7 NA B 2005 0.8597354
#> 8 NA B 2006 1.0000000
#> 9 NA B 2007 1.0065871
#> 10 NA B 2008 1.0055269
#> 11 NA C 2004 1.0000000
#> 12 NA C 2005 1.0000000
#> 13 NA C 2006 1.0000000
#> 14 NA C 2007 1.0000000
#> 15 NA C 2008 1.0039623Now the ahistorical mean matrix.
lambda3(cypr2mean)
#> pop patch lambda
#> 1 1 A 0.8992573
#> 2 1 B 0.9321193
#> 3 1 C 0.9577990
#> 4 1 All 0.9332170And we will finish by looking at the same numbers for the historical analyses. First the annual matrices.
lambda3(cypmatrix3r)
#> pop patch year2 lambda
#> 1 NA A 2005 0.7142857
#> 2 NA A 2006 1.0000000
#> 3 NA A 2007 1.0000000
#> 4 NA A 2008 1.0000000
#> 5 NA B 2005 0.6666667
#> 6 NA B 2006 1.0000000
#> 7 NA B 2007 1.0000000
#> 8 NA B 2008 1.0000000
#> 9 NA C 2005 1.0000000
#> 10 NA C 2006 1.0000000
#> 11 NA C 2007 1.0000000
#> 12 NA C 2008 1.0000000Now the mean.
lambda3(cypr3mean)
#> pop patch lambda
#> 1 1 A 0.7802134
#> 2 1 B 0.7601081
#> 3 1 C 1.0000000
#> 4 1 All 0.6905696Readers will likely observe both that there are fewer lambda estimates in the historical case, and that the mean lambda is typically lower. Both are expected. First, because there are 6 years of data, there are three ahistorical transitions possible for estimation: year 1 to 2, year 2 to 3, and year 3 to 4. However, in the historical case, only two are possible: from years 1 and 2 to years 2 and 3, and from years 2 and 3 to years 3 and 4. Second, historical matrices integrate temporal autocorrelation in vital rates in ways that ahistorical matrices generally cannot, and these autocorrelations are likely to be most strongly impacted by trade-offs operating across years (Shefferson and Roach 2010). One particularly common such trade-off is the cost of growth: an individual that grows a great deal in one time step due to great environmental conditions in that year might pay a large cost of survival, growth, or reproduction in the next if those environmental conditions deteriorate (Shefferson, Warren II, and Pulliam 2014). While we do not argue that the drop in lambda is due to this specific trade-off, we do argue that this lambda is likely to be more realistic than the higher lambda estimated in the ahistorical case.
We can also take a peek at the stable stage distributions, as follows for the ahistorical case:
stablestage3(cypr2mean)
#> matrix new_stage_id orig_stage_id original_size ss_prop
#> 1 1 1 SD 0.0 0.090037090
#> 2 1 2 P1 0.0 0.100049436
#> 3 1 3 P2 0.0 0.022251576
#> 4 1 4 P3 0.0 0.004948866
#> 5 1 5 SL 0.0 0.002478141
#> 6 1 6 D 0.0 0.030513881
#> 7 1 7 XSm 1.0 0.249352147
#> 8 1 8 Sm 2.5 0.290121522
#> 9 1 9 Md 4.5 0.101212097
#> 10 1 10 Lg 8.0 0.109035244
#> 11 1 11 XLg 17.5 0.000000000
#> 12 2 1 SD 0.0 0.147677868
#> 13 2 2 P1 0.0 0.163521115
#> 14 2 3 P2 0.0 0.035085868
#> 15 2 4 P3 0.0 0.007528201
#> 16 2 5 SL 0.0 0.003536880
#> 17 2 6 D 0.0 0.027764530
#> 18 2 7 XSm 1.0 0.172012760
#> 19 2 8 Sm 2.5 0.199776799
#> 20 2 9 Md 4.5 0.093236400
#> 21 2 10 Lg 8.0 0.109216888
#> 22 2 11 XLg 17.5 0.040642691
#> 23 3 1 SD 0.0 0.085711562
#> 24 3 2 P1 0.0 0.094660376
#> 25 3 3 P2 0.0 0.019766226
#> 26 3 4 P3 0.0 0.004127438
#> 27 3 5 SL 0.0 0.001849871
#> 28 3 6 D 0.0 0.020105859
#> 29 3 7 XSm 1.0 0.178308076
#> 30 3 8 Sm 2.5 0.182244061
#> 31 3 9 Md 4.5 0.221104786
#> 32 3 10 Lg 8.0 0.192121744
#> 33 3 11 XLg 17.5 0.000000000
#> 34 4 1 SD 0.0 0.114482899
#> 35 4 2 P1 0.0 0.126750486
#> 36 4 3 P2 0.0 0.027164217
#> 37 4 4 P3 0.0 0.005821612
#> 38 4 5 SL 0.0 0.002729480
#> 39 4 6 D 0.0 0.022921317
#> 40 4 7 XSm 1.0 0.177582482
#> 41 4 8 Sm 2.5 0.203366188
#> 42 4 9 Md 4.5 0.158605194
#> 43 4 10 Lg 8.0 0.146358608
#> 44 4 11 XLg 17.5 0.014217515We can do this for the historical case as well, but need to bear in mind that the distribution is in paired stages, and so can take some effort to interpret properly. We can look first at the paired stage distribution. Let’s see the summary and first 6 entries.
cypr3ss <- stablestage3(cypr3mean)
summary(cypr3ss$hist)
#> matrix orig_stage_id_2 orig_stage_id_1 new_stage_id_2 new_stage_id_1 ss_prop
#> Min. :1.00 Length:484 Length:484 Min. : 1 Min. : 1 Min. :0.0000000
#> 1st Qu.:1.75 Class :character Class :character 1st Qu.: 3 1st Qu.: 3 1st Qu.:0.0000000
#> Median :2.50 Mode :character Mode :character Median : 6 Median : 6 Median :0.0000000
#> Mean :2.50 Mean : 6 Mean : 6 Mean :0.0082645
#> 3rd Qu.:3.25 3rd Qu.: 9 3rd Qu.: 9 3rd Qu.:0.0000079
#> Max. :4.00 Max. :11 Max. :11 Max. :0.5273180
head(cypr3ss$hist)
#> matrix orig_stage_id_2 orig_stage_id_1 new_stage_id_2 new_stage_id_1 ss_prop
#> 1 1 SD SD 1 1 0.003083154
#> 2 1 P1 SD 2 1 0.006166308
#> 3 1 P2 SD 3 1 0.000000000
#> 4 1 P3 SD 4 1 0.000000000
#> 5 1 SL SD 5 1 0.000000000
#> 6 1 D SD 6 1 0.000000000And we can look at the resulting life history stable stage distribution associated with the historical case.
cypr3ss$ahist
#> matrix new_stage_id ss_prop
#> 1 1 1 0.0240552215
#> 2 1 2 0.0271383756
#> 3 1 3 0.0015806534
#> 4 1 4 0.0004052074
#> 5 1 5 0.0002664555
#> 6 1 6 0.0578366170
#> 7 1 7 0.4955790205
#> 8 1 8 0.3584457196
#> 9 1 9 0.0317567933
#> 10 1 10 0.0029359364
#> 11 1 11 0.0000000000
#> 12 2 1 0.0931549377
#> 13 2 2 0.1054104247
#> 14 2 3 0.0064493456
#> 15 2 4 0.0016969480
#> 16 2 5 0.0011780848
#> 17 2 6 0.0807390019
#> 18 2 7 0.2526776338
#> 19 2 8 0.3393131233
#> 20 2 9 0.0370767622
#> 21 2 10 0.0518189154
#> 22 2 11 0.0304848226
#> 23 3 1 0.2197158418
#> 24 3 2 0.2416874140
#> 25 3 3 0.0087886169
#> 26 3 4 0.0017577354
#> 27 3 5 0.0007323598
#> 28 3 6 0.0000000000
#> 29 3 7 0.0000000000
#> 30 3 8 0.0000000000
#> 31 3 9 0.0000000000
#> 32 3 10 0.0000000000
#> 33 3 11 0.5273180322
#> 34 4 1 0.0662983940
#> 35 4 2 0.0758989111
#> 36 4 3 0.0055609200
#> 37 4 4 0.0016105376
#> 38 4 5 0.0013856524
#> 39 4 6 0.0674429087
#> 40 4 7 0.4016556123
#> 41 4 8 0.2403823061
#> 42 4 9 0.0930652020
#> 43 4 10 0.0453846335
#> 44 4 11 0.0013149224These values are considerably different from those of the purely ahistorical stable stage distribution. Indeed, here seeds and protocorms are a much smaller part of the stable stage structure, while dormant, extra small, and small adults are a much larger component. It seems that individual history has a strong impact here.
Let’s look at the reproductive values next. First, the ahistorical case.
repvalue3(cypr2mean)
#> matrix new_stage_id orig_stage_id original_size left_vector rep_value
#> 1 1 1 SD 0.0 0.0014030 1.000000
#> 2 1 2 P1 0.0 0.0056066 3.996151
#> 3 1 3 P2 0.0 0.0252089 17.967855
#> 4 1 4 P3 0.0 0.1133465 80.788667
#> 5 1 5 SL 0.0 0.0000000 0.000000
#> 6 1 6 D 0.0 0.1533835 109.325374
#> 7 1 7 XSm 1.0 0.2536043 180.758589
#> 8 1 8 Sm 2.5 0.3188277 227.247113
#> 9 1 9 Md 4.5 0.5400528 384.927156
#> 10 1 10 Lg 8.0 0.7108684 506.677406
#> 11 1 11 XLg 17.5 0.0000000 0.000000
#> 12 2 1 SD 0.0 -0.0030046 1.000000
#> 13 2 2 P1 0.0 -0.0125010 4.160620
#> 14 2 3 P2 0.0 -0.0582622 19.391000
#> 15 2 4 P3 0.0 -0.2715368 90.373694
#> 16 2 5 SL 0.0 0.0000000 0.000000
#> 17 2 6 D 0.0 -0.2952139 98.253977
#> 18 2 7 XSm 1.0 -0.3418877 113.788092
#> 19 2 8 Sm 2.5 -0.3693806 122.938361
#> 20 2 9 Md 4.5 -0.4824601 160.573820
#> 21 2 10 Lg 8.0 -0.4897617 163.003961
#> 22 2 11 XLg 17.5 -0.3310488 110.180656
#> 23 3 1 SD 0.0 -0.0003103 1.000000
#> 24 3 2 P1 0.0 -0.0013308 4.288753
#> 25 3 3 P2 0.0 -0.0063730 20.538189
#> 26 3 4 P3 0.0 -0.0305204 98.357718
#> 27 3 5 SL 0.0 0.0000000 0.000000
#> 28 3 6 D 0.0 -0.0713462 229.926523
#> 29 3 7 XSm 1.0 -0.2923241 942.069288
#> 30 3 8 Sm 2.5 -0.3880920 1250.699323
#> 31 3 9 Md 4.5 -0.3910292 1260.165002
#> 32 3 10 Lg 8.0 -0.3836712 1236.452465
#> 33 3 11 XLg 17.5 -0.6765793 2180.403803
#> 34 4 1 SD 0.0 -0.0018065 1.000000
#> 35 4 2 P1 0.0 -0.0075259 4.166012
#> 36 4 3 P2 0.0 -0.0351166 19.439026
#> 37 4 4 P3 0.0 -0.1638568 90.704013
#> 38 4 5 SL 0.0 0.0000000 0.000000
#> 39 4 6 D 0.0 -0.1992797 110.312593
#> 40 4 7 XSm 1.0 -0.3691083 204.322336
#> 41 4 8 Sm 2.5 -0.4536172 251.102795
#> 42 4 9 Md 4.5 -0.5397128 298.761583
#> 43 4 10 Lg 8.0 -0.5288866 292.768669
#> 44 4 11 XLg 17.5 -0.1382934 76.553224Now the historical case, beginning with the stage pair distribution. We will first look at the summary.
cypr3rv <- repvalue3(cypr3mean)
summary(cypr3rv$hist)
#> matrix orig_stage_id_2 orig_stage_id_1 new_stage_id_2 new_stage_id_1 left_vector rep_value
#> Min. :1.00 Length:484 Length:484 Min. : 1 Min. : 1 Min. :-0.561342 Min. : 0.0000
#> 1st Qu.:1.75 Class :character Class :character 1st Qu.: 3 1st Qu.: 3 1st Qu.: 0.000000 1st Qu.: 0.0000
#> Median :2.50 Mode :character Mode :character Median : 6 Median : 6 Median : 0.000000 Median : 0.0000
#> Mean :2.50 Mean : 6 Mean : 6 Mean :-0.006077 Mean : 0.1994
#> 3rd Qu.:3.25 3rd Qu.: 9 3rd Qu.: 9 3rd Qu.: 0.000000 3rd Qu.: 0.0000
#> Max. :4.00 Max. :11 Max. :11 Max. : 1.000000 Max. :27.5077Now the first 6 rows.
head(cypr3rv$hist)
#> matrix orig_stage_id_2 orig_stage_id_1 new_stage_id_2 new_stage_id_1 left_vector rep_value
#> 1 1 SD SD 1 1 0 0
#> 2 1 P1 SD 2 1 0 0
#> 3 1 P2 SD 3 1 0 0
#> 4 1 P3 SD 4 1 0 0
#> 5 1 SL SD 5 1 0 0
#> 6 1 D SD 6 1 0 0And now the life history stage distribution associated with the historical case.
cypr3rv$ahist
#> matrix new_stage_id rep_value
#> 1 1 1 0.00000000
#> 2 1 2 0.00000000
#> 3 1 3 0.00000000
#> 4 1 4 0.00000000
#> 5 1 5 0.00000000
#> 6 1 6 1.00000000
#> 7 1 7 4.50823457
#> 8 1 8 9.00497739
#> 9 1 9 0.62246947
#> 10 1 10 0.06189564
#> 11 1 11 0.00000000
#> 12 2 1 0.00000000
#> 13 2 2 0.00000000
#> 14 2 3 0.00000000
#> 15 2 4 0.00000000
#> 16 2 5 0.00000000
#> 17 2 6 1.00000000
#> 18 2 7 3.00926000
#> 19 2 8 3.88439088
#> 20 2 9 0.88071161
#> 21 2 10 1.12439876
#> 22 2 11 0.02567988
#> 23 3 1 0.00000000
#> 24 3 2 0.00000000
#> 25 3 3 0.00000000
#> 26 3 4 0.00000000
#> 27 3 5 0.00000000
#> 28 3 6 0.00000000
#> 29 3 7 0.00000000
#> 30 3 8 0.00000000
#> 31 3 9 0.00000000
#> 32 3 10 0.00000000
#> 33 3 11 1.00000000
#> 34 4 1 0.00000000
#> 35 4 2 0.00000000
#> 36 4 3 0.00000000
#> 37 4 4 0.00000000
#> 38 4 5 0.00000000
#> 39 4 6 1.00000000
#> 40 4 7 6.64583147
#> 41 4 8 4.53914861
#> 42 4 9 1.35352344
#> 43 4 10 0.77936516
#> 44 4 11 0.01444315
The reproductive values of juveniles stages have essentially dropped to 0, leaving only the adult stages, starting from the dormant stage. So, indeed, we see dramatically different reproductive values here.
Now on to function-based matrices!
Analysis 2. Function-based matrix estimation
Step 1. Data characterization and reorganization
Before we go further, we need to describe the life history characterizing dataset properly with a stageframe for our Cypripedium candidum dataset, which will be different from the one that we created for the raw matrix example. Since this analysis will be function-based, we will include all possible size classes here. If constructing raw matrices, as in the previous example, all sizes that occur in the dataset need to be accounted for in a way that is both natural and parsimonious with respect to transition estimation. If constructing function-based matrices, such as IPMs, then representative sizes at systematic increments will be satisfactory. Since size is count-based in the Cypripedium candidum case, we will use all numbers of stems that might occur from 0 to the maximum in the dataset, representing the life history diagram shown in the beginning of this chapter.
sizevector <- c(0, 0, 0, 0, 0, seq(from = 0, t = 24), seq(from = 1, to = 24))
stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "V1", "V2", "V3", "V4", "V5",
"V6", "V7", "V8", "V9", "V10", "V11", "V12", "V13", "V14", "V15",
"V16", "V17", "V18", "V19", "V20", "V21", "V22", "V23", "V24",
"F1", "F2", "F3", "F4", "F5", "F6", "F7", "F8", "F9", "F10",
"F11", "F12", "F13", "F14", "F15", "F16", "F17", "F18", "F19",
"F20", "F21", "F22", "F23", "F24")
repvector <- c(0, 0, 0, 0, 0, rep(0, 25), rep(1, 24))
obsvector <- c(0, 0, 0, 0, 0, 0, rep(1, 48))
matvector <- c(0, 0, 0, 0, 0, rep(1, 49))
immvector <- c(0, 1, 1, 1, 1, rep(0, 49))
propvector <- c(1, rep(0, 53))
indataset <- c(0, 0, 0, 0, 0, rep(1, 49))
cypframe <- sf_create(sizes = sizevector, stagenames = stagevector,
repstatus = repvector, obsstatus = obsvector,
matstatus = matvector, propstatus = propvector,
immstatus = immvector, indataset = indataset)
A close look at the output object, cypframe, shows a data frame that includes in order for each stage: the stage’s name, the associated size, its reproductive status, its status as an observable stage, its status as a propagule stage, its status as an immature stage, its status as a mature stage, whether it occurs in the dataset, the half-width of a size class bin, the minima and maxima of size class bins, the centroid of the size class bin, the full size class bin width, and comments. Stage names and combinations of characteristics must be unique to prevent estimation errors, and the comments field may be edited to include any information deemed pertinent. We may edit the comments field as below.
cypframe$comments[(cypframe$stagenames == "SD")] <- "Dormant seed"
cypframe$comments[(cypframe$stagenames == "P1")] <- "1st yr protocorm"
cypframe$comments[(cypframe$stagenames == "P2")] <- "2nd yr protocorm"
cypframe$comments[(cypframe$stagenames == "P3")] <- "3rd yr protocorm"
cypframe$comments[(cypframe$stagenames == "SL")] <- "Seedling"
cypframe$comments[(cypframe$stagenames == "D")] <- "Dormant mature"
cypframe$comments[(cypframe$stagenames == "V1")] <- "Non-reproductive mature individual with 1 stem"
cypframe$comments[(cypframe$stagenames == "V2")] <- "Non-reproductive mature individual with 2 stems"
cypframe$comments[(cypframe$stagenames == "V3")] <- "Non-reproductive mature individual with 3 stems"
cypframe$comments[(cypframe$stagenames == "V4")] <- "Non-reproductive mature individual with 4 stems"
cypframe$comments[(cypframe$stagenames == "V5")] <- "Non-reproductive mature individual with 5 stems"
cypframe$comments[(cypframe$stagenames == "V6")] <- "Non-reproductive mature individual with 6 stems"
cypframe$comments[(cypframe$stagenames == "V7")] <- "Non-reproductive mature individual with 7 stems"
cypframe$comments[(cypframe$stagenames == "V8")] <- "Non-reproductive mature individual with 8 stems"
cypframe$comments[(cypframe$stagenames == "V9")] <- "Non-reproductive mature individual with 9 stems"
cypframe$comments[(cypframe$stagenames == "V10")] <- "Non-reproductive mature individual with 10 stems"
cypframe$comments[(cypframe$stagenames == "V11")] <- "Non-reproductive mature individual with 11 stems"
cypframe$comments[(cypframe$stagenames == "V12")] <- "Non-reproductive mature individual with 12 stems"
cypframe$comments[(cypframe$stagenames == "V13")] <- "Non-reproductive mature individual with 13 stems"
cypframe$comments[(cypframe$stagenames == "V14")] <- "Non-reproductive mature individual with 14 stems"
cypframe$comments[(cypframe$stagenames == "V15")] <- "Non-reproductive mature individual with 15 stems"
cypframe$comments[(cypframe$stagenames == "V16")] <- "Non-reproductive mature individual with 16 stems"
cypframe$comments[(cypframe$stagenames == "V17")] <- "Non-reproductive mature individual with 17 stems"
cypframe$comments[(cypframe$stagenames == "V18")] <- "Non-reproductive mature individual with 18 stems"
cypframe$comments[(cypframe$stagenames == "V19")] <- "Non-reproductive mature individual with 19 stems"
cypframe$comments[(cypframe$stagenames == "V20")] <- "Non-reproductive mature individual with 20 stems"
cypframe$comments[(cypframe$stagenames == "V21")] <- "Non-reproductive mature individual with 21 stems"
cypframe$comments[(cypframe$stagenames == "V22")] <- "Non-reproductive mature individual with 22 stems"
cypframe$comments[(cypframe$stagenames == "V23")] <- "Non-reproductive mature individual with 23 stems"
cypframe$comments[(cypframe$stagenames == "V24")] <- "Non-reproductive mature individual with 24 stems"
cypframe$comments[(cypframe$stagenames == "F1")] <- "Flowering mature individual with 1 stem"
cypframe$comments[(cypframe$stagenames == "F2")] <- "Flowering mature individual with 2 stems"
cypframe$comments[(cypframe$stagenames == "F3")] <- "Flowering mature individual with 3 stems"
cypframe$comments[(cypframe$stagenames == "F4")] <- "Flowering mature individual with 4 stems"
cypframe$comments[(cypframe$stagenames == "F5")] <- "Flowering mature individual with 5 stems"
cypframe$comments[(cypframe$stagenames == "F6")] <- "Flowering mature individual with 6 stems"
cypframe$comments[(cypframe$stagenames == "F7")] <- "Flowering mature individual with 7 stems"
cypframe$comments[(cypframe$stagenames == "F8")] <- "Flowering mature individual with 8 stems"
cypframe$comments[(cypframe$stagenames == "F9")] <- "Flowering mature individual with 9 stems"
cypframe$comments[(cypframe$stagenames == "F10")] <- "Flowering mature individual with 10 stems"
cypframe$comments[(cypframe$stagenames == "F11")] <- "Flowering mature individual with 11 stems"
cypframe$comments[(cypframe$stagenames == "F12")] <- "Flowering mature individual with 12 stems"
cypframe$comments[(cypframe$stagenames == "F13")] <- "Flowering mature individual with 13 stems"
cypframe$comments[(cypframe$stagenames == "F14")] <- "Flowering mature individual with 14 stems"
cypframe$comments[(cypframe$stagenames == "F15")] <- "Flowering mature individual with 15 stems"
cypframe$comments[(cypframe$stagenames == "F16")] <- "Flowering mature individual with 16 stems"
cypframe$comments[(cypframe$stagenames == "F17")] <- "Flowering mature individual with 17 stems"
cypframe$comments[(cypframe$stagenames == "F18")] <- "Flowering mature individual with 18 stems"
cypframe$comments[(cypframe$stagenames == "F19")] <- "Flowering mature individual with 19 stems"
cypframe$comments[(cypframe$stagenames == "F20")] <- "Flowering mature individual with 20 stems"
cypframe$comments[(cypframe$stagenames == "F21")] <- "Flowering mature individual with 21 stems"
cypframe$comments[(cypframe$stagenames == "F22")] <- "Flowering mature individual with 22 stems"
cypframe$comments[(cypframe$stagenames == "F23")] <- "Flowering mature individual with 23 stems"
cypframe$comments[(cypframe$stagenames == "F24")] <- "Flowering mature individual with 24 stems"
Now we may transform our vertical dataset into a historically-formatted vertical file. The resulting dataset will have each individual’s observed life history broken up into states corresponding to three consecutive years per row, with plant identity marked in each row. To handle this, we use the verticalize3() function, as below.
vertdata <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004,
patchidcol = "patch", individcol = "plantid",
blocksize = 4, size1col = "Inf2.04", size2col = "Inf.04",
size3col = "Veg.04", repstr1col = "Inf.04",
fec1col = "Pod.04", stageassign = cypframe,
stagesize = "sizeadded", NAas0 = TRUE)
summary(vertdata)
#> rowid popid patchid individ year2 sizea1 sizea2 sizea3
#> Min. : 1.00 Mode:logical A: 93 Min. : 164.0 Min. :2004 Min. :0.000000 Min. :0.000000 Min. :0.000000
#> 1st Qu.:21.00 NA's:320 B:154 1st Qu.: 391.0 1st Qu.:2005 1st Qu.:0.000000 1st Qu.:0.000000 1st Qu.:0.000000
#> Median :37.50 C: 73 Median : 453.0 Median :2006 Median :0.000000 Median :0.000000 Median :0.000000
#> Mean :38.45 Mean : 651.5 Mean :2006 Mean :0.009375 Mean :0.009375 Mean :0.009375
#> 3rd Qu.:56.00 3rd Qu.: 476.0 3rd Qu.:2007 3rd Qu.:0.000000 3rd Qu.:0.000000 3rd Qu.:0.000000
#> Max. :77.00 Max. :1560.0 Max. :2008 Max. :1.000000 Max. :1.000000 Max. :1.000000
#> repstra1 repstra2 repstra3 feca1 feca2 feca3 sizeb1
#> Min. : 0.0000 Min. : 0.0000 Min. : 0.000 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. : 0.0000
#> 1st Qu.: 0.0000 1st Qu.: 0.0000 1st Qu.: 0.000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.: 0.0000
#> Median : 0.0000 Median : 0.0000 Median : 0.000 Median :0.0000 Median :0.0000 Median :0.0000 Median : 0.0000
#> Mean : 0.7469 Mean : 0.8969 Mean : 1.069 Mean :0.2656 Mean :0.2906 Mean :0.4562 Mean : 0.7469
#> 3rd Qu.: 1.0000 3rd Qu.: 1.0000 3rd Qu.: 1.000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.: 1.0000
#> Max. :18.0000 Max. :18.0000 Max. :18.000 Max. :7.0000 Max. :7.0000 Max. :8.0000 Max. :18.0000
#> sizeb2 sizeb3 sizec1 sizec2 sizec3 size1added size2added size3added
#> Min. : 0.0000 Min. : 0.000 Min. : 0.0 Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. : 0.000
#> 1st Qu.: 0.0000 1st Qu.: 0.000 1st Qu.: 0.0 1st Qu.: 1.000 1st Qu.: 1.000 1st Qu.: 0.000 1st Qu.: 1.000 1st Qu.: 1.000
#> Median : 0.0000 Median : 0.000 Median : 1.0 Median : 2.000 Median : 1.000 Median : 2.000 Median : 2.000 Median : 2.000
#> Mean : 0.8969 Mean : 1.069 Mean : 1.9 Mean : 2.416 Mean : 2.209 Mean : 2.656 Mean : 3.322 Mean : 3.288
#> 3rd Qu.: 1.0000 3rd Qu.: 1.000 3rd Qu.: 3.0 3rd Qu.: 3.000 3rd Qu.: 3.000 3rd Qu.: 4.000 3rd Qu.: 4.000 3rd Qu.: 4.000
#> Max. :18.0000 Max. :18.000 Max. :13.0 Max. :13.000 Max. :13.000 Max. :21.000 Max. :24.000 Max. :24.000
#> obsstatus1 obsstatus2 obsstatus3 repstatus1 repstatus2 repstatus3 fecstatus1 fecstatus2
#> Min. :0.0000 Min. :0.0000 Min. :0.0 Min. :0.0000 Min. :0.0000 Min. :0.0 Min. :0.0000 Min. :0.0000
#> 1st Qu.:0.0000 1st Qu.:1.0000 1st Qu.:1.0 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0 1st Qu.:0.0000 1st Qu.:0.0000
#> Median :1.0000 Median :1.0000 Median :1.0 Median :0.0000 Median :0.0000 Median :0.0 Median :0.0000 Median :0.0000
#> Mean :0.7469 Mean :0.9531 Mean :0.9 Mean :0.2875 Mean :0.3688 Mean :0.4 Mean :0.1344 Mean :0.1562
#> 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0 3rd Qu.:0.0000 3rd Qu.:0.0000
#> Max. :1.0000 Max. :1.0000 Max. :1.0 Max. :1.0000 Max. :1.0000 Max. :1.0 Max. :1.0000 Max. :1.0000
#> fecstatus3 firstseen lastseen alive1 alive2 alive3 obsage obslifespan
#> Min. :0.0000 Min. :2004 Min. :2004 Min. :0.0000 Min. :1 Min. :0.0000 Min. :0.000 Min. :0.000
#> 1st Qu.:0.0000 1st Qu.:2004 1st Qu.:2009 1st Qu.:1.0000 1st Qu.:1 1st Qu.:1.0000 1st Qu.:1.000 1st Qu.:5.000
#> Median :0.0000 Median :2004 Median :2009 Median :1.0000 Median :1 Median :1.0000 Median :2.000 Median :5.000
#> Mean :0.2219 Mean :2004 Mean :2009 Mean :0.7688 Mean :1 Mean :0.9469 Mean :1.853 Mean :4.556
#> 3rd Qu.:0.0000 3rd Qu.:2004 3rd Qu.:2009 3rd Qu.:1.0000 3rd Qu.:1 3rd Qu.:1.0000 3rd Qu.:3.000 3rd Qu.:5.000
#> Max. :1.0000 Max. :2008 Max. :2009 Max. :1.0000 Max. :1 Max. :1.0000 Max. :4.000 Max. :5.000
#> stage1 stage2 stage3 matstatus1 matstatus2 matstatus3
#> Length:320 Length:320 Length:320 Min. :1 Min. :1 Min. :1
#> Class :character Class :character Class :character 1st Qu.:1 1st Qu.:1 1st Qu.:1
#> Mode :character Mode :character Mode :character Median :1 Median :1 Median :1
#> Mean :1 Mean :1 Mean :1
#> 3rd Qu.:1 3rd Qu.:1 3rd Qu.:1
#> Max. :1 Max. :1 Max. :1
dim(vertdata)
#> [1] 320 45
In the above code, we described the input dataset in a way that allows R to reorganize it appropriately. For this to work properly, the input dataset needs to be arranged in blocks of columns for each year, with variables in the same order every year. The output dataset includes a number of summary variables, but the data is essentially broken down into groups of three censecutive monitoring occasions each (time t+1, t, and t-1, corresponding to year3, year2, and year1 in the output, respectively), with individuals spread across multiple rows. The output dataset is further limited to those entries in which the individual is alive in time t (year2), meaning that all rows in which an individual is dead or not yet recruited in time t are dropped. Thus, we have 320 rows of data and 45 variables.
This reorganized dataset includes a set of interesting terms, the sizeadded group of three variables. These are sums of the size variables for each time. So, size1added is calculated as sizea1 + sizeb1 + sizec1. This may or may not be sensical depending on the dataset. In this particular dataset, the full size of the individual in each time step is this sum, because Veg gives the number of non-reproductive stems, Inf gives the number of single-flowered inflorescences, and Inf2 gives the number of double-flowered inflorescences per plant per time-step (an inflorescence takes a single stem, and no inflorescence has more than two flowers). So, we will use the sizeadded group of variables to code individual size in our analyses.
Step 2. Provide supplemental information for matrix estimation
The next steps involve the creation of a reproduction matrix and an overwrite data frame. These are optional, and only to be set if the life history of the organism calls for it. The reproduction matrix tells where where fecundity rates need to be set, and at what level. Cypripedium candidum can create seed that germinates by the following growing season (stage P1, or a first year protocorm), or that remains dormant for the next year (stage SD). In the following matrix, we detail that the fecundity of each reproductive stage needs to be split into two between each of these output stages. The actual split places 50% of the fecundity of a stage into each category of recruit, where the full fecundity is estimated by the linear models that we have already created.
Next we create a data frame that outlines transitions that cannot be estimated from the data set and need to be set by other means. For this task, we will use the overwrite function. The function will handle two kinds of given transitions:
transitions that will be set to specific probabilities or rates that we specify, and
transitions that will be set to the values of other transitions that are to be estimated.
Here is an example for the Cypripedium candidum analysis. Note that each row refers to a specific transition, and there are codes for 13 given transitions. Most of these transitions are set to specific probabilities, but four are transitions that will be set to other, estimated transitions (these are the non-NA transitions in eststage parameter set below).
cypover <- overwrite(stage3 = c("SD", "SD", "P1", "P1", "P2", "P3", "D", "V1",
"V2", "V3", "SL", "SL", "SL"), stage2 = c("SD", "SD", "SD",
"SD", "P1", "P2", "P3", "P3", "P3", "P3", "P3", "SL", "SL"),
stage1 = c("SD", "rep", "SD", "rep", "SD", "P1", "P2", "P2",
"P2", "P2", "P2", "P3", "SL"), eststage3 = c(NA, NA, NA, NA,
NA, NA, "D", "V1", "V2", "V3", NA, NA, NA), eststage2 = c(NA,
NA, NA, NA, NA, NA, "D", "D", "D", "D", NA, NA, NA),
eststage1 = c(NA, NA, NA, NA, NA, NA, "D", "D", "D", "D", NA,
NA, NA), givenrate = c(0.1, 0.1, 0.2, 0.2, 0.2, 0.2, NA, NA,
NA, NA, 0.25, 0.4, 0.4), type = c("S", "S", "S", "S", "S",
"S", "S", "S", "S", "S", "S", "S", "S"))
cypover
#> stage3 stage2 stage1 eststage3 eststage2 eststage1 givenrate convtype
#> 1 SD SD SD <NA> <NA> <NA> 0.10 1
#> 2 SD SD rep <NA> <NA> <NA> 0.10 1
#> 3 P1 SD SD <NA> <NA> <NA> 0.20 1
#> 4 P1 SD rep <NA> <NA> <NA> 0.20 1
#> 5 P2 P1 SD <NA> <NA> <NA> 0.20 1
#> 6 P3 P2 P1 <NA> <NA> <NA> 0.20 1
#> 7 D P3 P2 D D D NA 1
#> 8 V1 P3 P2 V1 D D NA 1
#> 9 V2 P3 P2 V2 D D NA 1
#> 10 V3 P3 P2 V3 D D NA 1
#> 11 SL P3 P2 <NA> <NA> <NA> 0.25 1
#> 12 SL SL P3 <NA> <NA> <NA> 0.40 1
#> 13 SL SL SL <NA> <NA> <NA> 0.40 1We can now proceed with matrix estimation!
Step 3. Develop models of demographic parameters
Matrix creation can proceed either as raw matrix creation, as initially outlined in Ehrlén (2000), or via the creation of function-based matrices, in many ways equivalent to complex integral projection models per Ellner and Rees (2006) and as further described in the non-Gaussian case in Shefferson et al. (2014). In the raw matrix case, no vital rate estimation is conducted prior to the construction of a projection matrix. In function-based matrix case, vital rates are first analyzed via linear or non-linear models of the raw demographic dataset, and functions are created that estimate these vital rates according the inputs given. Matrices are then estimated using these functions, as opposed to the raw data.
Package lefko3 estimates up to nine different sets of vital rates:
Survival probability, which is the probability of surviving in time step t to time step t+1, given that the individual is in some stage j. In
lefko3, this parameter may be modeled as a function of size, reproductive status, patch, year, and individual identity. This parameter is required in all function-based matrices.Observation probability, which is an optional parameter denoting the probability of observation in time step t+1 of an individual in stage k given survival from time t to time t+1. This parameter is only used when at least one stage is technically unobservable. For example, some plants are capable of vegetative dormancy, in which case they are alive but do not necessarily sprout in all years. In these cases, the probability of sprouting can be estimated as the observation probability. Note that this probability does NOT refer to observer effort, and so should ONLY be used when some stages are completely unobservable or when size is a count variable and a size of 0 is possible. In
lefko3, this parameter may be modeled as a function of size, reproductive status, patch, year, and individual identity.Size transition probability, which is the probability of becoming size k in time step t+1 assuming survival from time t to time step t+1 and observation in that time. In
lefko3, this parameter may be modeled as a function of size, reproductive status, patch, year, and individual identity. This parameter is required in all function-based size-classified matrices.Reproduction probability, which is an optional parameter denoting the probability of reproducing in time step t+1 given survival from time t to time t+1, and observation in that time. Note that this should be used to separate breeding from non-breeding mature stages. If all adult stages are potentially reproductive, then this parameter is not needed. In
lefko3, this parameter may be modeled as a function of size, reproductive status, patch, year, and individual identity.Fecundity rate, which is the rate of successful production of offspring in time t by individuals alive, observable, and reproductive in that time, and the survival of those offspring into time t+1 in whatever juvenile class is possible. Thus, fecundity rate would be in terms of seedlings produced if seeds are not capable of dormancy. In
lefko3, this parameter may be modeled as a function of size, reproductive status, patch, year, and individual identity. This parameter is required in all function-based matrices.Juvenile survival probability, which is the probability of surviving from juvenile stage j in time step t to a mature stage in time step t+1. In
lefko3, this parameter may be modeled as a function of patch, year, and individual identity.Juvenile observation probability, which is a parameter denoting the probability of observation in time step t+1 of an individual in mature stage k given survival from a juvenile stage in time t to time t+1. In
lefko3, this parameter may be modeled as a function of patch, year, and individual identity, and all other caveats noted in (2) above apply.Juvenile size transition probability, which is the probability of becoming mature size k in time step t+1 assuming survival from juvenile stage j in time t to time step t+1 and observation in that time. In
lefko3, this parameter may be modeled as a function of patch, year, and individual identity.Juvenile reproduction probability, which is a parameter denoting the probability of reproducing in mature stage k in time step t+1 given survival from juvenile stage j in time t to time t+1, and observation in that time. In
lefko3, this parameter may be modeled as a function of patch, year, and individual identity, and all other caveats in (4) apply.
All models require data, and so analyses of datasets that do not include the fates of juveniles should not attempt to estimate juvenile-associated rates and probabilities. Please see chapter 1 for further information.
Prior to vital rate estimation, a number of key decisions need to be made regarding the assumptions underlying the vital rates, and their relationships with the factors under investigation. These decisions include the general modeling strategy, and the size and fecundity distributions.
Step 3a. General modeling strategy
Most function-based matrices, whether integral projection models or otherwise, use either a generalized linear modeling (GLM) strategy, or a generalized linear mixed modeling (GLMM) strategy. The former is more common because of its simplicity, but the latter is theoretically more correct because it allows a correction for pseudoreplication of individuals and a more realistic approach handling time. The difference between the two with regards to vital rate modeling is strongly related to assumptions regarding the individual and the nature of temporal variation in vital rates.
In both cases, the underlying dataset utilized is a vertical dataset. Here, each row of data gives the state of the individual in either two consecutive time steps (the ahistorical case), or three consecutive time steps (the historical case). Under a GLM framework, time step is a fixed categorical variable, and individual identity is ignored. Using time step as a fixed categorical variable implies that the monitoring occasions worked with are the only time steps for which inference is wanted. Thus, it would not be correct to infer future population dynamics after 2009 for a dataset collected between 2004 and 2009, if year is treated as fixed. Ignoring individual identity treats all data points as independent, even though data originating from the same sampled individual is clearly not independent. This leads to pseudoreplication, which is observable in demographic modeling as higher significance of terms, and more significant terms, in modeling.
Under a GLMM framework, both time step and individual identity can be treated as random categorical terms. This has two major implications. First, both time step and individuals can be assumed to be random samples from a broader population of time steps and individuals for which we want to make inferences. Thus, sampled years represent a greater universe of years for which inference can be made, and so their associated coefficients are assumed to come from a normal distribution with mean = 0. Second, treating individual as a random categorical term eliminates the pseudoreplication that is inherent in the GLM approach to vital rate estimation when individuals are monitored potentially many times. We encourage researchers to use the GLMM approach in their demographic work, but we have also included easy-to-use GLM functionality, since many will find the GLM approach particularly useful in cases where mixed modeling breaks down.
Step 3b. Size and fecundity distributions
Once a general approach is decided upon, the next step is to decide the underlying distributions. The probabilities of survival, observation, and reproductive status are automatically set to the binomial distribution, and this cannot be altered. However, the probability of size transition and fecundity rate can be set to the Gaussian, Poisson, or negative binomial distributions. In general, if size or fecundity rate is a continuous variable (i.e., not an integer or count variable), then it should be set to the Gaussian distribution. In contrast, if size or fecundity rate is a count, then it should be set to the Poisson distribution. The negative binomial distribution is also provided in cases where the Poisson distribution’s assumptions, such as the mean equaling the variance, are clearly broken. We do not encourage the use of the negative binomial except in such extreme cases, as the extra parameters estimated for the negative binomial distribution reduce the power of the modeling exercises conducted.
Neither the Poisson nor the negative binomial distribution handle 0s. We do not include zero-inflated versions of these distributions because such distributions are technically incorrect, or theoretically sloppy. If 0s occur within the dataset, then the best option is to utilize the probability of observation or the probability of reproduction as vital rates that will also be modeled. This will allow the 0s to be analyzed within these binomial models, and leave the non-zero values for rate estimation under the Poisson or negative binomial distributions. If 0s still exist after doing this, then it is likely that some individuals assigned to visible classes are likely not observable, size class might be determined logarithmically, or reproductive classes are not really reproductive. Further exploration of the dataset should yield the correct answer.
Step 3c. Model building and selection
In lefko3, the modelsearch function is the main workhorse function that conducts vital rate estimation. Here, we will create a full suite of vital rate models for the Cypripedium candidum dataset. Before proceeding, we need to decide on the linear model building strategy, the correct vital rates to model, the proper statistical distributions for estimated vital rates, the proper parameterizations for each vital rate, and the strategy for determination of the best-fit models.
First, the building strategy. In most cases, the best procedure will be through linear mixed models in which monitoring occasion and individual identity are random terms. We set monitoring occasion as random because we wish to make inferences for the population as a whole and not restrict ourselves to inference only for the years monitored (i.e. our distribution of years sampled is itself a sample of the population in time). We set individual identity to random because we do not wish to pseudoreplicate, since not doing so will result in the assumption that each row in our vertical dataset is statistically independent of all other rows. Thus, we will use approach = "lme4" in the parameterization for modelsearch, and keep the defaults for year.as.random, indiv, and year, which are set to the default output for whether monitoring occasion is a random or fixed term (random by default), which variable corresponds to individual identity (individ by default), and which variable corresponds to time t(year2 by default).
The mixed modeling approach is usually better, particularly because of its elimination of pseudoreplication. However, a mixed modeling strategy results in lower statistical power and a greater time used in estimating models. lefko3 users wishing to use a standard generalized linear modeling strategy can set approach = "glm". In this case, individual identity is not used and all observed transitions are treated as independent.
Second, the correct vital rates to model. The function modelsearch() estimates up to 5 sets of linear models:
survival probability from time t to time t+1 (juvenile and mature),
observation probability in time t+1 assuming survival until that time (juvenile and mature),
size in time t+1 assuming survival and observation in that time (juvenile and mature),
reproduction status in time t+1 assuming survival and observation until that time (individuals leaving the juvenile stages, and those already mature), and
fecundity rate assuming survival in time t and observation and reproduction in time t+1 (mature only).
The default settings for modelsearch involve the estimation of 1), survival probability, 3) size distribution, and 5) fecundity, which are the minimum required for a full projection matrix. Observation probability (option obs in vitalrates) should only be included when a life history stage or size exists that cannot be observed. For example, in the case of a plant with vegetative dormancy, the observation probability can be thought of as the sprouting probability, which is a biologically meaningful vital rate (Shefferson et al. 2001). Further, reproduction status (option repst in vitalrates) should only be modeled if size classification should be stratified by the ability to reproduce, as when 0 fecundity occurs within the dataset. In this latter case, we can imagine that reproductive and non-reproductive individuals of each size class might theoretically exist, and we wish to parameterize transitions allowing individuals to be reproductive or non-reproductive. Since Cypripedium candidum is capable of longs bouts of vegetative dormancy, and since we wish to stratify the population into reproductive and non-reproductive adults, we will set vitalrates = c("surv", "obs", "size", "repst", "fec").
Third, we need to set the proper statistical distributions for each parameter. Survival probability, observation probability, and reproductive status (a probability) are all modeled as binomial variables, and this cannot be changed. In the case of Cypripedium, size was measured as the number of stems and so is a count variable. Likewise, fecundity is actually estimated as the number of fruits produced per plant, and so is also a count variable. Thus, we will use the Poisson distribution for both variables.
Users of lefko3 who wish to use a Poisson or negative binomial distribution should note that 0s are not allowed, because zero-inflated models are essentially conditional models with different statistical properties at 0 than elsewhere. The modelsearch function will subset the dataset used to parameterize size and fecundity to only rows without zeroes in time t+1. If size is potentially 0, then consider using observation probability to catch the probability of becoming size 0 vs. observable sizes.
Fourth, we need the proper model parameterizations for each vital rate, using the suite option. The default, suite = "main", under the mixed model setting (approach = "lme4") is for modelsearch to estimate a global model that includes sizes in times t and t-1, and reproductive status in times t and t-1, as fixed factors, with individual identity and time t set as random terms in a mixed model framework using R package lme4 (Bates et al. 2015). However, setting suite = "full" will yield a global model that also includes all two-way interactions. We will set to the latter. The default under the GLM setting (approach = "glm") makes time t a fixed term and drops individual identity from consideration. The global model under suite = "full" then includes all fixed factors noted before, plus time t and all two-way interactions with it. To eliminate all interactions from the model and only analyze main effects, use suite = "main". If the population is not stratified by reproductive status, then suite = "size" will eliminate reproductive status terms and use all others in the global model. If size is not important, then suite = "rep" will eiminate size but keep reproductive status and all other terms. Finally, suite = "cons" will result in a global model in which both reproductive status and size are not considered.
Fifth, and finally, we need to determine the proper strategy for determination of the best-fit model. Model building proceeds through the dredge function in package MuMIn (Bartoń 2014). Each model has a resulting AICc value. The default setting in lefko3 (bestfit = "AICc&k") will compare all models within 2.0 AICc units of the model with &DeltaAICc=0, and choose the one with the lowest degrees of freedom. This approach is generally better than the alternative, which simply uses the model with &DeltaAICc=0 (bestfit = "AICc"), as all models within 2.0 AICc units of that model are equally parsimonious and so fewer degrees of freedom result from fewer parameters estimated (Burnham and Anderson 2002).
In the model building exercise below, we will use the suite = "full" option to run all main effects and their two-way interactions.
cypmodels3 <- modelsearch(vertdata, historical = TRUE, approach = "lme4",
vitalrates = c("surv", "obs", "size", "repst", "fec"),
sizedist = "poisson", fecdist = "poisson", suite = "full",
size = c("size3added", "size2added", "size1added"),
quiet = TRUE)
#> Warning in modelsearch(vertdata, historical = TRUE, approach = "lme4", vitalrates = c("surv", : WARNING: Fecundity in time t cannot be Poisson-distributed and include 0s. Will develop fecundity models excluding all 0s. Consider adding a reproductive status variable to absorb 0 values.
#>
#> Developing global model of survival probability...
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge: degenerate Hessian with 2
#> negative eigenvalues
#>
#> Developing global model of observation probability...
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0217509 (tol =
#> 0.002, component 1)
#>
#> Developing global model of size (Poisson)...
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0181958 (tol =
#> 0.002, component 1)
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
#> - Rescale variables?
#>
#> Developing global model of the probability of reproduction...
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00594956 (tol =
#> 0.002, component 1)
#>
#> Developing global model of fecundity (Poisson)...
#> fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients
#> boundary (singular) fit: see ?isSingular
#>
#> All global models developed.
#>
#>
#> Commencing dredge of survival probability...
#>
#>
#> Commencing dredge of observation probability...
#>
#>
#> Commencing dredge of size...
#>
#>
#> Commencing dredge of reproduction probability...
#>
#>
#> Commencing dredge of fecundity...
#>
#>
#> Finished dredging all vital rates.
As modelsearch works, it produces a good amount of text to allow the user to understand what is going on. It is entirely possible, and actually quite likely, that many warning messages will appear, and these may be of use to users in understanding their data and how well it conforms to their analytical assumptions. We have silenced this with the quiet = TRUE option, but we encourage users to allow the function to run unsilenced, in case of large problems in modeling.
Once done, we can summarize the output with the summary() function.
summary(cypmodels3)
#> This LefkoMod object includes 5 linear models.
#> Best-fit model criterion used: AICc&k
#>
#> ────────────────────────────────────────────────────────────────────────────────
#> Survival model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#> Family: binomial ( logit )
#> Formula: alive3 ~ size2added + (1 | individ) + (1 | year2)
#> Data: surv.data
#> AIC BIC logLik deviance df.resid
#> 128.1324 143.2057 -60.0662 120.1324 316
#> Random effects:
#> Groups Name Std.Dev.
#> individ (Intercept) 1.198378
#> year2 (Intercept) 0.008826
#> Number of obs: 320, groups: individ, 74; year2, 5
#> Fixed Effects:
#> (Intercept) size2added
#> 2.0352 0.6344
#> convergence code 0; 0 optimizer warnings; 1 lme4 warnings
#>
#> ────────────────────────────────────────
#>
#> Observation model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#> Family: binomial ( logit )
#> Formula: obsstatus3 ~ size2added + (1 | individ) + (1 | year2)
#> Data: obs.data
#> AIC BIC logLik deviance df.resid
#> 118.2567 133.1117 -55.1284 110.2567 299
#> Random effects:
#> Groups Name Std.Dev.
#> individ (Intercept) 1.078e-05
#> year2 (Intercept) 8.776e-01
#> Number of obs: 303, groups: individ, 70; year2, 5
#> Fixed Effects:
#> (Intercept) size2added
#> 2.4904 0.3134
#> convergence code 0; 0 optimizer warnings; 1 lme4 warnings
#>
#> ────────────────────────────────────────
#>
#> Size model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#> Family: poisson ( log )
#> Formula: size3added ~ repstatus1 + repstatus2 + size1added + size2added +
#> (1 | individ) + (1 | year2) + repstatus1:repstatus2 + size1added:size2added
#> Data: size.data
#> AIC BIC logLik deviance df.resid
#> 1081.6019 1114.5685 -531.8009 1063.6019 279
#> Random effects:
#> Groups Name Std.Dev.
#> individ (Intercept) 0.06104
#> year2 (Intercept) 0.19627
#> Number of obs: 288, groups: individ, 70; year2, 5
#> Fixed Effects:
#> (Intercept) repstatus1 repstatus2 size1added size2added repstatus1:repstatus2
#> 0.213173 0.278913 0.255861 0.131367 0.153437 -0.368510
#> size1added:size2added
#> -0.008696
#> convergence code 0; 0 optimizer warnings; 2 lme4 warnings
#>
#> ────────────────────────────────────────
#>
#> Reproductive status model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#> Family: binomial ( logit )
#> Formula: repstatus3 ~ repstatus2 + size2added + (1 | individ) + (1 | year2)
#> Data: repst.data
#> AIC BIC logLik deviance df.resid
#> 333.6176 351.9324 -161.8088 323.6176 283
#> Random effects:
#> Groups Name Std.Dev.
#> individ (Intercept) 0.1829
#> year2 (Intercept) 0.6250
#> Number of obs: 288, groups: individ, 70; year2, 5
#> Fixed Effects:
#> (Intercept) repstatus2 size2added
#> -1.4630 1.6457 0.1715
#>
#> ────────────────────────────────────────
#>
#> Fecundity model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#> Family: poisson ( log )
#> Formula: feca2 ~ size2added + (1 | individ) + (1 | year2)
#> Data: fec.data
#> AIC BIC logLik deviance df.resid
#> 156.8531 164.5012 -74.4266 148.8531 46
#> Random effects:
#> Groups Name Std.Dev.
#> individ (Intercept) 0.000
#> year2 (Intercept) 0.172
#> Number of obs: 50, groups: individ, 26; year2, 5
#> Fixed Effects:
#> (Intercept) size2added
#> 0.21931 0.04429
#> convergence code 0; 0 optimizer warnings; 1 lme4 warnings
#>
#> ────────────────────────────────────────────────────────────────────────────────
#> Juvenile survival model:
#> [1] 1
#>
#> ────────────────────────────────────────
#>
#> Juvenile observation model:
#> [1] 1
#>
#> ────────────────────────────────────────
#>
#> Juvenile size model:
#> [1] 1
#>
#> ────────────────────────────────────────
#>
#> Juvenile reproduction model:
#> [1] 1
#>
#>
#> ────────────────────────────────────────────────────────────────────────────────
#>
#> Number of models in survival table:113
#>
#> Number of models in observation table:113
#>
#> Number of models in size table:113
#>
#> Number of models in reproduction status table:113
#>
#> Number of models in fecundity table:113
#>
#> Number of models in juvenile survival table: 1
#>
#> Number of models in juvenile observation table: 1
#>
#> Number of models in juvenile size table: 1
#>
#> Number of models in juvenile reproduction table: 1
#>
#>
#> ────────────────────────────────────────────────────────────────────────────────
#>
#> General model parameter names (column 1), and specific names used in these models (column 2):
#> mainparams modelparams
#> 1 year2 year2
#> 2 individ individ
#> 3 patch <NA>
#> 4 surv3 alive3
#> 5 obs3 obsstatus3
#> 6 size3 size3added
#> 7 repst3 repstatus3
#> 8 fec3 feca3
#> 9 fec2 feca2
#> 10 size2 size2added
#> 11 size1 size1added
#> 12 repst2 repstatus2
#> 13 repst1 repstatus1
#> 14 age <NA>
#>
#>
#> ────────────────────────────────────────────────────────────────────────────────
#>
#> Quality control:
#>
#> Survival estimated with 74 individuals and 320 individual transitions.
#> Observation estimated with 70 individuals and 303 individual transitions.
#> Size estimated with 70 individuals and 288 individual transitions.
#> Reproductive status estimated with 70 individuals and 288 individual transitions.
#> Fecundity estimated with 26 individuals and 50 individual transitions.
#> Juvenile survival not estimated.
#> Juvenile observation probability not estimated.
#> Juvenile size transition not estimated.
#> Juvenile reproduction probability not estimated.
#> NULL
The summary function gives us a great deal of information about all of the models, but also hides a number of more technical details. For example, the complete model selection tables are actually included in the modelsearch output, and users can see these by calling them directly, as below.
cypmodels3$survival_table
#> Global model call: lme4::glmer(formula = alive3 ~ size1added + size2added + repstatus1 +
#> repstatus2 + size1added:size2added + repstatus1:repstatus2 +
#> size1added:repstatus1 + size2added:repstatus2 + size1added:repstatus2 +
#> size2added:repstatus1 + (1 | individ) + (1 | year2), data = surv.data,
#> family = "binomial")
#> ---
#> Model selection table
#> (Int) rp1 rp2 sz1 sz2 rp1:rp2 rp1:sz1 rp1:sz2 rp2:sz1 rp2:sz2 sz1:sz2 df logLik used.criterion delta weight
#> 9 2.035 0.6344 4 -60.066 128.3 0.00 0.134
#> 11 2.006 0.54530 0.5899 5 -59.773 129.7 1.48 0.064
#> 13 1.538 0.18620 0.5185 5 -59.781 129.8 1.49 0.063
#> 10 1.984 0.31570 0.6240 5 -59.967 130.1 1.87 0.053
#> 28 2.063 -0.42760 -0.13110 0.5712 46.86 7 -58.039 130.4 2.18 0.045
#> 74 2.149 -0.77080 0.5004 0.81440 6 -59.407 131.1 2.82 0.033
#> 15 1.505 0.41530 0.18280 0.4824 6 -59.578 131.4 3.17 0.027
#> 46 2.042 -1.70300 0.03281 0.5542 1.0750 7 -58.546 131.5 3.19 0.027
#> 12 1.972 0.25460 0.51650 0.5866 6 -59.712 131.7 3.43 0.024
#> 525 1.453 0.23930 0.5678 -0.0242500 6 -59.745 131.8 3.50 0.023
#> 267 1.994 0.78800 0.6309 -1.287e-01 6 -59.748 131.8 3.50 0.023
#> 14 1.535 0.04577 0.18200 0.5183 6 -59.779 131.8 3.57 0.022
#> 64 2.193 -2.84700 -0.13040 0.03241 0.5207 177.70 1.2890 9 -56.623 131.8 3.57 0.022
#> 92 2.161 -0.97690 -0.09133 0.5073 66.09 0.43230 8 -57.864 132.2 3.93 0.019
#> 32 1.642 -0.62310 -0.13560 0.16720 0.4722 32.09 8 -57.910 132.3 4.02 0.018
#> 284 1.999 -0.42880 0.16980 0.6194 78.38 -1.665e-01 8 -57.991 132.4 4.18 0.017
#> 160 1.446 -0.80970 0.49640 0.35650 0.4784 29.16 -0.5133 9 -57.022 132.6 4.37 0.015
#> 76 2.126 -0.74500 0.43350 0.4772 0.75490 7 -59.228 132.8 4.56 0.014
#> 143 1.381 0.81620 0.27590 0.4952 -0.2702 7 -59.333 133.0 4.76 0.012
#> 78 1.696 -0.78060 0.15700 0.4225 0.66390 7 -59.333 133.0 4.77 0.012
#> 48 2.022 -1.76700 0.47620 0.03453 0.5164 1.0790 8 -58.336 133.1 4.88 0.012
#> 110 2.112 -2.41400 0.04005 0.5037 1.0360 0.64910 8 -58.346 133.2 4.89 0.012
#> 192 1.900 -2.52000 0.32410 0.18860 0.5129 50.38 1.0710 -0.3320 10 -56.289 133.3 5.03 0.011
#> 527 1.425 0.41160 0.23320 0.5288 -0.0227300 7 -59.546 133.5 5.19 0.010
#> 271 1.467 0.58300 0.18440 0.5106 -9.563e-02 7 -59.561 133.5 5.22 0.010
#> 16 1.507 -0.05451 0.42620 0.18760 0.4815 7 -59.575 133.5 5.25 0.010
#> 558 1.937 -1.68200 0.08943 0.6026 1.0560 -0.0230700 8 -58.524 133.5 5.25 0.010
#> 268 1.950 0.27090 0.78600 0.6328 -1.436e-01 7 -59.680 133.7 5.46 0.009
#> 526 1.451 0.04333 0.23500 0.5673 -0.0240900 7 -59.743 133.8 5.59 0.008
#> 96 1.988 -1.04200 -0.09154 0.13020 0.4754 67.44 0.37210 9 -57.649 133.9 5.62 0.008
#> 576 2.101 -2.79200 -0.13630 0.07267 0.5484 26.46 1.2560 -0.0149300 10 -56.614 133.9 5.68 0.008
#> 128 2.161 -2.75000 -0.13640 0.03317 0.5258 42.13 1.3190 -0.12630 10 -56.618 133.9 5.69 0.008
#> 320 2.156 -2.81200 -0.07952 0.03480 0.5263 16.11 1.2750 -2.692e-02 10 -56.623 134.0 5.70 0.008
#> 544 1.524 -0.63140 -0.14680 0.24430 0.5409 32.12 -0.0345200 9 -57.852 134.3 6.02 0.007
#> 348 2.123 -0.93650 0.04906 0.5338 48.18 0.40450 -7.679e-02 9 -57.855 134.3 6.03 0.007
#> 288 1.596 -0.62030 0.08667 0.16430 0.5111 102.40 -1.275e-01 9 -57.877 134.3 6.08 0.006
#> 80 1.951 -0.83530 0.42180 0.13510 0.4415 0.70600 8 -59.007 134.5 6.22 0.006
#> 5 2.203 0.42320 4 -63.249 134.6 6.36 0.006
#> 224 1.483 -0.98140 0.49340 0.34520 0.4569 37.83 0.16060 -0.4970 10 -56.997 134.7 6.45 0.005
#> 672 1.444 -0.80970 0.49510 0.35750 0.4797 33.19 -0.5125 -0.0005682 10 -57.022 134.8 6.50 0.005
#> 416 1.445 -0.80970 0.49780 0.35650 0.4788 30.05 -0.5131 -9.738e-04 10 -57.022 134.8 6.50 0.005
#> 332 2.060 -0.73950 0.75970 0.5324 0.75250 -1.818e-01 8 -59.169 134.8 6.54 0.005
#> 176 1.523 -1.75600 0.97530 0.20890 0.4554 1.0470 -0.4025 9 -58.111 134.8 6.54 0.005
#> 590 1.546 -0.84910 0.26430 0.5043 0.69230 -0.0454200 8 -59.216 134.9 6.64 0.005
#> 56 2.765 -3.14600 0.26410 0.22520 20.36 1.3760 8 -59.221 134.9 6.65 0.005
#> 112 2.082 -2.34100 0.40920 0.04015 0.4795 1.0320 0.55050 9 -58.191 135.0 6.70 0.005
#> 7 2.034 0.80290 0.40000 5 -62.393 135.0 6.72 0.005
#> 655 1.362 0.80260 0.28700 0.5080 -0.2618 -0.0063460 8 -59.330 135.1 6.86 0.004
#> 144 1.382 -0.02074 0.82000 0.27750 0.4949 -0.2698 8 -59.332 135.1 6.87 0.004
#> 399 1.376 0.83470 0.27570 0.4990 -0.2679 -1.284e-02 8 -59.332 135.1 6.87 0.004
#> 304 1.959 -1.75900 0.76360 0.03397 0.5650 1.0790 -1.612e-01 9 -58.292 135.2 6.90 0.004
#> 560 1.938 -1.75500 0.46630 0.07946 0.5562 1.0650 -0.0183300 9 -58.319 135.2 6.96 0.004
#> 622 2.018 -2.36400 0.09076 0.5452 1.0040 0.63730 -0.0198700 9 -58.327 135.2 6.97 0.004
#> 448 1.937 -2.55200 0.20510 0.18990 0.4885 21.69 1.0830 -0.3451 7.605e-02 11 -56.280 135.4 7.16 0.004
#> 256 1.890 -2.44700 0.31900 0.18750 0.5216 30.99 1.1140 -0.13040 -0.3329 11 -56.282 135.4 7.16 0.004
#> 704 1.921 -2.52600 0.33430 0.17940 0.5018 74.15 1.0740 -0.3394 0.0051310 11 -56.288 135.4 7.17 0.004
#> 783 1.401 0.55540 0.22940 0.5477 -8.203e-02 -0.0200800 8 -59.533 135.5 7.27 0.004
#> 528 1.427 -0.05600 0.42280 0.23850 0.5282 -0.0228700 8 -59.543 135.5 7.29 0.003
#> 272 1.470 -0.04976 0.59090 0.18880 0.5094 -9.440e-02 8 -59.558 135.6 7.32 0.003
#> 152 1.913 -0.93410 0.95620 0.65090 20.14 -0.6027 8 -59.560 135.6 7.32 0.003
#> 24 2.154 -0.69120 0.25110 0.39910 18.80 7 -60.686 135.7 7.47 0.003
#> 608 1.834 -1.11300 -0.09708 0.24060 0.5457 23.46 0.40290 -0.0405500 10 -57.568 135.8 7.59 0.003
#> [ reached getOption("max.print") -- omitted 51 rows ]
#> Models ranked by used.criterion(x)
#> Random terms (all models):
#> '1 | individ', '1 | year2'Before moving on, we note that the models created above are actually only usable for the construction of historical models. For comparison, we may wish to estimate ahistorical models. In that case, we also need linear models in which the global models tested do not include state at time t-1. Here, we produce these models.
cypmodels2 <- modelsearch(vertdata, historical = FALSE, approach = "lme4",
vitalrates = c("surv", "obs", "size", "repst", "fec"),
sizedist = "poisson", fecdist = "poisson", suite = "full",
size = c("size3added", "size2added"))
#> Warning in modelsearch(vertdata, historical = FALSE, approach = "lme4", : WARNING: Fecundity in time t cannot be Poisson-distributed and include 0s. Will develop fecundity models excluding all 0s. Consider adding a reproductive status variable to absorb 0 values.
#>
#> Developing global model of survival probability...
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0330463 (tol =
#> 0.002, component 1)
#>
#> Developing global model of observation probability...
#>
#>
#> Developing global model of size (Poisson)...
#>
#>
#> Developing global model of the probability of reproduction...
#>
#>
#> Developing global model of fecundity (Poisson)...
#> fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
#> boundary (singular) fit: see ?isSingular
#>
#> All global models developed.
#>
#>
#> Commencing dredge of survival probability...
#> Fixed term is "(Intercept)"
#> boundary (singular) fit: see ?isSingular
#> boundary (singular) fit: see ?isSingular
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0268727 (tol =
#> 0.002, component 1)
#> boundary (singular) fit: see ?isSingular
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge: degenerate Hessian with 1
#> negative eigenvalues
#> Warning in vcov.merMod(object, use.hessian = use.hessian): variance-covariance matrix computed from finite-difference Hessian is
#> not positive definite or contains NA values: falling back to var-cov estimated from RX
#> Warning in vcov.merMod(object, correlation = correlation, sigm = sig): variance-covariance matrix computed from finite-difference Hessian is
#> not positive definite or contains NA values: falling back to var-cov estimated from RX
#>
#> Commencing dredge of observation probability...
#> Fixed term is "(Intercept)"
#> boundary (singular) fit: see ?isSingular
#>
#> Commencing dredge of size...
#> Fixed term is "(Intercept)"
#>
#> Commencing dredge of reproduction probability...
#> Fixed term is "(Intercept)"
#>
#> Commencing dredge of fecundity...
#> Fixed term is "(Intercept)"
#> boundary (singular) fit: see ?isSingular
#> fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
#> boundary (singular) fit: see ?isSingular
#> boundary (singular) fit: see ?isSingular
#> fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
#> boundary (singular) fit: see ?isSingular
#> fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
#> boundary (singular) fit: see ?isSingular
#>
#> Finished dredging all vital rates.
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0268727 (tol =
#> 0.002, component 1)
#> boundary (singular) fit: see ?isSingular
#> boundary (singular) fit: see ?isSingular
summary(cypmodels2)
#> This LefkoMod object includes 5 linear models.
#> Best-fit model criterion used: AICc&k
#>
#> ────────────────────────────────────────────────────────────────────────────────
#> Survival model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#> Family: binomial ( logit )
#> Formula: alive3 ~ size2added + (1 | individ) + (1 | year2)
#> Data: surv.data
#> AIC BIC logLik deviance df.resid
#> 128.1324 143.2057 -60.0662 120.1324 316
#> Random effects:
#> Groups Name Std.Dev.
#> individ (Intercept) 1.198378
#> year2 (Intercept) 0.008826
#> Number of obs: 320, groups: individ, 74; year2, 5
#> Fixed Effects:
#> (Intercept) size2added
#> 2.0352 0.6344
#> convergence code 0; 0 optimizer warnings; 1 lme4 warnings
#>
#> ────────────────────────────────────────
#>
#> Observation model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#> Family: binomial ( logit )
#> Formula: obsstatus3 ~ size2added + (1 | individ) + (1 | year2)
#> Data: obs.data
#> AIC BIC logLik deviance df.resid
#> 118.2567 133.1117 -55.1284 110.2567 299
#> Random effects:
#> Groups Name Std.Dev.
#> individ (Intercept) 1.078e-05
#> year2 (Intercept) 8.776e-01
#> Number of obs: 303, groups: individ, 70; year2, 5
#> Fixed Effects:
#> (Intercept) size2added
#> 2.4904 0.3134
#> convergence code 0; 0 optimizer warnings; 1 lme4 warnings
#>
#> ────────────────────────────────────────
#>
#> Size model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#> Family: poisson ( log )
#> Formula: size3added ~ size2added + (1 | individ) + (1 | year2)
#> Data: size.data
#> AIC BIC logLik deviance df.resid
#> 1115.3617 1130.0135 -553.6808 1107.3617 284
#> Random effects:
#> Groups Name Std.Dev.
#> individ (Intercept) 0.55311
#> year2 (Intercept) 0.09895
#> Number of obs: 288, groups: individ, 70; year2, 5
#> Fixed Effects:
#> (Intercept) size2added
#> 0.83893 0.04244
#>
#> ────────────────────────────────────────
#>
#> Reproductive status model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#> Family: binomial ( logit )
#> Formula: repstatus3 ~ repstatus2 + size2added + (1 | individ) + (1 | year2)
#> Data: repst.data
#> AIC BIC logLik deviance df.resid
#> 333.6176 351.9324 -161.8088 323.6176 283
#> Random effects:
#> Groups Name Std.Dev.
#> individ (Intercept) 0.1829
#> year2 (Intercept) 0.6250
#> Number of obs: 288, groups: individ, 70; year2, 5
#> Fixed Effects:
#> (Intercept) repstatus2 size2added
#> -1.4630 1.6457 0.1715
#>
#> ────────────────────────────────────────
#>
#> Fecundity model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#> Family: poisson ( log )
#> Formula: feca2 ~ size2added + (1 | individ) + (1 | year2)
#> Data: fec.data
#> AIC BIC logLik deviance df.resid
#> 156.8531 164.5012 -74.4266 148.8531 46
#> Random effects:
#> Groups Name Std.Dev.
#> individ (Intercept) 0.000
#> year2 (Intercept) 0.172
#> Number of obs: 50, groups: individ, 26; year2, 5
#> Fixed Effects:
#> (Intercept) size2added
#> 0.21931 0.04429
#> convergence code 0; 0 optimizer warnings; 1 lme4 warnings
#>
#> ────────────────────────────────────────────────────────────────────────────────
#> Juvenile survival model:
#> [1] 1
#>
#> ────────────────────────────────────────
#>
#> Juvenile observation model:
#> [1] 1
#>
#> ────────────────────────────────────────
#>
#> Juvenile size model:
#> [1] 1
#>
#> ────────────────────────────────────────
#>
#> Juvenile reproduction model:
#> [1] 1
#>
#>
#> ────────────────────────────────────────────────────────────────────────────────
#>
#> Number of models in survival table:5
#>
#> Number of models in observation table:5
#>
#> Number of models in size table:5
#>
#> Number of models in reproduction status table:5
#>
#> Number of models in fecundity table:5
#>
#> Number of models in juvenile survival table: 1
#>
#> Number of models in juvenile observation table: 1
#>
#> Number of models in juvenile size table: 1
#>
#> Number of models in juvenile reproduction table: 1
#>
#>
#> ────────────────────────────────────────────────────────────────────────────────
#>
#> General model parameter names (column 1), and specific names used in these models (column 2):
#> mainparams modelparams
#> 1 year2 year2
#> 2 individ individ
#> 3 patch <NA>
#> 4 surv3 alive3
#> 5 obs3 obsstatus3
#> 6 size3 size3added
#> 7 repst3 repstatus3
#> 8 fec3 feca3
#> 9 fec2 feca2
#> 10 size2 size2added
#> 11 size1 <NA>
#> 12 repst2 repstatus2
#> 13 repst1 <NA>
#> 14 age <NA>
#>
#>
#> ────────────────────────────────────────────────────────────────────────────────
#>
#> Quality control:
#>
#> Survival estimated with 74 individuals and 320 individual transitions.
#> Observation estimated with 70 individuals and 303 individual transitions.
#> Size estimated with 70 individuals and 288 individual transitions.
#> Reproductive status estimated with 70 individuals and 288 individual transitions.
#> Fecundity estimated with 26 individuals and 50 individual transitions.
#> Juvenile survival not estimated.
#> Juvenile observation probability not estimated.
#> Juvenile size transition not estimated.
#> Juvenile reproduction probability not estimated.
#> NULLFewer models were estimated per dredge, since fewer parameters were tested in the global models (size and reproductive status in time t-1 were not included). So, the best-fit models should look a little bit different. However, a more thorough comparison will show that many of the best-fit models are similar between historical and ahistorical analysis - indeed, size appears to be the only major vital rate to differ. This is not guaranteed - in this case, it may be that the relatively small number of years and small overall sample size leaves too little power to find an impact of historical status on most vital rates.
Step 4. Projection matrix creation
Now we will create population function-based matrices for each year of the study. Function-based matrices have quickly been taking over in population ecology, probably because of their ability to parse out interesting trends and influential factors picked up by the linear modeling of vital rates. Let’s first create a set of ahistorical matrices.
cypmatrix2 <- flefko2(stageframe = cypframe, repmatrix = rep.assumptions,
modelsuite = cypmodels2, overwrite = cypover,
data = vertdata, year.as.random = TRUE)
summary(cypmatrix2)
#>
#> This lefkoMat object contains 5 matrices.
#>
#> Each matrix is a square matrix with 54 rows and columns, and a total of 2916 elements.
#> A total of 12035 survival-transitions were estimated, with 2407 per matrix.
#> A total of 240 fecundity transitions were estimated, with 48 per matrix.
#>
#> Vital rate modeling quality control:
#>
#> Survival estimated with 74 individuals and 320 individual transitions.
#> Observation estimated with 70 individuals and 303 individual transitions.
#> Size estimated with 70 individuals and 288 individual transitions.
#> Reproductive status estimated with 70 individuals and 288 individual transitions.
#> Fecundity estimated with 26 individuals and 50 individual transitions.
#> Juvenile survival estimated with 0 individuals and 0 individual transitions.
#> Juvenile observation estimated with 0 individuals and 0 individual transitions.
#> Juvenile size estimated with 0 individuals and 0 individual transitions.
#> Juvenile reproductive status estimated with 0 individuals and 0 individual transitions.
#> NULLA quick glance at the summary output will highlight that many more elements are estimated for function-based matrices. In raw matrices, elements associated with transitions from specific stages are only estimated when individuals actually exist within those particular stages. In function-based matrices, in contrast, the linear models estimated allow the estimation of all elements that are theoretically possible (i.e. only structural 0s are not estimated). Let’s take a look at an example matrix, but only on the top corner to deal with its size.
cypmatrix2$A[[1]][1:25, 1:8]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] 0.1 0.0 0.0 0.00 0.0 0.000000e+00 0.000000e+00 0.000000e+00
#> [2,] 0.2 0.0 0.0 0.00 0.0 0.000000e+00 0.000000e+00 0.000000e+00
#> [3,] 0.0 0.2 0.0 0.00 0.0 0.000000e+00 0.000000e+00 0.000000e+00
#> [4,] 0.0 0.0 0.2 0.00 0.0 0.000000e+00 0.000000e+00 0.000000e+00
#> [5,] 0.0 0.0 0.0 0.25 0.4 0.000000e+00 0.000000e+00 0.000000e+00
#> [6,] 0.0 0.0 0.0 0.00 0.0 4.701977e-02 3.686912e-02 2.809353e-02
#> [7,] 0.0 0.0 0.0 0.00 0.0 1.005860e-01 9.518443e-02 8.657188e-02
#> [8,] 0.0 0.0 0.0 0.00 0.0 1.491216e-01 1.472310e-01 1.397142e-01
#> [9,] 0.0 0.0 0.0 0.00 0.0 1.473847e-01 1.518244e-01 1.503187e-01
#> [10,] 0.0 0.0 0.0 0.00 0.0 1.092510e-01 1.174208e-01 1.212961e-01
#> [11,] 0.0 0.0 0.0 0.00 0.0 6.478709e-02 7.265046e-02 7.830161e-02
#> [12,] 0.0 0.0 0.0 0.00 0.0 3.201624e-02 3.745852e-02 4.212241e-02
#> [13,] 0.0 0.0 0.0 0.00 0.0 1.356143e-02 1.655450e-02 1.942267e-02
#> [14,] 0.0 0.0 0.0 0.00 0.0 5.026301e-03 6.401613e-03 7.836331e-03
#> [15,] 0.0 0.0 0.0 0.00 0.0 1.655919e-03 2.200444e-03 2.810373e-03
#> [16,] 0.0 0.0 0.0 0.00 0.0 4.909894e-04 6.807282e-04 9.071055e-04
#> [17,] 0.0 0.0 0.0 0.00 0.0 1.323465e-04 1.914452e-04 2.661699e-04
#> [18,] 0.0 0.0 0.0 0.00 0.0 3.270125e-05 4.935449e-05 7.159314e-05
#> [19,] 0.0 0.0 0.0 0.00 0.0 7.458545e-06 1.174483e-05 1.777550e-05
#> [20,] 0.0 0.0 0.0 0.00 0.0 1.579644e-06 2.595268e-06 4.098148e-06
#> [21,] 0.0 0.0 0.0 0.00 0.0 3.122490e-07 5.352471e-07 8.818407e-07
#> [22,] 0.0 0.0 0.0 0.00 0.0 5.786475e-08 1.034899e-07 1.778951e-07
#> [23,] 0.0 0.0 0.0 0.00 0.0 1.009249e-08 1.883269e-08 3.377604e-08
#> [24,] 0.0 0.0 0.0 0.00 0.0 1.662489e-09 3.236706e-09 6.056616e-09
#> [25,] 0.0 0.0 0.0 0.00 0.0 2.594408e-10 5.270030e-10 1.028893e-09The matrix is overwhelmingly composed of non-zero elements, unlike in the raw matrix case.
Next, we will create a set of historical Lefkovitch matrices.
cypmatrix3 <- flefko3(stageframe = cypframe, repmatrix = rep.assumptions,
modelsuite = cypmodels3, overwrite = cypover,
data = vertdata, yearcol = "year2", year.as.random = TRUE)
summary(cypmatrix3)
#>
#> This lefkoMat object contains 5 matrices.
#>
#> Each matrix is a square matrix with 2916 rows and columns, and a total of 8503056 elements.
#> A total of 588520 survival-transitions were estimated, with 117704 per matrix.
#> A total of 12960 fecundity transitions were estimated, with 2592 per matrix.
#>
#> Vital rate modeling quality control:
#>
#> Survival estimated with 74 individuals and 320 individual transitions.
#> Observation estimated with 70 individuals and 303 individual transitions.
#> Size estimated with 70 individuals and 288 individual transitions.
#> Reproductive status estimated with 70 individuals and 288 individual transitions.
#> Fecundity estimated with 26 individuals and 50 individual transitions.
#> Juvenile survival estimated with 0 individuals and 0 individual transitions.
#> Juvenile observation estimated with 0 individuals and 0 individual transitions.
#> Juvenile size estimated with 0 individuals and 0 individual transitions.
#> Juvenile reproductive status estimated with 0 individuals and 0 individual transitions.
#> NULLOnce again, we see many more elements estimated, and many more rows and columns (54 rows and columns in the ahistorical case, and 542 = 2916 rows and columns in the historical case). However, the dominance of structural 0s in historical matrices still yields matrices that are mostly composed of 0s. A quick glance at one matrix will show that. We will focus on only one small section of that matrix.
cypmatrix3$A[[1]][2001:2050,2036:2045]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 0.000000e+00 0 0 0 0 0 0 0 0 0
#> [2,] 0.000000e+00 0 0 0 0 0 0 0 0 0
#> [3,] 0.000000e+00 0 0 0 0 0 0 0 0 0
#> [4,] 4.551147e-03 0 0 0 0 0 0 0 0 0
#> [5,] 3.218013e-05 0 0 0 0 0 0 0 0 0
#> [6,] 1.743937e-04 0 0 0 0 0 0 0 0 0
#> [7,] 6.300612e-04 0 0 0 0 0 0 0 0 0
#> [8,] 1.707245e-03 0 0 0 0 0 0 0 0 0
#> [9,] 3.700829e-03 0 0 0 0 0 0 0 0 0
#> [10,] 6.685300e-03 0 0 0 0 0 0 0 0 0
#> [11,] 1.035133e-02 0 0 0 0 0 0 0 0 0
#> [12,] 1.402423e-02 0 0 0 0 0 0 0 0 0
#> [13,] 1.688923e-02 0 0 0 0 0 0 0 0 0
#> [14,] 1.830555e-02 0 0 0 0 0 0 0 0 0
#> [15,] 1.803696e-02 0 0 0 0 0 0 0 0 0
#> [16,] 1.629128e-02 0 0 0 0 0 0 0 0 0
#> [17,] 1.358267e-02 0 0 0 0 0 0 0 0 0
#> [18,] 1.051550e-02 0 0 0 0 0 0 0 0 0
#> [19,] 7.598222e-03 0 0 0 0 0 0 0 0 0
#> [20,] 5.147130e-03 0 0 0 0 0 0 0 0 0
#> [21,] 3.281628e-03 0 0 0 0 0 0 0 0 0
#> [22,] 1.976014e-03 0 0 0 0 0 0 0 0 0
#> [23,] 1.127222e-03 0 0 0 0 0 0 0 0 0
#> [24,] 6.108754e-04 0 0 0 0 0 0 0 0 0
#> [25,] 3.152873e-04 0 0 0 0 0 0 0 0 0
#> [26,] 1.553306e-04 0 0 0 0 0 0 0 0 0
#> [27,] 7.319850e-05 0 0 0 0 0 0 0 0 0
#> [28,] 3.305705e-05 0 0 0 0 0 0 0 0 0
#> [29,] 1.794085e-04 0 0 0 0 0 0 0 0 0
#> [30,] 9.722686e-04 0 0 0 0 0 0 0 0 0
#> [31,] 3.512676e-03 0 0 0 0 0 0 0 0 0
#> [32,] 9.518122e-03 0 0 0 0 0 0 0 0 0
#> [33,] 2.063262e-02 0 0 0 0 0 0 0 0 0
#> [34,] 3.727145e-02 0 0 0 0 0 0 0 0 0
#> [35,] 5.771004e-02 0 0 0 0 0 0 0 0 0
#> [36,] 7.818699e-02 0 0 0 0 0 0 0 0 0
#> [37,] 9.415971e-02 0 0 0 0 0 0 0 0 0
#> [38,] 1.020559e-01 0 0 0 0 0 0 0 0 0
#> [39,] 1.005585e-01 0 0 0 0 0 0 0 0 0
#> [40,] 9.082608e-02 0 0 0 0 0 0 0 0 0
#> [41,] 7.572519e-02 0 0 0 0 0 0 0 0 0
#> [42,] 5.862535e-02 0 0 0 0 0 0 0 0 0
#> [43,] 4.236111e-02 0 0 0 0 0 0 0 0 0
#> [44,] 2.869594e-02 0 0 0 0 0 0 0 0 0
#> [45,] 1.829552e-02 0 0 0 0 0 0 0 0 0
#> [46,] 1.101654e-02 0 0 0 0 0 0 0 0 0
#> [47,] 6.284415e-03 0 0 0 0 0 0 0 0 0
#> [48,] 3.405712e-03 0 0 0 0 0 0 0 0 0
#> [49,] 1.757769e-03 0 0 0 0 0 0 0 0 0
#> [50,] 8.659888e-04 0 0 0 0 0 0 0 0 0Now let’s estimate some mean matrices. In this case, with all relevant matrix elements holding non-zero values, we may use the geometric mean across time, as follows for the ahistorical matrices.
And now the historical mean. We will use this and the ahistorical mean in analyses in the next step.
Step 5. Analyses
Let’s now conduct some quick population analyses. First, we can estimate the deterministic growth rate for these means. First the ahistorical case.
Now the historical case.
We find that the historical lambda is lower than the ahistorical lambda, suggesting that trade-offs are working across time to offset survival and/or reproduction.
Now let’s take a look at the stable stage distributions. First the ahistorical case, with summary.
tm2ss <- stablestage3(tmeans2r)
summary(tm2ss)
#> matrix new_stage_id orig_stage_id original_size ss_prop
#> Min. :1.0 Min. : 1.0 Length:108 Min. : 0.00 Min. :0.000e+00
#> 1st Qu.:1.0 1st Qu.:14.0 Class :character 1st Qu.: 4.00 1st Qu.:0.000e+00
#> Median :1.5 Median :27.5 Mode :character Median :11.00 Median :9.451e-05
#> Mean :1.5 Mean :27.5 Mean :11.11 Mean :1.852e-02
#> 3rd Qu.:2.0 3rd Qu.:41.0 3rd Qu.:18.00 3rd Qu.:1.934e-02
#> Max. :2.0 Max. :54.0 Max. :24.00 Max. :2.283e-01Now the first 6 rows.
head(tm2ss)
#> matrix new_stage_id orig_stage_id original_size ss_prop
#> 1 1 1 SD 0 0.205168634
#> 2 1 2 P1 0 0.228302122
#> 3 1 3 P2 0 0.051483779
#> 4 1 4 P3 0 0.011609947
#> 5 1 5 SL 0 0.005961283
#> 6 1 6 D 0 0.019337881Next, the historical case. We will skip the historical stage-pair distributions and look directly at the predicted distributions of ahistorical stages.
tm3ss <- stablestage3(tmeans3r)
summary(tm3ss$ahist)
#> matrix new_stage_id ss_prop
#> Min. :1.0 Min. : 1.0 Min. :0.0000755
#> 1st Qu.:1.0 1st Qu.:14.0 1st Qu.:0.0014256
#> Median :1.5 Median :27.5 Median :0.0065553
#> Mean :1.5 Mean :27.5 Mean :0.0185185
#> 3rd Qu.:2.0 3rd Qu.:41.0 3rd Qu.:0.0127742
#> Max. :2.0 Max. :54.0 Max. :0.3245035And the first 6 rows.
head(tm3ss$ahist)
#> matrix new_stage_id ss_prop
#> 1 1 1 0.289074019
#> 2 1 2 0.324503508
#> 3 1 3 0.017369247
#> 4 1 4 0.004257621
#> 5 1 5 0.002559200
#> 6 1 6 0.009451046Overall, these look like fairly large shifts. For comparison, we can try graphing these together.
plot(tm2ss$ss_prop, ylab = "Stable Stage Proportion", xlab = "Stage Index", ylim = c(0, 0.35), type = "l", lty = 1, lwd = 2)
lines(tm3ss$ahist[,3], col = "red", lty = 2, lwd = 2)
legend("topright", c("ahistorical", "historical"), lty = c(1, 2), col = c("black", "red"), lwd = 2, bty = "n")
plot of chunk Ch3an2st5ln710
Although the two have much general agreement, there is nonetheless evidence of some disagreement. We leave it to the user to assess these differences.
Finally, let’s take a peek at the reproductive values associated with both ahistorical and historical approaches. First, the ahistorical set.
repvalue3(tmeans2r)
#> matrix new_stage_id orig_stage_id original_size left_vector rep_value
#> 1 1 1 SD 0 0.0000000 0.000000
#> 2 1 2 P1 0 0.0000000 0.000000
#> 3 1 3 P2 0 0.0000000 0.000000
#> 4 1 4 P3 0 0.0000000 0.000000
#> 5 1 5 SL 0 0.0000000 0.000000
#> 6 1 6 D 0 -0.1149892 1.000000
#> 7 1 7 V1 1 -0.1230093 1.069747
#> 8 1 8 V2 2 -0.1281946 1.114841
#> 9 1 9 V3 3 -0.1315767 1.144253
#> 10 1 10 V4 4 -0.1338763 1.164251
#> 11 1 11 V5 5 -0.1355452 1.178765
#> 12 1 12 V6 6 -0.1368542 1.190149
#> 13 1 13 V7 7 -0.1379633 1.199794
#> 14 1 14 V8 8 -0.1389671 1.208524
#> 15 1 15 V9 9 -0.1399212 1.216821
#> 16 1 16 V10 10 -0.1408574 1.224963
#> 17 1 17 V11 11 -0.1417924 1.233094
#> 18 1 18 V12 12 -0.1427327 1.241271
#> 19 1 19 V13 13 -0.1436783 1.249495
#> 20 1 20 V14 14 -0.1446251 1.257728
#> 21 1 21 V15 15 -0.1455662 1.265913
#> 22 1 22 V16 16 -0.1464943 1.273983
#> 23 1 23 V17 17 -0.1474018 1.281876
#> 24 1 24 V18 18 -0.1482822 1.289532
#> 25 1 25 V19 19 -0.1491304 1.296908
#> 26 1 26 V20 20 -0.1499427 1.303972
#> 27 1 27 V21 21 -0.1507169 1.310705
#> 28 1 28 V22 22 -0.1514522 1.317100
#> 29 1 29 V23 23 -0.1521491 1.323160
#> 30 1 30 V24 24 -0.1528086 1.328896
#> 31 1 31 F1 1 -0.1218207 1.059411
#> 32 1 32 F2 2 -0.1274049 1.107974
#> 33 1 33 F3 3 -0.1312803 1.141676
#> 34 1 34 F4 4 -0.1341254 1.166417
#> 35 1 35 F5 5 -0.1363550 1.185807
#> 36 1 36 F6 6 -0.1382082 1.201924
#> 37 1 37 F7 7 -0.1398184 1.215927
#> 38 1 38 F8 8 -0.1412592 1.228456
#> 39 1 39 F9 9 -0.1425710 1.239864
#> 40 1 40 F10 10 -0.1437769 1.250352
#> 41 1 41 F11 11 -0.1448909 1.260040
#> 42 1 42 F12 12 -0.1459225 1.269011
#> 43 1 43 F13 13 -0.1468789 1.277328
#> 44 1 44 F14 14 -0.1477666 1.285048
#> 45 1 45 F15 15 -0.1485914 1.292221
#> 46 1 46 F16 16 -0.1493591 1.298898
#> 47 1 47 F17 17 -0.1500753 1.305125
#> 48 1 48 F18 18 -0.1507451 1.310951
#> 49 1 49 F19 19 -0.1513738 1.316418
#> 50 1 50 F20 20 -0.1519658 1.321567
#> 51 1 51 F21 21 -0.1525257 1.326435
#> 52 1 52 F22 22 -0.1530571 1.331057
#> 53 1 53 F23 23 -0.1535637 1.335462
#> 54 1 54 F24 24 -0.1540484 1.339677
#> 55 2 1 SD 0 0.0000000 0.000000
#> 56 2 2 P1 0 0.0000000 0.000000
#> 57 2 3 P2 0 0.0000000 0.000000
#> 58 2 4 P3 0 0.0000000 0.000000
#> 59 2 5 SL 0 0.0000000 0.000000
#> 60 2 6 D 0 -0.1149892 1.000000
#> 61 2 7 V1 1 -0.1230093 1.069747
#> 62 2 8 V2 2 -0.1281946 1.114841
#> 63 2 9 V3 3 -0.1315767 1.144253
#> 64 2 10 V4 4 -0.1338763 1.164251
#> 65 2 11 V5 5 -0.1355452 1.178765
#> 66 2 12 V6 6 -0.1368542 1.190149
#> 67 2 13 V7 7 -0.1379633 1.199794
#> 68 2 14 V8 8 -0.1389671 1.208524
#> 69 2 15 V9 9 -0.1399212 1.216821
#> 70 2 16 V10 10 -0.1408574 1.224963
#> 71 2 17 V11 11 -0.1417924 1.233094
#> 72 2 18 V12 12 -0.1427327 1.241271
#> 73 2 19 V13 13 -0.1436783 1.249495
#> 74 2 20 V14 14 -0.1446251 1.257728
#> 75 2 21 V15 15 -0.1455662 1.265913
#> 76 2 22 V16 16 -0.1464943 1.273983
#> 77 2 23 V17 17 -0.1474018 1.281876
#> 78 2 24 V18 18 -0.1482822 1.289532
#> 79 2 25 V19 19 -0.1491304 1.296908
#> 80 2 26 V20 20 -0.1499427 1.303972
#> 81 2 27 V21 21 -0.1507169 1.310705
#> 82 2 28 V22 22 -0.1514522 1.317100
#> 83 2 29 V23 23 -0.1521491 1.323160
#> 84 2 30 V24 24 -0.1528086 1.328896
#> 85 2 31 F1 1 -0.1218207 1.059411
#> 86 2 32 F2 2 -0.1274049 1.107974
#> 87 2 33 F3 3 -0.1312803 1.141676
#> 88 2 34 F4 4 -0.1341254 1.166417
#> 89 2 35 F5 5 -0.1363550 1.185807
#> 90 2 36 F6 6 -0.1382082 1.201924
#> 91 2 37 F7 7 -0.1398184 1.215927
#> 92 2 38 F8 8 -0.1412592 1.228456
#> 93 2 39 F9 9 -0.1425710 1.239864
#> 94 2 40 F10 10 -0.1437769 1.250352
#> 95 2 41 F11 11 -0.1448909 1.260040
#> 96 2 42 F12 12 -0.1459225 1.269011
#> 97 2 43 F13 13 -0.1468789 1.277328
#> 98 2 44 F14 14 -0.1477666 1.285048
#> 99 2 45 F15 15 -0.1485914 1.292221
#> 100 2 46 F16 16 -0.1493591 1.298898
#> 101 2 47 F17 17 -0.1500753 1.305125
#> 102 2 48 F18 18 -0.1507451 1.310951
#> 103 2 49 F19 19 -0.1513738 1.316418
#> 104 2 50 F20 20 -0.1519658 1.321567
#> 105 2 51 F21 21 -0.1525257 1.326435
#> 106 2 52 F22 22 -0.1530571 1.331057
#> 107 2 53 F23 23 -0.1535637 1.335462
#> 108 2 54 F24 24 -0.1540484 1.339677And let’s compare these to the historical case.
repvalue3(tmeans3r)$ahist
#> matrix new_stage_id rep_value
#> 1 1 1 0.00000000
#> 2 1 2 0.00000000
#> 3 1 3 0.00000000
#> 4 1 4 0.00000000
#> 5 1 5 0.00000000
#> 6 1 6 1.00000000
#> 7 1 7 1.11337163
#> 8 1 8 1.25989440
#> 9 1 9 1.36130980
#> 10 1 10 1.40413383
#> 11 1 11 1.37804439
#> 12 1 12 1.28111047
#> 13 1 13 1.12415468
#> 14 1 14 0.93071691
#> 15 1 15 0.73091305
#> 16 1 16 0.55172715
#> 17 1 17 0.40917911
#> 18 1 18 0.30650513
#> 19 1 19 0.23796212
#> 20 1 20 0.19442578
#> 21 1 21 0.16729948
#> 22 1 22 0.15005974
#> 23 1 23 0.13827432
#> 24 1 24 0.12913590
#> 25 1 25 0.12102582
#> 26 1 26 0.11319538
#> 27 1 27 0.10544336
#> 28 1 28 0.09785310
#> 29 1 29 0.09055702
#> 30 1 30 0.08363655
#> 31 1 31 1.24482371
#> 32 1 32 1.35373424
#> 33 1 33 1.40137667
#> 34 1 34 1.37966649
#> 35 1 35 1.28654989
#> 36 1 36 1.13098199
#> 37 1 37 0.93537863
#> 38 1 38 0.73167196
#> 39 1 39 0.55062241
#> 40 1 40 0.41041837
#> 41 1 41 0.31324037
#> 42 1 42 0.25087072
#> 43 1 43 0.21246065
#> 44 1 44 0.18901643
#> 45 1 45 0.17439549
#> 46 1 46 0.16473622
#> 47 1 47 0.15766164
#> 48 1 48 0.15170294
#> 49 1 49 0.14596621
#> 50 1 50 0.13995298
#> 51 1 51 0.13345689
#> 52 1 52 0.12647446
#> 53 1 53 0.11910600
#> 54 1 54 0.11147341
#> 55 2 1 0.00000000
#> 56 2 2 0.00000000
#> 57 2 3 0.00000000
#> 58 2 4 0.00000000
#> 59 2 5 0.00000000
#> 60 2 6 1.00000000
#> 61 2 7 1.11337163
#> 62 2 8 1.25989440
#> 63 2 9 1.36130980
#> 64 2 10 1.40413383
#> 65 2 11 1.37804439
#> 66 2 12 1.28111047
#> 67 2 13 1.12415468
#> 68 2 14 0.93071691
#> 69 2 15 0.73091305
#> 70 2 16 0.55172715
#> 71 2 17 0.40917911
#> 72 2 18 0.30650513
#> 73 2 19 0.23796212
#> 74 2 20 0.19442578
#> 75 2 21 0.16729948
#> 76 2 22 0.15005974
#> 77 2 23 0.13827432
#> 78 2 24 0.12913590
#> 79 2 25 0.12102582
#> 80 2 26 0.11319538
#> 81 2 27 0.10544336
#> 82 2 28 0.09785310
#> 83 2 29 0.09055702
#> 84 2 30 0.08363655
#> 85 2 31 1.24482371
#> 86 2 32 1.35373424
#> 87 2 33 1.40137667
#> 88 2 34 1.37966649
#> 89 2 35 1.28654989
#> 90 2 36 1.13098199
#> 91 2 37 0.93537863
#> 92 2 38 0.73167196
#> 93 2 39 0.55062241
#> 94 2 40 0.41041837
#> 95 2 41 0.31324037
#> 96 2 42 0.25087072
#> 97 2 43 0.21246065
#> 98 2 44 0.18901643
#> 99 2 45 0.17439549
#> 100 2 46 0.16473622
#> 101 2 47 0.15766164
#> 102 2 48 0.15170294
#> 103 2 49 0.14596621
#> 104 2 50 0.13995298
#> 105 2 51 0.13345689
#> 106 2 52 0.12647446
#> 107 2 53 0.11910600
#> 108 2 54 0.11147341Note that the historical case predicts a greater importance of small adults and lesser importance of large adults than the ahistorical matrix does.