Chapter 6 Matrix Models III: Age-by-Stage and Leslie MPMs
“The really frightening thing about middle age is that you know you’ll grow out of it.”
The first MPMs were age-classified, and are also known as Leslie MPMs in honor of their creator (Leslie 1945). However, in many organisms age is seen as less important in determining demography than other variables. Stage-classified, or Lefkovitch, MPMs were developed with this in mind, as they allow the life history of an organism to be stratified by life history stages (Lefkovitch 1965). Life history stages can be defined as necessary by any status variable or combination of status variables (Caswell 2001). So stage-based MPMs provide a great deal of flexibility and power to analyze population dynamics.
For much of the history of the MPM, the dichotomy between age-classified and stage-classified MPMs was driven by the lack of computing power to analyze complex relationships in large matrices. With the advent of powerful home computers, more and more ecologists have addressed questions using age-by-stage MPMs, which are MPMs that incorporate both and age and stage (Caswell & Salguero-Gómez 2013). Here, we provide a short introduction to the core theory before building an age-by-stage MPM. We strongly advise readers interested in this style of MPM to read Caswell et al. (2018), which is the most comprehensive overview of the theory behind age-by-stage MPMs that we are aware of.
6.1 What is an age-by-stage MPM?
To answer the question of what an age-by-stage MPM is, let us first address the question of what an age-classified MPM is. Age-classified, or Leslie, MPMs treat the life history of an organism as a series of ages rather than stages. Each age is of equal duration, and the number of ages to include is decided through demographic studies that attempt to determine which groups of ages have unique demographic characteristics. The only transitions that are possible are one-way survival transitions from one age to the next, survival transitions between the final age and itself if the lifespan has no hard limit, and fecundity from any reproductive age to age 0 (if post-breeding) or 1 (if pre-breeding).
Let’s look at an example. Let us suppose that we are studying an organism that can live for potentially many years, such as a long-lived seabird. We conduct several years of monitoring, and find that the seabird at our monitored population seems to follow a basic pattern in which survival is relatively low in the first year, higher in the second, and from the third year on survival is typically quite high and stable. Although the maximum longevity is unknown, nonetheless birds seem to live not much longer than 20 years. Fecundity occurs from the third year on, with the first two years spent as a juvenile. In this circumstance, we might develop a Leslie matrix that looks like this:
\[\begin{equation} \tag{6.1} \left[\begin{array} {rrr} 0 & 0 & F _{1,3} \\ S _{2,1} & 0 & 0 \\ 0 & S _{3,2} & S _{3,3} \\ \end{array}\right] \end{equation}\]
This matrix includes only three rows and columns. Survival-transition probabilities are typically just below the diagonal, and represent the probability of surviving from age \(i\) (column \(i\)) to age \(i+1\) (row \(i+1\)). There is a notable exception in element \(S _{3,3}\), which is the survival-transition probability for an organism to stay alive within the final age. This allows us to model the tendency of organisms to live on without an explicit maximum longevity in a state in which the probability of survival is not expected to change. Fecundity is in the top row (\(F_{1,3}\)), since all offspring begin in the first age.
Age-by-stage MPMs differ from the above because they incorporate both age and stage, and so the matrix needs to be expanded to allow both characteristics. The standard approach is to create block matrices in which each age has potentially several stages, and survival-transition probabilities must describe the probability of transitioning from each stage at one age to each stage in the next.
As an example, let’s take the three stage life history above, which has a single newborn stage and two adult stages, only one of which is reproductive. Let’s further assume that age impacts survival and fecundity only for the first four years, beyond which survival and fecundity transitions remain essentially the same. This leads to the following age-by-stage matrix:
\[\begin{equation} \tag{6.2} \left[\begin{array} {rrrrrrr} 0 & 0 & F _{1A,2C} & 0 & F _{1A,3C} & 0 & F _{1A,4+C} \\ S _{2B,1A} & 0 & 0 & 0 & 0 & 0 & 0 \\ S _{2C,1A} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & S _{3B,2B} & S _{3B,2C} & 0 & 0 & 0 & 0 \\ 0 & S _{3C,2B} & S _{3C,2C} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & S _{4B,3B} & S _{4B,3C} & S _{5+B,4B} & S _{5+B,4+C} \\ 0 & 0 & 0 & S _{4C,3B} & S _{4C,3C} & S _{5+C,4B} & S _{5+C,4+C} \\ \end{array}\right] \end{equation}\]
Here, we have only one possible initial stage, followed by two possible stages in all later ages. Looking only at the columns, column 1 refers to age 1, which only includes the seedling stage; columns 2 and 3 refer to age 2, which includes non-reproductive and reproductive adult stages; columns 4 and 5 refer to age 3, which also includes non-reproductive and reproductive adult stages; and columns 6 and 7 refer to ages 4 and higher, once again also including non-reproductive and reproductive adult stages. While the columns refer to the From ages and stages, the rows refer to the To ages and stages and are in the same order as the columns. If we ignore stage and consider only the placement of age-related transitions, then we see that we are still generally following the Leslie MPM pattern of placing survival-transition probabilities below the diagonal for all but the final transitions. We also have more fecundity terms because although only one stage is reproductive, multiple reproductive ages occur, and all of these stay within the first row (more entry rows are possible only if the life history of the organism can start from multiple different stages).
The above matrix is ahistorical. However, unlike a pure Leslie MPM, the possibility of different stages within each age allows for the estimation of historical age-by-stage MPMs. Package lefko3 does not currently support historical age-by-stage matrices, but this capability will eventually be added.
6.2 Developing function-based age-by-stage MPMs in lefko3
Package lefko3 can produce both raw (empirical) and function-based age-by-stage MPMs. Let’s start with the latter. For our example, we will use the lathyrus dataset that comes with lefko3 to illustrate the estimation of age-by-stage function-based MPMs. First, let’s load the dataset, and then look at its dimensions.
data(lathyrus)
dim(lathyrus)
> [1] 1119 38
This dataset includes information on 1,119 individuals, so there are 1,119 rows with data (not counting the header). There are 38 columns. The first two columns are variables identifying each individual (SUBPLOT refers to the patch, and GENET refers to individual identity), with each individual’s data entirely restricted to one row. This is followed by four sets of nine columns, each named VolumeXX, lnVolXX, FCODEXX, FlowXX, IntactseedXX, Dead19XX, DormantXX, Missing19XX, and SeedlingXX, where XX corresponds to the year of observation and with years organized consecutively. Thus, columns 3-11 refer to year 1988, columns 12-20 refer to year 1989, etc.
To begin, we will create a stageframe that describes the organism’s life cycle for this dataset. In this case, the life history model is a life cycle graph (figure 6.1).
Figure 6.1: Life history model of Lathyrus vernus using log leaf volume as the size classification metric
This model is based on the life history model provided in Ehrlén (2000), but we utilize a different size classification based on the log leaf volume to make the size distribution more closely match a symmetric and somewhat normal distribution.
Our stageframe needs to include complete descriptions of all stages that occur in the dataset, with each stage defined uniquely, and also needs to describe the ages for each stage portrayed in our life history model. Since this object can be used for automated classification of individuals, all sizes, reproductive states, and other characteristics defining each proper stage in the dataset need to be accounted for explicitly. This can be difficult if a few data points exist outside the range of sizes specified in the stageframe. Such points can cause problems, because rare stages can cause an overestimation of survival for existing stages under some circumstances, and can also yield spurious values for survival-transition probabilities and fecundity rates. The final description of each stage occurring in the dataset must also avoid complete overlap with any other stage also found in the dataset, although partial overlap is allowed and expected.
Before creating the stage frame, let’s explore the possible size variables. We will particularly look at summaries of the distribution of original and log sizes.
summary(c(lathyrus$Volume88, lathyrus$Volume89, lathyrus$Volume90,
lathyrus$Volume91))
> Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
> 1.8 14.7 123.0 484.2 732.5 7032.0 1248
summary(c(lathyrus$lnVol88, lathyrus$lnVol89, lathyrus$lnVol90,
lathyrus$lnVol91))
> Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
> 0.600 2.700 4.800 4.777 6.600 8.900 1248The upper summary shows the original size, while the lower line shows the size given in logarithmic terms. We should note the size minima and maxima, because we have been using 0 as the size of vegetatively dormant individuals. The lowest uncorrected size is 1.8, with a maximum of 7032.0. The minimum corrected size is 0.600, and the maximum corrected size is 8.900. Since the minimum corrected size is above 0 (i.e. all log sizes should be positive), and since the number of NAs has not increased (increased NAs would suggest that some unusable log sizes occur in the dataset), we are still able to use the log size value 0 as an indicator of vegetative dormancy. Note, however, that vegetative dormancy is also currently included in the many NAs that occur in size variables in this dataset.
It can also help to take a look at plots of these distributions. We will plot raw and log volume (figure 6.2).
par(mfrow=c(1,2))
plot(density(c(lathyrus$Volume88, lathyrus$Volume89, lathyrus$Volume90,
lathyrus$Volume91), na.rm = TRUE), main = "", xlab = "Volume", bty = "n")
title("a)", adj = 0)
plot(density(c(lathyrus$lnVol88, lathyrus$lnVol89, lathyrus$lnVol90,
lathyrus$lnVol91), na.rm = TRUE), main = "", xlab = "Log volume", bty = "n")
title("b)", adj = 0)
Figure 6.2: Density plot of aboveground plant volume (a) and log volume (b).
Note how highly skewed the raw volume distribution is. This might cause some difficulty if we use raw size untransformed and with a Gaussian distribution. Certainly, a gamma distribution would be perfectly justified, and users are urged to try that approach. We will use the log volume here, which looks ‘better’ than the raw volume distribution in the sense that it is closer to some semblance of a Gaussian distribution, mostly through an increased level of symmetry. We can then assume that log volume is Gaussian-distributed and that the mean bears no relationship to the variance.
We will now develop a stageframe that incorporates the log volume of size. We will build this by creating vectors of the values describing each stage, always in the same order. Because we wish to build an age-by-stage MPM, we will also incorporate age information for each stage. Here, we include minimum and maximum ages for each stage via the vectors minima and maxima (NA in the maximum age vector is interpreted as meaning that there is no maximum).
sizevector <- c(0, 4.6, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9)
stagevector <- c("Sd", "Sdl", "Dorm", "Sz1nr", "Sz2nr", "Sz3nr", "Sz4nr",
"Sz5nr","Sz6nr", "Sz7nr", "Sz8nr", "Sz9nr", "Sz1r", "Sz2r", "Sz3r", "Sz4r",
"Sz5r", "Sz6r", "Sz7r", "Sz8r", "Sz9r")
repvector <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1)
obsvector <- c(0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
matvector <- c(0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
immvector <- c(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
indataset <- c(0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
binvec <- c(0, 4.6, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5)
minima <- c(1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)
maxima <- c(NA, 1, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA)
comments <- c("Dormant seed", "Seedling", "Dormant", "Size 1 Veg", "Size 2 Veg",
"Size 3 Veg", "Size 4 Veg", "Size 5 Veg", "Size 6 Veg", "Size 7 Veg",
"Size 8 Veg", "Size 9 Veg", "Size 1 Flo", "Size 2 Flo", "Size 3 Flo",
"Size 4 Flo", "Size 5 Flo", "Size 6 Flo", "Size 7 Flo", "Size 8 Flo",
"Size 9 Flo")
lathframeln <- sf_create(sizes = sizevector, stagenames = stagevector,
repstatus = repvector, obsstatus = obsvector, propstatus = propvector,
immstatus = immvector, matstatus = matvector, indataset = indataset,
binhalfwidth = binvec, minage = minima, maxage = maxima, comments = comments)
lathframeln
> stage size size_b size_c min_age max_age repstatus obsstatus propstatus
> 1 Sd 0.0 NA NA 1 NA 0 0 1
> 2 Sdl 4.6 NA NA 1 1 0 1 0
> 3 Dorm 0.0 NA NA 2 NA 0 0 0
> 4 Sz1nr 1.0 NA NA 2 NA 0 1 0
> 5 Sz2nr 2.0 NA NA 2 NA 0 1 0
> 6 Sz3nr 3.0 NA NA 2 NA 0 1 0
> 7 Sz4nr 4.0 NA NA 2 NA 0 1 0
> 8 Sz5nr 5.0 NA NA 2 NA 0 1 0
> 9 Sz6nr 6.0 NA NA 2 NA 0 1 0
> 10 Sz7nr 7.0 NA NA 2 NA 0 1 0
> 11 Sz8nr 8.0 NA NA 2 NA 0 1 0
> 12 Sz9nr 9.0 NA NA 2 NA 0 1 0
> 13 Sz1r 1.0 NA NA 2 NA 1 1 0
> 14 Sz2r 2.0 NA NA 2 NA 1 1 0
> 15 Sz3r 3.0 NA NA 2 NA 1 1 0
> 16 Sz4r 4.0 NA NA 2 NA 1 1 0
> 17 Sz5r 5.0 NA NA 2 NA 1 1 0
> 18 Sz6r 6.0 NA NA 2 NA 1 1 0
> 19 Sz7r 7.0 NA NA 2 NA 1 1 0
> 20 Sz8r 8.0 NA NA 2 NA 1 1 0
> 21 Sz9r 9.0 NA NA 2 NA 1 1 0
> immstatus matstatus indataset binhalfwidth_raw sizebin_min sizebin_max
> 1 1 0 0 0.0 0.0 0.0
> 2 1 0 1 4.6 0.0 9.2
> 3 0 1 1 0.5 -0.5 0.5
> 4 0 1 1 0.5 0.5 1.5
> 5 0 1 1 0.5 1.5 2.5
> 6 0 1 1 0.5 2.5 3.5
> 7 0 1 1 0.5 3.5 4.5
> 8 0 1 1 0.5 4.5 5.5
> 9 0 1 1 0.5 5.5 6.5
> 10 0 1 1 0.5 6.5 7.5
> 11 0 1 1 0.5 7.5 8.5
> 12 0 1 1 0.5 8.5 9.5
> 13 0 1 1 0.5 0.5 1.5
> 14 0 1 1 0.5 1.5 2.5
> 15 0 1 1 0.5 2.5 3.5
> 16 0 1 1 0.5 3.5 4.5
> 17 0 1 1 0.5 4.5 5.5
> 18 0 1 1 0.5 5.5 6.5
> 19 0 1 1 0.5 6.5 7.5
> 20 0 1 1 0.5 7.5 8.5
> 21 0 1 1 0.5 8.5 9.5
> sizebin_center sizebin_width binhalfwidthb_raw sizebinb_min sizebinb_max
> 1 0.0 0.0 NA NA NA
> 2 4.6 9.2 NA NA NA
> 3 0.0 1.0 NA NA NA
> 4 1.0 1.0 NA NA NA
> 5 2.0 1.0 NA NA NA
> 6 3.0 1.0 NA NA NA
> 7 4.0 1.0 NA NA NA
> 8 5.0 1.0 NA NA NA
> 9 6.0 1.0 NA NA NA
> 10 7.0 1.0 NA NA NA
> 11 8.0 1.0 NA NA NA
> 12 9.0 1.0 NA NA NA
> 13 1.0 1.0 NA NA NA
> 14 2.0 1.0 NA NA NA
> 15 3.0 1.0 NA NA NA
> 16 4.0 1.0 NA NA NA
> 17 5.0 1.0 NA NA NA
> 18 6.0 1.0 NA NA NA
> 19 7.0 1.0 NA NA NA
> 20 8.0 1.0 NA NA NA
> 21 9.0 1.0 NA NA NA
> sizebinb_center sizebinb_width binhalfwidthc_raw sizebinc_min sizebinc_max
> 1 NA NA NA NA NA
> 2 NA NA NA NA NA
> 3 NA NA NA NA NA
> 4 NA NA NA NA NA
> 5 NA NA NA NA NA
> 6 NA NA NA NA NA
> 7 NA NA NA NA NA
> 8 NA NA NA NA NA
> 9 NA NA NA NA NA
> 10 NA NA NA NA NA
> 11 NA NA NA NA NA
> 12 NA NA NA NA NA
> 13 NA NA NA NA NA
> 14 NA NA NA NA NA
> 15 NA NA NA NA NA
> 16 NA NA NA NA NA
> 17 NA NA NA NA NA
> 18 NA NA NA NA NA
> 19 NA NA NA NA NA
> 20 NA NA NA NA NA
> 21 NA NA NA NA NA
> sizebinc_center sizebinc_width group comments
> 1 NA NA 0 Dormant seed
> 2 NA NA 0 Seedling
> 3 NA NA 0 Dormant
> 4 NA NA 0 Size 1 Veg
> 5 NA NA 0 Size 2 Veg
> 6 NA NA 0 Size 3 Veg
> 7 NA NA 0 Size 4 Veg
> 8 NA NA 0 Size 5 Veg
> 9 NA NA 0 Size 6 Veg
> 10 NA NA 0 Size 7 Veg
> 11 NA NA 0 Size 8 Veg
> 12 NA NA 0 Size 9 Veg
> 13 NA NA 0 Size 1 Flo
> 14 NA NA 0 Size 2 Flo
> 15 NA NA 0 Size 3 Flo
> 16 NA NA 0 Size 4 Flo
> 17 NA NA 0 Size 5 Flo
> 18 NA NA 0 Size 6 Flo
> 19 NA NA 0 Size 7 Flo
> 20 NA NA 0 Size 8 Flo
> 21 NA NA 0 Size 9 Flo
Once the stageframe is created, we can reorganize the dataset into historically-formatted vertical (hfv) format. To handle this, we will use the verticalize3() function, which creates historically-formatted vertical datasets, as below. We will also replace NAs in size and fecundity variables with zeros for modelsearch to work properly when we build models of vital rates, so we will now set NAas0 = TRUE. Some care needs to be taken with this last step, since some authors give missing values extra meaning not present in a value of zero. In our case, a missing value indicates that a plant was dead (both size and fecundity are missing), was alive but not sprouting (size was missing), or was alive but did not produce seed (fecundity was missing). In all cases, these NA values may be replaced by zero, because other variables indicate those conditions.
We also have two choices for use as our reproductive status and fecundity variables. The first choice, FCODE88 indicates whether a plant flowered. The second choice, Intactseed88, indicates the number of seed produced. The choice of which to use depends strongly on the aims of the study. In our case, we would like to treat all plants that flowered as reproductive, but treat fecundity in terms of real seed produced. The reason is that we believe that flowering plants have a different demography than non-flowering plants, either reflecting reproductive costs, or, conversely, because flowering plants might have more resources and hence higher survival than non-flowering plants, and so we wish to separate transitions among these two groups. So, let’s use FCODE88 to indicate reproductive status, and Intactseed88 to indicate fecundity. Once complete, we will look at a summary.
Finally, note that in the input to the following function, we utilize a strictly repeating pattern of variable names arranged in the same order for each monitoring occasion. This arrangement allows us to enter only the first variable in each set, as long as noyears and blocksize are set properly and no gaps or shuffles appear in the dataset. The data management functions that we have created for lefko3 do not require such repeating patterns, but they do make the required input in the function much shorter and more succinct.
lathvertln <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
juvcol = "Seedling1988", sizeacol = "lnVol88", repstracol = "FCODE88",
fecacol = "Intactseed88", deadacol = "Dead1988", nonobsacol = "Dormant1988",
stageassign = lathframeln, stagesize = "sizea", censorcol = "Missing1988",
censorkeep = NA, NAas0 = TRUE, censor = TRUE)
dim(lathvertln)
> [1] 2552 42
summary(lathvertln)
> rowid popid patchid individ year2
> Min. : 1.0 :2552 1:674 Length:2552 Min. :1988
> 1st Qu.: 239.8 2:206 Class :character 1st Qu.:1988
> Median : 523.0 3:560 Mode :character Median :1989
> Mean : 538.7 4:478 Mean :1989
> 3rd Qu.: 823.2 5:312 3rd Qu.:1990
> Max. :1119.0 6:322 Max. :1990
> firstseen lastseen obsage obslifespan sizea1
> Min. :1988 Min. :1988 Min. :1.000 Min. :0.00 Min. :0.000
> 1st Qu.:1988 1st Qu.:1991 1st Qu.:1.000 1st Qu.:2.00 1st Qu.:0.000
> Median :1988 Median :1991 Median :2.000 Median :3.00 Median :2.200
> Mean :1988 Mean :1991 Mean :1.829 Mean :2.45 Mean :2.947
> 3rd Qu.:1988 3rd Qu.:1991 3rd Qu.:2.000 3rd Qu.:3.00 3rd Qu.:6.400
> Max. :1990 Max. :1991 Max. :3.000 Max. :3.00 Max. :8.900
> repstra1 feca1 juvgiven1 obsstatus1
> Min. :0.0000 Min. : 0.0000 Min. :0.00000 Min. :0.0000
> 1st Qu.:0.0000 1st Qu.: 0.0000 1st Qu.:0.00000 1st Qu.:0.0000
> Median :0.0000 Median : 0.0000 Median :0.00000 Median :1.0000
> Mean :0.1795 Mean : 0.9793 Mean :0.06309 Mean :0.5537
> 3rd Qu.:0.0000 3rd Qu.: 0.0000 3rd Qu.:0.00000 3rd Qu.:1.0000
> Max. :1.0000 Max. :66.0000 Max. :1.00000 Max. :1.0000
> repstatus1 fecstatus1 matstatus1 alive1
> Min. :0.0000 Min. :0.00000 Min. :0.0000 Min. :0.0000
> 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.0000
> Median :0.0000 Median :0.00000 Median :1.0000 Median :1.0000
> Mean :0.1795 Mean :0.08895 Mean :0.5192 Mean :0.5823
> 3rd Qu.:0.0000 3rd Qu.:0.00000 3rd Qu.:1.0000 3rd Qu.:1.0000
> Max. :1.0000 Max. :1.00000 Max. :1.0000 Max. :1.0000
> stage1 stage1index sizea2 repstra2
> Length:2552 Min. : 0.000 Min. :0.000 Min. :0.0000
> Class :character 1st Qu.: 0.000 1st Qu.:2.500 1st Qu.:0.0000
> Mode :character Median : 3.000 Median :4.700 Median :0.0000
> Mean : 6.115 Mean :4.562 Mean :0.2355
> 3rd Qu.:10.000 3rd Qu.:6.600 3rd Qu.:0.0000
> Max. :21.000 Max. :8.900 Max. :1.0000
> feca2 juvgiven2 obsstatus2 repstatus2
> Min. : 0.000 Min. :0.0000 Min. :0.0000 Min. :0.0000
> 1st Qu.: 0.000 1st Qu.:0.0000 1st Qu.:1.0000 1st Qu.:0.0000
> Median : 0.000 Median :0.0000 Median :1.0000 Median :0.0000
> Mean : 1.136 Mean :0.1117 Mean :0.9459 Mean :0.2355
> 3rd Qu.: 0.000 3rd Qu.:0.0000 3rd Qu.:1.0000 3rd Qu.:0.0000
> Max. :66.000 Max. :1.0000 Max. :1.0000 Max. :1.0000
> fecstatus2 matstatus2 alive2 stage2
> Min. :0.0000 Min. :0.0000 Min. :1 Length:2552
> 1st Qu.:0.0000 1st Qu.:1.0000 1st Qu.:1 Class :character
> Median :0.0000 Median :1.0000 Median :1 Mode :character
> Mean :0.1046 Mean :0.8883 Mean :1
> 3rd Qu.:0.0000 3rd Qu.:1.0000 3rd Qu.:1
> Max. :1.0000 Max. :1.0000 Max. :1
> stage2index sizea3 repstra3 feca3 juvgiven3
> Min. : 2.000 Min. :0.000 Min. :0.0000 Min. : 0.000 Min. :0
> 1st Qu.: 5.000 1st Qu.:2.300 1st Qu.:0.0000 1st Qu.: 0.000 1st Qu.:0
> Median : 8.000 Median :4.300 Median :0.0000 Median : 0.000 Median :0
> Mean : 9.305 Mean :4.003 Mean :0.2159 Mean : 1.277 Mean :0
> 3rd Qu.:10.000 3rd Qu.:6.300 3rd Qu.:0.0000 3rd Qu.: 0.000 3rd Qu.:0
> Max. :21.000 Max. :8.800 Max. :1.0000 Max. :66.000 Max. :0
> obsstatus3 repstatus3 fecstatus3 matstatus3
> Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :1
> 1st Qu.:1.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:1
> Median :1.0000 Median :0.0000 Median :0.0000 Median :1
> Mean :0.8264 Mean :0.2159 Mean :0.1105 Mean :1
> 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:1
> Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1
> alive3 stage3 stage3index
> Min. :0.0000 Length:2552 Min. : 0.000
> 1st Qu.:1.0000 Class :character 1st Qu.: 5.000
> Median :1.0000 Mode :character Median : 7.000
> Mean :0.9138 Mean : 8.632
> 3rd Qu.:1.0000 3rd Qu.:10.000
> Max. :1.0000 Max. :21.000This dataset has 2,552 rows representing our original dataset of over 1100 individuals. Ordinarily, we would now go on to produce either vital rate models to create function-based MPMs. However, the fact that we are incorporating age in our analysis leads to the problem that there are many individuals in our dataset that are of unknown age. The hfv dataset that we have created includes an estimated age for each individual in each year, but this is estimated from the time in which the individual is first seen. Many individuals are first seen in the first year of the study, by which time many could have already been alive for years. So, we will subset our data to eliminate individuals first seen in the very first year of the study.
lathvertln_small <- subset(lathvertln, firstseen > 1988)
dim(lathvertln_small)
> [1] 542 42
summary(lathvertln_small)
> rowid popid patchid individ year2
> Min. : 30.0 :542 1:131 Length:542 Min. :1989
> 1st Qu.: 302.5 2: 72 Class :character 1st Qu.:1989
> Median : 583.0 3:101 Mode :character Median :1990
> Mean : 586.8 4:132 Mean :1990
> 3rd Qu.: 801.0 5: 70 3rd Qu.:1990
> Max. :1097.0 6: 36 Max. :1990
> firstseen lastseen obsage obslifespan sizea1
> Min. :1989 Min. :1989 Min. :1.000 Min. :0.00 Min. :0.000
> 1st Qu.:1989 1st Qu.:1990 1st Qu.:1.000 1st Qu.:1.00 1st Qu.:0.000
> Median :1989 Median :1991 Median :1.000 Median :2.00 Median :0.000
> Mean :1989 Mean :1991 Mean :1.352 Mean :1.36 Mean :1.186
> 3rd Qu.:1989 3rd Qu.:1991 3rd Qu.:2.000 3rd Qu.:2.00 3rd Qu.:2.200
> Max. :1990 Max. :1991 Max. :2.000 Max. :2.00 Max. :8.400
> repstra1 feca1 juvgiven1 obsstatus1
> Min. :0.0000 Min. : 0.0000 Min. :0.0000 Min. :0.0000
> 1st Qu.:0.0000 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.:0.0000
> Median :0.0000 Median : 0.0000 Median :0.0000 Median :0.0000
> Mean :0.0203 Mean : 0.1218 Mean :0.1624 Mean :0.3524
> 3rd Qu.:0.0000 3rd Qu.: 0.0000 3rd Qu.:0.0000 3rd Qu.:1.0000
> Max. :1.0000 Max. :34.0000 Max. :1.0000 Max. :1.0000
> repstatus1 fecstatus1 matstatus1 alive1
> Min. :0.0000 Min. :0.00000 Min. :0.00 Min. :0.0000
> 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.00 1st Qu.:0.0000
> Median :0.0000 Median :0.00000 Median :0.00 Median :0.0000
> Mean :0.0203 Mean :0.01292 Mean :0.19 Mean :0.3524
> 3rd Qu.:0.0000 3rd Qu.:0.00000 3rd Qu.:0.00 3rd Qu.:1.0000
> Max. :1.0000 Max. :1.00000 Max. :1.00 Max. :1.0000
> stage1 stage1index sizea2 repstra2
> Length:542 Min. : 0.000 Min. :0.000 Min. :0.00000
> Class :character 1st Qu.: 0.000 1st Qu.:2.200 1st Qu.:0.00000
> Mode :character Median : 0.000 Median :2.500 Median :0.00000
> Mean : 1.893 Mean :3.059 Mean :0.03506
> 3rd Qu.: 2.000 3rd Qu.:3.600 3rd Qu.:0.00000
> Max. :20.000 Max. :8.400 Max. :1.00000
> feca2 juvgiven2 obsstatus2 repstatus2
> Min. : 0.0000 Min. :0.0000 Min. :0.0000 Min. :0.00000
> 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.:1.0000 1st Qu.:0.00000
> Median : 0.0000 Median :0.0000 Median :1.0000 Median :0.00000
> Mean : 0.1421 Mean :0.3469 Mean :0.9815 Mean :0.03506
> 3rd Qu.: 0.0000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:0.00000
> Max. :34.0000 Max. :1.0000 Max. :1.0000 Max. :1.00000
> fecstatus2 matstatus2 alive2 stage2
> Min. :0.00000 Min. :0.0000 Min. :1 Length:542
> 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:1 Class :character
> Median :0.00000 Median :1.0000 Median :1 Mode :character
> Mean :0.01661 Mean :0.6531 Mean :1
> 3rd Qu.:0.00000 3rd Qu.:1.0000 3rd Qu.:1
> Max. :1.00000 Max. :1.0000 Max. :1
> stage2index sizea3 repstra3 feca3
> Min. : 2.00 Min. :0.000 Min. :0.00000 Min. : 0.0000
> 1st Qu.: 2.00 1st Qu.:0.000 1st Qu.:0.00000 1st Qu.: 0.0000
> Median : 5.00 Median :2.400 Median :0.00000 Median : 0.0000
> Mean : 5.19 Mean :2.303 Mean :0.04059 Mean : 0.2454
> 3rd Qu.: 7.00 3rd Qu.:3.200 3rd Qu.:0.00000 3rd Qu.: 0.0000
> Max. :20.00 Max. :8.800 Max. :1.00000 Max. :48.0000
> juvgiven3 obsstatus3 repstatus3 fecstatus3 matstatus3
> Min. :0 Min. :0.0000 Min. :0.00000 Min. :0.00000 Min. :1
> 1st Qu.:0 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:1
> Median :0 Median :1.0000 Median :0.00000 Median :0.00000 Median :1
> Mean :0 Mean :0.7306 Mean :0.04059 Mean :0.01292 Mean :1
> 3rd Qu.:0 3rd Qu.:1.0000 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:1
> Max. :0 Max. :1.0000 Max. :1.00000 Max. :1.00000 Max. :1
> alive3 stage3 stage3index
> Min. :0.0000 Length:542 Min. : 0.000
> 1st Qu.:1.0000 Class :character 1st Qu.: 3.000
> Median :1.0000 Mode :character Median : 5.000
> Mean :0.7989 Mean : 4.969
> 3rd Qu.:1.0000 3rd Qu.: 6.000
> Max. :1.0000 Max. :21.000We are now down to only 542 rows, so we have lost close to 80% of the transition data here. That is a steep price to pay for accurate inference, but it is necessary in this case.
Let’s look at fecundity. Fecundity is integer-based, suggesting that it can be treated as a count variable. This package currently allows eight choices of count distributions: Gaussian, gamma, Poisson, negative binomial, zero-inflated Poisson, zero-inflated negative binomial, zero-truncated Poisson, and zero-truncated negative binomial. To assess which to use, we should first assess whether the mean and variance of the count are equal using a dispersion test. This test allows us to test whether the variance is greater than that expected under our mean value for fecundity using a chi-squared test. If it is not significantly different, then we may use some variant of the Poisson distribution. If the data are overdispersed, then we should use the negative binomial distribution. We should also test whether the number of zeroes is more than expected under these distributions, and make the distribution zero-inflated if so. Note that, because we have not excluded 0s from fecundity using reproductive status, we should not use a zero-truncated distribution.
Let’s start off by looking at a bar plot of the distribution of fecundity (figure 6.3).
hist(subset(lathvertln_small, repstatus2 == 1)$feca2, main = "Fecundity",
xlab = "Intact seeds produced in occasion t")
Figure 6.3: Histogram of fecundity in occasion t
We see that the distribution conforms to a classic count variable with a very low mean value. The first bar suggests that there may be too many zeroes, as well. Let’s now formally test our assumptions, that the mean equals the variance and that the number of zeroes meets expectation. Both tests use chi-squared distribution-based approaches, with the zero-inflation test in particular based on Broek (1995).
sf_distrib(lathvertln_small, sizea = c("sizea3", "sizea2"),
fec = c("feca3", "feca2"), obs3 = "obsstatus3",
repst = c("repstatus3", "repstatus2"))
> Non-integer values detected, so will not test for overdispersion and zero-inflation in sizea
> Mean fec is 4.6
>
> The variance in fec is 73.54
>
> The probability of this dispersion level by chance assuming that
> the true mean fec = variance in fec,
> and an alternative hypothesis of overdispersion, is 3.787e-145
>
> Fecundity is significantly overdispersed.
>
>
> Mean lambda in fec is 0.01005
> The actual number of 0s in fec is 7
> The expected number of 0s in fec under the null hypothesis is 0.1508
> The probability of this deviation in 0s from expectation by chance is 1.122e-73
>
> Fecundity is significantly zero-inflated.We see that fecundity is significantly overdispersed, and has a significant excess of zeros. So, we should use a zero-inflated negative binomial distribution here.
Next, let’s create an ahistorical supplement table organizing the extra data that we need to incorporate into our matrices. Each row refers to a specific transition. The first two of these transitions are set to specific probabilities, which are the probabilities of germination and seed dormancy, estimated from a separate study. The final two terms are fecundity multipliers, which mark which transitions correspond to fecundity and provide information on what multiple of fecundity estimated via linear modeling applies to each case. Note that we can also include proxy transitions, in which we define a specific transition as being equal to another in the matrix. The latter approach is useful when some transitions cannot be estimated given a particular dataset, and so need to be set to other, proxy values that are estimable.
lathsupp2 <- supplemental(stage3 = c("Sd", "Sdl", "Sd", "Sdl"),
stage2 = c("Sd", "Sd", "rep", "rep"), givenrate = c(0.345, 0.054, NA, NA),
multiplier = c(NA, NA, 0.345, 0.054), type = c(1, 1, 3, 3),
stageframe = lathframeln, historical = FALSE)
lathsupp2
> stage3 stage2 stage1 eststage3 eststage2 eststage1 givenrate multiplier
> 1 Sd Sd <NA> <NA> <NA> <NA> 0.345 1.000
> 2 Sdl Sd <NA> <NA> <NA> <NA> 0.054 1.000
> 3 Sd rep <NA> <NA> <NA> <NA> NA 0.345
> 4 Sdl rep <NA> <NA> <NA> <NA> NA 0.054
> convtype convtype_t12
> 1 1 1
> 2 1 1
> 3 3 1
> 4 3 1
Next we will run the modelsearch function with the new vertical dataset. This function will develop our best-fit vital rate models for us. This function looks simple, but it automates several crucial and complex tasks in MPM analysis. Specifically, it automates 1) the building of global models for each vital rate requested, 2) the exhaustive construction of all reduced models, and 3) the selection of the best-fit models. In relation to our previous uses of this function in chapter 5, the most noteworthy difference is the inclusion of an age term (age = "obsage", which we know from looking at the summary of the vertical dataset). We will only develop the full effects models here, which include main effects and all two-way interactions between age and other terms, including size and reproductive status. Let’s start off by generating a set of vital rate models that covers the entire population.
lathmodelsln2 <- modelsearch(lathvertln_small, historical = FALSE,
approach = "mixed", suite = "full", bestfit = "AICc&k", juvestimate = "Sdl",
vitalrates = c("surv", "obs", "size", "repst", "fec"), sizedist = "gaussian",
fecdist = "negbin", indiv = "individ", year = "year2", age = "obsage",
year.as.random = TRUE, patch.as.random = TRUE, show.model.tables = TRUE,
fec.zero = TRUE, quiet = TRUE)Let’s see a summary of the lefkoMod object that we have created.
summary(lathmodelsln2)
> This LefkoMod object includes 8 linear models.
> Best-fit model criterion used: aicc&k
>
>
>
> Survival model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: alive3 ~ obsage + repstatus2 + (1 | year2) + (1 | individ)
> Data: subdata
> AIC BIC logLik deviance df.resid
> 314.3559 333.7024 -152.1779 304.3559 349
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 0
> year2 (Intercept) 0
> Number of obs: 354, groups: individ, 190; year2, 2
> Fixed Effects:
> (Intercept) obsage repstatus2
> 2.5636 -0.6252 29.6185
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
>
>
>
> Observation model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: obsstatus3 ~ obsage + sizea2 + (1 | year2) + (1 | individ) +
> obsage:sizea2
> Data: subdata
> AIC BIC logLik deviance df.resid
> 150.6174 172.7595 -69.3087 138.6174 290
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 1.741e+01
> year2 (Intercept) 1.408e-04
> Number of obs: 296, groups: individ, 162; year2, 2
> Fixed Effects:
> (Intercept) obsage sizea2 obsage:sizea2
> -9.892 15.058 6.489 -4.448
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 2 lme4 warnings
>
>
>
> Size model:
> Linear mixed model fit by REML ['lmerMod']
> Formula: sizea3 ~ sizea2 + (1 | year2) + (1 | individ)
> Data: subdata
> REML criterion at convergence: 700.2483
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 0.0000
> year2 (Intercept) 0.4826
> Residual 0.8850
> Number of obs: 266, groups: individ, 158; year2, 2
> Fixed Effects:
> (Intercept) sizea2
> 0.7088 0.7777
> optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
>
>
>
> Secondary size model:
> [1] 1
>
>
>
> Tertiary size model:
> [1] 1
>
>
>
> Reproductive status model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: repstatus3 ~ repstatus2 + sizea2 + (1 | year2) + (1 | individ)
> Data: subdata
> AIC BIC logLik deviance df.resid
> 74.2007 92.1181 -32.1003 64.2007 261
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 101.78
> year2 (Intercept) 52.35
> Number of obs: 266, groups: individ, 158; year2, 2
> Fixed Effects:
> (Intercept) repstatus2 sizea2
> -184.21 -13.28 19.93
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
>
>
>
> Fecundity model:
> Formula: feca2 ~ obsage + sizea2 + (1 | year2) + (1 | individ)
> Zero inflation: ~(1 | year2) + (1 | individ)
> Data: subdata
> AIC BIC logLik df.resid
> 80.13600 88.63595 -31.06800 10
> Random-effects (co)variances:
>
> Conditional model:
> Groups Name Std.Dev.
> year2 (Intercept) 2.487e-04
> individ (Intercept) 2.897e-05
>
> Zero-inflation model:
> Groups Name Std.Dev.
> year2 (Intercept) 0.001679
> individ (Intercept) 2.021180
>
> Number of obs: 19 / Conditional model: year2, 2; individ, 16 / Zero-inflation model: year2, 2; individ, 16
>
> Dispersion parameter for nbinom2 family (): 2.45e+08
>
> Fixed Effects:
>
> Conditional model:
> (Intercept) obsage sizea2
> -7.1322 0.4786 1.2036
>
> Zero-inflation model:
> (Intercept)
> -0.1414
>
>
> Juvenile survival model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: alive3 ~ (1 | year2) + (1 | individ)
> Data: subdata
> AIC BIC logLik deviance df.resid
> 225.7193 235.4286 -109.8596 219.7193 185
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 0.3755
> year2 (Intercept) 0.0000
> Number of obs: 188, groups: individ, 155; year2, 2
> Fixed Effects:
> (Intercept)
> 1.024
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
>
>
>
> Juvenile observation model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: obsstatus3 ~ (1 | year2) + (1 | individ)
> Data: subdata
> AIC BIC logLik deviance df.resid
> 28.1909 36.9509 -11.0955 22.1909 134
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 5.638e+01
> year2 (Intercept) 1.312e-05
> Number of obs: 137, groups: individ, 122; year2, 2
> Fixed Effects:
> (Intercept)
> 13.8
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
>
>
>
> Juvenile size model:
> Linear mixed model fit by REML ['lmerMod']
> Formula: sizea3 ~ (1 | year2) + (1 | individ)
> Data: subdata
> REML criterion at convergence: 145.4426
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 4.504e-06
> year2 (Intercept) 1.430e-01
> Residual 4.139e-01
> Number of obs: 130, groups: individ, 115; year2, 2
> Fixed Effects:
> (Intercept)
> 2.288
> optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
>
>
>
> Juvenile secondary size model:
> [1] 1
>
>
>
> Juvenile tertiary size model:
> [1] 1
>
>
>
> Juvenile reproduction model:
> [1] 0
>
>
>
> Juvenile maturity model:
> [1] 1
>
>
>
>
>
> Number of models in survival table: 18
>
> Number of models in observation table: 18
>
> Number of models in size table: 18
>
> Number of models in secondary size table: 1
>
> Number of models in tertiary size table: 1
>
> Number of models in reproduction status table: 17
>
> Number of models in fecundity table: 285
>
> Number of models in juvenile survival table: 1
>
> Number of models in juvenile observation table: 1
>
> Number of models in juvenile size table: 1
>
> Number of models in juvenile secondary size table: 1
>
> Number of models in juvenile tertiary size table: 1
>
> Number of models in juvenile reproduction table: 1
>
> Number of models in juvenile maturity table: 1
>
>
>
>
>
> General model parameter names (column 1), and
> specific names used in these models (column 2):
> parameter_names mainparams
> 1 time t year2
> 2 individual individ
> 3 patch patch
> 4 alive in time t+1 surv3
> 5 observed in time t+1 obs3
> 6 sizea in time t+1 size3
> 7 sizeb in time t+1 sizeb3
> 8 sizec in time t+1 sizec3
> 9 reproductive status in time t+1 repst3
> 10 fecundity in time t+1 fec3
> 11 fecundity in time t fec2
> 12 sizea in time t size2
> 13 sizea in time t-1 size1
> 14 sizeb in time t sizeb2
> 15 sizeb in time t-1 sizeb1
> 16 sizec in time t sizec2
> 17 sizec in time t-1 sizec1
> 18 reproductive status in time t repst2
> 19 reproductive status in time t-1 repst1
> 20 maturity status in time t+1 matst3
> 21 maturity status in time t matst2
> 22 age in time t age
> 23 density in time t density
> 24 individual covariate a in time t indcova2
> 25 individual covariate a in time t-1 indcova1
> 26 individual covariate b in time t indcovb2
> 27 individual covariate b in time t-1 indcovb1
> 28 individual covariate c in time t indcovc2
> 29 individual covariate c in time t-1 indcovc1
> 30 stage group in time t group2
> 31 stage group in time t-1 group1
>
>
>
>
>
> Quality control:
>
> Survival estimated with 190 individuals and 354 individual transitions.
> Survival accuracy is 0.836.
> Observation estimated with 162 individuals and 296 individual transitions.
> Observation accuracy is 0.973.
> Primary size estimated with 158 individuals and 266 individual transitions.
> Primary size pseudo R-squared is 0.692.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 158 individuals and 266 individual transitions.
> Reproductive status accuracy is 1.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Fecundity pseudo R-squared is 0.863.
> Juvenile survival estimated with 155 individuals and 188 individual transitions.
> Juvenile survival accuracy is 0.729.
> Juvenile observation estimated with 122 individuals and 137 individual transitions.
> Juvenile observation accuracy is 1.
> Juvenile primary size estimated with 115 individuals and 130 individual transitions.
> Juvenile primary size pseudo R-squared is 0.107.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity probability not estimated.We see that age is influential in survival, observation, and fecundity, though not in other vital rates. Accuracy in our binomial models is high though not fabulous.
Next, we will create a second set that includes patch as a random factor. This model set will allow us to explore patch dynamics in addition to the population dynamics of the previous set.
lathmodelsln2p <- modelsearch(lathvertln_small, historical = FALSE,
approach = "mixed", suite = "full", bestfit = "AICc&k", juvestimate = "Sdl",
vitalrates = c("surv", "obs", "size", "repst", "fec"), sizedist = "gaussian",
fecdist = "negbin", indiv = "individ", patch = "patchid", year = "year2",
age = "obsage", year.as.random = TRUE, patch.as.random = TRUE,
show.model.tables = TRUE, fec.zero = TRUE, quiet = TRUE)And a summary.
summary(lathmodelsln2p)
> This LefkoMod object includes 8 linear models.
> Best-fit model criterion used: aicc&k
>
>
>
> Survival model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: alive3 ~ obsage + repstatus2 + (1 | year2) + (1 | patchid) +
> (1 | individ)
> Data: subdata
> AIC BIC logLik deviance df.resid
> 307.2136 330.4294 -147.6068 295.2136 348
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 7.285e-07
> patchid (Intercept) 5.828e-01
> year2 (Intercept) 0.000e+00
> Number of obs: 354, groups: individ, 190; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept) obsage repstatus2
> 2.7389 -0.6484 21.7793
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
>
>
>
> Observation model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: obsstatus3 ~ obsage + sizea2 + (1 | year2) + (1 | patchid) +
> (1 | individ) + obsage:sizea2
> Data: subdata
> AIC BIC logLik deviance df.resid
> 153.0167 178.8493 -69.5084 139.0167 289
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 17.8608
> patchid (Intercept) 0.1586
> year2 (Intercept) 0.5055
> Number of obs: 296, groups: individ, 162; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept) obsage sizea2 obsage:sizea2
> -12.146 16.752 7.071 -4.800
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 2 lme4 warnings
>
>
>
> Size model:
> Linear mixed model fit by REML ['lmerMod']
> Formula: sizea3 ~ sizea2 + (1 | year2) + (1 | patchid) + (1 | individ)
> Data: subdata
> REML criterion at convergence: 700.2483
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 0.0000
> patchid (Intercept) 0.0000
> year2 (Intercept) 0.4826
> Residual 0.8850
> Number of obs: 266, groups: individ, 158; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept) sizea2
> 0.7088 0.7777
> optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
>
>
>
> Secondary size model:
> [1] 1
>
>
>
> Tertiary size model:
> [1] 1
>
>
>
> Reproductive status model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: repstatus3 ~ obsage + sizea2 + (1 | year2) + (1 | patchid) +
> (1 | individ)
> Data: subdata
> AIC BIC logLik deviance df.resid
> 75.5686 97.0696 -31.7843 63.5686 260
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 105.936
> patchid (Intercept) 7.812
> year2 (Intercept) 1.106
> Number of obs: 266, groups: individ, 158; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept) obsage sizea2
> -236.07 66.23 18.97
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 2 lme4 warnings
>
>
>
> Fecundity model:
> Formula: feca2 ~ sizea2 + (1 | year2) + (1 | patchid) + (1 | individ)
> Zero inflation: ~(1 | year2) + (1 | patchid) + (1 | individ)
> Data: subdata
> AIC BIC logLik df.resid
> 83.29320 92.73759 -31.64660 9
> Random-effects (co)variances:
>
> Conditional model:
> Groups Name Std.Dev.
> year2 (Intercept) 1.936e-04
> patchid (Intercept) 6.955e-05
> individ (Intercept) 2.645e-04
>
> Zero-inflation model:
> Groups Name Std.Dev.
> year2 (Intercept) 6.761e-03
> patchid (Intercept) 8.938e-06
> individ (Intercept) 1.777e+00
>
> Number of obs: 19 / Conditional model: year2, 2; patchid, 6; individ, 16 / Zero-inflation model: year2, 2; patchid, 6; individ, 16
>
> Dispersion parameter for nbinom2 family (): 2.61e+09
>
> Fixed Effects:
>
> Conditional model:
> (Intercept) sizea2
> -6.265 1.157
>
> Zero-inflation model:
> (Intercept)
> -0.1917
>
>
> Juvenile survival model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: alive3 ~ (1 | year2) + (1 | patchid) + (1 | individ)
> Data: subdata
> AIC BIC logLik deviance df.resid
> 227.5465 240.4922 -109.7732 219.5465 184
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 0.3386
> patchid (Intercept) 0.2282
> year2 (Intercept) 0.0000
> Number of obs: 188, groups: individ, 155; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept)
> 1.015
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
>
>
>
> Juvenile observation model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: obsstatus3 ~ (1 | year2) + (1 | patchid) + (1 | individ)
> Data: subdata
> AIC BIC logLik deviance df.resid
> 30.1909 41.8709 -11.0955 22.1909 133
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 5.638e+01
> patchid (Intercept) 0.000e+00
> year2 (Intercept) 1.494e-04
> Number of obs: 137, groups: individ, 122; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept)
> 13.8
> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
>
>
>
> Juvenile size model:
> Linear mixed model fit by REML ['lmerMod']
> Formula: sizea3 ~ (1 | year2) + (1 | patchid) + (1 | individ)
> Data: subdata
> REML criterion at convergence: 145.4426
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 0.0000
> patchid (Intercept) 0.0000
> year2 (Intercept) 0.1430
> Residual 0.4139
> Number of obs: 130, groups: individ, 115; patchid, 6; year2, 2
> Fixed Effects:
> (Intercept)
> 2.288
> optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings
>
>
>
> Juvenile secondary size model:
> [1] 1
>
>
>
> Juvenile tertiary size model:
> [1] 1
>
>
>
> Juvenile reproduction model:
> [1] 0
>
>
>
> Juvenile maturity model:
> [1] 1
>
>
>
>
>
> Number of models in survival table: 18
>
> Number of models in observation table: 18
>
> Number of models in size table: 18
>
> Number of models in secondary size table: 1
>
> Number of models in tertiary size table: 1
>
> Number of models in reproduction status table: 7
>
> Number of models in fecundity table: 265
>
> Number of models in juvenile survival table: 1
>
> Number of models in juvenile observation table: 1
>
> Number of models in juvenile size table: 1
>
> Number of models in juvenile secondary size table: 1
>
> Number of models in juvenile tertiary size table: 1
>
> Number of models in juvenile reproduction table: 1
>
> Number of models in juvenile maturity table: 1
>
>
>
>
>
> General model parameter names (column 1), and
> specific names used in these models (column 2):
> parameter_names mainparams
> 1 time t year2
> 2 individual individ
> 3 patch patch
> 4 alive in time t+1 surv3
> 5 observed in time t+1 obs3
> 6 sizea in time t+1 size3
> 7 sizeb in time t+1 sizeb3
> 8 sizec in time t+1 sizec3
> 9 reproductive status in time t+1 repst3
> 10 fecundity in time t+1 fec3
> 11 fecundity in time t fec2
> 12 sizea in time t size2
> 13 sizea in time t-1 size1
> 14 sizeb in time t sizeb2
> 15 sizeb in time t-1 sizeb1
> 16 sizec in time t sizec2
> 17 sizec in time t-1 sizec1
> 18 reproductive status in time t repst2
> 19 reproductive status in time t-1 repst1
> 20 maturity status in time t+1 matst3
> 21 maturity status in time t matst2
> 22 age in time t age
> 23 density in time t density
> 24 individual covariate a in time t indcova2
> 25 individual covariate a in time t-1 indcova1
> 26 individual covariate b in time t indcovb2
> 27 individual covariate b in time t-1 indcovb1
> 28 individual covariate c in time t indcovc2
> 29 individual covariate c in time t-1 indcovc1
> 30 stage group in time t group2
> 31 stage group in time t-1 group1
>
>
>
>
>
> Quality control:
>
> Survival estimated with 190 individuals and 354 individual transitions.
> Survival accuracy is 0.836.
> Observation estimated with 162 individuals and 296 individual transitions.
> Observation accuracy is 0.973.
> Primary size estimated with 158 individuals and 266 individual transitions.
> Primary size pseudo R-squared is 0.692.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 158 individuals and 266 individual transitions.
> Reproductive status accuracy is 1.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Fecundity pseudo R-squared is 0.877.
> Juvenile survival estimated with 155 individuals and 188 individual transitions.
> Juvenile survival accuracy is 0.729.
> Juvenile observation estimated with 122 individuals and 137 individual transitions.
> Juvenile observation accuracy is 1.
> Juvenile primary size estimated with 115 individuals and 130 individual transitions.
> Juvenile primary size pseudo R-squared is 0.107.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity probability not estimated.Incorporating patch has led to some changes in the models, most notably that now reproductive status is a function of age while fecundity is not. Interesting patterns worth exploring….
Next, we will estimate the ahistorical sets of matrices. We will match the ahistorical age-by-stage matrix estimation function, aflefko2(), with the appropriate ahistorical input, including the ahistorical lefkoMod objects lathmodelsln2 and lathmodelsln2p. Model sets that include historical terms should not be used to create ahistorical matrices, since the coefficients in the best-fit models are estimated assuming a specific model structure that either relies on historical terms or does not. Historical vital rate models may yield biased results if used to construct ahistorical matrices. Also note that lefko3 does not currently allow the construction of historical age-by-stage MPMs. Let’s start off by developing the population-only MPM.
lathmat2age <- aflefko2(year = "all", stageframe = lathframeln,
modelsuite = lathmodelsln2, data = lathvertln_small, supplement = lathsupp2,
continue = TRUE, reduce = FALSE)
summary(lathmat2age)
>
> This ahistorical lefkoMat object contains 2 matrices.
>
> Each matrix is square with 42 rows and columns, and a total of 1764 elements.
> A total of 682 survival transitions were estimated, with 341 per matrix.
> A total of 36 fecundity transitions were estimated, with 18 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 2 time steps.
>
> Vital rate modeling quality control:
>
> Survival estimated with 190 individuals and 354 individual transitions.
> Observation estimated with 162 individuals and 296 individual transitions.
> Primary size estimated with 158 individuals and 266 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 158 individuals and 266 individual transitions.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival estimated with 155 individuals and 188 individual transitions.
> Juvenile observation estimated with 122 individuals and 137 individual transitions.
> Juvenile primary size estimated with 115 individuals and 130 individual transitions.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
>
> Survival probability sum check (each matrix represented by column in order):
> [,1] [,2]
> Min. 0.000 0.000
> 1st Qu. 0.000 0.000
> Median 0.172 0.172
> Mean 0.395 0.409
> 3rd Qu. 0.788 0.788
> Max. 1.000 1.000
This first model set led to the development of two matrices, because although there are four years of data, we have limited our use of the data to individuals first seen from the second year on. Function aflefko2() has assessed patterns in the data set and has found that the maximum age is 2 years (it would have been 3 years if we allowed the first year). We have informed the function that this age is not terminal and that the demography of that age should continue onward (continue = TRUE). This has resulted in a block matrix with two ages and 21 stages, and so 42 age-stage combinations and 42 rows and columns. Of course, this particular plant species is actually long-lived, and so there may very well be further vital rate variability in across the lifespan of the plant. However, as our dataset only includes four years of study and we do not have absolute ages for any plant, we can only incur four years of relative age at best. Thus, our settings are actually the most parsimonious possible under the circumstances.
The quality control section gives us a sense of the amount of data used to model each vital rate, and also shows us that the survival-transition (U) matrices are composed entirely of proper probabilities yielding stage survival probability falling between 0.0 and 1.0. Matrix estimation can sometimes create spurious values, such as stage survival greater than 1.0. Such values can occur for a variety of reasons, but the most common is the inclusion through a supplemental table or overwrite table of externally-determined survival probabilities that are too high. Make sure to check your matrix column sums each time you estimate MPMs to prevent this problem. Survival greater than 1.0 can lead to strange effects on metrics of population dynamics.
Let’s now develop the patch-level MPMs.
lathmat2agep <- aflefko2(year = "all", patch = "all", stageframe = lathframeln,
modelsuite = lathmodelsln2p, data = lathvertln_small, supplement = lathsupp2,
continue = TRUE, reduce = FALSE)
summary(lathmat2agep)
>
> This ahistorical lefkoMat object contains 12 matrices.
>
> Each matrix is square with 42 rows and columns, and a total of 1764 elements.
> A total of 3876 survival transitions were estimated, with 323 per matrix.
> A total of 216 fecundity transitions were estimated, with 18 per matrix.
> This lefkoMat object covers 1 population, 6 patches, and 2 time steps.
>
> Vital rate modeling quality control:
>
> Survival estimated with 190 individuals and 354 individual transitions.
> Observation estimated with 162 individuals and 296 individual transitions.
> Primary size estimated with 158 individuals and 266 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 158 individuals and 266 individual transitions.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival estimated with 155 individuals and 188 individual transitions.
> Juvenile observation estimated with 122 individuals and 137 individual transitions.
> Juvenile primary size estimated with 115 individuals and 130 individual transitions.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
>
> Survival probability sum check (each matrix represented by column in order):
> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
> Min. 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
> 1st Qu. 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
> Median 0.172 0.172 0.168 0.172 0.172 0.172 0.148 0.172 0.166 0.172 0.172 0.172
> Mean 0.415 0.430 0.389 0.402 0.412 0.427 0.369 0.382 0.387 0.401 0.409 0.424
> 3rd Qu. 0.879 0.879 0.757 0.757 0.865 0.865 0.667 0.667 0.750 0.750 0.851 0.851
> Max. 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000This second model set led to the development of 12 matrices, reflecting our subset number of three years and six patches. The rest of the output looks quite similar and even the survival-transition matrix column sum summaries look extremely similar, suggesting little impact of patch.
We can get a sense of what these matrices look like by visualizing them. Let’s use the image3() function to look at just one (figure 6.4).
image3(lathmat2age, used = 1)
Figure 6.4: Visualization of 1st A Matrix. Red area corresponds to non-zero elements
> [[1]]
> NULLThe clear squares refer to zero elements, and the red elements refer to non-zero values corresponding to survival transitions and fecundity. The vast number of zeros may be surprising, but this matrix is a supermatrix organized by age first, with stage organizing within-age blocks. The first age is age 1, which cannot be adult, and so we find zeros in adult stages at age 1. The adult block occurs from age 2, and this block can perpetuate indefinitely. The number of elements estimated is greater now than in the typical ahistorical MPM, because now we have added age as a major factor for analysis. This matrix is overwhelmingly composed of elements that must be zero, and so it is a rather sparse matrix ((341 + 18) / 1764 = 20.4% of elements).
Given the amount of white space, we might prefer to remove impossible age-stage combinations and reduce the size of the matrices. We can do this by recreating our matrices with reduce = TRUE. Let’s try that with the main population set.
lathmat2age_red <- aflefko2(year = "all", stageframe = lathframeln,
modelsuite = lathmodelsln2, data = lathvertln_small, supplement = lathsupp2,
continue = TRUE, reduce = TRUE)
summary(lathmat2age_red)
>
> This ahistorical lefkoMat object contains 2 matrices.
>
> Each matrix is square with 22 rows and columns, and a total of 484 elements.
> A total of 682 survival transitions were estimated, with 341 per matrix.
> A total of 36 fecundity transitions were estimated, with 18 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 2 time steps.
>
> Vital rate modeling quality control:
>
> Survival estimated with 190 individuals and 354 individual transitions.
> Observation estimated with 162 individuals and 296 individual transitions.
> Primary size estimated with 158 individuals and 266 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 158 individuals and 266 individual transitions.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival estimated with 155 individuals and 188 individual transitions.
> Juvenile observation estimated with 122 individuals and 137 individual transitions.
> Juvenile primary size estimated with 115 individuals and 130 individual transitions.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
>
> Survival probability sum check (each matrix represented by column in order):
> [,1] [,2]
> Min. 0.000 0.000
> 1st Qu. 0.752 0.781
> Median 0.788 0.788
> Mean 0.755 0.782
> 3rd Qu. 0.997 0.995
> Max. 1.000 1.000This exercise has eliminated 20 rows and columns, yielding a matrix with 22 rows and 22 columns. The total number of estimated elements has not changed, meaning that our matrices are now much denser ((341 + 18) / 484 = 74.2%). Let’s take a look at an image of the first matrix (figure 6.5).
image3(lathmat2age_red, used = 1)
Figure 6.5: Visualization of 1st A Matrix. Red area corresponds to non-zero elements
> [[1]]
> NULLSome of the white space has been reduced, and we see greater coverage of the matrix by non-zero elements.
We can see the order of ages and stages using the agestages element of the lefkoMat object we produced, as below. Note that our matrix is 22 rows by 22 columns, and this object gives us the exact order used.
lathmat2age_red$agestages
> stage_id stage age
> 1 1 Sd 1
> 2 2 Sdl 1
> 1.1 1 Sd 2
> 3.1 3 Dorm 2
> 4.1 4 Sz1nr 2
> 5.1 5 Sz2nr 2
> 6.1 6 Sz3nr 2
> 7.1 7 Sz4nr 2
> 8.1 8 Sz5nr 2
> 9.1 9 Sz6nr 2
> 10.1 10 Sz7nr 2
> 11.1 11 Sz8nr 2
> 12.1 12 Sz9nr 2
> 13.1 13 Sz1r 2
> 14.1 14 Sz2r 2
> 15.1 15 Sz3r 2
> 16.1 16 Sz4r 2
> 17.1 17 Sz5r 2
> 18.1 18 Sz6r 2
> 19.1 19 Sz7r 2
> 20.1 20 Sz8r 2
> 21.1 21 Sz9r 2
Now let’s estimate the element-wise arithmetic mean matrices. The first lefkoMat object created will include a single mean matrix, while the second will include six mean matrices for the patches followed by a grand mean matrix, yielding a total of seven matrices.
lathmat2mean <- lmean(lathmat2age)
lathmat2pmean <- lmean(lathmat2agep)
summary(lathmat2mean)
>
> This ahistorical lefkoMat object contains 1 matrix.
>
> Each matrix is square with 42 rows and columns, and a total of 1764 elements.
> A total of 361 survival transitions were estimated, with 361 per matrix.
> A total of 18 fecundity transitions were estimated, with 18 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 0 time steps.
>
> Vital rate modeling quality control:
>
> Survival estimated with 190 individuals and 354 individual transitions.
> Observation estimated with 162 individuals and 296 individual transitions.
> Primary size estimated with 158 individuals and 266 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 158 individuals and 266 individual transitions.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival estimated with 155 individuals and 188 individual transitions.
> Juvenile observation estimated with 122 individuals and 137 individual transitions.
> Juvenile primary size estimated with 115 individuals and 130 individual transitions.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
>
> Survival probability sum check (each matrix represented by column in order):
> [,1]
> Min. 0.000
> 1st Qu. 0.000
> Median 0.172
> Mean 0.402
> 3rd Qu. 0.788
> Max. 1.000
summary(lathmat2pmean)
>
> This ahistorical lefkoMat object contains 7 matrices.
>
> Each matrix is square with 42 rows and columns, and a total of 1764 elements.
> A total of 2275 survival transitions were estimated, with 325 per matrix.
> A total of 126 fecundity transitions were estimated, with 18 per matrix.
> This lefkoMat object covers 1 population, 7 patches, and 0 time steps.
>
> Vital rate modeling quality control:
>
> Survival estimated with 190 individuals and 354 individual transitions.
> Observation estimated with 162 individuals and 296 individual transitions.
> Primary size estimated with 158 individuals and 266 individual transitions.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproductive status estimated with 158 individuals and 266 individual transitions.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival estimated with 155 individuals and 188 individual transitions.
> Juvenile observation estimated with 122 individuals and 137 individual transitions.
> Juvenile primary size estimated with 115 individuals and 130 individual transitions.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
>
> Survival probability sum check (each matrix represented by column in order):
> [,1] [,2] [,3] [,4] [,5] [,6] [,7]
> Min. 0.000 0.000 0.000 0.000 0.000 0.000 0.000
> 1st Qu. 0.000 0.000 0.000 0.000 0.000 0.000 0.000
> Median 0.172 0.172 0.172 0.172 0.172 0.172 0.172
> Mean 0.423 0.395 0.420 0.375 0.394 0.417 0.404
> 3rd Qu. 0.879 0.757 0.865 0.667 0.750 0.851 0.795
> Max. 1.000 1.000 1.000 1.000 1.000 1.000 1.000We see one overall population mean in the first case, and a set of six patch-level means and one population mean in the second case. Note that the population mean in each case should be a little different, because the population mean in the second set weights each patch equally, while the overall population mean in the first case weights individual transitions equally regardless of patch of origin.
6.3 Developing raw (empirical) age-by-stage MPMs in lefko3
Now let’s focus our attention on raw (empirical) age-by-stage MPMs. For this, we will use the function arlefko2(). Raw MPMs generally require a great deal more data to properly parameterize. We will also need a stageframe with fewer stages, in order to limit the number of artificial zeros entering our matrices. Let’s start with the following life history model, which assumes a grand total of seven stages.
Figure 6.6: Life history model of Lathyrus vernus using log leaf volume as the size classification metric
Let’s now look at our stageframe for this model. The actual stage definitions in terms of size bins are no longer based on the log leaf volume. Instead, they are based on the raw leaf volume. In order to deal with the fact that small plants are far more common than large plants, the bins have increasing widths with increasing leaf size. We will also lump all reproduction into a single flowering class that can include individuals of any size, provided that they have sprouted.
sizevector <- c(0, 100, 13, 127, 3730, 3800, 0)
stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm")
repvector <- c(0, 0, 0, 0, 0, 1, 0)
obsvector <- c(0, 1, 1, 1, 1, 1, 0)
matvector <- c(0, 0, 1, 1, 1, 1, 1)
immvector <- c(1, 1, 0, 0, 0, 0, 0)
propvector <- c(1, 0, 0, 0, 0, 0, 0)
indataset <- c(0, 1, 1, 1, 1, 1, 1)
binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5)
lathframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector,
repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
immstatus = immvector, indataset = indataset, binhalfwidth = binvec,
propstatus = propvector)
lathframe_raw
> stage size size_b size_c min_age max_age repstatus obsstatus propstatus
> 1 Sd 0 NA NA NA NA 0 0 1
> 2 Sdl 100 NA NA NA NA 0 1 0
> 3 VSm 13 NA NA NA NA 0 1 0
> 4 Sm 127 NA NA NA NA 0 1 0
> 5 VLa 3730 NA NA NA NA 0 1 0
> 6 Flo 3800 NA NA NA NA 1 1 0
> 7 Dorm 0 NA NA NA NA 0 0 0
> immstatus matstatus indataset binhalfwidth_raw sizebin_min sizebin_max
> 1 1 0 0 0.0 0.0 0.0
> 2 1 0 1 100.0 0.0 200.0
> 3 0 1 1 11.0 2.0 24.0
> 4 0 1 1 103.0 24.0 230.0
> 5 0 1 1 3500.0 230.0 7230.0
> 6 0 1 1 3800.0 0.0 7600.0
> 7 0 1 1 0.5 -0.5 0.5
> sizebin_center sizebin_width binhalfwidthb_raw sizebinb_min sizebinb_max
> 1 0 0 NA NA NA
> 2 100 200 NA NA NA
> 3 13 22 NA NA NA
> 4 127 206 NA NA NA
> 5 3730 7000 NA NA NA
> 6 3800 7600 NA NA NA
> 7 0 1 NA NA NA
> sizebinb_center sizebinb_width binhalfwidthc_raw sizebinc_min sizebinc_max
> 1 NA NA NA NA NA
> 2 NA NA NA NA NA
> 3 NA NA NA NA NA
> 4 NA NA NA NA NA
> 5 NA NA NA NA NA
> 6 NA NA NA NA NA
> 7 NA NA NA NA NA
> sizebinc_center sizebinc_width group comments
> 1 NA NA 0 No description
> 2 NA NA 0 No description
> 3 NA NA 0 No description
> 4 NA NA 0 No description
> 5 NA NA 0 No description
> 6 NA NA 0 No description
> 7 NA NA 0 No descriptionLet’s now restandardize our core dataset, using the new stageframe for stage assignment and the original size variable for size classification. We will also cut sightings of individuals in the first year, as before, to result in a dataset in which individuals are ages that we can be reasonably sure about.
lathvert_raw <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88",
fecacol = "Intactseed88", deadacol = "Dead1988", nonobsacol = "Dormant1988",
stageassign = lathframe_raw, stagesize = "sizea", censorcol = "Missing1988",
censorkeep = NA, censor = TRUE)
lathvert_raw_small <- subset(lathvert_raw, firstseen > 1988)
dim(lathvert_raw_small)
> [1] 542 45
summary(lathvert_raw_small)
> rowid popid patchid individ year2
> Min. : 30.0 :542 1:131 Length:542 Min. :1989
> 1st Qu.: 302.5 2: 72 Class :character 1st Qu.:1989
> Median : 583.0 3:101 Mode :character Median :1990
> Mean : 586.8 4:132 Mean :1990
> 3rd Qu.: 801.0 5: 70 3rd Qu.:1990
> Max. :1097.0 6: 36 Max. :1990
>
> firstseen lastseen obsage obslifespan sizea1
> Min. :1989 Min. :1989 Min. :1.000 Min. :0.00 Min. : 0.00
> 1st Qu.:1989 1st Qu.:1990 1st Qu.:1.000 1st Qu.:1.00 1st Qu.: 0.00
> Median :1989 Median :1991 Median :1.000 Median :2.00 Median : 0.00
> Mean :1989 Mean :1991 Mean :1.352 Mean :1.36 Mean : 59.41
> 3rd Qu.:1989 3rd Qu.:1991 3rd Qu.:2.000 3rd Qu.:2.00 3rd Qu.: 9.00
> Max. :1990 Max. :1991 Max. :2.000 Max. :2.00 Max. :4394.20
>
> repstra1 feca1 fec1added juvgiven1
> Min. :0.0000 Min. : 0 Min. : 0.0000 Min. :0.0000
> 1st Qu.:0.0000 1st Qu.: 0 1st Qu.: 0.0000 1st Qu.:0.0000
> Median :0.0000 Median : 4 Median : 0.0000 Median :0.0000
> Mean :0.0576 Mean : 6 Mean : 0.1218 Mean :0.1624
> 3rd Qu.:0.0000 3rd Qu.: 6 3rd Qu.: 0.0000 3rd Qu.:0.0000
> Max. :1.0000 Max. :34 Max. :34.0000 Max. :1.0000
> NA's :351 NA's :531
> obsstatus1 repstatus1 fecstatus1 matstatus1
> Min. :0.0000 Min. :0.0000 Min. :0.00000 Min. :0.00
> 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.00
> Median :0.0000 Median :0.0000 Median :0.00000 Median :0.00
> Mean :0.3524 Mean :0.0203 Mean :0.01292 Mean :0.19
> 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:0.00000 3rd Qu.:0.00
> Max. :1.0000 Max. :1.0000 Max. :1.00000 Max. :1.00
>
> alive1 stage1 stage1index sizea2
> Min. :0.0000 Length:542 Min. :0.0 Min. : 0.00
> 1st Qu.:0.0000 Class :character 1st Qu.:0.0 1st Qu.: 9.00
> Median :0.0000 Mode :character Median :0.0 Median : 12.60
> Mean :0.3524 Mean :1.1 Mean : 97.24
> 3rd Qu.:1.0000 3rd Qu.:2.0 3rd Qu.: 37.40
> Max. :1.0000 Max. :6.0 Max. :4394.20
>
> repstra2 feca2 fec2added juvgiven2
> Min. :0.00000 Min. : 0.000 Min. : 0.0000 Min. :0.0000
> 1st Qu.:0.00000 1st Qu.: 0.000 1st Qu.: 0.0000 1st Qu.:0.0000
> Median :0.00000 Median : 0.000 Median : 0.0000 Median :0.0000
> Mean :0.03571 Mean : 4.053 Mean : 0.1421 Mean :0.3469
> 3rd Qu.:0.00000 3rd Qu.: 6.000 3rd Qu.: 0.0000 3rd Qu.:1.0000
> Max. :1.00000 Max. :34.000 Max. :34.0000 Max. :1.0000
> NA's :10 NA's :523
> obsstatus2 repstatus2 fecstatus2 matstatus2
> Min. :0.0000 Min. :0.00000 Min. :0.00000 Min. :0.0000
> 1st Qu.:1.0000 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.0000
> Median :1.0000 Median :0.00000 Median :0.00000 Median :1.0000
> Mean :0.9815 Mean :0.03506 Mean :0.01661 Mean :0.6531
> 3rd Qu.:1.0000 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:1.0000
> Max. :1.0000 Max. :1.00000 Max. :1.00000 Max. :1.0000
>
> alive2 stage2 stage2index sizea3
> Min. :1 Length:542 Min. :2.000 Min. : 0.00
> 1st Qu.:1 Class :character 1st Qu.:2.000 1st Qu.: 0.00
> Median :1 Mode :character Median :3.000 Median : 10.50
> Mean :1 Mean :3.159 Mean : 76.21
> 3rd Qu.:1 3rd Qu.:4.000 3rd Qu.: 24.50
> Max. :1 Max. :7.000 Max. :6645.80
>
> repstra3 feca3 fec3added juvgiven3
> Min. :0.00000 Min. : 0.000 Min. : 0.0000 Min. :0
> 1st Qu.:0.00000 1st Qu.: 0.000 1st Qu.: 0.0000 1st Qu.:0
> Median :0.00000 Median : 0.000 Median : 0.0000 Median :0
> Mean :0.05556 Mean : 6.045 Mean : 0.2454 Mean :0
> 3rd Qu.:0.00000 3rd Qu.: 7.500 3rd Qu.: 0.0000 3rd Qu.:0
> Max. :1.00000 Max. :48.000 Max. :48.0000 Max. :0
> NA's :146 NA's :520
> obsstatus3 repstatus3 fecstatus3 matstatus3
> Min. :0.0000 Min. :0.00000 Min. :0.00000 Min. :1
> 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:1
> Median :1.0000 Median :0.00000 Median :0.00000 Median :1
> Mean :0.7306 Mean :0.04059 Mean :0.01292 Mean :1
> 3rd Qu.:1.0000 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:1
> Max. :1.0000 Max. :1.00000 Max. :1.00000 Max. :1
>
> alive3 stage3 stage3index
> Min. :0.0000 Length:542 Min. :0.000
> 1st Qu.:1.0000 Class :character 1st Qu.:3.000
> Median :1.0000 Mode :character Median :3.000
> Mean :0.7989 Mean :3.037
> 3rd Qu.:1.0000 3rd Qu.:4.000
> Max. :1.0000 Max. :7.000
> Once again we have a dataset with 542 rows.
Now let’s build our matrices. We will use the same supplement table as before in ahistorical analysis, since there are no stages used in that supplement table unique to the function-based analysis.
lathmat2p_raw <- arlefko2(data = lathvert_raw_small, stageframe = lathframe_raw,
supplement = lathsupp2, stages = c("stage3", "stage2", "stage1"),
patch = "all", patchcol = "patchid", yearcol = "year2", agecol = "obsage",
indivcol = "individ")
summary(lathmat2p_raw)
>
> This ahistorical lefkoMat object contains 12 matrices.
>
> Each matrix is square with 14 rows and columns, and a total of 196 elements.
> A total of 138 survival transitions were estimated, with 11.5 per matrix.
> A total of 12 fecundity transitions were estimated, with 1 per matrix.
> This lefkoMat object covers 1 population, 6 patches, and 2 time steps.
>
> The dataset contains a total of 232 unique individuals and 542 unique transitions.
>
> Survival probability sum check (each matrix represented by column in order):
> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
> Min. 0.000 0.000 0.000 0.0000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
> 1st Qu. 0.000 0.000 0.000 0.0000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
> Median 0.199 0.000 0.199 0.0625 0.199 0.000 0.199 0.000 0.000 0.000 0.000
> Mean 0.393 0.220 0.350 0.2802 0.383 0.320 0.374 0.293 0.306 0.211 0.319
> 3rd Qu. 0.893 0.399 0.688 0.4748 0.842 0.725 0.781 0.577 0.635 0.383 0.600
> Max. 1.000 1.000 1.000 1.0000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
> [,12]
> Min. 0.000
> 1st Qu. 0.000
> Median 0.000
> Mean 0.184
> 3rd Qu. 0.299
> Max. 1.000Here we have 12 matrices, each with 14 rows and columns. There are a total of 196 elements per matrix, but only 12.5 are estimated per matrix due to the sparsity of data relative to the size of the matrix. Let’s take a look at the first matrix.
lathmat2p_raw$A[[1]]
> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
> [1,] 0.000 0.0000000 0.0000000 0.0000 0.000 0 1.035 0.000 0 0 0
> [2,] 0.000 0.0000000 0.0000000 0.0000 0.000 0 0.162 0.000 0 0 0
> [3,] 0.000 0.0000000 0.0000000 0.0000 0.000 0 0.000 0.000 0 0 0
> [4,] 0.000 0.0000000 0.0000000 0.0000 0.000 0 0.000 0.000 0 0 0
> [5,] 0.000 0.0000000 0.0000000 0.0000 0.000 0 0.000 0.000 0 0 0
> [6,] 0.000 0.0000000 0.0000000 0.0000 0.000 0 0.000 0.000 0 0 0
> [7,] 0.000 0.0000000 0.0000000 0.0000 0.000 0 0.000 0.000 0 0 0
> [8,] 0.345 0.0000000 0.0000000 0.0000 0.000 0 0.000 0.345 0 0 0
> [9,] 0.054 0.0000000 0.0000000 0.0000 0.000 0 0.000 0.054 0 0 0
> [10,] 0.000 0.9047619 0.7142857 0.1875 0.250 0 0.000 0.000 0 0 0
> [11,] 0.000 0.0000000 0.1428571 0.6875 0.500 0 1.000 0.000 0 0 0
> [12,] 0.000 0.0000000 0.0000000 0.0000 0.000 0 0.000 0.000 0 0 0
> [13,] 0.000 0.0000000 0.0000000 0.0625 0.125 0 0.000 0.000 0 0 0
> [14,] 0.000 0.0000000 0.0000000 0.0000 0.125 0 0.000 0.000 0 0 0
> [,12] [,13] [,14]
> [1,] 0 0 0
> [2,] 0 0 0
> [3,] 0 0 0
> [4,] 0 0 0
> [5,] 0 0 0
> [6,] 0 0 0
> [7,] 0 0 0
> [8,] 0 0 0
> [9,] 0 0 0
> [10,] 0 0 0
> [11,] 0 0 0
> [12,] 0 0 0
> [13,] 0 0 0
> [14,] 0 0 0One of the more interesting problems that we can see here is the fact that since this matrix reflects a very small dataset, both in terms of numbers of observed transitions and in terms of years of observation, we cannot really estimate transitions staying within the final age. To yield such transitions, we would typically need a large and long enough dataset that we could deliberately reduce the number of ages to a number less than the longest seen, perhaps reflecting stable demographic patterns beyond some age. Clearly, here is one area in which function-based MPMs have a clear advantage. To help increase the numbers of estimated elements, let’s develop an MPM that does not discriminate patches.
lathmat2_raw <- arlefko2(data = lathvert_raw_small, stageframe = lathframe_raw,
supplement = lathsupp2, stages = c("stage3", "stage2", "stage1"),
patchcol = "patchid", yearcol = "year2", agecol = "obsage",
indivcol = "individ")
summary(lathmat2_raw)
>
> This ahistorical lefkoMat object contains 2 matrices.
>
> Each matrix is square with 14 rows and columns, and a total of 196 elements.
> A total of 40 survival transitions were estimated, with 20 per matrix.
> A total of 6 fecundity transitions were estimated, with 3 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 2 time steps.
>
> The dataset contains a total of 232 unique individuals and 542 unique transitions.
>
> Survival probability sum check (each matrix represented by column in order):
> [,1] [,2]
> Min. 0.000 0.000
> 1st Qu. 0.000 0.000
> Median 0.199 0.199
> Mean 0.374 0.356
> 3rd Qu. 0.833 0.758
> Max. 1.000 1.000Note the increase in estimated elements - roughly twice the number per matrix as last time. Let’s take a look at the first matrix here.
lathmat2_raw$A[[1]]
> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
> [1,] 0.000 0.00000000 0.00000000 0.00000000 0.0000000 0 2.07000000 0.000
> [2,] 0.000 0.00000000 0.00000000 0.00000000 0.0000000 0 0.32400000 0.000
> [3,] 0.000 0.00000000 0.00000000 0.00000000 0.0000000 0 0.00000000 0.000
> [4,] 0.000 0.00000000 0.00000000 0.00000000 0.0000000 0 0.00000000 0.000
> [5,] 0.000 0.00000000 0.00000000 0.00000000 0.0000000 0 0.00000000 0.000
> [6,] 0.000 0.00000000 0.00000000 0.00000000 0.0000000 0 0.00000000 0.000
> [7,] 0.000 0.00000000 0.00000000 0.00000000 0.0000000 0 0.00000000 0.000
> [8,] 0.345 0.00000000 0.00000000 0.00000000 0.0000000 0 0.00000000 0.345
> [9,] 0.054 0.00000000 0.00000000 0.00000000 0.0000000 0 0.00000000 0.054
> [10,] 0.000 0.69491525 0.78787879 0.25490196 0.1052632 0 0.00000000 0.000
> [11,] 0.000 0.02542373 0.12121212 0.52941176 0.3157895 0 0.36363636 0.000
> [12,] 0.000 0.00000000 0.00000000 0.00000000 0.2105263 0 0.27272727 0.000
> [13,] 0.000 0.02542373 0.03030303 0.05882353 0.1052632 0 0.09090909 0.000
> [14,] 0.000 0.00000000 0.00000000 0.01960784 0.1578947 0 0.27272727 0.000
> [,9] [,10] [,11] [,12] [,13] [,14]
> [1,] 0 0 0 0 0 0
> [2,] 0 0 0 0 0 0
> [3,] 0 0 0 0 0 0
> [4,] 0 0 0 0 0 0
> [5,] 0 0 0 0 0 0
> [6,] 0 0 0 0 0 0
> [7,] 0 0 0 0 0 0
> [8,] 0 0 0 0 0 0
> [9,] 0 0 0 0 0 0
> [10,] 0 0 0 0 0 0
> [11,] 0 0 0 0 0 0
> [12,] 0 0 0 0 0 0
> [13,] 0 0 0 0 0 0
> [14,] 0 0 0 0 0 0Let’s now move on to Leslie MPMs.
6.4 Age-classified (Leslie) MPMs
Package lefko3 can also be used to estimate Leslie MPMs, which are purely age-based and so include no individual history. Here, we will illustrate how to create both raw and function-based Leslie MPMs using the rleslie() and fleslie() functions, respectively. We will ignore the dormant seed stage, as including dormant seeds would require using an age-by-stage approach. Our final matrices will take the following form.
\[\begin{equation} \tag{6.3} \left[\begin{array} {rrrr} F _{0,0} & F _{0,1} & F _{0,2} & F _{0,3} \\ S _{1,0} & 0 & 0 & 0 \\ 0 & S _{2,1} & 0 & 0 \\ 0 & 0 & S _{3,2} & S_{3,3} \\ \end{array}\right] \end{equation}\]
Normally, we would need to develop a stageframe. However, in the purely age-based case, a stageframe is unnecessary and instead all we need is the standardized hfv dataset. So, let’s start by creating that. Most of the settings will be as before. However, we need to include both NRasRep = TRUE and NOasObs = TRUE to make sure that stage classification ignores whether the individual was actually reproductive and actually observed at each time (the latter has to do with the fact that individuals can be alive with a size of zero if they are vegetatively dormant). We will also subset our data to only those individuals whose age we are reasonably sure of, by eliminating those individuals first observed in the first monitoring session.
lathvert_base <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
sizeacol = "Volume88", repstracol = "FCODE88", fecacol = "Intactseed88",
deadacol = "Dead1988", censorcol = "Missing1988", censorkeep = NA,
censor = TRUE, NRasRep = TRUE, NOasObs = TRUE)
lathvert_age <- subset(lathvert_base, firstseen > 1988)
summary(lathvert_age)
> rowid popid patchid individ year2
> Min. : 30.0 :542 1:131 Length:542 Min. :1989
> 1st Qu.: 302.5 2: 72 Class :character 1st Qu.:1989
> Median : 583.0 3:101 Mode :character Median :1990
> Mean : 586.8 4:132 Mean :1990
> 3rd Qu.: 801.0 5: 70 3rd Qu.:1990
> Max. :1097.0 6: 36 Max. :1990
>
> firstseen lastseen obsage obslifespan sizea1
> Min. :1989 Min. :1989 Min. :1.000 Min. :0.00 Min. : 0.00
> 1st Qu.:1989 1st Qu.:1990 1st Qu.:1.000 1st Qu.:1.00 1st Qu.: 0.00
> Median :1989 Median :1991 Median :1.000 Median :2.00 Median : 0.00
> Mean :1989 Mean :1991 Mean :1.352 Mean :1.36 Mean : 59.41
> 3rd Qu.:1989 3rd Qu.:1991 3rd Qu.:2.000 3rd Qu.:2.00 3rd Qu.: 9.00
> Max. :1990 Max. :1991 Max. :2.000 Max. :2.00 Max. :4394.20
>
> repstra1 feca1 fec1added juvgiven1 obsstatus1
> Min. :0.0000 Min. : 0 Min. : 0.0000 Min. :0 Min. :0.0000
> 1st Qu.:0.0000 1st Qu.: 0 1st Qu.: 0.0000 1st Qu.:0 1st Qu.:0.0000
> Median :0.0000 Median : 4 Median : 0.0000 Median :0 Median :0.0000
> Mean :0.0576 Mean : 6 Mean : 0.1218 Mean :0 Mean :0.3524
> 3rd Qu.:0.0000 3rd Qu.: 6 3rd Qu.: 0.0000 3rd Qu.:0 3rd Qu.:1.0000
> Max. :1.0000 Max. :34 Max. :34.0000 Max. :0 Max. :1.0000
> NA's :351 NA's :531
> repstatus1 fecstatus1 matstatus1 alive1
> Min. :0.0000 Min. :0.00000 Min. :1 Min. :0.0000
> 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:1 1st Qu.:0.0000
> Median :0.0000 Median :0.00000 Median :1 Median :0.0000
> Mean :0.0203 Mean :0.01292 Mean :1 Mean :0.3524
> 3rd Qu.:0.0000 3rd Qu.:0.00000 3rd Qu.:1 3rd Qu.:1.0000
> Max. :1.0000 Max. :1.00000 Max. :1 Max. :1.0000
>
> stage1 stage1index sizea2 repstra2
> Length:542 Min. :0 Min. : 0.00 Min. :0.00000
> Class :character 1st Qu.:0 1st Qu.: 9.00 1st Qu.:0.00000
> Mode :character Median :0 Median : 12.60 Median :0.00000
> Mean :0 Mean : 97.24 Mean :0.03571
> 3rd Qu.:0 3rd Qu.: 37.40 3rd Qu.:0.00000
> Max. :0 Max. :4394.20 Max. :1.00000
> NA's :10
> feca2 fec2added juvgiven2 obsstatus2
> Min. : 0.000 Min. : 0.0000 Min. :0 Min. :0.0000
> 1st Qu.: 0.000 1st Qu.: 0.0000 1st Qu.:0 1st Qu.:1.0000
> Median : 0.000 Median : 0.0000 Median :0 Median :1.0000
> Mean : 4.053 Mean : 0.1421 Mean :0 Mean :0.9815
> 3rd Qu.: 6.000 3rd Qu.: 0.0000 3rd Qu.:0 3rd Qu.:1.0000
> Max. :34.000 Max. :34.0000 Max. :0 Max. :1.0000
> NA's :523
> repstatus2 fecstatus2 matstatus2 alive2
> Min. :0.00000 Min. :0.00000 Min. :1 Min. :1
> 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:1 1st Qu.:1
> Median :0.00000 Median :0.00000 Median :1 Median :1
> Mean :0.03506 Mean :0.01661 Mean :1 Mean :1
> 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:1 3rd Qu.:1
> Max. :1.00000 Max. :1.00000 Max. :1 Max. :1
>
> stage2 stage2index sizea3 repstra3
> Length:542 Min. :0 Min. : 0.00 Min. :0.00000
> Class :character 1st Qu.:0 1st Qu.: 0.00 1st Qu.:0.00000
> Mode :character Median :0 Median : 10.50 Median :0.00000
> Mean :0 Mean : 76.21 Mean :0.05556
> 3rd Qu.:0 3rd Qu.: 24.50 3rd Qu.:0.00000
> Max. :0 Max. :6645.80 Max. :1.00000
> NA's :146
> feca3 fec3added juvgiven3 obsstatus3
> Min. : 0.000 Min. : 0.0000 Min. :0 Min. :0.0000
> 1st Qu.: 0.000 1st Qu.: 0.0000 1st Qu.:0 1st Qu.:0.0000
> Median : 0.000 Median : 0.0000 Median :0 Median :1.0000
> Mean : 6.045 Mean : 0.2454 Mean :0 Mean :0.7306
> 3rd Qu.: 7.500 3rd Qu.: 0.0000 3rd Qu.:0 3rd Qu.:1.0000
> Max. :48.000 Max. :48.0000 Max. :0 Max. :1.0000
> NA's :520
> repstatus3 fecstatus3 matstatus3 alive3
> Min. :0.00000 Min. :0.00000 Min. :1 Min. :0.0000
> 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:1 1st Qu.:1.0000
> Median :0.00000 Median :0.00000 Median :1 Median :1.0000
> Mean :0.04059 Mean :0.01292 Mean :1 Mean :0.7528
> 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:1 3rd Qu.:1.0000
> Max. :1.00000 Max. :1.00000 Max. :1 Max. :1.0000
>
> stage3 stage3index
> Length:542 Min. :0
> Class :character 1st Qu.:0
> Mode :character Median :0
> Mean :0
> 3rd Qu.:0
> Max. :0
>
So far so good! Now let’s create the vital rate models. This will go much more quickly than last time, because we no longer care about size and reproductive status as factors determining vital rates. So, we will set suite = "cons" to prevent these factors from being tested, and set age = "obsage" to incorporate our age at time t variable into all models. We will also use global.only = TRUE to prevent the age term from being dropped (normally we would not set this, but in our case the dataset is so small that age will drop out in the best-fit models, so we will force the models to include age in this case).
lathmodels2_age <- modelsearch(lathvert_age, historical = FALSE,
approach = "mixed", suite = "cons", bestfit = "AICc&k", age = "obsage",
vitalrates = c("surv", "fec"), fecdist = "negbin", indiv = "individ",
year = "year2", year.as.random = TRUE, patch.as.random = TRUE,
show.model.tables = TRUE, fec.zero = TRUE, global.only = TRUE, quiet = TRUE)Let’s see a summary of the models.
summary(lathmodels2_age)
> This LefkoMod object includes 2 linear models.
> Best-fit model criterion used: global model only
>
>
>
> Survival model:
> Generalized linear mixed model fit by maximum likelihood (Laplace
> Approximation) [glmerMod]
> Family: binomial ( logit )
> Formula: alive3 ~ obsage + (1 | year2) + (1 | individ)
> Data: subdata
> AIC BIC logLik deviance df.resid
> 608.1500 625.3311 -300.0750 600.1500 538
> Random effects:
> Groups Name Std.Dev.
> individ (Intercept) 0.0002353
> year2 (Intercept) 0.3738235
> Number of obs: 542, groups: individ, 232; year2, 2
> Fixed Effects:
> (Intercept) obsage
> 0.9445 0.1817
>
>
>
> Observation model:
> [1] 1
>
>
>
> Size model:
> [1] 1
>
>
>
> Secondary size model:
> [1] 1
>
>
>
> Tertiary size model:
> [1] 1
>
>
>
> Reproductive status model:
> [1] 1
>
>
>
> Fecundity model:
> Formula: feca2 ~ obsage + (1 | year2) + (1 | individ)
> Zero inflation: ~.
> Data: subdata
> AIC BIC logLik df.resid
> NA NA NA 10
> Random-effects (co)variances:
>
> Conditional model:
> Groups Name Std.Dev.
> year2 (Intercept) 0.02271
> individ (Intercept) 1.19287
>
> Zero-inflation model:
> Groups Name Std.Dev.
> year2 (Intercept) 0.09092
> individ (Intercept) 0.17023
>
> Number of obs: 19 / Conditional model: year2, 2; individ, 16 / Zero-inflation model: year2, 2; individ, 16
>
> Dispersion parameter for nbinom2 family (): 88.6
>
> Fixed Effects:
>
> Conditional model:
> (Intercept) obsage
> 1.3220 0.3938
>
> Zero-inflation model:
> (Intercept) obsage
> -3.842 2.635
>
>
> Juvenile survival model:
> [1] 1
>
>
>
> Juvenile observation model:
> [1] 1
>
>
>
> Juvenile size model:
> [1] 1
>
>
>
> Juvenile secondary size model:
> [1] 1
>
>
>
> Juvenile tertiary size model:
> [1] 1
>
>
>
> Juvenile reproduction model:
> [1] 1
>
>
>
> Juvenile maturity model:
> [1] 1
>
>
>
>
>
> Number of models in survival table: 1
>
> Number of models in observation table: 1
>
> Number of models in size table: 1
>
> Number of models in secondary size table: 1
>
> Number of models in tertiary size table: 1
>
> Number of models in reproduction status table: 1
>
> Number of models in fecundity table: 1
>
> Number of models in juvenile survival table: 1
>
> Number of models in juvenile observation table: 1
>
> Number of models in juvenile size table: 1
>
> Number of models in juvenile secondary size table: 1
>
> Number of models in juvenile tertiary size table: 1
>
> Number of models in juvenile reproduction table: 1
>
> Number of models in juvenile maturity table: 1
>
>
>
>
>
> General model parameter names (column 1), and
> specific names used in these models (column 2):
> parameter_names mainparams
> 1 time t year2
> 2 individual individ
> 3 patch patch
> 4 alive in time t+1 surv3
> 5 observed in time t+1 obs3
> 6 sizea in time t+1 size3
> 7 sizeb in time t+1 sizeb3
> 8 sizec in time t+1 sizec3
> 9 reproductive status in time t+1 repst3
> 10 fecundity in time t+1 fec3
> 11 fecundity in time t fec2
> 12 sizea in time t size2
> 13 sizea in time t-1 size1
> 14 sizeb in time t sizeb2
> 15 sizeb in time t-1 sizeb1
> 16 sizec in time t sizec2
> 17 sizec in time t-1 sizec1
> 18 reproductive status in time t repst2
> 19 reproductive status in time t-1 repst1
> 20 maturity status in time t+1 matst3
> 21 maturity status in time t matst2
> 22 age in time t age
> 23 density in time t density
> 24 individual covariate a in time t indcova2
> 25 individual covariate a in time t-1 indcova1
> 26 individual covariate b in time t indcovb2
> 27 individual covariate b in time t-1 indcovb1
> 28 individual covariate c in time t indcovc2
> 29 individual covariate c in time t-1 indcovc1
> 30 stage group in time t group2
> 31 stage group in time t-1 group1
>
>
>
>
>
> Quality control:
>
> Survival estimated with 232 individuals and 542 individual transitions.
> Survival accuracy is 0.753.
> Observation probability not estimated.
> Primary size transition not estimated.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproduction probability not estimated.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Fecundity pseudo R-squared is NA.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity probability not estimated.The two models all include observed age, as well as our random terms. Now we will create the raw matrices.
lathmat2ageonly_raw <- rleslie(data = lathvert_age, year = "all",
yearcol = "year2", indivcol = "individ")
summary(lathmat2ageonly_raw)
>
> This ahistorical lefkoMat object contains 1 matrix.
>
> Each matrix is square with 3 rows and columns, and a total of 9 elements.
> A total of 2 survival transitions were estimated, with 2 per matrix.
> A total of 2 fecundity transitions were estimated, with 2 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 1 time step.
>
> The dataset contains a total of 232 unique individuals and 542 unique transitions.
>
> Survival probability sum check (each matrix represented by column in order):
> [,1]
> Min. 0.000
> 1st Qu. 0.332
> Median 0.664
> Mean 0.462
> 3rd Qu. 0.693
> Max. 0.723We have just created the raw Leslie MPM. We have a single matrix, because we needed to eliminate the first year of data to estimate age properly, leaving only three observation periods. The quality control looks reasonable in terms of the survival-transition probabilities falling within the range of 0.0 to 1.0. Let’s also take a look at our one matrix, to get a handle on its structure.
lathmat2ageonly_raw$A[[1]]
> [,1] [,2] [,3]
> [1,] 0.02521008 0.04188482 0
> [2,] 0.66386555 0.00000000 0
> [3,] 0.00000000 0.72251309 0Our matrix certainly looks like a valid Leslie matrix, with survival transitions in the subdiagonal and fecundity at the top. The last column is entirely composed of zeros, and that is likely due to a combination of a small dataset and a short dataset in which we cannot find individuals whom we know are older than three years old.
Let’s now build the function-based versions and see what happens.
lathmat2ageonly_func <- fleslie(year = "all", data = lathvert_age,
modelsuite = lathmodels2_age)
summary(lathmat2ageonly_func)
>
> This ahistorical lefkoMat object contains 2 matrices.
>
> Each matrix is square with 3 rows and columns, and a total of 9 elements.
> A total of 6 survival transitions were estimated, with 3 per matrix.
> A total of 6 fecundity transitions were estimated, with 3 per matrix.
> This lefkoMat object covers 1 population, 1 patch, and 2 time steps.
>
> Vital rate modeling quality control:
>
> Survival estimated with 232 individuals and 542 individual transitions.
> Observation probability not estimated.
> Primary size transition not estimated.
> Secondary size transition not estimated.
> Tertiary size transition not estimated.
> Reproduction probability not estimated.
> Fecundity estimated with 16 individuals and 19 individual transitions.
> Juvenile survival not estimated.
> Juvenile observation probability not estimated.
> Juvenile primary size transition not estimated.
> Juvenile secondary size transition not estimated.
> Juvenile tertiary size transition not estimated.
> Juvenile reproduction probability not estimated.
> Juvenile maturity transition probability not estimated.
>
> Survival probability sum check (each matrix represented by column in order):
> [,1] [,2]
> Min. 0.813 0.685
> 1st Qu. 0.826 0.704
> Median 0.839 0.723
> Mean 0.838 0.722
> 3rd Qu. 0.850 0.740
> Max. 0.862 0.758
lathmat2ageonly_func$A[[1]]
> [,1] [,2] [,3]
> [1,] 0.02118612 0.2318637 0.8080455
> [2,] 0.81273691 0.0000000 0.0000000
> [3,] 0.00000000 0.8388330 0.8619104Notice that we have one more matrix here (reflecting the two estimable years), and we have estimates survival transition probabilities staying within the last age. Clearly, function-based approaches have an advantage in this small, short dataset.
Everything looks quite good! Congratulations on creating your Leslie MPMs!
6.5 Points to remember
- Age-based (Leslie) MPMs can be be built in raw or function-based form using functions
rleslie()andfleslie(), respectively. Age-by-stage MPMs can be built ion raw or function-based forms using functionsarlefko2()andaflefko2(), respectively. - Age-by-stage MPMs are larger than either age-based or stage-based MPMs, and so generally require more data to parameterize properly.
- Dot all possible ages need to be modeled in age-based and age-by-stage MPMs. Often, the final age in these MPMs represents a long stretch of the adult span of life during which vital rates are not expected to change dramatically, and so is capable of self-transition while other ages are not.