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 adapt3 will be made available over
time. Please check available resources on
the
projects page of the {r}evolutionary demography
website.
Pairwise invasibility analysis focused on matrix element
manipulation
Our first example will focus on running pairwise invasibility analyses
in which matrix elements are manipulated to create our pairs of
residents and mutants. Because we use existing MPMs, the process
generally proceeds quite quickly, and only slows when using giant
matrices such as might occur in historical MPM projection. To make this
process easier, we will develop a relatively small, stage-based
empirical MPM.
Per lefko3 terminology, we will develop a stageframe that
encapsulates our life history model by noting all of the relevant
characteristics of each stage, as below.
sizevector <- c(0, 0, 0, 0, 1, 2.5, 4.5, 8, 17.5)
stagevector <- c("SD", "P1", "SL", "D", "XSm", "Sm", "Md", "Lg", "XLg")
repvector <- c(0, 0, 0, 0, 1, 1, 1, 1, 1)
obsvector <- c(0, 0, 0, 0, 1, 1, 1, 1, 1)
matvector <- c(0, 0, 0, 1, 1, 1, 1, 1, 1)
immvector <- c(0, 1, 1, 0, 0, 0, 0, 0, 0)
propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0)
indataset <- c(0, 0, 0, 1, 1, 1, 1, 1, 1)
binvec <- c(0, 0, 0, 0.5, 0.5, 1, 1, 2.5, 7)
cypframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector,
repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
propstatus = propvector, immstatus = immvector, indataset = indataset,
binhalfwidth = binvec)
cypframe_raw
#> stage size size_b size_c min_age max_age repstatus obsstatus propstatus immstatus matstatus indataset binhalfwidth_raw sizebin_min
#> 1 SD 0.0 NA NA NA NA 0 0 1 0 0 0 0.0 0.0
#> 2 P1 0.0 NA NA NA NA 0 0 0 1 0 0 0.0 0.0
#> 3 SL 0.0 NA NA NA NA 0 0 0 1 0 0 0.0 0.0
#> 4 D 0.0 NA NA NA NA 0 0 0 0 1 1 0.5 -0.5
#> 5 XSm 1.0 NA NA NA NA 1 1 0 0 1 1 0.5 0.5
#> 6 Sm 2.5 NA NA NA NA 1 1 0 0 1 1 1.0 1.5
#> 7 Md 4.5 NA NA NA NA 1 1 0 0 1 1 1.0 3.5
#> 8 Lg 8.0 NA NA NA NA 1 1 0 0 1 1 2.5 5.5
#> 9 XLg 17.5 NA NA NA NA 1 1 0 0 1 1 7.0 10.5
#> sizebin_max sizebin_center sizebin_width binhalfwidthb_raw sizebinb_min sizebinb_max sizebinb_center sizebinb_width
#> 1 0.0 0.0 0 NA NA NA NA NA
#> 2 0.0 0.0 0 NA NA NA NA NA
#> 3 0.0 0.0 0 NA NA NA NA NA
#> 4 0.5 0.0 1 NA NA NA NA NA
#> 5 1.5 1.0 1 NA NA NA NA NA
#> 6 3.5 2.5 2 NA NA NA NA NA
#> 7 5.5 4.5 2 NA NA NA NA NA
#> 8 10.5 8.0 5 NA NA NA NA NA
#> 9 24.5 17.5 14 NA NA NA NA NA
#> binhalfwidthc_raw sizebinc_min sizebinc_max sizebinc_center sizebinc_width group comments
#> 1 NA NA NA NA NA 0 No description
#> 2 NA NA NA NA NA 0 No description
#> 3 NA NA NA NA NA 0 No description
#> 4 NA NA NA NA NA 0 No description
#> 5 NA NA NA NA NA 0 No description
#> 6 NA NA NA NA NA 0 No description
#> 7 NA NA NA NA NA 0 No description
#> 8 NA NA NA NA NA 0 No description
#> 9 NA NA NA NA NA 0 No description
Next we will create our verticalized dataset, in which our horizontal
dataset is restructured to have all data for each individual split into
blocks of three consecutive monitoring occasions, as below. Note that,
because we are using the stageassign argument, we can have
R assign stages to each individual in each occasion. Note
also our use of the summary_hfv() function, which provides
specialized, useful output for this style of data frame.
cypraw_v1 <- verticalize3(data = cypa_data, noyears = 18, firstyear = 1994,
individcol = "plant_id", blocksize = 2, sizeacol = "Inf.94",
sizebcol = "Veg.94", repstracol = "Inf.94", fecacol = "Inf.94",
stageassign = cypframe_raw, stagesize = "sizeadded", NAas0 = TRUE,
NRasRep = TRUE)
summary_hfv(cypraw_v1)
#>
#> This hfv dataset contains 6871 rows, 51 variables, 1 population,
#> 1 patch, 1025 individuals, and 17 time steps.
Now we will create our MPM. Later on, this will be projected through our
invasibility analyses, so we need to be sure that it is constructed
correctly. See our online book,
lefko3:
a gentle introduction, for further material on how to construct MPMs
properly.
cypsupp2r <- supplemental(stage3 = c("SD", "P1", "SL", "D", "XSm", "Sm", "SD",
"P1"),
stage2 = c("SD", "SD", "P1", "SL", "SL", "SL", "rep", "rep"),
eststage3 = c(NA, NA, NA, "D", "XSm", "Sm", NA, NA),
eststage2 = c(NA, NA, NA, "XSm", "XSm", "XSm", NA, NA),
givenrate = c(0.10, 0.20, 0.25, NA, NA, NA, NA, NA),
multiplier = c(NA, NA, NA, NA, NA, NA, 5000, 5000),
type =c(1, 1, 1, 1, 1, 1, 3, 3),
stageframe = cypframe_raw, historical = FALSE)
#> Warning: NA values in argument multiplier will be treated as 1 values.
cypmatrix2r <- rlefko2(data = cypraw_v1, stageframe = cypframe_raw,
year = "all", patch = "all", stages = c("stage3", "stage2", "stage1"),
size = c("size3added", "size2added"), supplement = cypsupp2r,
yearcol = "year2", patchcol = "patchid", indivcol = "individ")
cypmean <- lmean(cypmatrix2r)
summary(cypmean)
#>
#> This ahistorical lefkoMat object contains 1 matrix.
#>
#> Each matrix is square with 9 rows and columns, and a total of 81 elements.
#> A total of 31 survival transitions were estimated, with 31 per matrix.
#> A total of 8 fecundity transitions were estimated, with 8 per matrix.
#> This lefkoMat object covers 1 population, 1 patch, and 0 time steps.
#>
#> The dataset contains a total of 1025 unique individuals and 6871 unique transitions.
#>
#> Survival probability sum check (each matrix represented by column in order):
#> [,1]
#> Min. 0.000
#> 1st Qu. 0.300
#> Median 0.783
#> Mean 0.594
#> 3rd Qu. 0.858
#> Max. 0.941
#>
#>
#> Stages without connections leading to the rest of the life cycle found in matrices: 1
Now that we have our core MPM, we will need to clarify some of the
assumptions of the projection. First, let’s set up a starting vector
residents and mutants. We will use the typical setting in this case,
assuming that both variants will start with a single individual. We will
make that individual a protocorm, since that is the youngest
post-germination individual possible.
cyp_start <- start_input(cypmean, stage2 = "P1", value = 1)
cyp_start
#> stage2 stage_id_2 stage1 stage_id_1 age2 row_num value
#> 1 P1 2 NA NA NA 2 1
Next we need to clarify how density dependence will operate on the
demography of this population. In herbaceous perennial plants, the
germination stage is most commonly negatively density dependent. In
orchids, germination may be from a fresh seed, or from a previously
dormant seed, and we will assume density dependence in both. Determining
the exact nature and proper parameterization of density dependence is
too great a topic for this vignette, and so we utilize the following
two-parameter Ricker function-based approach, which we know works in
this case based on our own previous studies.
c2d_4 <- density_input(cypmean, stage3 = c("P1", "P1"), stage2= c("SD", "rep"),
style = 2, time_delay = 1, alpha = 0.005, beta = 0.000005, type = c(2, 2))
c2d_4
#> stage3 stage2 stage1 age2 style alpha beta gamma time_delay type type_t12
#> 1 P1 SD NA NA 2 0.005 5e-06 0 1 2 1
#> 2 P1 XSm NA NA 2 0.005 5e-06 0 1 2 1
#> 3 P1 Sm NA NA 2 0.005 5e-06 0 1 2 1
#> 4 P1 Md NA NA 2 0.005 5e-06 0 1 2 1
#> 5 P1 Lg NA NA 2 0.005 5e-06 0 1 2 1
#> 6 P1 XLg NA NA 2 0.005 5e-06 0 1 2 1
Next, we need to propose the variants to test. The information regarding
the variants is probably the most important information to provide to R,
because it clarifies the trait(s) under investigation. In the following
data frame, which we call a trait axis, we define 15 genetic variants,
each of which differs in fecundity multipliers. In a pairwise
invasibility analysis, each possible pair permutation will be projected,
with one serving as the resident population and the other as the mutant
invader.
cyp_ta <- trait_axis(stageframe = cypframe_raw, stage3 = rep("P1", 15),
stage2 = rep("rep", 15), multiplier = seq(from=1, to=3, length.out = 15),
type = rep(2, 15))
cyp_ta
#> variant stage3 stage2 stage1 age3 age2 eststage3 eststage2 eststage1 estage3 estage2 givenrate offset multiplier convtype
#> 1 1 P1 rep <NA> NA NA <NA> <NA> <NA> NA NA NA NA 1.000000 2
#> 2 2 P1 rep <NA> NA NA <NA> <NA> <NA> NA NA NA NA 1.142857 2
#> 3 3 P1 rep <NA> NA NA <NA> <NA> <NA> NA NA NA NA 1.285714 2
#> 4 4 P1 rep <NA> NA NA <NA> <NA> <NA> NA NA NA NA 1.428571 2
#> 5 5 P1 rep <NA> NA NA <NA> <NA> <NA> NA NA NA NA 1.571429 2
#> 6 6 P1 rep <NA> NA NA <NA> <NA> <NA> NA NA NA NA 1.714286 2
#> 7 7 P1 rep <NA> NA NA <NA> <NA> <NA> NA NA NA NA 1.857143 2
#> 8 8 P1 rep <NA> NA NA <NA> <NA> <NA> NA NA NA NA 2.000000 2
#> 9 9 P1 rep <NA> NA NA <NA> <NA> <NA> NA NA NA NA 2.142857 2
#> 10 10 P1 rep <NA> NA NA <NA> <NA> <NA> NA NA NA NA 2.285714 2
#> 11 11 P1 rep <NA> NA NA <NA> <NA> <NA> NA NA NA NA 2.428571 2
#> 12 12 P1 rep <NA> NA NA <NA> <NA> <NA> NA NA NA NA 2.571429 2
#> 13 13 P1 rep <NA> NA NA <NA> <NA> <NA> NA NA NA NA 2.714286 2
#> 14 14 P1 rep <NA> NA NA <NA> <NA> <NA> NA NA NA NA 2.857143 2
#> 15 15 P1 rep <NA> NA NA <NA> <NA> <NA> NA NA NA NA 3.000000 2
#> convtype_t12 surv_dev obs_dev size_dev sizeb_dev sizec_dev repst_dev fec_dev jsurv_dev jobs_dev jsize_dev jsizeb_dev jsizec_dev
#> 1 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> 2 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> 3 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> 4 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> 5 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> 6 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> 7 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> 8 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> 9 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> 10 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> 11 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> 12 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> 13 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> 14 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> 15 NA NA NA NA NA NA NA NA NA NA NA NA NA
#> jrepst_dev jmatst_dev indcova indcovb indcovc
#> 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
Finally, we will run our invasibility analysis.
cyp_inv <- invade3(axis = cyp_ta, mpm = cypmean, density = c2d_4, times = 350,
starts = cyp_start, entry_time = c(0, 250), fitness_times = 30,
var_per_run = 2)
Let’s see some diagnostics, in particular by plotting the sizes of the
resident and invader subpopulations in each time step. We will look at
two separate runs.
par(mfrow = c(2, 1))
plot(cyp_inv, pip = FALSE, run = 25)
plot(cyp_inv, pip = FALSE, run = 50)
Figure I1. Population growth of resident (solid
black) and invader (dashed red) variants in two runs
In the plots above, the resident variant is given in solid black, while
the invader variant, which is introduced at time 250, is given as a
dashed red line. It appears that the invader never gained a foothold in
the population, likely going extinct while the resident grew to a stable
population size. However, we can also take a peek at the fitness values
in the pairwise simulations, as below.
summary(cyp_inv$fitness)
#> simulation_num rep variant1 variant2 entrytime1 entrytime2 fitness_variant1 fitness_variant2
#> Min. : 1 Min. :1 Min. : 1 Min. : 1 Min. :0 Min. :250 Min. :-64933.135 Min. : -1.343
#> 1st Qu.: 57 1st Qu.:1 1st Qu.: 4 1st Qu.: 4 1st Qu.:0 1st Qu.:250 1st Qu.: -118.362 1st Qu.: -0.660
#> Median :113 Median :1 Median : 8 Median : 8 Median :0 Median :250 Median : 0.000 Median : 0.000
#> Mean :113 Mean :1 Mean : 8 Mean : 8 Mean :0 Mean :250 Mean : -2442.037 Mean : 5880.032
#> 3rd Qu.:169 3rd Qu.:1 3rd Qu.:12 3rd Qu.:12 3rd Qu.:0 3rd Qu.:250 3rd Qu.: 0.928 3rd Qu.: 162.300
#> Max. :225 Max. :1 Max. :15 Max. :15 Max. :0 Max. :250 Max. : 1.531 Max. :180963.650
The data frame that we see the summary of above shows the fitness for
both the resident (variant1) and the invader
(variant2). Note that the invaders’ fitness, assessed as
the Lyapunov coefficient for the final 30 times in each simulation,
shows that some invaders actually had very high growth. These are likely
cases in which low fecundity multiplier residents were run with high
fecundity multiplier invaders. Let’s view the pairwise invasibility plot
(PIP) to see if an ESS resulted from this. This will give us a sense of
whether we have any evolutionarily stable equilibria, or whether natural
selection appears to completely directional.
Figure I2. Pairwise invasibility plot showing
directional selection for increased reproduction. Black and white
regions show positive and negative invader fitness, respectively
It appears that natural selection is completely directional here,
favoring increasingly large fruiting value multipliers.
Pairwise invasibility analysis focused on vital rate model
manipulation
Let’s now move on to looking at analyzing function-based MPMs. This
approach works best when vital rates are given as linear or non-linear
functions of size, reproductive status, and other variables. Here, the
linear or non-linear models can be manipulated by altering their
y-intercepts.
Our example will once again focus on Cypripedium parviflorum.
However, this time we will attempt to replicate the pairwise
invasibility analysis results from the ahistorical, deterministic
analysis of the same population published in Shefferson, Warren II, and Pulliam (2014). That
analysis was used to explore the adaptive context of vegetative
dormancy, and hence imperfect sprouting, and found an evolutionarily
stable strategy at an intermediate level of sprouting. The dataset that
we will use is not exactly the same as in the paper, nor is the
methodology exactly the same, but our results should be similar. Because
we will use a function-based approach, which is very similar to building
an integral projection model (IPM), we will use a larger resolution life
history model than used in that paper, allowing for all observed sizes
of both vegetative and flowering adults, where size is defined by the
number of sprouts.
stagevector <- c("SD", "P1", "SL", "D", "V1", "V2", "V3", "V4", "V5", "V6",
"V7", "F1", "F2", "F3", "F4", "F5", "F6", "F7")
indataset <- c(0, 0, 0, 1, rep(1, 7), rep(1, 7))
sizevector <- c(0, 0, 0, 0, c(1:7), c(1:7))
repvector <- c(0, 0, 0, 0, rep(0, 7), rep(1, 7))
obsvector <- c(0, 0, 0, 0, rep(1, 7), rep(1, 7))
matvector <- c(0, 0, 0, 1, rep(1, 7), rep(1, 7))
immvector <- c(0, 1, 1, 0, rep(0, 7), rep(0, 7))
propvector <- c(1, 0, 0, 0, rep(0, 7), rep(0, 7))
comments <- c("Dormant seed", "Protocorm", "Seedling", "Veg dorm",
"Veg adult 1 stem", "Veg adult 2 stems", "Veg adult 3 stems",
"Veg adult 4 stems", "Veg adult 5 stems", "Veg adult 6 stems",
"Veg adult 7 stems", "Flo adult 1 stem", "Flo adult 2 stems",
"Flo adult 3 stems", "Flo adult 4 stems", "Flo adult 5 stems",
"Flo adult 6 stems", "Flo adult 7 stems")
cypframe_fb <- sf_create(sizes = sizevector, stagenames = stagevector,
repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
propstatus = propvector, immstatus = immvector, indataset = indataset,
comments = comments)
cypframe_fb
#> stage size size_b size_c min_age max_age repstatus obsstatus propstatus immstatus matstatus indataset binhalfwidth_raw sizebin_min
#> 1 SD 0 NA NA NA NA 0 0 1 0 0 0 0.5 -0.5
#> 2 P1 0 NA NA NA NA 0 0 0 1 0 0 0.5 -0.5
#> 3 SL 0 NA NA NA NA 0 0 0 1 0 0 0.5 -0.5
#> 4 D 0 NA NA NA NA 0 0 0 0 1 1 0.5 -0.5
#> 5 V1 1 NA NA NA NA 0 1 0 0 1 1 0.5 0.5
#> 6 V2 2 NA NA NA NA 0 1 0 0 1 1 0.5 1.5
#> 7 V3 3 NA NA NA NA 0 1 0 0 1 1 0.5 2.5
#> 8 V4 4 NA NA NA NA 0 1 0 0 1 1 0.5 3.5
#> 9 V5 5 NA NA NA NA 0 1 0 0 1 1 0.5 4.5
#> 10 V6 6 NA NA NA NA 0 1 0 0 1 1 0.5 5.5
#> 11 V7 7 NA NA NA NA 0 1 0 0 1 1 0.5 6.5
#> 12 F1 1 NA NA NA NA 1 1 0 0 1 1 0.5 0.5
#> 13 F2 2 NA NA NA NA 1 1 0 0 1 1 0.5 1.5
#> 14 F3 3 NA NA NA NA 1 1 0 0 1 1 0.5 2.5
#> 15 F4 4 NA NA NA NA 1 1 0 0 1 1 0.5 3.5
#> 16 F5 5 NA NA NA NA 1 1 0 0 1 1 0.5 4.5
#> 17 F6 6 NA NA NA NA 1 1 0 0 1 1 0.5 5.5
#> 18 F7 7 NA NA NA NA 1 1 0 0 1 1 0.5 6.5
#> sizebin_max sizebin_center sizebin_width binhalfwidthb_raw sizebinb_min sizebinb_max sizebinb_center sizebinb_width
#> 1 0.5 0 1 NA NA NA NA NA
#> 2 0.5 0 1 NA NA NA NA NA
#> 3 0.5 0 1 NA NA NA NA NA
#> 4 0.5 0 1 NA NA NA NA NA
#> 5 1.5 1 1 NA NA NA NA NA
#> 6 2.5 2 1 NA NA NA NA NA
#> 7 3.5 3 1 NA NA NA NA NA
#> 8 4.5 4 1 NA NA NA NA NA
#> 9 5.5 5 1 NA NA NA NA NA
#> 10 6.5 6 1 NA NA NA NA NA
#> 11 7.5 7 1 NA NA NA NA NA
#> 12 1.5 1 1 NA NA NA NA NA
#> 13 2.5 2 1 NA NA NA NA NA
#> 14 3.5 3 1 NA NA NA NA NA
#> 15 4.5 4 1 NA NA NA NA NA
#> 16 5.5 5 1 NA NA NA NA NA
#> 17 6.5 6 1 NA NA NA NA NA
#> 18 7.5 7 1 NA NA NA NA NA
#> binhalfwidthc_raw sizebinc_min sizebinc_max sizebinc_center sizebinc_width group comments
#> 1 NA NA NA NA NA 0 Dormant seed
#> 2 NA NA NA NA NA 0 Protocorm
#> 3 NA NA NA NA NA 0 Seedling
#> 4 NA NA NA NA NA 0 Veg dorm
#> 5 NA NA NA NA NA 0 Veg adult 1 stem
#> 6 NA NA NA NA NA 0 Veg adult 2 stems
#> 7 NA NA NA NA NA 0 Veg adult 3 stems
#> 8 NA NA NA NA NA 0 Veg adult 4 stems
#> 9 NA NA NA NA NA 0 Veg adult 5 stems
#> 10 NA NA NA NA NA 0 Veg adult 6 stems
#> 11 NA NA NA NA NA 0 Veg adult 7 stems
#> 12 NA NA NA NA NA 0 Flo adult 1 stem
#> 13 NA NA NA NA NA 0 Flo adult 2 stems
#> 14 NA NA NA NA NA 0 Flo adult 3 stems
#> 15 NA NA NA NA NA 0 Flo adult 4 stems
#> 16 NA NA NA NA NA 0 Flo adult 5 stems
#> 17 NA NA NA NA NA 0 Flo adult 6 stems
#> 18 NA NA NA NA NA 0 Flo adult 7 stems
Now we will standardize our dataset using the new life history model,
allowing lefko3 to develop new stage calls in all cases.
cypraw_v2 <- verticalize3(data = cypa_data, noyears = 18, firstyear = 1994,
individcol = "plant_id", blocksize = 2, sizeacol = "Inf.94",
sizebcol = "Veg.94", repstracol = "Inf.94", fecacol = "Inf.94",
stageassign = cypframe_fb, stagesize = "sizeadded", NAas0 = TRUE)
summary_hfv(cypraw_v2)
#>
#> This hfv dataset contains 6871 rows, 51 variables, 1 population,
#> 1 patch, 1025 individuals, and 17 time steps.
Now we will develop the vital rate models. These will be developed as
mixed, non-linear models. Survival probability, observation probability
(equivalent to sprouting probability, which is assumed to be less than
1.0 due to vegetative dormancy), and reproduction probability will all
be assessed assuming a binomial response under a logit link. Size and
fecundity will be assessed assuming a Poisson distribution using a log
link, though fecundity will be zero-inflated.
cypmodels2 <- modelsearch(cypraw_v2, historical = FALSE, approach = "mixed",
vitalrates = c("surv", "obs", "size", "repst", "fec"),
sizedist = "poisson", fecdist = "poisson", fec.zero = TRUE,
suite = "main", size = c("size3added", "size2added", "size1added"),
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...
#>
#> Finished selecting best-fit models.
summary(cypmodels2)
#> This LefkoMod object includes 5 linear models.
#> Best-fit model criterion used: aicc&k
#>
#>
#>
#> Survival model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#> Family: binomial ( logit )
#> Formula: alive3 ~ repstatus2 + size2added + (1 | year2) + (1 | individ)
#> Data: surv.data
#> AIC BIC logLik -2*log(L) df.resid
#> 5064.423 5098.598 -2527.211 5054.423 6866
#> Random effects:
#> Groups Name Std.Dev.
#> individ (Intercept) 1.1961
#> year2 (Intercept) 0.4063
#> Number of obs: 6871, groups: individ, 1025; year2, 17
#> Fixed Effects:
#> (Intercept) repstatus2 size2added
#> 2.1917 0.2536 -0.4617
#> optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 2 lme4 warnings
#>
#>
#>
#> Observation model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#> Family: binomial ( logit )
#> Formula: obsstatus3 ~ size2added + (1 | year2) + (1 | individ)
#> Data: obs.data
#> AIC BIC logLik -2*log(L) df.resid
#> 7246.503 7273.305 -3619.252 7238.503 6001
#> Random effects:
#> Groups Name Std.Dev.
#> individ (Intercept) 0.6777
#> year2 (Intercept) 1.0130
#> Number of obs: 6005, groups: individ, 779; year2, 17
#> Fixed Effects:
#> (Intercept) size2added
#> 0.6961 0.1791
#>
#>
#>
#> Size model:
#> Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
#> Family: poisson ( log )
#> Formula: size3added ~ size2added + (1 | year2) + (1 | individ)
#> Data: size.data
#> AIC BIC logLik -2*log(L) df.resid
#> 10296.756 10321.854 -5144.378 10288.756 3919
#> Random effects:
#> Groups Name Std.Dev.
#> individ (Intercept) 0.00000
#> year2 (Intercept) 0.06962
#> Number of obs: 3923, groups: individ, 779; year2, 17
#> Fixed Effects:
#> (Intercept) size2added
#> 0.1682 0.2038
#> optimizer (Nelder_Mead) 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 + size2added + (1 | year2) + (1 | individ)
#> Data: repst.data
#> AIC BIC logLik -2*log(L) df.resid
#> 4635.477 4666.850 -2312.738 4625.477 3918
#> Random effects:
#> Groups Name Std.Dev.
#> individ (Intercept) 0.7218
#> year2 (Intercept) 0.7104
#> Number of obs: 3923, groups: individ, 779; year2, 17
#> Fixed Effects:
#> (Intercept) repstatus2 size2added
#> -1.0444 1.0124 0.2542
#>
#>
#>
#> Fecundity model:
#> Formula: feca2 ~ (1 | year2) + (1 | individ)
#> Zero inflation: ~(1 | year2) + (1 | individ)
#> Data: fec.data
#> AIC BIC logLik -2*log(L) df.resid
#> 5199.186 5232.910 -2593.593 5187.186 2034
#> Random-effects (co)variances:
#>
#> Conditional model:
#> Groups Name Std.Dev.
#> year2 (Intercept) 0.0004197
#> individ (Intercept) 0.0748143
#>
#> Zero-inflation model:
#> Groups Name Std.Dev.
#> year2 (Intercept) 5.0750
#> individ (Intercept) 0.3871
#>
#> Number of obs: 2040 / Conditional model: year2, 17; individ, 699 / Zero-inflation model: year2, 17; individ, 699
#>
#> Fixed Effects:
#>
#> Conditional model:
#> (Intercept)
#> 0.3264
#>
#> Zero-inflation model:
#> (Intercept)
#> -20.05
#>
#>
#> 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: 4
#>
#> Number of models in observation table: 4
#>
#> Number of models in size table: 4
#>
#> Number of models in secondary size table: 1
#>
#> Number of models in tertiary size table: 1
#>
#> Number of models in reproduction status table: 4
#>
#> Number of models in fecundity table: 16
#>
#> 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
#> 32 annual covariate a in time t annucova2
#> 33 annual covariate a in time t-1 annucova1
#> 34 annual covariate b in time t annucovb2
#> 35 annual covariate b in time t-1 annucovb1
#> 36 annual covariate c in time t annucovc2
#> 37 annual covariate c in time t-1 annucovc1
#>
#>
#>
#>
#>
#> Quality control:
#>
#> Survival model estimated with 1025 individuals and 6871 individual transitions.
#> Survival model accuracy is 0.875.
#> Observation status model estimated with 779 individuals and 6005 individual transitions.
#> Observation status model accuracy is 0.726.
#> Primary size model estimated with 779 individuals and 3923 individual transitions.
#> Primary size model R-squared is 0.271.
#> Secondary size model not estimated.
#> Tertiary size model not estimated.
#> Reproductive status model estimated with 779 individuals and 3923 individual transitions.
#> Reproductive status model accuracy is 0.737.
#> Fecundity model estimated with 699 individuals and 2040 individual transitions.
#> Fecundity model R-squared is 0.042.
#> Juvenile survival model not estimated.
#> Juvenile observation status model not estimated.
#> Juvenile primary size model not estimated.
#> Juvenile secondary size model not estimated.
#> Juvenile tertiary size model not estimated.
#> Juvenile reproductive status model not estimated.
#> Juvenile maturity status model not estimated.
The output of the vital rate modeling exercise is worth looking over.
Note that our best-fit models vary in interesting ways. For example,
increasing size has a negative impact on survival, while reproductive
status has a positive impact. These patterns suggest interesting
trade-offs that might constrain natural selection. Also note the model
accuracy in the best-fit models, which includes some models with high
accuracy (e.g., survival, at around 88%), and others with fairly low
accuracy (e.g., primary size, with roughly 27%, and fecundity, with a
pseudo-R2 of 4.2%).
Our next step will be to develop an MPM. Although we will not project
this MPM through our invasibility analysis, this is our chance to see
what the matrices will look like that will be built in the course of
analysis. So, we will set up the proper supplemental parameters, and
then create MPMs to look at more carefully.
germ <- 0.2
sl_mult <- 0.2
sl_surv <- 0.2
seeds_per_fruit <- 10000
cypsupp2f <- supplemental(stage3 = c("SD", "P1", "SL", "SL", "D", "V1", "mat",
"mat", "mat", "SD", "P1"),
stage2 = c("SD", "SD", "P1", "SL", "SL", "SL", "D", "V1", "V2", "rep", "rep"),
eststage3 = c(NA, NA, NA, NA, "D", "V1", "mat", "mat", "mat", NA, NA),
eststage2 = c(NA, NA, NA, NA, "V1", "V1", "D", "V1", "V2", NA, NA),
givenrate = c(0.1, germ, sl_surv, sl_surv, NA, NA, NA, NA, NA, NA, NA),
multiplier = c(NA, NA, NA, NA, sl_mult, sl_mult, 0.4, 0.4, 0.4,
0.5 * seeds_per_fruit * 0.1, 0.5 * seeds_per_fruit * germ),
type =c(1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3),
stageframe = cypframe_fb, historical = FALSE)
#> Warning: NA values in argument multiplier will be treated as 1 values.
cypmatrix2f <- flefko2(stageframe = cypframe_fb, supplement = cypsupp2f,
modelsuite = cypmodels2, data = cypraw_v2)
summary(cypmatrix2f)
#>
#> This ahistorical lefkoMat object contains 17 matrices.
#>
#> Each matrix is square with 18 rows and columns, and a total of 324 elements.
#> A total of 3927 survival transitions were estimated, with 231 per matrix.
#> A total of 238 fecundity transitions were estimated, with 14 per matrix.
#> This lefkoMat object covers 1 population, 1 patch, and 17 time steps.
#>
#> The dataset contains a total of 1025 unique individuals and 6871 unique transitions.
#>
#> Vital rate modeling quality control:
#>
#> Survival estimated with 1025 individuals and 6871 individual transitions.
#> Observation estimated with 779 individuals and 6005 individual transitions.
#> Primary size estimated with 779 individuals and 3923 individual transitions.
#> Secondary size transition not estimated.
#> Tertiary size transition not estimated.
#> Reproductive status estimated with 779 individuals and 3923 individual transitions.
#> Fecundity estimated with 699 individuals and 2040 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] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17]
#> Min. 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.200 0.178 0.160 0.136 0.151 0.200 0.109 0.126
#> 1st Qu. 0.291 0.295 0.294 0.300 0.305 0.281 0.282 0.296 0.280 0.281 0.272 0.248 0.232 0.250 0.280 0.202 0.216
#> Median 0.337 0.430 0.395 0.450 0.466 0.335 0.374 0.408 0.334 0.356 0.312 0.299 0.291 0.296 0.329 0.253 0.264
#> Mean 0.420 0.469 0.446 0.481 0.491 0.416 0.437 0.458 0.415 0.421 0.389 0.365 0.337 0.363 0.410 0.297 0.318
#> 3rd Qu. 0.567 0.659 0.612 0.674 0.686 0.568 0.605 0.634 0.568 0.572 0.526 0.475 0.430 0.471 0.559 0.373 0.422
#> Max. 0.759 0.761 0.746 0.773 0.782 0.740 0.743 0.762 0.738 0.723 0.711 0.703 0.672 0.711 0.730 0.602 0.604
#>
#>
#> All matrices have no stage discontinuities.
Since everything looks alright, we will now need to translate the
lefkoMod object that we created to hold our vital rate
models into an object accessible by adapt3. We will do that
with the miniMod() function in package lefko3.
cyp_vrm <- miniMod(cypmodels2, hfv_data = cypraw_v2)
Next, as before, we need to define a starting vector for the simulation.
We will use the same starting vector as before, using a single protocorm
individual.
cyp_start <- start_input(cypmatrix2f, stage2 = "P1", value = 1)
cyp_start
#> stage2 stage_id_2 stage1 stage_id_1 age2 row_num value
#> 1 P1 2 NA NA NA 2 1
Next, we will set up the proper density dependence to use in
simulations. As before, we advise that users explore the proper
development of density dependence elsewhere.
c2d_1 <- density_input(cypmatrix2f, stage3 = c("P1", "P1"), stage2= c("SD", "rep"),
style = 1, time_delay = 1, alpha = 1.1, beta = 0.0000005, type = c(1, 2))
c2d_1
#> stage3 stage2 stage1 age2 style alpha beta gamma time_delay type type_t12
#> 1 P1 SD NA NA 1 1.1 5e-07 0 1 1 1
#> 2 P1 F1 NA NA 1 1.1 5e-07 0 1 2 1
#> 3 P1 F2 NA NA 1 1.1 5e-07 0 1 2 1
#> 4 P1 F3 NA NA 1 1.1 5e-07 0 1 2 1
#> 5 P1 F4 NA NA 1 1.1 5e-07 0 1 2 1
#> 6 P1 F5 NA NA 1 1.1 5e-07 0 1 2 1
#> 7 P1 F6 NA NA 1 1.1 5e-07 0 1 2 1
#> 8 P1 F7 NA NA 1 1.1 5e-07 0 1 2 1
Next, we will quickly run a projection of the sample MPM under the
density dependence that we set above, using lefko3’s
f_projection3() function. This will give us a sense of
whether everything appears to be working.
cyp_proj1 <- f_projection3(format = 3, stageframe = cypframe_fb,
supplement = cypsupp2f, modelsuite = cypmodels2, times = 1000,
start_frame = cyp_start, data = cypraw_v2, density = c2d_1,
integeronly = FALSE, year = 2008)
#> Warning: Option patch not set, so will set to first patch/population.
plot(cyp_proj1)
Figure I3. Density dependent projected
population growth under density frame c2d_1
Next, we will set a trait axis defining all variants to test. Here, we
will use 10 variants total, and please note that the resulting analysis
may take a while to run. Note that these will be defined solely by
additive deviation to the y-intercept of the sprouting model
(observation probability), which is a binomial mixed effects model.
Normally, we would advise using many more variants, such as 50 or 100,
in order to get a good resolution of trait variation for the pairwise
invasibility plot, in case the actual adaptive dynamics are complex, and
to estimate the trait optima properly.
current_traits <- trait_axis(cypframe_fb, stagebased = TRUE, agebased = FALSE,
historical = FALSE, obs_dev = seq(from = -0.5, to = 2.5, length.out = 10))
current_traits
#> variant stage3 stage2 stage1 age3 age2 eststage3 eststage2 eststage1 estage3 estage2 givenrate offset multiplier convtype
#> 1 1 <NA> <NA> <NA> NA NA <NA> <NA> <NA> NA NA NA NA NA NA
#> 2 2 <NA> <NA> <NA> NA NA <NA> <NA> <NA> NA NA NA NA NA NA
#> 3 3 <NA> <NA> <NA> NA NA <NA> <NA> <NA> NA NA NA NA NA NA
#> 4 4 <NA> <NA> <NA> NA NA <NA> <NA> <NA> NA NA NA NA NA NA
#> 5 5 <NA> <NA> <NA> NA NA <NA> <NA> <NA> NA NA NA NA NA NA
#> 6 6 <NA> <NA> <NA> NA NA <NA> <NA> <NA> NA NA NA NA NA NA
#> 7 7 <NA> <NA> <NA> NA NA <NA> <NA> <NA> NA NA NA NA NA NA
#> 8 8 <NA> <NA> <NA> NA NA <NA> <NA> <NA> NA NA NA NA NA NA
#> 9 9 <NA> <NA> <NA> NA NA <NA> <NA> <NA> NA NA NA NA NA NA
#> 10 10 <NA> <NA> <NA> NA NA <NA> <NA> <NA> NA NA NA NA NA NA
#> convtype_t12 surv_dev obs_dev size_dev sizeb_dev sizec_dev repst_dev fec_dev jsurv_dev jobs_dev jsize_dev jsizeb_dev jsizec_dev
#> 1 NA NA -0.5000000 NA NA NA NA NA NA NA NA NA NA
#> 2 NA NA -0.1666667 NA NA NA NA NA NA NA NA NA NA
#> 3 NA NA 0.1666667 NA NA NA NA NA NA NA NA NA NA
#> 4 NA NA 0.5000000 NA NA NA NA NA NA NA NA NA NA
#> 5 NA NA 0.8333333 NA NA NA NA NA NA NA NA NA NA
#> 6 NA NA 1.1666667 NA NA NA NA NA NA NA NA NA NA
#> 7 NA NA 1.5000000 NA NA NA NA NA NA NA NA NA NA
#> 8 NA NA 1.8333333 NA NA NA NA NA NA NA NA NA NA
#> 9 NA NA 2.1666667 NA NA NA NA NA NA NA NA NA NA
#> 10 NA NA 2.5000000 NA NA NA NA NA NA NA NA NA NA
#> jrepst_dev jmatst_dev indcova indcovb indcovc
#> 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
Finally, we get to our invasibility analysis. As before, this is set up
as a pairwise invasibility analysis (hence the
var_per_run = 2 setting). As this will take a while, grab a
cappuccino in the meantime. To speed things up, change the number of
variants in the trait axis from 10 to something lower, like 3 or 4
(though this will result in much poorer resolution to see any potential
ESS values). Note that we are also setting
trait_optima = TRUE, which allows us to estimate the ESS
value of the additive deviation to the sprouting model, if one exists.
cyp_inv3 <- invade3(axis = current_traits, vrm = cyp_vrm, format = 3,
starts = cyp_start, years = 2008, stageframe = cypframe_fb,
supplement = cypsupp2f, entry_time = c(0, 800), integeronly = FALSE,
times = 1000, fitness_times = 100, var_per_run = 2, stochastic = FALSE,
density = c2d_1, trait_optima = TRUE)
Let’s see a summary of the resulting object.
summary(cyp_inv3)
#>
#> The input adaptInv object covers 10 variants, 1001 projected steps, 2 variants per run, and 1 projected replicate.
#> It includes optimization data suggesting 1 purported ESS optimum.
The resulting output suggests that we have an ESS value of our trait. As
before, we can take a peek at what some of these simulations look like.
We will try two different runs. Note that if the number of trait
variants is set to less than 20, then the run numbers might need to be
adjusted.
par(mfrow = c(2, 1))
plot(cyp_inv3, pip = FALSE, run = 1)
plot(cyp_inv3, pip = FALSE, run = 11)
Figure I4. Population growth of resident (solid
black) and invader (dashed red) variants in two runs
At first glance the invader seems to be so low as to be extinct in both
cases. However, we can take a look at the actual growth of the invader
in run 11 as follows.
cyp_inv3$N_out[[1]][2,,11]
#> [1] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [11] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [21] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [31] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [41] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [51] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [61] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [71] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [81] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [91] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [101] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [111] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [121] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [131] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [141] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [151] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [161] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [171] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [181] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [191] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [201] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [211] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [221] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [231] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [241] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [251] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [261] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [271] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [281] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [291] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [301] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [311] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [321] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [331] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [341] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [351] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [361] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [371] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [381] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [391] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [401] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [411] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [421] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [431] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [441] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [451] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [461] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [471] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [481] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [491] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [501] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [511] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [521] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [531] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [541] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [551] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [561] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [571] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [581] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [591] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [601] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [611] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [621] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [631] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [641] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [651] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [661] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [671] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [681] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [691] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [701] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [711] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [721] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [731] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [741] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [751] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [761] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [771] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [781] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [791] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
#> [801] 1.00000000 0.20000000 0.05689588 0.01625894 0.85675865 0.84117873 0.53852514 0.28750896 0.45963811 0.66678898
#> [811] 0.68243668 0.55158649 0.50181915 0.57147537 0.64877794 0.64806505 0.60995111 0.60484331 0.63969541 0.67070287
#> [821] 0.67441691 0.66852374 0.67637508 0.69677699 0.71407013 0.72200993 0.72802137 0.73947314 0.75483310 0.76862311
#> [831] 0.77940077 0.79008904 0.80293534 0.81704440 0.83055851 0.84323485 0.85617777 0.87004217 0.88440595 0.89865881
#> [841] 0.91280727 0.92725562 0.94218445 0.95741579 0.97275724 0.98824334 1.00401935 1.02013240 1.03651383 1.05310726
#> [851] 1.06993717 1.08705624 1.10447721 1.12217671 1.14014135 1.15838575 1.17693063 1.19578131 1.21493234 1.23438357
#> [861] 1.25414436 1.27422487 1.29462964 1.31536017 1.33642043 1.35781755 1.37955853 1.40164843 1.42409155 1.44689330
#> [871] 1.47006012 1.49359834 1.51751367 1.54181173 1.56649860 1.59158078 1.61706473 1.64295677 1.66926331 1.69599099
#> [881] 1.72314665 1.75073717 1.77876947 1.80725059 1.83618771 1.86558817 1.89545941 1.92580894 1.95664440 1.98797358
#> [891] 2.01980440 2.05214489 2.08500321 2.11838764 2.15230661 2.18676868 2.22178255 2.25735704 2.29350115 2.33022398
#> [901] 2.36753480 2.40544304 2.44395825 2.48309015 2.52284862 2.56324369 2.60428555 2.64598456 2.68835124 2.73139628
#> [911] 2.77513054 2.81956507 2.86471106 2.91057992 2.95718321 3.00453270 3.05264033 3.10151825 3.15117878 3.20163447
#> [921] 3.25289803 3.30498240 3.35790073 3.41166638 3.46629290 3.52179408 3.57818393 3.63547667 3.69368676 3.75282889
#> [931] 3.81291799 3.87396921 3.93599797 3.99901990 4.06305092 4.12810719 4.19420511 4.26136137 4.32959291 4.39891695
#> [941] 4.46935098 4.54091277 4.61362039 4.68749218 4.76254677 4.83880311 4.91628045 4.99499831 5.07497659 5.15623544
#> [951] 5.23879538 5.32267724 5.40790219 5.49449172 5.58246770 5.67185231 5.76266812 5.85493803 5.94868534 6.04393369
#> [961] 6.14070713 6.23903006 6.33892730 6.44042406 6.54354595 6.64831898 6.75476960 6.86292466 6.97281146 7.08445773
#> [971] 7.19789163 7.31314179 7.43023729 7.54920767 7.67008296 7.79289366 7.91767075 8.04444573 8.17325057 8.30411778
#> [981] 8.43708038 8.57217193 8.70942651 8.84887875 8.99056384 9.13451753 9.28077614 9.42937659 9.58035636 9.73375356
#> [991] 9.88960689 10.04795567 10.20883986 10.37230006 10.53837752 10.70711413 10.87855247 11.05273581 11.22970810 11.40951398
#> [1001] 11.59219884
We can see that the invader does grow in at least some of the
simulations. Now let’s see the pairwise invasibility plot (PIP).
plot(cyp_inv3, auto_axis = TRUE)
Figure I5. Pairwise invasibility plot showing an
intermediate level of sprouting as the predicted ESS. Black and white
regions show positive and negative invader fitness, respectively
A fairly beautiful PIP, if I say so myself! We find that there appears
to be an ESS value for sprouting. Let’s now take a look at the
associated elasticity plot.
plot(cyp_inv3, pip = FALSE, elast = TRUE, auto_axis = TRUE)
Figure I6. Elasticity plot showing the ESS
deviation to the y-intercept of sprouting as a value somewhere around
the 7th or 8th variant in the trait frame
Elasticity analyses are performed using a variant of the procedure
outlined in Benton and Grant (1999), and
on the version of it that appears in Roff
(2010). In this case, each original variant in the supplied trait
axis is run as the resident, against an invader that has the key
variable trait reduced by 0.5%. Here, we see a potential optimum near 0,
corresponding roughly to variant 7 in our trait axis. Let’s take a look
at the predicted ESS value of the trait.
cyp_inv3$optim$ESS_values
#> variant stage3 stage2 stage1 age3 age2 eststage3 eststage2 eststage1 estage3 estage2 givenrate offset multiplier convtype
#> 1 1 NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> convtype_t12 surv_dev obs_dev size_dev sizeb_dev sizec_dev repst_dev fec_dev jsurv_dev jobs_dev jsize_dev jsizeb_dev jsizec_dev
#> 1 NA NA 0.1666667 NA NA NA NA NA NA NA NA NA NA
#> jrepst_dev jmatst_dev indcova indcovb indcovc year2 mpm_altered vrm_altered fitness converged
#> 1 NA NA NA NA NA <NA> 0 1 -8.917901e-06 TRUE
Our converged ESS value is at a slightly positive additive deviation to
our sprouting probability y-intercept, meaning that our observed value
appears to be an ESS or very close to it. This is similar to what was
found in Shefferson, Warren II, and Pulliam
(2014), so we’ve done reasonably well.
Let’s now estimate the standard error for the ESS value. We will use the
demographic shuffling approach used in Shefferson, Warren II, and Pulliam (2014).
library(lme4)
library(glmmTMB)
sample_traits <- trait_axis(cypframe_fb, stagebased = TRUE, agebased = FALSE,
historical = FALSE, obs_dev = seq(from = -1.5, to = 3.0, length.out = 10))
num_runs <- 10
output_vec <- rep(NA, num_runs)
sample_size <- dim(cypa_data)[1]
for (i in c(1:num_runs)) {
set.seed(i)
sample_order <- sample(1:sample_size, sample_size, replace = TRUE)
sample_data <- cypa_data[sample_order,]
sample_data$plant_id <- c(1:sample_size)
sampleraw_v2 <- verticalize3(data = sample_data, noyears = 18, firstyear = 1994,
individcol = "plant_id", blocksize = 2, sizeacol = "Inf.94",
sizebcol = "Veg.94", repstracol = "Inf.94", fecacol = "Inf.94",
stageassign = cypframe_fb, stagesize = "sizeadded", NAas0 = TRUE)
sample_surv <- try(glmer(alive3 ~ repstatus2 + size2added + (1 | year2) +
(1 | individ), data = sampleraw_v2, family = binomial))
sample_obsdata <- subset(sampleraw_v2, alive3 == 1)
sample_obs <- try(glmer(obsstatus3 ~ size2added + (1 | year2) + (1 | individ),
data = sample_obsdata, family = binomial))
sample_sizedata <- subset(sample_obsdata, obsstatus3 == 1)
sample_sizem <- try(glmer(size3added ~ size2added + (1 | year2) + (1 | individ),
data = sample_sizedata, family = poisson))
sample_repst <- try(glmer(repstatus3 ~ repstatus2 + size2added + (1 | year2) +
(1 | individ), data = sample_sizedata, family = poisson))
sample_fecdata <- subset(sampleraw_v2, repstatus2 == 1)
sample_fec <- try(glmmTMB(feca2 ~ (1 | year2) + (1 | individ),
zi = feca2 ~ (1 | year2) + (1 | individ), data = sample_fecdata,
family = poisson))
samplemodels2 <- cypmodels2
samplemodels2$survival_model <- sample_surv
samplemodels2$observation_model <- sample_obs
samplemodels2$size_model <- sample_sizem
samplemodels2$repstatus_model <- sample_repst
samplemodels2$fecundity_model <- sample_fec
sample_vrm <- try(miniMod(samplemodels2, hfv_data = sampleraw_v2))
sample_inv3 <- try(invade3(axis = sample_traits, vrm = sample_vrm, format = 3,
starts = cyp_start, years = 2008, stageframe = cypframe_fb,
supplement = cypsupp2f, entry_time = c(0, 800), integeronly = FALSE,
times = 1000, fitness_times = 190, var_per_run = 2, stochastic = FALSE,
density = c2d_1, trait_optima = TRUE))
try(output_vec[i] <- sample_inv3$optim$ESS_values$obs_dev[1])
try(writeLines(paste0("ESS: ", sample_inv3$optim$ESS_values$obs_dev[1], "\n")))
}
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.224395 (tol = 0.002, component 1)
#> See ?lme4::convergence and ?lme4::troubleshooting.
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
#> - Rescale variables?
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> Warning in f(init, x[[i]]): discarding LHS of second argument
#> Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence problem; non-positive-definite Hessian matrix. See
#> vignette('troubleshooting')
#> ESS: -0.169397333333333
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.227223 (tol = 0.002, component 1)
#> See ?lme4::convergence and ?lme4::troubleshooting.
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
#> - Rescale variables?
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> Warning in f(init, x[[i]]): discarding LHS of second argument
#> Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence problem; non-positive-definite Hessian matrix. See
#> vignette('troubleshooting')
#> Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence problem; singular convergence (7). See
#> vignette('troubleshooting'), help('diagnose')
#> Error in eval(expr, envir) :
#> Data frames must have the same number of columns.
#> Error in sample_inv3$optim : $ operator is invalid for atomic vectors
#> Error in sample_inv3$optim : $ operator is invalid for atomic vectors
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> Warning in f(init, x[[i]]): discarding LHS of second argument
#> ESS: -0.230833333333333
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.208192 (tol = 0.002, component 1)
#> See ?lme4::convergence and ?lme4::troubleshooting.
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
#> - Rescale variables?
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> Warning in f(init, x[[i]]): discarding LHS of second argument
#> Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence problem; non-positive-definite Hessian matrix. See
#> vignette('troubleshooting')
#> ESS: -0.161197916666667
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : unable to evaluate scaled gradient
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge: degenerate Hessian with 1 negative eigenvalues
#> See ?lme4::convergence and ?lme4::troubleshooting.
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> Warning in f(init, x[[i]]): discarding LHS of second argument
#> ESS: -0.145182291666667
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> Warning in f(init, x[[i]]): discarding LHS of second argument
#> Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence problem; non-positive-definite Hessian matrix. See
#> vignette('troubleshooting')
#> ESS: -0.19
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.216018 (tol = 0.002, component 1)
#> See ?lme4::convergence and ?lme4::troubleshooting.
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
#> - Rescale variables?
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> Warning in f(init, x[[i]]): discarding LHS of second argument
#> Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence problem; non-positive-definite Hessian matrix. See
#> vignette('troubleshooting')
#> Error in eval(expr, envir) :
#> Data frames must have the same number of columns.
#> Error in sample_inv3$optim : $ operator is invalid for atomic vectors
#> Error in sample_inv3$optim : $ operator is invalid for atomic vectors
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.202177 (tol = 0.002, component 1)
#> See ?lme4::convergence and ?lme4::troubleshooting.
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
#> - Rescale variables?
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> Warning in f(init, x[[i]]): discarding LHS of second argument
#> Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence problem; non-positive-definite Hessian matrix. See
#> vignette('troubleshooting')
#> ESS: -0.148177083333333
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.228734 (tol = 0.002, component 1)
#> See ?lme4::convergence and ?lme4::troubleshooting.
#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
#> - Rescale variables?
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> Warning in f(init, x[[i]]): discarding LHS of second argument
#> Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): failed to invert Hessian from numDeriv::jacobian(), falling back to
#> internal vcov estimate
#> ESS: -0.168064768
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> Warning in f(init, x[[i]]): discarding LHS of second argument
#> Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence problem; non-positive-definite Hessian matrix. See
#> vignette('troubleshooting')
#> ESS: -0.170762666666667
Let’s see the mean and SD, as well as whether the ESS value is
significantly different from 0.
mean_sim <- mean(output_vec, na.rm = TRUE)
se_sim <- sd(output_vec, na.rm = TRUE)
true_ESS <- cyp_inv3$optim$ESS_values$obs_dev[1]
z_value <- abs(true_ESS) / se_sim
p_value <- pnorm(z_value, mean = 0, sd = 1, lower.tail = FALSE)
writeLines(paste0("True ESS: ", true_ESS))
#> True ESS: 0.166666666666667
writeLines(paste0("Simulated ESS: ", mean_sim))
#> Simulated ESS: -0.172951924125
writeLines(paste0("Simulated SE: ", se_sim))
#> Simulated SE: 0.0272693145084921
writeLines(paste0("Z score: ", z_value))
#> Z score: 6.11187591880112
writeLines(paste0("p-value: ", p_value))
#> p-value: 4.9233366559445e-10
Now let’s work out a figure showing the predicted sprouting probability.
cyp_pred_df <- data.frame(size2added = c(0:7), year2 = 2004, individ = "A 1")
new_model <- cypmodels2$observation_model
new_model@beta[1] <- cypmodels2$observation_model@beta[1] + true_ESS
cyp_pred_df$pred_y <- predict(cypmodels2$observation_model,
newdata = cyp_pred_df, type = "response")
cyp_pred_df$new_pred_y <- predict(new_model, newdata = cyp_pred_df,
type = "response")
plot(pred_y ~ size2added, data = cyp_pred_df, ylab = "P(sprouting)",
xlab = "# of sprouts", ylim = c(0.60, 0.90), lwd = 2, type = "l", bty = "n")
lines(new_pred_y ~ size2added, data = cyp_pred_df, lty = 2, lwd = 2,
col = "red")
legend("bottomright", c("Observed", "ESS"), lty = c(1, 2), lwd = 2,
col = c("black", "red"), bty = "n")
plot of chunk ch2-2.18