Lathyrus vernus raw MPMs
In this vignette, we will use the lathyrus dataset to
illustrate the estimation of empirical or
raw MPMs. We will produce matrices similar to those
published in Ehrlén (2000), though there
will be some differences because our dataset includes data for more
individuals as well as an extra year of monitoring. It also includes
differences in classification due to different assumptions regarding
transitions to and from the vegetative dormancy life stage.
To reduce vignette size, we have prevented some statements from running
if they produce long stretches of output. Examples include most
summary() calls. In these cases, we include hashtagged
versions of these calls, and we encourage the user to run these
statements without hashtags.
This vignette is only a sample analysis. Detailed information and
instructions on using lefko3 are available through a free
online e-book called lefko3: a gentle introduction, available
on the projects
page of the Shefferson lab website.
ORGANISM AND POPULATION
Lathyrus vernus (family Fabaceae) is a long-lived forest herb, native to Europe and large parts of northern Asia. Individuals increase slowly in size and usually flower only after 10-15 years of vegetative growth. Flowering individuals have an average conditional lifespan of 44.3 years (Ehrlén and Lehtila 2002). L. vernus lacks organs for vegetative spread and individuals are well delimited (Ehrlén 2002). One or several erect shoots of up to 40 cm height emerge from a subterranean rhizome in March and April. Flowering occurs about four weeks after shoot emergence. Shoot growth is determinate, and the number of flowers is determined in the previous year (Ehrlén and Van Groenendael 2001). Individuals do not necessarily produce aboveground structures every year, and instead can remain vegetatively dormant in one or more seasons. L. vernus is self-compatible but requires visits from bumble-bees to produce seeds. Individuals produce few, large seeds and establishment from seeds is relatively frequent (Ehrlén and Eriksson 1996). The pre-dispersal seed predator Bruchus atomarius often consumes a large fraction of developing seeds, and roe deer (Capreolus capreolus) sometimes consume the shoots (Ehrlén and Munzbergova 2009).
Data for this study were collected from six permanent plots in a population of L. vernus located in a deciduous forest in the Tullgarn area, SE Sweden, from 1988 to 1991 (Ehrlén 1995). The plots had similar soil type, elevation, slope, and canopy cover. Within each plot, all individuals were marked with numbered tags that remained over the study period, and their locations were mapped. At the time of shoot emergence, we recorded whether individuals were alive and produced above-ground shoots, and if shoots had been grazed. During flowering, we recorded flower number and the height and diameter of all shoots. At fruit maturation, we counted the number of intact and damaged seeds. To derive a measure of aboveground size for each individual, we calculated the volume of each shoot as $\pi × (\frac{1}{2} diameter)^2 × height$, and summed the volumes of all shoots. This measure is strongly correlated with the dry mass of aboveground tissues (R2 = 0.924, P < 0.001, n = 50, log-transformed values; Ehrlén 1995). Size of grazed individuals was estimated based on measures of shoot diameter in grazed shoots, and the relationship between shoot diameter and shoot height in non-grazed individuals. Individuals that lacked aboveground structures in one season but reappeared in the following year were considered dormant. Individuals that lacked aboveground structures in two subsequent seasons were considered dead from the year in which they first lacked aboveground structures. Probabilities of seeds surviving to the next year, and of being present as seedlings or seeds in the soil seed bank, were derived from separate yearly sowing experiments in separate plots adjacent to each subplot (Ehrlén and Eriksson 1996).
BASIC WORKFLOW
Our dataset is organized in horizontal format, with rows corresponding to unique individuals and columns corresponding to individual condition in particular observation occasions (which we also refer to as years here). The original spreadsheet file used to keep the dataset has a repeating pattern to these columns, with each year having a similarly arranged group of variables. Let’s take a look at this dataset.
This dataset includes information on 1,119 individuals arranged
horizontally, by row. There are 38 variables, by column. The first two
columns are variables giving identifying information about each
individual. 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. This
strictly repeating pattern allows us to manipulate the original dataset
quickly and efficiently via lefko3. There are four years of
data, from 1988 to 1991. Ideally, we should also have arranged the
columns in the same order for each year, with years in consecutive order
with no extra columns between them. This order is not required, provided
that we input all variable names in correct order when transforming the
dataset later.
Step 1. Life history model development
We will now create a stageframe, which is a data frame that describes all stages in the life history of the organism, in a way usable by the functions in this package and using stage names and descriptions that completely match those used in the dataset. It links the dataset to the life cycle graph used to model the organism’s life history. It should include complete descriptions of all stages that occur in the dataset, with each stage defined uniquely. Since this object can be used for automated classification, all sizes, reproductive states, and other characteristics defining each stage in the dataset need to be accounted for explicitly. The final description of each stage occurring in the dataset must not completely overlap with any other stage also found in the dataset, although partial overlap is allowed and expected.
Here, we create a stageframe named lathframe based on the
classification used in Ehrlén (2000). We
build this by creating vectors of the characteristics describing each
stage, with each element always in the same order within the vector,
using the sf_create() function. The most important settings
have to do with the size bins for the stages. Two methods are currently
allowed: 1) a vector of representative sizes in the option called
sizes along with the associated minimum and maximum sizes
for each stage in the options called sizemin and
sizemax, or 2) a vector of stage size bin midpoints in the
sizes option along with the half-width of each size bin in
the binhalfwidth option. Note that the representative
sizes, whether midpoint or otherwise, are not actually used in any
calculations - the size bin minima and maxima are the most important
data used in calculations. Here, we will use approach 2), focused on the
use of size bin half-widths.
If size values are not to be binned, then narrow half bin-widths can be
used (the default is 0.5, and all entries must be positive). For
example, in this dataset, vegetatively dormant individuals necessarily
have a size of zero, and so we can set the halfbinwidth for
this stage to 0.5 provided that the resulting size bin does not overlap
with any other size bin matching the other characteristics of vegetative
dormancy. We can also set this manually with the sizemin
and sizemax options, if we prefer.
Here we use a single size variable, but up to three size variables may
be used. Vector stagenames must include unique names only.
Vectors repstatus, obsstatus,
matstatus, immstatus, propstatus,
and indataset are binomial vectors referring to status as a
reproductive stage, status as an observed stage, status as a mature
stage, status as an immature stage, status as a propagule stage, and
status as a stage occurring within the user-supplied dataset,
respectively. The combination of these characteristics must be
completely unique for each stage. The final vector, called
comments, holds strings of stage descriptions.
sizevector <- c(0, 100, 13, 127, 3730, 3800, 0)
stagevector <- c("Sd", "Sdl", "Tm", "Sm", "La", "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)
comments <- c("Dormant seed", "Seedling", "Tiny vegetative", "Small vegetative",
"Large vegetative", "Flowering", "Vegetatively dormant")
lathframe <- sf_create(sizes = sizevector, stagenames = stagevector,
repstatus = repvector, obsstatus = obsvector, propstatus = propvector,
immstatus = immvector, matstatus = matvector, indataset = indataset,
binhalfwidth = binvec, comments = comments)
#lathframe
Care should be taken in assigning sizes to stages, particularly when
stages occur with size = 0. In most cases, a size of zero
will mean that the individual is alive but not observable, such as in
the case of vegetative dormancy. However, a size of zero may have other
meanings. For example, if the size metric used is a logarithm of the
measured size, then observable sizes of zero and lower may be possible.
These situations may impact matrix construction and analysis,
particularly when dealing with function-based MPMs, such as IPMs.
Step 2a. Data standardization
Next, we will standardize the dataset into vertical format.
Vertically formatted datasets are structured such that each row
corresponds to the state of a single individual in two (if ahistorical)
or three (if historical) consecutive monitoring occasions. Function
verticalize3() will create a historically formatted
vertical dataset, or hfv dataset. We can use one of two
approaches for this task, both using this function. In the most general
case, we may input all column names from the original dataset in each
input option in the same order, corresponding to monitoring occasion.
This allows any format of dataset to be used. An easier approach may be
used in datasets with strictly repeating sets of columns. In the latter
case, most of the inputs should constitute the names of the first
variables coding for particular states. For example,
Seedling1988 is the first variable coding for status as a
juvenile, Volume88 is the first variable coding for the
main size metric, and Intactseed88 is the first variable
coding for fecundity. Function verticalize3() can determine
the locations of variables for later monitoring occasions assuming a
repeating pattern of such variables organized as blocks for each year
(noted as blocksize, which here equals 9 because there are
nine such variables per year in the same order). There are four
monitoring occasions (noted as noyears). We also have a
repeated censor variable, the first of which is
Missing1988, and we note here that we wish to censor the
data and to keep data points with NA values in the censor
term. The patchidcol and individcol terms are
variables denoting which patch/subpopulation and individual each row of
data belongs to, respectively. The setting NAas0 = TRUE
tells R that missing values in size and fecundity should be interpreted
as zeros, which allows us to infer incidents of vegetative dormancy as
cases where size is equal to zero. Finally, stageassign
ties the dataset to the correct stageframe. Note that prior to
standardizing this dataset, we need to create a new individual identity
variable, indiv_id, composed of the identity of the patch
that the individual is in as well as its patch-specific ID, because
otherwise individual identity will be shared by plants in different
patches.
lathyrus$indiv_id <- paste(lathyrus$SUBPLOT, lathyrus$GENET)
lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
patchidcol = "SUBPLOT", individcol = "indiv_id", blocksize = 9,
juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88",
fecacol = "Intactseed88", deadacol = "Dead1988", nonobsacol = "Dormant1988",
stageassign = lathframe, stagesize = "sizea", censorcol = "Missing1988",
censorkeep = NA, censorRepeat = TRUE, censor = TRUE, NAas0 = TRUE)
summary_hfv(lathvert, full = FALSE)
#>
#> This hfv dataset contains 2527 rows, 54 variables, 1 population,
#> 6 patches, 1053 individuals, and 3 time steps.
The resulting reorganization has dramatically changed the dimensions of
the dataset, which started with 1119 rows and 38 variables, and now has
2527 rows and 54 variables. The verticalize3() function
includes error-checking measures designed to find instances where
individual characteristics do not match those assigned to stages in the
associated stageframe. In those instances, warning messages will be
displayed and the instances will be marked NoMatch in
stage1, stage2, or stage3. The
summary_hfv() function allows us to quickly summarize the
main characteristics of our resulting hfv dataset, and we can
see a longer summary by removing full = FALSE.
If we remove full = FALSE from the
summary_hfv() call, then our summary will reveal that the
verticalize3() function has automatically subset the data
to only those instances in which the individual is alive in occasion
t (please look at the variable alive2). Knowing
this can help us interpret other variables. For example, the mean value
for alive3 suggests very high survival to occasion
t+1 (92.2%). Further, the minimum values for all size variables
are zero, suggesting that unobservable stages occur within the dataset
(these are instances of vegetative dormancy).
The verticalize3() function has made stage assignments for
all individuals at each time. This can be seen in the subset summary in
the stage2index column, which shows all individuals that
are alive and have a size of 0 in occasion t to be in the
seventh stage. Type lathframe and enter at the prompt to
check that the seventh stage in the stageframe is really vegetative
dormancy.
A data subset summary can teach us how reproduction is handled here. The
reproductive status of flowering adults is set as reproductive in
lathframe, and a subset summary of just those individuals
flowering in occasion t shows all of those individuals set to
reproductive (see the distribution of values for repstatus2
in the output below). However, fecundity ranges from 0 to 66 (see
feca2, which codes for fecundity in occasion t),
meaning that some flowering adults do not actually produce any
offspring. In fact, only 44.2% of flowering plants produced seed in
occasion t. This happened because our reproductive status
variable, FCODE88, notes whether these individuals flowered
but not whether they produced seed. Since this plant must be pollinated
by an insect vector, some flowers should yield no seed. This issue does
not cause problems in the creation of raw matrices, but it might cause
difficulties in the creation of function-based matrices under some
conditions. It helps to consider whether the definitions used for stages
are appropriate, and so whether reproductive status must necessarily be
associated with successful reproduction or merely the attempt.
Here, we associate it with the latter, but in other vignettes we will
reconsider this assumption.
summary(lathvert[which(lathvert$stage2 == "Flo"),c(22:27)])
#> alive1 stage1 stage1index sizea2 size2added repstra2
#> Min. :0.000 Length:599 Min. :0.000 Min. : 98.4 Min. : 98.4 Min. :1
#> 1st Qu.:0.000 Class :character 1st Qu.:0.000 1st Qu.: 732.5 1st Qu.: 732.5 1st Qu.:1
#> Median :1.000 Mode :character Median :5.000 Median :1141.8 Median :1141.8 Median :1
#> Mean :0.576 Mean :3.297 Mean :1388.5 Mean :1388.5 Mean :1
#> 3rd Qu.:1.000 3rd Qu.:6.000 3rd Qu.:1758.0 3rd Qu.:1758.0 3rd Qu.:1
#> Max. :1.000 Max. :7.000 Max. :7032.0 Max. :7032.0 Max. :1
summary(lathvert[which(lathvert$stage2 == "Flo"),c(28:33)])
#> repstr2added feca2 fec2added censor2 juvgiven2 obsstatus2
#> Min. :1 Min. : 0.000 Min. : 0.000 Min. :0 Min. :0 Min. :1
#> 1st Qu.:1 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.:0 1st Qu.:0 1st Qu.:1
#> Median :1 Median : 0.000 Median : 0.000 Median :0 Median :0 Median :1
#> Mean :1 Mean : 4.793 Mean : 4.793 Mean :0 Mean :0 Mean :1
#> 3rd Qu.:1 3rd Qu.: 6.000 3rd Qu.: 6.000 3rd Qu.:0 3rd Qu.:0 3rd Qu.:1
#> Max. :1 Max. :66.000 Max. :66.000 Max. :0 Max. :0 Max. :1Step 2b. Provide supplemental information for matrix estimation
Now we will create supplement tables, which provide external data for matrix estimation not included in the main demographic dataset. Specifically, we will provide the seed dormancy probability and germination rate, which are given as transitions from the dormant seed stage to another year of seed dormancy or to the germinated seedling stage, respectively. We assume that the germination rate is the same regardless of whether the seed was produced in the previous year or has been in the seedbank for longer. We will incorporate both terms as constants for specific transitions within our matrices, and as constant multipliers for fecundity, since fecundity will be estimated as the product of seed produced and either the seed germination rate or the seed dormancy/survival rate. The fecundity multipliers will also serve to tell R which transitions are the fecundity transitions.
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 = lathframe, historical = FALSE)
#lathsupp2
The supplement table above will only work with ahistorical MPMs, while
the next supplement table will work for historical MPMs. The primary
difference is the incorporation of stage in time t-1. Going
back a further step to look at time t-1 shows us that there are
two ways that a dormant seed in time t can be a dormant seed in
time t+1 - it could have been a dormant seed in time
t-1, or it could have been produced by a flowering adult in
time t-1. Likewise, a seed in occasion t could have
become a seedling in occasion t+1 after being produced by a
flowering adult in occasion t-1, or remaining in the seed bank
in occasion t-1. So, we will enter four given transitions in
the historical case, rather than two transitions as in the ahistorical
case. We will also designate that transitions from seed in occasion
t-1 and seedling in occasion t to mature stages in
occasion t+1 should originate from the NotAlive
stage in occasion t-1, because the seed stage does not actually
occur in the dataset and so function rlefko3() will not
estimate these transitions without knowing the appropriate proxy set.
lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"),
stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"),
stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"),
eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"),
eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"),
eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"),
givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA),
multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA),
type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1),
stageframe = lathframe, historical = TRUE)
#lathsupp3
These two supplement tables show us that we have survival-transition
probabilities (type = 1, whereas fecundity rates would be
given as type = 2) and fecundity multipliers
(type = 3), that the second and fourth transitions involve
reproduction from occasion t-1 to occasion t while the
others involve survival (type_t12 = 1 for survival between
occasions t-1 and t, and type_t12 = 2 for
fecundity), that the given transitions originate from the dormant seed
stage (Sd) in occasion t (and seeds or
reproductive stages in occasion t-1 in the historical case),
and the specific values to be used in overwriting and to multiply
fecundity values by: 0.345 and 0.054. If we
wished, we could have used the values of transitions to be estimated
within this matrix as proxies for these values, in which case we would
enter the stages corresponding to the correct transitions in the
eststageX columns, and the givenrate column
would be blank. This is precisely what we did for one set of transitions
- namely transitions from the seedling class to mature stages. Finally,
we have also included fecundity multipliers for newly produced seed in
the bottom two rows.
Note that the supplemental() function allowed us to isolate
specific transitions to alter, and to use shorthand to identify large
groups of transitions (e.g. using mat, immat,
rep, nrep, prop,
npr, obs, nobs,
groupX, or all to signify all mature stages,
immature stages, reproductive stages, non-reproductive stages, propagule
stages, non-propagule stages, observable stages, non-observable stages,
stages within group X, or simply all stages, respectively).
Function supplemental() also allowed us to tell R which
transitions to treat as primary reproductive transitions, which was
accomplished by identifying these transitions with reproductive
multipliers.
Step 3. Tests of history
We have chosen to build both ahistorical and historical MPMs in this vignette. However, in a typical analysis, it is most parsimonious to test whether history influences the demography of the population significantly first, and only use historical MPMs if the test supports the hypothesis that it does. A number of methods exist to conduct these tests, and we recommend Brownie et al. (1993), Pradel, Wintrebert, and Gimenez (2003), Pradel (2005), and Cole et al. (2014) for good discussions and tools to help with this.
In lefko3, we provide a function that can be used to test
for the effects of history directly from the historical vertical
dataset: modelsearch(). Function modelsearch()
estimates best-fit linear models of the key vital rates used to
propagate elements in function-based MPMs. There are up to 14 different
vital rates possible to test, of which seven are adult vital rates and
seven are juvenile vital rates. The standard vital rates that we may
wish to test are survival (marked as surv in
vitalrates), primary size (size), and
fecundity (fec), which are the default vital rates assumed
by the function. Here, we also test observation status
(obs), which can serve as a proxy for sprouting probability
in cases where plants do not necessarily sprout, and reproductive status
(repst), which assesses the probability of reproduction in
cases where reproduction is not certain. This dataset also includes
juveniles whose vital rates we wish to estimate. We designate this by
setting juvestimate to the correct juvenile stage in the
dataset. Because size varies in juveniles, we also set
juvsize = TRUE.
To test history with function modelsearch: 1) use the
historical vertical dataset as input, 2) set
historical = TRUE, 3) input the relevant vital rates to
estimate, 4) set the suite of independent factors to test (size and
reproductive status in occasions t and t-1 and all
interactions, or some subset thereof), 5) set the name of the juvenile
stage (if juvenile vital rates are to be estimated and such stages occur
in the dataset), 6) set the proper distributions to use for size and
fecundity, and 7) note which variables code for individual identity
(used to treat identity as a random factor in mixed linear models),
patch identity (if multiple patches occur in the dataset and vital rates
should be estimated with patch as a factor), and observation occasion.
We also set quiet = "partial" to limit the amount of text
output while the function runs.
histtest <- modelsearch(lathvert, historical = TRUE, suite = "main",
vitalrates = c("surv", "obs", "size", "repst", "fec"), juvestimate = "Sdl",
sizedist = "gaussian", fecdist = "gaussian", indiv = "individ",
year = "year2", juvsize = TRUE, quiet = "partial")
#>
#> Developing global model of survival probability...
#>
#> Global model of survival probability developed. Proceeding with model dredge...
#>
#> Developing global model of observation probability...
#>
#> Global model of observation probability developed. Proceeding with model dredge...
#>
#> Developing global model of primary size...
#>
#> Global model of primary size developed. Proceeding with model dredge...
#>
#> Developing global model of reproduction probability...
#>
#> Global model of reproduction probability developed. Proceeding with model dredge...
#>
#> Developing global model of fecundity...
#>
#> Global model of fecundity developed. Proceeding with model dredge...
#>
#> Developing global model of juvenile survival probability...
#>
#> Global model of juvenile survival probability developed. Proceeding with model dredge...
#> Warning: Juvenile maturity status in time t+1 appears to be constant (1). Setting to constant.
#>
#> Developing global model of juvenile observation probability...
#>
#> Global model of juvenile observation probability developed. Proceeding with model dredge...
#>
#> Developing global model of juvenile primary size...
#>
#>
#> Developing global model of juvenile primary size...
#>
#> Global model of juvenile primary size developed. Proceeding with model dredge...
#> Warning: Juvenile reproductive status in time t+1 appears to be constant, and so will be set to constant.
#>
#> Finished selecting best-fit models.
#summary(histtest)
The summary generated is quite long, and likely resulted in a series of
warnings from the model-building functions utilized. For our purposes,
we need only look at elements corresponding to the best-fit models for
our tested vital rates, which are the elements marked
survival_model, observation_model,
size_model, repstatus_model,
fecundity_model, juv_survival_model,
juv_observation _model, juv_size_model,
juv_reproduction_model, and
juv_maturity_model. The line beginning
Formula: in each of these sections shows the best-fit
model, in standard R formula notation (i.e. y = ax + b
is given as y ~ x). The independent terms tested include
variables coding for size in occasions t and t-1,
given as sizea2 and sizea1, and reproductive
status in occasions t and t-1, given as
repstatus2 and repstatus1, respectively. Since
several vital rates show sizea1 or repstatus1
as term in the best-fit models, most notably for adult survival, primary
size, reproduction probability, and fecundity, we see that individual
history has a significant impact on the demography of this population.
The quality control information in element qc is also of
interest, showing that some of our models are quite strong, with
accuracy of 0.982 for adult survival and 0.903 for observation, while
others are less accurate, most notably primary size and fecundity.
Step 4. MPM estimation
Now let’s create some raw Lefkovitch MPMs based on Ehrlén (2000). We have seen that history should
be included in these analyses, which justifies creating only historical
matrices. However, to introduce these functions in greater depth and
detail, we will first create ahistorical MPMs. The functions we will use
to build these MPMs include rlefko2() for raw ahistorical
MPMs, and rlefko3() for raw historical MPMs.
Ehrlén (2000) shows a mean matrix covering
years 1989 and 1990 as occasion t. We will utilize the entire
dataset instead, covering 1988 to 1991. Note that we will not create
matrices for subpopulations in this case (to include them, add the
options patch = "all", patchcol = "patchid" to the input
below).
ehrlen2 <- rlefko2(data = lathvert, stageframe = lathframe, year = "all",
stages = c("stage3", "stage2"), supplement = lathsupp2, yearcol = "year2",
indivcol = "individ")
#ehrlen2
The output from this analysis is a lefkoMat object, which
is list with the following elements:
A: a list of full population projection matrices, in
order of population, patch, and year (order given in
labels)
U: a list of matrices where non-zero entries are
limited to survival-transition elements, in the same order as
A
F: a list of matrices where non-zero entries are
limited to fecundity elements, in the same order as A
hstages: a data frame showing the order of paired
stages (given if matrices are historical; otherwise NA)
agestages: a data frame showing the order of
age-stage pairs (given if matrices are age-by-stage ahistorical;
otherwise NA)
ahstages: the stageframe used in analysis, with stages potentially reordered and edited as they occur in the matrix
labels: a table showing the order of matrices by population, patch, and year
matrixqc: a short vector used in
summary statements to describe the overall quality of each
matrix (used in summary() calls)
dataqc: a short vector used in summary
statements to describe key sampling aspects of the dataset (used in
summary() calls)
Input options for function rlefko2() include
year = "all", which can be changed to
year = c(1989, 1990) to focus just on years 1989 and 1990,
as in Ehrlén (2000), or
year = 1989 to focus exclusively on the transition from
1989 to 1990 (the year entered is the year in occasion t).
Matrix-estimating functions in lefko3 have a default
behavior of creating matrices for each year in the dataset except the
final year, rather than lumping all years together to produce a single
matrix. However, patches or subpopulations will only be separated if a
patch ID variable is provided as input.
We can understand lefkoMat objects in greater detail
through the summary() function.
We learn that three full A matrices and their
U and F decompositions were estimated, and
that they are ahistorical. This is expected given that there are four
consecutive years of data, yielding three time steps, and an ahistorical
matrix requires two consecutive years to estimate transitions. The
following line notes the dimensions of those matrices. The third,
fourth, fifth, and sixth lines of the summary show how many survival
transition and fecundity elements were actually estimated, both overall
and per matrix, the number of populations, patches, and time steps
covered by the MPM, and the number of individuals and transitions the
matrices are based on (these last numbers can be used to understand
matrix quality). Finally, the last section shows summaries of the column
sums from the survival-transition (U) matrices in this
lefkoMat object. These column sums correspond to the
survival probabilities of the different life stages, and so the
summaries must show numbers ranging from 0.0 to 1.0.
The matrix creation functions in this package sort the stages provided
in the stageframe according to a standardized rubric, and so the order
of stages in the ahstages element of a
lefkoMat object may differ from the original input.
Particularly, they sort propagules first, followed by immature stages,
followed by non-reproductive adult stages, followed by reproductive
adult stages. Let’s see the order in our matrices. We see that the key
difference from the original input is that the order of the flowering
stage and the vegetative dormancy stage have been flipped.
Now we’ll estimate historical matrices. Because of the size of these
matrices, we will only show the lefkoMat summary. We will
only create Ehrlén-format matrices - users wishing to create
deVries-format hMPMs should add format = "deVries" to the
input options (the resulting matrices will be bigger, but contain the
same number of estimated elements).
ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all",
stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3,
yearcol = "year2", indivcol = "individ")
#summary(ehrlen3)The summary output shows several differences from the ahistorical case. First, there is one less matrix estimated in the historical case than in the ahistorical case, because raw historical matrices require three consecutive occasions to estimate each transition instead of two. Second, these matrices are larger than ahistorical matrices, with the numbers of rows and columns generally equaling the number of ahistorical rows and columns squared (although deVries-format increases the dimensions through the addition of a prior stage for newborns, and in both Ehrlén and deVries formats the numbers of rows and columns may be reduced under some circumstances as shown in the next block). Finally, a much greater proportion of each matrix is composed of zeros in the historical case than in the ahistorical case, although there are certainly more non-zero elements as well. This sparseness results from historical matrices being composed primarily of structural zeros. As a result, the historical matrices in this example have non-zero entries in (79 + 7) / 2401 = 3.5% of matrix elements, while the equivalent ahistorical matrices have non-zero entries in (24.67 + 2) / 49 = 54.4% of matrix elements.
We can see the impact of structural zeros by eliminating some of them in
the process of matrix estimation. This can be done by setting
reduce = TRUE, which tells rlefko3() to
eliminate stage pairs in which both column and row are zero vectors.
Here, we now have matrices with 19 fewer rows and columns, and (79 + 7)
/ 900 = 9.6% of elements as potentially non-zero.
ehrlen3red <- rlefko3(data = lathvert, stageframe = lathframe,
year = "all", stages = c("stage3", "stage2", "stage1"),
supplement = lathsupp3, yearcol = "year2", indivcol = "individ",
reduce = TRUE)
#summary(ehrlen3red)
Next we will create the element-wise mean ahistorical matrix. Function
lmean() creates a lefkoMat object retaining
all descriptive information from the original lefkoMat
object. The full output not in any order includes the composite mean
matrix (shown in element A), the mean survival-transition
matrix (U), the mean fecundity matrix (F), a
data frame outlining the definitions and order of historical paired
stages (hstages, shown as NA in this case
because the matrices are ahistorical), a data frame outlining the
age-by-stage combinations used (agestages, shown as
NA because these are not age-by-stage matrices), a data
frame outlining the actual stages as outlined in the
stageframe object used to create these matrices
(ahstages), a data frame outlining the definitions and
order of the matrices (labels), and two quality control
vectors used in output for the summary() function
(matrixqc and dataqc).
Now we will estimate the historical element-wise mean matrix. We will show only the top-left corner of the rather large matrix (a section comprised of the first twenty rows and eight columns of the 30 x 30 matrix).
The prevalence of zeros in this matrix is normal because most elements
are structural zeros and so cannot equal anything else. This matrix is
also a raw matrix, meaning that transitions that do not actually occur
in the dataset cannot equal anything other than zero. Now let’s take a
look at the hstages object associated with this mean
matrix.
There are 49 pairs of life history stages, corresponding to the rows and columns of the historical matrices. The pairs are interpreted so that matrix columns represent stage pairs in occasions t-1 and t, and rows represent stage pairs in occasions t and t+1. For an element in the matrix to contain a number other than zero, it must represent the same stage at occasion t in both the column stage pairs and the row stage pairs. For example, The element [1, 1] represents the transition probability from dormant seed at occasions t-1 and t (column pair), to dormant seed at occasions t and t+1 (row pair) - the occasion t stages match, and so this element is possible. However, element [1, 2] represents the transition probability from seedling in occasion t-1 and very small adult in occasion t (column pair), to dormant seed in occasion t and in occasion t+1 (row pair). Clearly [1, 2] is a structural zero because it is impossible for an individual to be both a dormant seed and a very small adult at the same time.
Error-checking is more difficult with historical matrices because they
are typically one or two orders of magnitude bigger than their
ahistorical counterparts, but the same basic strategy can be used here
as with ahistorical matrices. In these cases we can use
summary() function to assess the key quality control
characteristics of the mean hMPM, such as the distribution of survival
probability estimates for historical stages.
Let’s try another approach, looking at some conditional historical
matrices. Conditional matrices are matrices of ahistorical dimension
that show transitions from stage at occasion t to occasion
t+1, conditional on all individuals having been in the same
stage in occasion t-1. They are calculated from historical
MPMs, and the output below shows all conditional matrices developed from
the first A matrix of ehrlen3mean. The first
matrix, for example, shows all transitions involving individuals that
had been in the dormant seed stage Sd in occasion
t-1, while the last matrix shows transitions involving
individuals that had been vegetatively dormant in occasion t-1.
Quick scans will show many transitions missing, because each stage has only certain stages that can transition from it, and to which it can transition. Further transitions are missing because the MPM is raw, and some transitions are not parameterized because no individuals made them.
One last error-checking technique before we analyze our MPMs: matrix
visualization plots. These plots provide a relatively easy way to
understand the “spatial spread” of values throughout a matrix. Package
lefko3 includes function image3(), which
provides an easy way to make these images. Let’s start off by looking at
the ahistorical mean matrix.
The resulting image shows non-zero elements as red spaces, and zero elements as white spaces. Rows and columns are numbered, and we can see that this matrix is reasonably dense.
Now let’s take a look at the mean historical matrix. The historical mean matrix has many more rows and columns, and has more of both zero and non-zero elements. However, it has become a sparse matrix, and it turns out that increasing the numbers of life stages will increase the sparsity of historical projection matrices.
Step 5. MPM analysis
Package lefko3 includes functions to conduct some analyses
of population dynamics. We will start by estimating the asymptotic
population growth rate (λ) and
the stochastic population growth rate (a = logλS)
from the ahistorical MPMs, including both the annual MPM and the mean.
For the stochastic case, we will set the seed for R’s random number
generation to make our output reproducible. Note that each λ estimate includes a data frame
describing the matrices in order (given as the labels
object within the output list). Here is the set of ahistorical annual
λ estimates, followed by λ for the mean matrix, and the
stochastic population growth rate (a = logλS).
# Deterministic ahistorical
lambda3(ehrlen2)
#> pop patch year2 lambda
#> 1 1 1 1988 0.8952585
#> 2 1 1 1989 0.9235493
#> 3 1 1 1990 1.0096490
# Deterministic mean ahistorical
lambda3(ehrlen2mean)
#> pop patch lambda
#> 1 1 1 0.9574162
# Stochastic ahistorical
set.seed(42)
slambda3(ehrlen2)
#> pop patch a var sd se
#> 1 1 1 -0.04490197 0.03154986 0.1776228 0.001776228We will now look at the same numbers for the historical analyses. There are several differences in the output in addition to the lower growth rate estimates. First, because there are four years of data, there are three ahistorical transitions possible for estimation: year 1 to 2, year 2 to 3, and year 3 to 4. However, in the historical case, only two are possible: from years 1 and 2 to 3 (technically, from years 1 [t-1] and 2 [t] to years 2 [t] and 3 [t+1]), and from years 2 and 3 to 4. Second, historical matrices cover more of the individual heterogeneity in a population by splitting ahistorical transitions by stage in occasion t-1. This heterogeneity may reflect many sources, for example the impacts of trade-offs operating across years (Shefferson and Roach 2010). One particularly common trade-off is the cost of growth: an individual that grows a great deal in one time step due to great environmental conditions in that year might pay a large cost of survival, growth, or reproduction in the next if those environmental conditions deteriorate (Shefferson, Warren II, and Pulliam 2014; Shefferson et al. 2018).
# Deterministic historical
lambda3(ehrlen3)
#> pop patch year2 lambda
#> 1 1 1 1989 0.8863920
#> 2 1 1 1990 0.9855435
# Deterministic mean historical
lambda3(ehrlen3mean)
#> pop patch lambda
#> 1 1 1 0.9182809
# Stochastic historical
set.seed(42)
slambda3(ehrlen3)
#> pop patch a var sd se
#> 1 1 1 -0.08872396 0.0177556 0.1332501 0.001332501We can also examine the stable stage distributions, as follows for the ahistorical case.
ehrlen2mss <- stablestage3(ehrlen2mean)
ehrlen2mss
#> matrix stage_id stage ss_prop
#> 1 1 1 Sd 0.29261073
#> 2 1 2 Sdl 0.04579994
#> 3 1 3 Tm 0.22875637
#> 4 1 4 Sm 0.18625613
#> 5 1 5 La 0.07716891
#> 6 1 7 Dorm 0.05588715
#> 7 1 6 Flo 0.11352077
The data frame output shows us the stages themselves
(stage, and associated number in stage_id),
which matrix they refer to (matrix), and the stable stage
distribution (ss_prop). Interpreting these values, we find
that the mean matrix suggests that, if we project the population forward
indefinitely assuming the population dynamics are static and represented
by this matrix, we will find that approximately 29% of individuals
should be dormant seeds (suggesting a large seedbank). A further 23% and
19% should be very small and small adults, respectively, and 11% should
be flowering adults. Almost 6% of the population should eventually be
composed of vegetatively dormant adults.
We can estimate the stable stage distribution for the historical case,
as well. Because the historical output for the
stablestage3() function is a list with two data frames,
let’s take a look at each of these data frames in turn. The first will
be the stage-pair output.
This data frame is structured in historical format, and so shows the
stable stage distribution of stage pairs. We may wish to see which stage
pair dominates, in which case we might look at the row with the maximum
ss_prop value.
ehrlen3mss$hist[which(ehrlen3mss$hist$ss_prop == max(ehrlen3mss$hist$ss_prop)),]
#> matrix stage_id_2 stage_id_1 stage_2 stage_1 ss_prop
#> 17 1 3 3 Tm Tm 0.2404459Here we see that about 24% of the population is expected to be composed of tiny adults maintaining themselves as tiny adults.
The longer format of the historical stable stage output makes it harder
to read. However, historical values can also be combined by stage at
occasion t (stage_2) to estimate the
historically-corrected stable stage distribution shown
in the ahist element, which allows comparison to a stable
stage distribution estimated from a purely ahistorical MPM. Notice below
that the ss_prop column shows values that are different
from the purely ahistorical case, suggesting the influence of individual
history.
ehrlen3mss$ahist
#> matrix stage_id stage ss_prop
#> 1 1 1 Sd 0.26275496
#> 2 1 2 Sdl 0.04112686
#> 3 1 3 Tm 0.29476860
#> 4 1 4 Sm 0.19794066
#> 5 1 5 La 0.06553412
#> 6 1 7 Dorm 0.04645369
#> 7 1 6 Flo 0.09142110
To see the impact of history on the stable stage distribution, let’s
plot the ahistorical and historically-corrected stable stage
distributions together. We will also include the stochastic long-run
stage distribution in our output, which is estimated with the same
function but using the stochastic = TRUE option. This will
allow us to see the impact of random temporal variation.
ehrlen2mss_s <- stablestage3(ehrlen2, stochastic = TRUE, seed = 42)
ehrlen3mss_s <- stablestage3(ehrlen3, stochastic = TRUE, seed = 42)
ss_put_together <- cbind.data.frame(ehrlen2mss$ss_prop, ehrlen3mss$ahist$ss_prop,
ehrlen2mss_s$ss_prop, ehrlen3mss_s$ahist$ss_prop)
names(ss_put_together) <- c("d_ahist", "d_hist", "s_ahist", "s_hist")
rownames(ss_put_together) <- ehrlen2mss$stage
barplot(t(ss_put_together), beside=T, ylab = "Proportion", xlab = "Stage",
ylim = c(0, 0.35), col = c("black", "orangered", "grey", "darkred"), bty = "n")
legend("topright", c("det ahistorical", "det historical", "sto ahistorical",
"sto historical"), pch = 22, col = "black",
pt.bg = c("black", "orangered", "grey", "darkred"), bty = "n")
Let’s take the deterministic portion first. Accounting for individual history increased the prevalence of tiny and small adults, but decreased the prevalence of dormant seeds, and dormant, large, and flowering adults. Now when we also take in the impact of temporal stochasticity, we can see differences in the proportion of dormant seeds, seedlings, and all adult stages, with the greatest differences in dormant seeds and tiny adults.
Let’s take a look at the reproductive values now, in similar order to the stable stage distribution case. Initially, we will create all sets of reproductive value objects, and then we will plot them. The structure of these objects is the same as that of the stable stage structure outputs. Because the four vectors are all standardized such that the first non-zero reproductive value is set to 1.0, they are on different scales, and so we will make them comparable for plotting purposes by standardizing them against their vector sums.
ehrlen2mrv <- repvalue3(ehrlen2mean)
ehrlen3mrv <- repvalue3(ehrlen3mean)
ehrlen2mrv_s <- repvalue3(ehrlen2, stochastic = TRUE, seed = 42)
ehrlen3mrv_s <- repvalue3(ehrlen3red, stochastic = TRUE, seed = 42)
rv_put_together <- cbind.data.frame((ehrlen2mrv$rep_value / sum(ehrlen2mrv$rep_value)),
(ehrlen3mrv$ahist$rep_value / sum(ehrlen3mrv$ahist$rep_value)),
(ehrlen2mrv_s$rep_value / sum(ehrlen2mrv_s$rep_value)),
(ehrlen3mrv_s$ahist$rep_value / sum(ehrlen3mrv_s$ahist$rep_value)))
names(rv_put_together) <- c("det ahist", "det hist", "sto ahist", "sto hist")
rownames(rv_put_together) <- ehrlen2mrv$stage
barplot(t(rv_put_together), beside=T, ylab = "Relative reproductive value",
ylim = c(0, 0.4), xlab = "Stage", col = c("black", "orangered", "grey", "darkred"),
bty = "n")
legend("topleft", c("det ahistorical", "det historical", "sto ahistorical",
"sto historical"), pch = 22, col = "black",
pt.bg = c("black", "orangered", "grey", "darkred"), bty = "n")
Both deterministic and stochastic analyses show that flowering adults have the greatest reproductive value in both ahistorical and historical analysis, while dormant seeds have the least. However, the historical MPMs suggest lower contributions of seedlings and vegetative dormancy, but larger contributions of flowering adults.
Now we will look at the deterministic sensitivities of λ to the ahistorical mean matrix elements.
ehrlen2sens <- sensitivity3(ehrlen2mean)
#> Running deterministic analysis...
print(ehrlen2sens$ah_sensmats[[1]], digits = 3)
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#> [1,] 0.0182 0.00284 0.0142 0.0116 0.00479 0.00347 0.00705
#> [2,] 0.2061 0.03225 0.1611 0.1312 0.05434 0.03936 0.07994
#> [3,] 0.2565 0.04016 0.2006 0.1633 0.06766 0.04900 0.09953
#> [4,] 0.4110 0.06432 0.3213 0.2616 0.10838 0.07849 0.15943
#> [5,] 0.5639 0.08826 0.4408 0.3589 0.14871 0.10770 0.21877
#> [6,] 0.3067 0.04801 0.2398 0.1952 0.08089 0.05858 0.11899
#> [7,] 0.7221 0.11302 0.5645 0.4596 0.19043 0.13791 0.28014
# The highest sensitivity value is:
max(ehrlen2sens$ah_sensmats[[1]])
#> [1] 0.7220817
# This occurs in element:
which(ehrlen2sens$ah_sensmats[[1]] == max(ehrlen2sens$ah_sensmats[[1]]))
#> [1] 7
# The highest sensitivity value among biologically plausible elements is:
max(ehrlen2sens$ah_sensmats[[1]][which(ehrlen2mean$A[[1]] > 0)])
#> [1] 0.4596282
# This occurs in element:
which(ehrlen2sens$ah_sensmats[[1]] ==
max(ehrlen2sens$ah_sensmats[[1]][which(ehrlen2mean$A[[1]] > 0)]))
#> [1] 28The highest sensitivity value is associated with a biologically impossible transition - dormant seeds (stage/column 1) cannot transition to flowering (stage/row 7) (element 7). Among biologically plausible elements, the highest sensitivity is associated with element 28, which is the transition from small adult (stage/column 4) to flowering (stage/row 7).
We will now look at the sensitivity of λ to elements in the historical mean MPM.
The first element produced in this analysis is h_sensmats,
which is a list composed of sensitivity matrices of the historical
matrices, in order (we have only one in our mean matrix object). These
matrices are the same dimensions as the historical matrices used as
input, and so can be quite huge. This is followed by
ah_sensmats, which is a list composed of
historically-corrected sensitivity matrices of corresponding ahistorical
matrix elements (calculated using the historically-corrected stable
stage distribution and reproductive value vector produced in
ehrlen3mss and ehrlen3mrv, respectively). So,
these are different than the sensitivities estimated from the
ahistorical matrices themselves, but have the same dimensions. Next,
h_stages and ah_stages give the order of
paired stages and life history stages used in the historical and
historically-corrected sensitivity matrices, respectively. Finally, we
have the original A, U, and F
matrices used as input.
Our historical matrices are large and full of zeros. So, we will look for the highest sensitivity associated with a biologically plausible element (i.e. non-zero matrix elements). Then, we will assess the highest biologically plausible sensitivity in the historically-corrected sensitivity matrices, to compare against the ahistorical sensitivity analysis.
# The highest sensitivity value among biologically plausible elements:
max(ehrlen3sens$h_sensmats[[1]][which(ehrlen3mean$A[[1]] > 0)])
#> [1] 0.3148525
# This value is associated with element:
which(ehrlen3sens$h_sensmats[[1]] ==
max(ehrlen3sens$h_sensmats[[1]][which(ehrlen3mean$A[[1]] > 0)]))
#> [1] 802
# Highest historically-corrected sensitivity value among plausible elements is:
max(ehrlen3sens$ah_sensmats[[1]][which(ehrlen2mean$A[[1]] > 0)])
#> [1] 0.8591236
# This occurs in element:
which(ehrlen3sens$ah_sensmats[[1]] ==
max(ehrlen3sens$ah_sensmats[[1]][which(ehrlen2mean$A[[1]] > 0)]))
#> [1] 20The maximum biologically plausible sensitivity value in the historical matrix is element 802, which is associated with column 17 (tiny adult in occasions t-1 and t) and row 18 (tiny adult in occasion t to small adult in occasion t+1). This transition is from tiny adult in occasions t-1 and t to small adult in t+1. The historically-corrected sensitivity analysis finds that λ is most sensitive to element 20, which is the transition from tiny adult to dormant. This is a little different from the ahistorical MPM, which suggested element 28 (small adult to flowering adult).
We can perform the same analyses as above with stochastic sensitivities.
In that circumstance, simply run the sensitivity3()
function with the stochastic = TRUE argument set, as below.
ehrlen2sens_s <- sensitivity3(ehrlen2, stochastic = TRUE)
#> Running stochastic analysis...
ehrlen3sens_s <- sensitivity3(ehrlen3, stochastic = TRUE)
#> Running stochastic analysis...A complementary approach to sensitivity analysis is elasticity analysis. Elasticities are easier to interpret because projection matrix elements valued at 0 produce elasticity values also equal to 0, thus eliminating biologically impossible transitions from consideration. Elasticity values are also scaled to sum to 1.0, making elasticities of survival transitions easier to compare to those of fecundity. However, they are also interpreted differently, because while sensitivity analysis shows the impact of a tiny but absolute change to a matrix element on λ, elasticity analysis shows the impact of a tiny but proportional change to a matrix element on λ. In fact, both sensitivities and elasticities are essentially local slopes, and so are not unit free. It is therefore not unusual for sensitivity and elasticity analysis to yield different inferences.
Let’s look at the elasticity of λ to matrix elements in the ahistorical mean matrix.
ehrlen2elas <- elasticity3(ehrlen2mean)
#> Running deterministic analysis...
print(ehrlen2elas$ah_elasmats, digits = 3)
#> [[1]]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#> [1,] 0.00655 0.00000 0.000000 0.0000 0.00000 0.00000 0.0116
#> [2,] 0.01162 0.00000 0.000000 0.0000 0.00000 0.00000 0.0206
#> [3,] 0.00000 0.02784 0.151678 0.0164 0.00052 0.00415 0.0000
#> [4,] 0.00000 0.00127 0.041102 0.1586 0.02210 0.02184 0.0167
#> [5,] 0.00000 0.00000 0.000604 0.0290 0.04851 0.01744 0.0532
#> [6,] 0.00000 0.00315 0.007180 0.0223 0.00896 0.00628 0.0107
#> [7,] 0.00000 0.00000 0.000000 0.0353 0.06862 0.00887 0.1673
# The maximum elasticity value:
max(ehrlen2elas$ah_elasmats[[1]])
#> [1] 0.167339
# This value is associated with element:
which(ehrlen2elas$ah_elasmats[[1]] == max(ehrlen2elas$ah_elasmats[[1]]))
#> [1] 49Elasticity analysis exhibits strong differences from sensitivity analysis. In particular, we find that λ is most strongly elastic in response to changes in element 49, which is the stasis transition for flowering adults (stage 7). We can sum the columns of the elasticity matrix to see which stages λ is most and least elastic in response to, as below.
print(colSums(ehrlen2elas$ah_elasmats[[1]]), digits = 3)
#> [1] 0.0182 0.0323 0.2006 0.2616 0.1487 0.0586 0.2801Here we see that λ is most strongly elastic in response to changes in transitions associated with flowering adults, followed by transitions involving small adults. Dormant seeds and seedlings have the smallest impact on λ, and the impacts of fecundity (shown in the top-right corner of the elasticity matrix) appear quite small.
Now on to elasticity analysis of the historical MPMs. Once again, we
will not output the matrices. Type ehrlen3elas at the
prompt to see these matrices.
ehrlen3elas <- elasticity3(ehrlen3mean)
#> Running deterministic analysis...
# The highest deterministic elasticity value:
max(ehrlen3elas$h_elasmats[[1]])
#> [1] 0.1718096
# This value is associated with element:
which(ehrlen3elas$h_elasmats[[1]] == max(ehrlen3elas$h_elasmats[[1]]))
#> [1] 801The highest elasticity appears to be associated with the element 801, which is at row 17, column 17. This corresponds to the stasis transition of tiny adults (tiny in t-1 to tiny in t to tiny in t+1).
Elasticities are often treated as additive, making the calculation of
historically-corrected elasticity matrices easy. These are stored in the
ah_elasmats element of elasticity3() output
originating from a historical MPM. Eyeballing this matrix, it appears
that the historically-corrected elasticity matrix supports stasis in the
tiny adult stage as the transition that λ is most elastic to, with stasis in
vegetative dormancy a close second.
print(ehrlen3elas$ah_elasmats, digits = 3)
#> [[1]]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7]
#> [1,] 0.00699 0.000000 0.00000 0.0000 0.00e+00 0.01161 0.00000
#> [2,] 0.01161 0.000000 0.00000 0.0000 0.00e+00 0.00000 0.00000
#> [3,] 0.00000 0.010736 0.19877 0.0344 8.26e-05 0.00000 0.00219
#> [4,] 0.00000 0.000292 0.04500 0.1615 2.56e-02 0.02781 0.01907
#> [5,] 0.00000 0.000000 0.00000 0.0306 4.54e-02 0.05418 0.00920
#> [6,] 0.00000 0.000000 0.00000 0.0366 6.09e-02 0.16395 0.00497
#> [7,] 0.00000 0.000583 0.00238 0.0162 7.42e-03 0.00889 0.00315Next, we will create a barplot of the elasticities of life history stages from ahistorical vs. historically-corrected analyses. We will also incorporate stochastic elasticity analysis here to assess the importance of temporal environmental stochasticity on population growth.
ehrlen2elas_s <- elasticity3(ehrlen2, stochastic = TRUE)
#> Running stochastic analysis...
ehrlen3elas_s <- elasticity3(ehrlen3, stochastic = TRUE)
#> Running stochastic analysis...
elas_put_together <- cbind.data.frame(colSums(ehrlen2elas$ah_elasmats[[1]]),
colSums(ehrlen3elas$ah_elasmats[[1]]), colSums(ehrlen2elas_s$ah_elasmats[[1]]),
colSums(ehrlen3elas_s$ah_elasmats[[1]]))
names(elas_put_together) <- c("det ahist", "det hist", "sto ahist", "sto hist")
rownames(elas_put_together) <- ehrlen2elas$ahstages$stage
barplot(t(elas_put_together), beside=T, ylab = "Elasticity", xlab = "Stage",
col = c("black", "orangered", "grey", "darkred"), bty = "n")
legend("topleft", c("det ahistorical", "det historical", "sto ahistorical",
"sto historical"), pch = 22, col = "black",
pt.bg = c("black", "orangered", "grey", "darkred"), bty = "n")
Historical analyses generally find that population growth rate is less elastic in response to seedlings, and flowering adults, and more elastic to tiny and dormant adults than in ahistorical analyses. Stochastic and deterministic population growth rates are barely elastic in response to dormant seeds and seedlings. Note the dramatic impact of environmental stochasticity and individual history combined on the impact of tiny adults.
Let’s now look at the elasticities of different kinds of transitions. We
will use the summary() function, which outputs data frames
summarizing elasticity sums by the kind of transition. First, we will
compare ahistorical against historically-corrected transitions.
ehrlen2elas_sums <- summary(ehrlen2elas)
ehrlen3elas_sums <- summary(ehrlen3elas)
ehrlen2elas_s_sums <- summary(ehrlen2elas_s)
ehrlen3elas_s_sums <- summary(ehrlen3elas_s)
elas_sums_together <- cbind.data.frame(ehrlen2elas_sums$ahist[,2],
ehrlen3elas_sums$ahist[,2], ehrlen2elas_s_sums$ahist[,2],
ehrlen3elas_s_sums$ahist[,2])
names(elas_sums_together) <- c("det ahist", "det hist", "sto ahist", "sto hist")
rownames(elas_sums_together) <- ehrlen2elas_sums$ahist$category
barplot(t(elas_sums_together), beside=T, ylab = "Elasticity",
xlab = "Transition", col = c("black", "orangered", "grey", "darkred"), bty = "n")
legend("topright", c("det ahistorical", "det historical", "sto ahistorical",
"sto historical"), pch = 22, col = "black",
pt.bg = c("black", "orangered", "grey", "darkred"), bty = "n")
We see similar patterns across types of elasticity analysis. Particularly, population growth rate is most elastic in response to changes in stasis transitions, and least elastic to changes in fecundity. Environmental stochasticity appears to exacerbate these differences.
Package lefko3 also includes functions to conduct many
other analyses, including deterministic and stochastic life table
response experiments, and general projection including quasi-extinction
analysis and density dependent analysis. Users wishing to conduct these
analyses should see our free e-manual called lefko3: a
gentle introduction and other vignettes, including
long-format and video vignettes, on
the projects
page of the Shefferson lab website.
Acknowledgements
We are grateful to two anonymous reviewers whose scrutiny improved the quality of this vignette. The project resulting in this package and this tutorial was funded by Grant-In-Aid 19H03298 from the Japan Society for the Promotion of Science.