flefko3 {lefko3}

Create Function-based Historical Matrix Projection Model

Description

Function flefko3() returns function-based historical MPMs corresponding to the patches and occasions given, including the associated component transition and fecundity matrices, data frames detailing the characteristics of the ahistorical stages used and historical stage pairs created, and a data frame characterizing the patch and occasion combinations corresponding to these matrices.

Usage

flefko3(
  year = "all",
  patch = "all",
  stageframe,
  supplement = NULL,
  repmatrix = NULL,
  overwrite = NULL,
  data = NULL,
  modelsuite = NULL,
  surv_model = NULL,
  obs_model = NULL,
  size_model = NULL,
  sizeb_model = NULL,
  sizec_model = NULL,
  repst_model = NULL,
  fec_model = NULL,
  jsurv_model = NULL,
  jobs_model = NULL,
  jsize_model = NULL,
  jsizeb_model = NULL,
  jsizec_model = NULL,
  jrepst_model = NULL,
  jmatst_model = NULL,
  paramnames = NULL,
  inda = NULL,
  indb = NULL,
  indc = NULL,
  surv_dev = 0,
  obs_dev = 0,
  size_dev = 0,
  sizeb_dev = 0,
  sizec_dev = 0,
  repst_dev = 0,
  fec_dev = 0,
  jsurv_dev = 0,
  jobs_dev = 0,
  jsize_dev = 0,
  jsizeb_dev = 0,
  jsizec_dev = 0,
  jrepst_dev = 0,
  jmatst_dev = 0,
  density = NA,
  fecmod = 1,
  random.inda = FALSE,
  random.indb = FALSE,
  random.indc = FALSE,
  negfec = FALSE,
  format = "ehrlen",
  ipm_method = "CDF",
  reduce = FALSE,
  simple = FALSE,
  err_check = FALSE,
  exp_tol = 700,
  theta_tol = 1e+08,
  sparse_output = FALSE
)

Arguments

year

A variable corresponding to the observation occasion, or a set of such values, given in values associated with the year term used in linear model development. Defaults to "all", in which case matrices will be estimated for all occasions.

patch

A variable designating which patches or subpopulations will have matrices estimated. Defaults to "all", but can also be set to specific patch names or a vector thereof.

stageframe

An object of class stageframe. These objects are generated by function sf_create(), and include information on the size, observation status, propagule status, reproduction status, immaturity status, maturity status, stage group, size bin widths, and other key characteristics of each ahistorical stage.

supplement

An optional data frame of class lefkoSD that provides supplemental data that should be incorporated into the MPM. Three kinds of data may be integrated this way: transitions to be estimated via the use of proxy transitions, transition overwrites from the literature or supplemental studies, and transition multipliers for survival and fecundity. This data frame should be produced using the supplemental() function. Can be used in place of or in addition to an overwrite table (see overwrite below) and a reproduction matrix (see repmatrix below).

repmatrix

An optional reproduction matrix. This matrix is composed mostly of 0s, with non-zero entries acting as element identifiers and multipliers for fecundity (with 1 equaling full fecundity). If left blank, and no supplement is provided, then flefko3() will assume that all stages marked as reproductive produce offspring at 1x that of estimated fecundity, and that offspring production will yield the first stage noted as propagule or immature. May be the dimensions of either a historical or an ahistorical matrix. If the latter, then all stages will be used in occasion t-1 for each suggested ahistorical transition.

overwrite

An optional data frame developed with the overwrite() function describing transitions to be overwritten either with given values or with other estimated transitions. Note that this function supplements overwrite data provided in supplement.

data

The historical vertical demographic data frame used to estimate vital rates (class hfvdata), which is required to initialize times and patches properly. Variable names should correspond to the naming conventions in verticalize3() and historicalize3(). Not required if option modelsuite is set to a vrm_input object.

modelsuite

One of two kinds of lists. The first is a lefkoMod object holding the vital rate models and associated metadata. Alternatively, an object of class vrm_input may be provided. If given, then surv_model, obs_model, size_model, sizeb_model, sizec_model, repst_model, fec_model, jsurv_model, jobs_model, jsize_model, jsizeb_model, jsizec_model, jrepst_model, jmatst_model, and paramnames are not required. One or more of these models should include size or reproductive status in occasion t-1. Although this is optional input, it is recommended, and without it all vital rate model inputs (named XX_model) are required.

surv_model

A linear model predicting survival probability. This can be a model of class glm or glmer, and requires a predicted binomial variable under a logit link. Ignored if modelsuite is provided. This model must have been developed in a modeling exercise testing the impacts of occasions t and t-1.

obs_model

A linear model predicting sprouting or observation probability. This can be a model of class glm or glmer, and requires a predicted binomial variable under a logit link. Ignored if modelsuite is provided. This model must have been developed in a modeling exercise testing the impacts of occasions t and t-1.

size_model

A linear model predicting primary size. This can be a model of class glm, glmer, glmmTMB, zeroinfl, vglm, lm, or lmer. Ignored if modelsuite is provided. This model must have been developed in a modeling exercise testing the impacts of occasions t and t-1.

sizeb_model

A linear model predicting secondary size. This can be a model of class glm, glmer, glmmTMB, zeroinfl, vglm, lm, or lmer. Ignored if modelsuite is provided. This model must have been developed in a modeling exercise testing the impacts of occasions t and t-1.

sizec_model

A linear model predicting tertiary size. This can be a model of class glm, glmer, glmmTMB, zeroinfl, vglm, lm, or lmer. Ignored if modelsuite is provided. This model must have been developed in a modeling exercise testing the impacts of occasions t and t-1.

repst_model

A linear model predicting reproduction probability. This can be a model of class glm or glmer, and requires a predicted binomial variable under a logit link. Ignored if modelsuite is provided. This model must have been developed in a modeling exercise testing the impacts of occasions t and t-1.

fec_model

A linear model predicting fecundity. This can be a model of class glm, glmer, glmmTMB, zeroinfl, vglm, lm, or lmer. Ignored if modelsuite is provided. This model must have been developed in a modeling exercise testing the impacts of occasions t and t-1.

jsurv_model

A linear model predicting juvenile survival probability. This can be a model of class glm or glmer, and requires a predicted binomial variable under a logit link. Ignored if modelsuite is provided. This model must have been developed in a modeling exercise testing the impacts of occasions t and t-1.

jobs_model

A linear model predicting juvenile sprouting or observation probability. This can be a model of class glm or glmer, and requires a predicted binomial variable under a logit link. Ignored if modelsuite is provided. This model must have been developed in a modeling exercise testing the impacts of occasions t and t-1.

jsize_model

A linear model predicting juvenile primary size. This can be a model of class glm, glmer, glmmTMB, zeroinfl, vglm, lm, or lmer. Ignored if modelsuite is provided. This model must have been developed in a modeling exercise testing the impacts of occasions t and t-1.

jsizeb_model

A linear model predicting juvenile secondary size. This can be a model of class glm, glmer, glmmTMB, zeroinfl, vglm, lm, or lmer. Ignored if modelsuite is provided. This model must have been developed in a modeling exercise testing the impacts of occasions t and t-1.

jsizec_model

A linear model predicting juvenile tertiary size. This can be a model of class glm, glmer, glmmTMB, zeroinfl, vglm, lm, or lmer. Ignored if modelsuite is provided. This model must have been developed in a modeling exercise testing the impacts of occasions t and t-1.

jrepst_model

A linear model predicting reproduction probability of a mature individual that was immature in time t. This can be a model of class glm or glmer, and requires a predicted binomial variable under a logit link. Ignored if modelsuite is provided. This model must have been developed in a modeling exercise testing the impacts of occasions t and t-1.

jmatst_model

A linear model predicting maturity probability of an individual that was immature in time t. This can be a model of class glm or glmer, and requires a predicted binomial variable under a logit link. Ignored if modelsuite is provided. This model must have been developed in a modeling exercise testing the impacts of occasions t and t-1.

paramnames

A data frame with three columns, the first describing all terms used in linear modeling, the second (must be called mainparams) giving the general model terms that will be used in matrix creation, and the third showing the equivalent terms used in modeling (must be named modelparams). Function create_pm() can be used to create a skeleton paramnames object, which can then be edited. Only required if modelsuite is not supplied.

inda

Can be a single value to use for individual covariate a in all matrices, a pair of values to use for times t and t-1 in historical matrices, or a vector of such values corresponding to each occasion in the dataset. Defaults to NULL.

indb

Can be a single value to use for individual covariate b in all matrices, a pair of values to use for times t and t-1 in historical matrices, or a vector of such values corresponding to each occasion in the dataset. Defaults to NULL.

indc

Can be a single value to use for individual covariate c in all matrices, a pair of values to use for times t and t-1 in historical matrices, or a vector of such values corresponding to each occasion in the dataset. Defaults to NULL.

surv_dev

A numeric value to be added to the y-intercept in the linear model for survival probability. Defaults to 0.

obs_dev

A numeric value to be added to the y-intercept in the linear model for observation probability. Defaults to 0.

size_dev

A numeric value to be added to the y-intercept in the linear model for primary size. Defaults to 0.

sizeb_dev

A numeric value to be added to the y-intercept in the linear model for secondary size. Defaults to 0.

sizec_dev

A numeric value to be added to the y-intercept in the linear model for tertiary size. Defaults to 0.

repst_dev

A numeric value to be added to the y-intercept in the linear model for probability of reproduction. Defaults to 0.

fec_dev

A numeric value to be added to the y-intercept in the linear model for fecundity. Defaults to 0.

jsurv_dev

A numeric value to be added to the y-intercept in the linear model for juvenile survival probability. Defaults to 0.

jobs_dev

A numeric value to be added to the y-intercept in the linear model for juvenile observation probability. Defaults to 0.

jsize_dev

A numeric value to be added to the y-intercept in the linear model for juvenile primary size. Defaults to 0.

jsizeb_dev

A numeric value to be added to the y-intercept in the linear model for juvenile secondary size. Defaults to 0.

jsizec_dev

A numeric value to be added to the y-intercept in the linear model for juvenile tertiary size. Defaults to 0.

jrepst_dev

A numeric value to be added to the y-intercept in the linear model for juvenile reproduction probability. Defaults to 0.

jmatst_dev

A numeric value to be added to the y-intercept in the linear model for juvenile maturity probability. Defaults to 0.

density

A numeric value indicating density value to use to propagate matrices. Only needed if density is an explanatory term used in one or more vital rate models. Defaults to NA.

fecmod

A scalar multiplier of fecundity. Defaults to 1.0.

random.inda

A logical value denoting whether to treat individual covariate a as a random, categorical variable. Otherwise is treated as a fixed, numeric variable. Defaults to FALSE.

random.indb

A logical value denoting whether to treat individual covariate b as a random, categorical variable. Otherwise is treated as a fixed, numeric variable. Defaults to FALSE.

random.indc

A logical value denoting whether to treat individual covariate c as a random, categorical variable. Otherwise is treated as a fixed, numeric variable. Defaults to FALSE.

negfec

A logical value denoting whether fecundity values estimated to be negative should be reset to 0. Defaults to FALSE.

format

A string indicating whether to estimate matrices in ehrlen format or deVries format. The latter adds one extra prior stage to account for the prior state of newborns. Defaults to ehrlen format.

ipm_method

A string indicating what method to use to estimate size transition probabilities, if size is treated as continuous. Options include: "midpoint", which utilizes the midpoint method; and "CDF", which uses the cumulative distribution function. Defaults to "CDF".

reduce

A logical value denoting whether to remove historical stages associated solely with 0 transitions. These are only removed in cases where the associated row and column sums in ALL matrices estimated equal 0. Defaults to FALSE.

simple

A logical value indicating whether to produce A, U, and F matrices, or only the latter two. Defaults to FALSE, in which case all three are output.

err_check

A logical value indicating whether to append extra information used in matrix calculation within the output list. Defaults to FALSE.

exp_tol

A numeric value used to indicate a maximum value to set exponents to in the core kernel to prevent numerical overflow. Defaults to 700.

theta_tol

A numeric value used to indicate a maximum value to theta as used in the negative binomial probability density kernel. Defaults to 100000000, but can be reset to other values during error checking.

sparse_output

A logical value indicating whether to output matrices in sparse format. Defaults to FALSE, in which case all matrices are output in standard matrix format.

Value

If all inputs are properly formatted, then this function will return an object of class lefkoMat, which is a list that holds the matrix projection model and all of its metadata. The structure has the following elements:

A

A list of full projection matrices in order of sorted patches and occasion times. All matrices output in R's matrix class, or in the dgCMatrix class from the Matrix package if sparse.

U

A list of survival transition matrices sorted as in A. All matrices output in R's matrix class, or in the dgCMatrix class from the Matrix package if sparse.

F

A list of fecundity matrices sorted as in A. All matrices output in R's matrix class, or in the dgCMatrix class from the Matrix package if sparse.

hstages

A data frame matrix showing the pairing of ahistorical stages used to create historical stage pairs.

agestages

A data frame showing age-stage pairs. In this function, it is set to NA. Only used in output to function aflefko2().

ahstages

A data frame detailing the characteristics of associated ahistorical stages, in the form of a modified stageframe that includes status as an entry stage through reproduction.

labels

A data frame giving the population, patch, and year of each matrix in order. In flefko3(), only one population may be analyzed at once.

dataqc

A vector showing the numbers of individuals and rows in the vertical dataset used as input.

matrixqc

A short vector describing the number of non-zero elements in U and F matrices, and the number of annual matrices.

modelqc

This is the qc portion of the modelsuite input.

prob_out

An optional element only added if err_check = TRUE. This is a list of vital rate probability matrices, with 7 columns in the order of survival, observation probability, reproduction probability, primary size transition probability, secondary size transition probability, tertiary size transition probability, and probability of juvenile transition to maturity.

allstages

An optional element only added if err_check = TRUE. This is a data frame giving the values used to determine each matrix element capable of being estimated.

Notes

Unlike rlefko2(), rlefko3(), arlefko2(), and rleslie(), this function does not currently distinguish populations. Users wishing to use the same vital rate models across populations should label them as patches (though we do not advise this approach, as populations should typically be treated as statistically independent).

The default behavior of this function is to estimate fecundity with regards to transitions specified via associated fecundity multipliers in either supplement or repmatrix. If both of these fields are left empty, then fecundity will be estimated at full for all transitions leading from reproductive stages to immature and propagule stages. However, if a supplement is provided and a repmatrix is not, or if repmatrix is set to 0, then only fecundity transitions noted in the supplement will be set to non-zero values. To use the default behavior of setting all reproductive stages to reproduce at full fecundity into immature and propagule stages, but also incorporate given or proxy survival transitions, input those given and proxy transitions through the overwrite option.

If used, the reproduction matrix (field repmatrix) may be supplied as either historical or ahistorical. If provided as ahistorical, then flefko3() will assume that all historical transitions involving stages noted for occasions t and t+1 should be set to the respective fecundity multipliers noted.

Users may at times wish to estimate MPMs using a dataset incorporating multiple patches or subpopulations, but without discriminating between those patches or subpopulations. Should the aim of analysis be a general MPM that does not distinguish these patches or subpopulations, the modelsearch() run should not include patch terms.

Input options including multiple variable names must be entered in the order of variables in occasion t+1, t, and t-1. Rearranging the order will lead to erroneous calculations, and will may lead to fatal errors.

The ipm_method function gives the option of using two different means of estimating the probability of size transition. The midpoint method ("midpoint") refers to the method in which the probability is estimated by first estimating the probability associated with transition from the exact size at the midpoint of the size class using the corresponding probability density function, and then multiplying that value by the bin width of the size class. Doak et al. 2021 (Ecological Monographs) noted that this method can produce biased results, with total size transitions associated with a specific size not totaling to 1.0 and even specific size transition probabilities capable of being estimated at values greater than 1.0. The alternative and default method, "CDF", uses the corresponding cumulative density function to estimate the probability of size transition as the cumulative probability of size transition at the greater limit of the size class minus the cumulative probability of size transition at the lower limit of the size class. This latter method avoids this bias. Note, however, that both methods are exact and unbiased for negative binomial and Poisson distributions.

Under the Gaussian and gamma size distributions, the number of estimated parameters may differ between the two ipm_method settings. Because the midpoint method has a tendency to incorporate upward bias in the estimation of size transition probabilities, it is more likely to yield non- zero values when the true probability is extremely close to 0. This will result in the summary.lefkoMat function yielding higher numbers of estimated parameters than the ipm_method = "CDF" yields in some cases.

Using the err_check option will produce a matrix of 7 columns, each characterizing a different vital rate. The product of each row yields an element in the associated U matrix. The number and order of elements in each column of this matrix matches the associated matrix in column vector format. Use of this option is generally for the purposes of debugging code. ' Individual covariates are treated as categorical only if they are set as random terms. Fixed categorical individual covariates are currently not allowed. However, such terms may be supplied if the modelsuite option is set to a vrm_input object. In that case, the user should also set the logical random switch for the individual covariate to be used to TRUE (e.g., random.inda = TRUE).

See Also

mpm_create()

flefko2()

aflefko2()

arlefko2()

fleslie()

rlefko3()

rlefko2()

rleslie()

Examples

# Lathyrus example
data(lathyrus)

sizevector <- c(0, 4.6, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8,
  9)
stagevector <- c("Sd", "Sdl", "Dorm", "Sz1nr", "Sz2nr", "Sz3nr", "Sz4nr",
  "Sz5nr", "Sz6nr", "Sz7nr", "Sz8nr", "Sz9nr", "Sz1r", "Sz2r", "Sz3r", 
  "Sz4r", "Sz5r", "Sz6r", "Sz7r", "Sz8r", "Sz9r")
repvector <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1)
obsvector <- c(0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
matvector <- c(0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
immvector <- c(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0)
indataset <- c(0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
binvec <- c(0, 4.6, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 
  0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5)

lathframeln <- sf_create(sizes = sizevector, stagenames = stagevector, 
  repstatus = repvector, obsstatus = obsvector, matstatus = matvector, 
  immstatus = immvector, indataset = indataset, binhalfwidth = binvec, 
  propstatus = propvector)

lathvertln <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
  patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9, 
  juvcol = "Seedling1988", sizeacol = "lnVol88", repstracol = "Intactseed88",
  fecacol = "Intactseed88", deadacol = "Dead1988", 
  nonobsacol = "Dormant1988", stageassign = lathframeln, stagesize = "sizea",
  censorcol = "Missing1988", censorkeep = NA, NAas0 = TRUE, censor = TRUE)

lathvertln$feca2 <- round(lathvertln$feca2)
lathvertln$feca1 <- round(lathvertln$feca1)
lathvertln$feca3 <- round(lathvertln$feca3)

lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "mat", "Sd", "Sdl"), 
  stage2 = c("Sd", "Sd", "Sd", "Sd", "Sdl", "rep", "rep"),
  stage1 = c("Sd", "rep", "Sd", "rep", "Sd", "mat", "mat"),
  eststage3 = c(NA, NA, NA, NA, "mat", NA, NA),
  eststage2 = c(NA, NA, NA, NA, "Sdl", NA, NA),
  eststage1 = c(NA, NA, NA, NA, "Sdl", NA, NA),
  givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA),
  multiplier = c(NA, NA, NA, NA, NA, 0.345, 0.054),
  type = c(1, 1, 1, 1, 1, 3, 3), type_t12 = c(1, 2, 1, 2, 1, 1, 1),
  stageframe = lathframeln, historical = TRUE)

lathvertln_adults <- subset(lathvertln, stage2index > 2)
surv_model <- glm(alive3 ~ sizea2 + sizea1 + as.factor(patchid) +
  as.factor(year2), data = lathvertln_adults, family = "binomial")

obs_data <- subset(lathvertln_adults, alive3 == 1)
obs_model <- glm(obsstatus3 ~ as.factor(patchid), data = obs_data,
  family = "binomial")

size_data <- subset(obs_data, obsstatus3 == 1)
siz_model <- lm(sizea3 ~ sizea2 + sizea1 + repstatus1 + as.factor(patchid) +
  as.factor(year2), data = size_data)

reps_model <- glm(repstatus3 ~ sizea2 + sizea1 + as.factor(patchid) +
  as.factor(year2), data = size_data, family = "binomial")

fec_data <- subset(lathvertln_adults, repstatus2 == 1)
fec_model <- glm(feca2 ~ sizea2 + sizea1 + repstatus1 + as.factor(patchid),
  data = fec_data, family = "poisson")

lathvertln_juvs <- subset(lathvertln, stage2index < 3)
jsurv_model <- glm(alive3 ~ as.factor(patchid), data = lathvertln_juvs,
  family = "binomial")

jobs_data <- subset(lathvertln_juvs, alive3 == 1)
jobs_model <- glm(obsstatus3 ~ 1, family = "binomial", data = jobs_data)

jsize_data <- subset(jobs_data, obsstatus3 == 1)
jsiz_model <- lm(sizea3 ~ as.factor(year2), data = jsize_data)

jrepst_model <- 0
jmatst_model <- 1

mod_params <- create_pm(name_terms = TRUE)
mod_params$modelparams[3] <- "patchid"
mod_params$modelparams[4] <- "alive3"
mod_params$modelparams[5] <- "obsstatus3"
mod_params$modelparams[6] <- "sizea3"
mod_params$modelparams[9] <- "repstatus3"
mod_params$modelparams[11] <- "feca2"
mod_params$modelparams[12] <- "sizea2"
mod_params$modelparams[13] <- "sizea1"
mod_params$modelparams[18] <- "repstatus2"
mod_params$modelparams[19] <- "repstatus1"

lathmat3ln <- flefko3(year = "all", patch = "all", data = lathvertln,
  stageframe = lathframeln, supplement = lathsupp3, paramnames = mod_params,
  surv_model = surv_model, obs_model = obs_model, size_model = siz_model,
  repst_model = reps_model, fec_model = fec_model, jsurv_model = jsurv_model,
  jobs_model = jobs_model, jsize_model = jsiz_model,
  jrepst_model = jrepst_model, jmatst_model = jmatst_model, reduce = FALSE)


# Cypripedium example using three size metrics for classification
data(cypdata)

sizevector_f <- c(0, 0, 0, 0, 0, 0, seq(1, 12, by = 1), seq(0, 9, by = 1),
  seq(0, 8, by = 1), seq(0, 7, by = 1), seq(0, 6, by = 1), seq(0, 5, by = 1),
  seq(0, 4, by = 1), seq(0, 3, by = 1), 0, 1, 2, 0, 1, 0, 
  0, 0, 1, 0)
sizebvector_f <- c(0, 0, 0, 0, 0, 0, rep(0, 12), rep(1, 10), rep(2, 9),
  rep(3, 8), rep(4, 7), rep(5, 6), rep(6, 5), rep(7, 4), rep(8, 3), 9, 9, 10, 
  0, 1, 1, 2)
sizecvector_f <- c(0, 0, 0, 0, 0, 0, rep(0, 12), rep(0, 10), rep(0, 9),
  rep(0, 8), rep(0, 7), rep(0, 6), rep(0, 5), rep(0, 4), 0, 0, 0, 0, 0, 0, 
  1, 1, 1, 1)
stagevector_f <- c("DS", "P1", "P2", "P3", "Sdl", "Dorm", "V1 I0 D0",
  "V2 I0 D0", "V3 I0 D0", "V4 I0 D0", "V5 I0 D0", "V6 I0 D0", "V7 I0 D0",
  "V8 I0 D0", "V9 I0 D0", "V10 I0 D0", "V11 I0 D0", "V12 I0 D0", "V0 I1 D0",
  "V1 I1 D0", "V2 I1 D0", "V3 I1 D0", "V4 I1 D0", "V5 I1 D0", "V6 I1 D0",
  "V7 I1 D0", "V8 I1 D0", "V9 I1 D0", "V0 I2 D0", "V1 I2 D0", "V2 I2 D0",
  "V3 I2 D0", "V4 I2 D0", "V5 I2 D0", "V6 I2 D0", "V7 I2 D0", "V8 I2 D0",
  "V0 I3 D0", "V1 I3 D0", "V2 I3 D0", "V3 I3 D0", "V4 I3 D0", "V5 I3 D0",
  "V6 I3 D0", "V7 I3 D0", "V0 I4 D0", "V1 I4 D0", "V2 I4 D0", "V3 I4 D0",
  "V4 I4 D0", "V5 I4 D0", "V6 I4 D0", "V0 I5 D0", "V1 I5 D0", "V2 I5 D0",
  "V3 I5 D0", "V4 I5 D0", "V5 I5 D0", "V0 I6 D0", "V1 I6 D0", "V2 I6 D0",
  "V3 I6 D0", "V4 I6 D0", "V0 I7 D0", "V1 I7 D0", "V2 I7 D0", "V3 I7 D0",
  "V0 I8 D0", "V1 I8 D0", "V2 I8 D0", "V0 I9 D0", "V1 I9 D0", "V0 I10 D0",
  "V0 I0 D1", "V0 I1 D1", "V1 I1 D1", "V0 I2 D1")
repvector_f <- c(0, 0, 0, 0, 0, rep(0, 13), rep(1, 59))
obsvector_f <- c(0, 0, 0, 0, 0, 0, rep(1, 71))
matvector_f <- c(0, 0, 0, 0, 0, rep(1, 72))
immvector_f <- c(0, 1, 1, 1, 1, rep(0, 72))
propvector_f <- c(1, rep(0, 76))
indataset_f <- c(0, 0, 0, 0, 0, rep(1, 72))
binvec_f <- c(0, 0, 0, 0, 0, rep(0.5, 72))
binbvec_f <- c(0, 0, 0, 0, 0, rep(0.5, 72))
bincvec_f <- c(0, 0, 0, 0, 0, rep(0.5, 72))

vertframe_f <- sf_create(sizes = sizevector_f, sizesb = sizebvector_f,
  sizesc = sizecvector_f, stagenames = stagevector_f, repstatus = repvector_f,
  obsstatus = obsvector_f, propstatus = propvector_f, immstatus = immvector_f,
  matstatus = matvector_f, indataset = indataset_f, binhalfwidth = binvec_f,
  binhalfwidthb = binbvec_f, binhalfwidthc = bincvec_f)

vert_data_f <- verticalize3(cypdata, noyears = 6, firstyear = 2004,
  individcol = "plantid", blocksize = 4, sizeacol = "Veg.04",
  sizebcol = "Inf.04", sizeccol = "Inf2.04", repstracol = "Inf.04",
  repstrbcol = "Inf2.04", fecacol = "Pod.04", censorcol = "censor",
  censorkeep = 1, censorRepeat = FALSE, stageassign = vertframe_f,
  stagesize = "sizeabc", NAas0 = TRUE, censor = FALSE)

vertsupp3f <- supplemental(stage3 = c("DS", "P1", "DS", "P1", "P2", "P2", "P3",
    "Sdl", "Sdl", "Sdl", "Dorm", "V1 I0 D0", "V2 I0 D0", "V3 I0 D0", "Dorm",
    "V1 I0 D0", "V2 I0 D0", "V3 I0 D0", "mat", "mat", "mat", "mat", "DS", "P1"),
  stage2 = c("DS", "DS", "DS", "DS", "P1", "P1", "P2", "P3", "Sdl", "Sdl", "Sdl",
    "Sdl", "Sdl", "Sdl", "Sdl", "Sdl", "Sdl", "Sdl", "Dorm", "V1 I0 D0",
    "V2 I0 D0", "V3 I0 D0", "rep", "rep"),
  stage1 = c("DS", "DS", "rep", "rep", "DS", "rep", "P1", "P2", "P3", "Sdl",
    "Sdl", "Sdl", "Sdl", "Sdl", "P3", "P3", "P3", "P3", "Sdl", "Sdl", "Sdl",
    "Sdl", "mat", "mat"),
  eststage3 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "Dorm", "V1 I0 D0",
    "V2 I0 D0", "V3 I0 D0", "Dorm", "V1 I0 D0", "V2 I0 D0", "V3 I0 D0", "mat",
    "mat", "mat", "mat", NA, NA), 
  eststage2 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "V1 I0 D0", "V1 I0 D0",
    "V1 I0 D0", "V1 I0 D0", "V1 I0 D0", "V1 I0 D0", "V1 I0 D0", "V1 I0 D0",
    "Dorm", "V1 I0 D0", "V2 I0 D0", "V3 I0 D0", NA, NA),
  eststage1 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "V1 I0 D0", "V1 I0 D0",
    "V1 I0 D0", "V1 I0 D0", "V1 I0 D0", "V1 I0 D0", "V1 I0 D0", "V1 I0 D0",
    "V1 I0 D0", "V1 I0 D0", "V1 I0 D0", "V1 I0 D0", NA, NA),
  givenrate = c(0.10, 0.20, 0.10, 0.20, 0.20, 0.20, 0.20, 0.25, 0.40, 0.40, NA,
    NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA),
  multiplier = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
    NA, NA, NA, NA, NA, NA, NA, 0.5 * 5000, 0.5 * 5000),
  type =c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    3, 3),
  type_t12 = c(1, 1, 2, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 1),
  stageframe = vertframe_f, historical = TRUE)

surv_model <- glm(alive3 ~ sizea2 + sizeb2, data = vert_data_f,
  family = "binomial")

obs_data <- subset(vert_data_f, alive3 == 1)
obs_model <- glm(obsstatus3 ~ sizeb2 + sizec1 + as.factor(year2),
  data = obs_data, family = "binomial")

size_data <- subset(obs_data, obsstatus3 == 1)
siz_model <- MASS::glm.nb(sizea3 ~ sizea2 + sizea1 + sizeb1, data = size_data)
sizb_model <- glm(sizeb3 ~ sizea2 + sizeb2 + sizec1 + repstatus2 + repstatus1 +
  as.factor(year2), data = size_data, family = "poisson")
sizc_model <- glm(sizec3 ~ sizea1 + repstatus2, data = size_data,
  family = "poisson")

reps_model <- glm(repstatus3 ~ sizea2 + sizeb2 + repstatus2 + as.factor(year2),
  data = size_data, family = "binomial")

fec_data <- subset(vert_data_f, repstatus2 == 1)
fec_model <- glm(feca2 ~ sizeb2 + as.factor(year2), data = fec_data,
  family = "poisson")

mod_params <- create_pm(name_terms = TRUE)
mod_params$modelparams[3] <- "patchid"
mod_params$modelparams[4] <- "alive3"
mod_params$modelparams[5] <- "obsstatus3"
mod_params$modelparams[6] <- "sizea3"
mod_params$modelparams[9] <- "repstatus3"
mod_params$modelparams[11] <- "feca2"
mod_params$modelparams[12] <- "sizea2"
mod_params$modelparams[13] <- "sizea1"
mod_params$modelparams[18] <- "repstatus2"
mod_params$modelparams[19] <- "repstatus1"

vert_mats_f3 <- flefko3(stageframe = vertframe_f, supplement = vertsupp3f,
  data = vert_data_f, surv_model = surv_model, obs_model = obs_model,
  size_model = siz_model, sizeb_model = sizb_model, sizec_model = sizc_model,
  repst_model = reps_model, fec_model = fec_model, paramnames = mod_params,
  sparse_output = TRUE)