# R/popdyn.R In lefko3: Historical and Ahistorical Population Projection Matrix Analysis

#### Documented in elasticity3elasticity3.dgCMatrixelasticity3.lefkoMatelasticity3.listelasticity3.matrixltre3plot.lefkoProjrepvalue3repvalue3.dgCMatrixrepvalue3.lefkoMatrepvalue3.listrepvalue3.matrixsensitivity3sensitivity3.dgCMatrixsensitivity3.lefkoMatsensitivity3.listsensitivity3.matrixstablestage3stablestage3.dgCMatrixstablestage3.lefkoMatstablestage3.liststablestage3.matrixsummary.lefkoElassummary.lefkoLTREsummary.lefkoProj

#' Estimate Stable Stage Distribution
#'
#' \code{stablestage3()} is a generic function that returns the stable stage
#' distribution for a population projection matrix or set of matrices. This
#' function is made to handle very large and sparse matrices supplied as
#' \code{lefkoMat} objects or as individual matrices, and can be used with large
#' historical matrices, IPMs, age x stage matrices, as well as ahistorical
#' matrices.
#'
#' @name stablestage3
#'
#' @param mats A lefkoMat object, a population projection matrix, or a list of
#' population projection matrices for which the stable stage distribution is
#' desired.
#' @param ... Other parameters.
#'
#' @return The value returned depends on the class of the \code{mats} argument.
#' See related functions for details.
#'
#'
#' @examples
#' # Lathyrus deterministic example
#' data(lathyrus)
#'
#' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0)
#' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm")
#' repvector <- c(0, 0, 0, 0, 0, 1, 0)
#' obsvector <- c(0, 1, 1, 1, 1, 1, 0)
#' matvector <- c(0, 0, 1, 1, 1, 1, 1)
#' immvector <- c(1, 1, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5)
#'
#' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   immstatus = immvector, indataset = indataset, binhalfwidth = binvec,
#'   propstatus = propvector)
#'
#' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
#'   patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
#'   juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88",
#'   nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea",
#'   censorcol = "Missing1988", censorkeep = NA, censor = TRUE)
#'
#' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"),
#'   stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"),
#'   stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"),
#'   eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"),
#'   eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"),
#'   eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"),
#'   givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA),
#'   type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1),
#'   stageframe = lathframe, historical = TRUE)
#'
#' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all",
#'   stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3,
#'   yearcol = "year2", indivcol = "individ")
#'
#' ehrlen3mean <- lmean(ehrlen3)
#' stablestage3(ehrlen3mean)
#'
#' # Cypripedium stochastic example
#' data(cypdata)
#'
#' sizevector <- c(0, 0, 0, 0, 0, 0, 1, 2.5, 4.5, 8, 17.5)
#' stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "XSm", "Sm", "Md", "Lg",
#'   "XLg")
#' repvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
#' obsvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
#' matvector <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)
#' immvector <- c(0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 0, 0, 0, 0, 0.5, 0.5, 1, 1, 2.5, 7)
#'
#' cypframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   propstatus = propvector, immstatus = immvector, indataset = indataset,
#'   binhalfwidth = binvec)
#'
#' cypraw_v1 <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004,
#'   patchidcol = "patch", individcol = "plantid", blocksize = 4,
#'   sizeacol = "Inf2.04", sizebcol = "Inf.04", sizeccol = "Veg.04",
#'   repstracol = "Inf.04", repstrbcol = "Inf2.04", fecacol = "Pod.04",
#'   stageassign = cypframe_raw, stagesize = "sizeadded", NAas0 = TRUE,
#'   NRasRep = TRUE)
#'
#' # Here we use supplemental() to provide overwrite and reproductive info
#' cypsupp2r <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "D",
#'     "XSm", "Sm", "SD", "P1"),
#'   stage2 = c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "rep",
#'     "rep"),
#'   eststage3 = c(NA, NA, NA, NA, NA, "D", "XSm", "Sm", NA, NA),
#'   eststage2 = c(NA, NA, NA, NA, NA, "XSm", "XSm", "XSm", NA, NA),
#'   givenrate = c(0.10, 0.20, 0.20, 0.20, 0.25, NA, NA, NA, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, NA, NA, NA, NA, 0.5, 0.5),
#'   type =c(1, 1, 1, 1, 1, 1, 1, 1, 3, 3),
#'   stageframe = cypframe_raw, historical = FALSE)
#'
#' cypmatrix2r <- rlefko2(data = cypraw_v1, stageframe = cypframe_raw,
#'   year = "all", patch = "all", stages = c("stage3", "stage2", "stage1"),
#'   yearcol = "year2", patchcol = "patchid", indivcol = "individ")
#'
#' stablestage3(cypmatrix2r, stochastic = TRUE)
#'
#' @export
stablestage3 <- function(mats, ...) UseMethod("stablestage3")

#' Estimate Stable Stage Distribution of Matrices in lefkoMat Object
#'
#' \code{stablestage3.lefkoMat()} returns the deterministic stable stage
#' distributions of all \code{A} matrices in an object of class \code{lefkoMat},
#' as well as the long-run projected mean stage distribution in stochastic
#' analysis. This function can handle large and sparse matrices, and so can be
#' used with large historical matrices, IPMs, age x stage matrices, as well as
#' ahistorical matrices.
#'
#' @name stablestage3.lefkoMat
#'
#' @param mats An object of class \code{lefkoMat}.
#' @param stochastic A logical value indicating whether to use deterministic
#' (\code{FALSE}) or stochastic (\code{TRUE}) analysis. Defaults to
#' \code{FALSE}.
#' @param times An integer variable indicating number of occasions to project if
#' using stochastic analysis. Defaults to 10000.
#' @param tweights An optional numeric vector or matrix denoting the
#' probabilities of choosing each matrix in a stochastic projection. If a matrix
#' is input, then a first-order Markovian environment is assumed, in which the
#' probability of choosing a specific annual matrix depends on which annual
#' matrix is currently chosen. If a vector is input, then the choice of annual
#' matrix is assumed to be independent of the current matrix. Defaults to equal
#' weighting among matrices.
#' @param seed A number to use as a random number seed in stochastic projection.
#' @param force_sparse A text string indicating whether to use sparse matrix
#' encoding (\code{"yes"}) if standard matrices are provided. Defaults to
#' \code{"auto"}, in which case sparse matrix encoding is used with square
#' matrices with at least 50 rows and no more than 50\% of elements with values
#' greater than zero.
#' @param ... Other parameters.
#'
#' @return This function returns the stable stage distributions (and long-run
#' mean stage distributions in stochastic analysis) corresponding to the
#' matrices in a \code{lefkoMat} object.
#'
#' The output depends on whether the \code{lefkoMat} object used as input is
#' ahistorical or historical, and whether the analysis is deterministic or
#' stochastic. If deterministic and ahistorical, then a single data frame is
#' output, which includes the number of the matrix within the \code{A} element
#' of the input \code{lefkoMat} object, followed by the stage id (numeric and
#' assigned through \code{\link{sf_create}()}), the stage name, and the
#' estimated proportion of the stable stage distribution (\code{ss_prop}). If
#' stochastic and ahistorical, then a single data frame is output starting with
#' the number of the population-patch (\code{matrix_set}), a string
#' concatenating the names of the population and the patch (\code{poppatch}),
#' the assigned stage id number (\code{stage_id}), and the stage name
#' (\code{stage}), and the long-run average stage distribution (\code{ss_prop}).
#'
#' If a historical matrix is used as input, then two data frames are output
#' into a list object. The \code{hist} element describes the historical
#' stage-pair distribution, while the \code{ahist} element describes the stage
#' distribution. If deterministic, then \code{hist} contains a data frame
#' including the matrix number (\code{matrix}), the numeric stage designations for
#' stages in occasions \emph{t} and \emph{t}-1, (\code{stage_id_2} and
#' \code{stage_id_1}, respectively), followed by the respective stage names (
#' \code{stage_2} and \code{stage_1}), and ending with the estimated stable
#' stage-pair distribution. The associated \code{ahist} element is as before. If
#' stochastic, then the \code{hist} element contains a single data frame with
#' the number of the population-patch (\code{matrix_set}), a string
#' concatenating the names of the population and the patch (\code{poppatch}),
#' the assigned stage id numbers in times \emph{t} and \emph{t}-1 (
#' \code{stage_id_2} and \code{stage_id_2}, respectively), and the associated
#' stage names (\code{stage_2} and \code{stage_1}, respectively), and the
#' long-run average stage distribution (\code{ss_prop}). The associated
#' \code{ahist} element is as before in the ahistorical, stochastic case.
#'
#' In addition to the data frames noted above, stochastic analysis will result
#' in the additional output of a list of matrices containing the actual
#' projected stage distributions across all projected occasions, in the order of
#' population-patch combinations in the \code{lefkoMat} input.
#'
#' @section Notes:
#' In stochastic analysis, the projected mean distribution is the arithmetic
#' mean across the final 1000 projected occasions if the simulation is at least
#' 2000 projected occasions long. If between 500 and 2000 projected occasions
#' long, then only the final 200 are used, and if fewer than 500 occasions are
#' used, then all are used. Note that because stage distributions in stochastic
#' simulations can change greatly in the initial portion of the run, we
#' encourage a minimum of 2000 projected occasions per simulation, with 10000
#' preferred.
#'
#' Speed can sometimes be increased by shifting from automatic sparse matrix
#' determination to forced dense or sparse matrix projection. This will most
#' likely occur when matrices have between 30 and 300 rows and columns.
#' Defaults work best when matrices are very small and dense, or very large and
#' sparse.
#'
#'
#' @examples
#' # Lathyrus deterministic example
#' data(lathyrus)
#'
#' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0)
#' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm")
#' repvector <- c(0, 0, 0, 0, 0, 1, 0)
#' obsvector <- c(0, 1, 1, 1, 1, 1, 0)
#' matvector <- c(0, 0, 1, 1, 1, 1, 1)
#' immvector <- c(1, 1, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5)
#'
#' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   immstatus = immvector, indataset = indataset, binhalfwidth = binvec,
#'   propstatus = propvector)
#'
#' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
#'   patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
#'   juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88",
#'   nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea",
#'   censorcol = "Missing1988", censorkeep = NA, censor = TRUE)
#'
#' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"),
#'   stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"),
#'   stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"),
#'   eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"),
#'   eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"),
#'   eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"),
#'   givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA),
#'   type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1),
#'   stageframe = lathframe, historical = TRUE)
#'
#' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all",
#'   stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3,
#'   yearcol = "year2", indivcol = "individ")
#'
#' ehrlen3mean <- lmean(ehrlen3)
#' stablestage3(ehrlen3mean)
#'
#' # Cypripedium stochastic example
#' data(cypdata)
#'
#' sizevector <- c(0, 0, 0, 0, 0, 0, 1, 2.5, 4.5, 8, 17.5)
#' stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "XSm", "Sm", "Md", "Lg",
#'   "XLg")
#' repvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
#' obsvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
#' matvector <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)
#' immvector <- c(0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 0, 0, 0, 0, 0.5, 0.5, 1, 1, 2.5, 7)
#'
#' cypframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   propstatus = propvector, immstatus = immvector, indataset = indataset,
#'   binhalfwidth = binvec)
#'
#' cypraw_v1 <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004,
#'   patchidcol = "patch", individcol = "plantid", blocksize = 4,
#'   sizeacol = "Inf2.04", sizebcol = "Inf.04", sizeccol = "Veg.04",
#'   repstracol = "Inf.04", repstrbcol = "Inf2.04", fecacol = "Pod.04",
#'   stageassign = cypframe_raw, stagesize = "sizeadded", NAas0 = TRUE,
#'   NRasRep = TRUE)
#'
#' # Here we use supplemental() to provide overwrite and reproductive info
#' cypsupp2r <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "D",
#'     "XSm", "Sm", "SD", "P1"),
#'   stage2 = c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "rep",
#'     "rep"),
#'   eststage3 = c(NA, NA, NA, NA, NA, "D", "XSm", "Sm", NA, NA),
#'   eststage2 = c(NA, NA, NA, NA, NA, "XSm", "XSm", "XSm", NA, NA),
#'   givenrate = c(0.10, 0.20, 0.20, 0.20, 0.25, NA, NA, NA, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, NA, NA, NA, NA, 0.5, 0.5),
#'   type =c(1, 1, 1, 1, 1, 1, 1, 1, 3, 3),
#'   stageframe = cypframe_raw, historical = FALSE)
#'
#' cypmatrix2r <- rlefko2(data = cypraw_v1, stageframe = cypframe_raw,
#'   year = "all", patch = "all", stages = c("stage3", "stage2", "stage1"),
#'   yearcol = "year2", patchcol = "patchid", indivcol = "individ")
#'
#' stablestage3(cypmatrix2r, stochastic = TRUE)
#'
#' @export
stablestage3.lefkoMat <- function(mats, stochastic = FALSE, times = 10000,
tweights = NA, seed = NA, force_sparse = "auto", ...) {

matrix_set <- theprohecy <- NULL
sparsemethod <- sparse_auto <- sparse_input <- FALSE

if (is(mats$A[[1]], "dgCMatrix")) sparse_input <- TRUE if (!sparse_input) { if (is.logical(force_sparse)) { if (force_sparse) { sparsemethod <- TRUE } else sparsemethod <- FALSE } else if (is.element(tolower(force_sparse), c("y", "yes", "yea", "yeah", "t", "true", "ja", "tak"))) { sparsemethod <- TRUE } else if (is.element(tolower(force_sparse), c("n", "no", "non", "nah", "f", "false", "nein", "nie"))) { sparsemethod <- FALSE } else { if (is.element(tolower(force_sparse), c("au", "aut", "auto"))) sparse_auto <- TRUE elements_total <- length(mats$A[[1]])
dense_elements <- length(which(mats$A[[1]] != 0)) if ((dense_elements / elements_total) < 0.5 & elements_total > 2500) { sparsemethod <- 1 } else sparsemethod <- 0 } } if (!stochastic) { baldrick <- if (is.matrix(mats$A)) {
if (!sparse_input) {
.ss3matrix(mats$A, sparsemethod) } else .ss3matrix_sp(mats$A)

} else if (is.list(mats$A)) { if (!sparse_input) { unlist(lapply(mats$A, .ss3matrix, sparsemethod))
} else unlist(lapply(mats$A, .ss3matrix_sp)) } else { stop("Input not recognized.") } } else { if (!is.na(seed[1])) { set.seed(seed[1]) } mats$labels$poppatch <- paste(mats$labels$pop, mats$labels$patch) used_poppatches <- as.list(unique(mats$labels$poppatch)) #Here we get the full stage distribution series for all occasions, as a list princegeorge <- lapply(used_poppatches, function(X) { used_slots <- which(mats$labels$poppatch == X) if (length(used_slots) < 2) { warning("Only 1 annual matrix found for some population-patch combinations. Stochastic analysis requires multiple annual matrices per population-patch combination.", call. = FALSE) } if (!any(is.na(tweights))) { tw_error_1 <- "Option tweights must be NA, a numeric vector equal to the" tw_error_2 <- "number of years or matrices supplied, or a square matrix" tw_error_3 <- "with dimensions equal to the number of matrices per patch." tw_error <- paste(tw_error_1, tw_error_2, tw_error_3) if (is.matrix(tweights)) { theprophecy <- rep(0, times) used_weights <- tweights[, 1] for (i in c(1:times)) { used_weights <- used_weights / sum(used_weights) chosen_timepath <- sample(used_slots, 1, replace = TRUE, prob = used_weights) theprophecy[i] <- chosen_timepath[1] - 1 used_weights <- tweights[, chosen_timepath[1]] } } else if (is.numeric(tweights)) { if (length(tweights) != length(used_slots)) { if (length(tweights) == length(mats$A)) {
used_weights <- tweights[used_slots] / sum(tweights[used_slots])
} else {
stop(tw_error, call. = FALSE)
}
} else {
used_weights <- tweights / sum(tweights)
}

theprophecy <- sample(used_slots, times, replace = TRUE,
prob = used_weights) - 1
} else {
stop(tw_error, call. = FALSE)
}
} else {
used_weights <- rep(1, length(used_slots))
used_weights <- used_weights / sum(used_weights)

theprophecy <- sample(used_slots, times, replace = TRUE,
prob = used_weights) - 1
}

#theprophecy <- sample(used_slots, times, replace = TRUE, prob = used_weights) - 1
starter <- if (!sparse_input) {
.ss3matrix(mats$A[[used_slots[1]]], sparsemethod) } else .ss3matrix_sp(mats$A[[used_slots[1]]])

theseventhmatrix <- if (!sparse_input) {
.proj3(starter, mats$A, theprophecy, 1, 0, 0, sparse_auto, sparsemethod) } else .proj3sp(starter, mats$A, theprophecy, 1, 0, 0)
ssonly <- theseventhmatrix[((dim(mats$A[[1]])[1]) + 1):(2 *(dim(mats$A[[1]])[1])),]

return(ssonly)
})

# Create mean distributions
baldrick <- unlist(
lapply(princegeorge, function(X) {
if (times > 2000) {
usedX <- X[,(times - 999):(times)]
} else if (times > 500) {
usedX <- X[,(times-199):(times)]
} else {
usedX <- X
}
apply(usedX, 1, mean)
})
)
}

# The final bits sort everything and clean it up, and create the ahistorical
# version if a historical input was used
if (is.list(mats$A)) { if (!stochastic) { multiplier <- length(mats$A)
} else {
multiplier <- length(used_poppatches)
}
} else multiplier <- 1

if (all(is.na(mats$hstages))) { labels_orig <- mats$ahstages[,1:2]
mat_dims <- dim(mats$A[[1]])[1] if (dim(labels_orig)[1] == mat_dims) { labels <- labels_orig } else { newmult <- mat_dims / dim(labels_orig)[1] if (mat_dims %% dim(labels_orig)[1] != 0) { stop("Matrices do not appear to be ahistorical, historical, or age x stage. Cannot proceed. Please make sure that matrix dimensions match stage descriptions.", call. = FALSE) } if (!all(is.na(mats$agestages))) {
check_min <- min(mats$agestages$age)
}
age_bit <- c(apply(as.matrix(c((0 + check_min):(newmult + check_min - 1))), 1, rep, dim(labels_orig)[1]))
core_labels <- cbind.data.frame(age_bit, do.call("rbind.data.frame", replicate(newmult, labels_orig, simplify = FALSE)))
core_labels$agestage_id <- c(1:length(core_labels$stage_id))
core_labels$agestage <- apply(as.matrix(core_labels$agestage_id), 1, function(X) {
paste(core_labels$age_bit[X], core_labels$stage[X])
})
labels <- core_labels
names(labels) <- c("age", "stage_id", "stage", "agestage_id", "agestage")
}

if (!stochastic) {
modlabels <- cbind.data.frame(as.matrix(rep(c(1:multiplier), each = dim(labels)[1])),
do.call("rbind.data.frame", apply(as.matrix(c(1:multiplier)), 1, function(X) return(labels))))
output <- cbind.data.frame(modlabels, baldrick)
if (is.element("age", names(labels))) {
names(output) <- c("matrix", "age", "stage_id", "stage", "agestage_id", "agestage", "ss_prop")
} else {
names(output) <- c("matrix", "stage_id", "stage", "ss_prop")
}
} else {
modlabels <- cbind.data.frame(as.matrix(rep(unlist(used_poppatches), each = dim(labels)[1])),
do.call("rbind.data.frame", apply(as.matrix(c(1:multiplier)), 1, function(X) return(labels))))
unipoppatches <- unique(modlabels[,1])

modnums <- apply(as.matrix(modlabels[,1]), 1, function(X) {
return(which(unipoppatches == X))
})

output <- cbind.data.frame(modnums, modlabels, baldrick)
if (is.element("age", names(labels))) {
names(output) <- c("matrix_set", "poppatch", "age", "stage_id", "stage",
"agestage_id", "agestage", "ss_prop")
} else {
names(output) <- c("matrix_set", "poppatch", "stage_id", "stage",
"ss_prop")
}
}
rownames(output) <- c(1:dim(output)[1])
} else {
labels <- mats$hstages if (!stochastic) { modlabels <- cbind.data.frame(as.matrix(rep(c(1:multiplier), each = dim(labels)[1])), do.call("rbind.data.frame", apply(as.matrix(c(1:multiplier)), 1, function(X) return(labels)))) outputh <- cbind.data.frame(modlabels, baldrick) names(outputh) <- c("matrix", "stage_id_2", "stage_id_1", "stage_2", "stage_1", "ss_prop") } else { modlabels <- cbind.data.frame(as.matrix(rep(unlist(used_poppatches), each = dim(labels)[1])), do.call("rbind.data.frame", apply(as.matrix(c(1:multiplier)), 1, function(X) return(labels)))) unipoppatches <- unique(modlabels[,1]) modnums <- apply(as.matrix(modlabels[,1]), 1, function(X) { return(which(unipoppatches == X)) }) outputh <- cbind.data.frame(modnums, modlabels, baldrick) names(outputh) <- c("matrix_set", "poppatch", "stage_id_2", "stage_id_1", "stage_2", "stage_1", "ss_prop") } rownames(outputh) <- c(1:dim(outputh)[1]) ahlabels <- mats$ahstages[,c("stage_id", "stage")]
ss2 <- c(apply(as.matrix(c(1:multiplier)), 1, function(X) {
rightset <- if (!stochastic) {
subset(outputh, matrix == X)
} else {
subset(outputh, matrix_set == X)
}
apply(as.matrix(ahlabels[,1]), 1, function(Y) {
sum(rightset$ss_prop[which(rightset$stage_id_2 == Y)])
})
}))

if (!stochastic) {
outputah <- cbind.data.frame(rep(c(1:multiplier), each = length(ahlabels[,1])),
rep(ahlabels[,1], multiplier), rep(ahlabels[,2], multiplier), ss2)
names(outputah) <- c("matrix", "stage_id", "stage", "ss_prop")
} else {
modlabels <- cbind.data.frame(as.matrix(rep(unlist(used_poppatches), each = dim(ahlabels)[1])),
do.call("rbind.data.frame", apply(as.matrix(c(1:multiplier)), 1, function(X) return(ahlabels))))
modnums <- apply(as.matrix(modlabels[,1]), 1, function(X) {
return(which(unipoppatches == X))
})

outputah <- cbind.data.frame(modnums, modlabels, ss2)
names(outputah) <- c("matrix_set", "poppatch", "stage_id", "stage", "ss_prop")
}

rownames(outputah) <- c(1:dim(outputah)[1])

if (!stochastic) {
output <- list(hist = outputh, ahist = outputah)
} else {
output <- list(hist = outputh, ahist = outputah, projections = princegeorge)
}
}

return(output)
}

#' Estimate Stable Stage Distribution of a Single Population Projection Matrix
#'
#' \code{stablestage3.matrix()} returns the stable stage distribution for a
#' population projection matrix. This function can handle large and sparse
#' matrices, and so can be used with large historical matrices, IPMs,
#' age x stage matrices, as well as smaller ahistorical matrices.
#'
#' @name stablestage3.matrix
#'
#' @param mats A population projection matrix of class \code{matrix}.
#' @param force_sparse A text string indicating whether to use sparse matrix
#' encoding (\code{"yes"}) when supplied with standard matrices. Defaults to
#' \code{"auto"}, in which case sparse matrix encoding is used with square
#' matrices with at least 50 rows and no more than 50\% of elements with values
#' greater than zero.
#' @param ... Other parameters.
#'
#' @return This function returns the stable stage distribution corresponding to
#' the input matrix.
#'
#' @section Notes:
#' Speed can sometimes be increased by shifting from automatic sparse matrix
#' determination to forced dense or sparse matrix projection. This will most
#' likely occur when matrices have between 30 and 300 rows and columns.
#' Defaults work best when matrices are very small and dense, or very large and
#' sparse.
#'
#'
#' @examples
#' data(lathyrus)
#'
#' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0)
#' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm")
#' repvector <- c(0, 0, 0, 0, 0, 1, 0)
#' obsvector <- c(0, 1, 1, 1, 1, 1, 0)
#' matvector <- c(0, 0, 1, 1, 1, 1, 1)
#' immvector <- c(1, 1, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5)
#'
#' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   immstatus = immvector, indataset = indataset, binhalfwidth = binvec,
#'   propstatus = propvector)
#'
#' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
#'   patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
#'   juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88",
#'   nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea",
#'   censorcol = "Missing1988", censorkeep = NA, censor = TRUE)
#'
#' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"),
#'   stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"),
#'   stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"),
#'   eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"),
#'   eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"),
#'   eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"),
#'   givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA),
#'   type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1),
#'   stageframe = lathframe, historical = TRUE)
#'
#' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all",
#'   stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3,
#'   yearcol = "year2", indivcol = "individ")
#'
#' ehrlen3mean <- lmean(ehrlen3)
#' stablestage3(ehrlen3mean$A[[1]]) #' #' @export stablestage3.matrix <- function(mats, force_sparse = "auto", ...) { sparsemethod <- 0 sparse_input <- FALSE if (is(mats, "dgCMatrix")) sparse_input <- TRUE if (is.logical(force_sparse) & !sparse_input) { if (force_sparse) { sparsemethod <- 1 } else sparsemethod <- 0 } else if (is.element(tolower(force_sparse), c("y", "yes", "yea", "yeah", "t", "true", "ja", "tak"))) { sparsemethod <- 1 } else if (is.element(tolower(force_sparse), c("n", "no", "non", "nah", "f", "false", "nein", "nie"))) { sparsemethod <- 0 } else { elements_total <- length(mats) dense_elements <- length(which(mats != 0)) if ((dense_elements / elements_total) < 0.5 & elements_total > 2500) { sparsemethod <- 1 } else sparsemethod <- 0 } if (!sparse_input) { wcorr <- .ss3matrix(mats, sparsemethod) } else wcorr <- .ss3matrix_sp(mats) return(wcorr) } #' Estimate Stable Stage Distribution of a Single Population Projection Matrix #' #' \code{stablestage3.dgCMatrix()} returns the stable stage distribution for a #' sparse population projection matrix. #' #' @name stablestage3.dgCMatrix #' #' @param mats A population projection matrix of class \code{dgCMatrix}. #' @param ... Other parameters. #' #' @return This function returns the stable stage distribution corresponding to #' the input matrix. #' #' @seealso \code{\link{stablestage3}()} #' @seealso \code{\link{stablestage3.lefkoMat}()} #' @seealso \code{\link{stablestage3.list}()} #' @seealso \code{\link{stablestage3.matrix}()} #' #' @examples #' data(lathyrus) #' #' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0) #' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm") #' repvector <- c(0, 0, 0, 0, 0, 1, 0) #' obsvector <- c(0, 1, 1, 1, 1, 1, 0) #' matvector <- c(0, 0, 1, 1, 1, 1, 1) #' immvector <- c(1, 1, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5) #' #' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' immstatus = immvector, indataset = indataset, binhalfwidth = binvec, #' propstatus = propvector) #' #' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988, #' patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9, #' juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88", #' fecacol = "Intactseed88", deadacol = "Dead1988", #' nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea", #' censorcol = "Missing1988", censorkeep = NA, censor = TRUE) #' #' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"), #' stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"), #' stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"), #' eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"), #' eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"), #' eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"), #' givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA), #' multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA), #' type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1), #' stageframe = lathframe, historical = TRUE) #' #' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all", #' stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3, #' yearcol = "year2", indivcol = "individ", sparse_output = TRUE) #' #' stablestage3(ehrlen3$A[[1]])
#'
#' @export
stablestage3.dgCMatrix <- function(mats, ...)
{

wcorr <- .ss3matrix_sp(mats)

return(wcorr)
}

#' Estimate Stable Stage Distribution of a List of Projection Matrices
#'
#' \code{stablestage3.list()} returns the stable stage distributions for stages
#' in population projection matrices arranged in a general list. The function
#' makes no assumptions about whether the matrix is ahistorical and simply
#' provides stable stage distribution values corresponding to each row, meaning
#' that the overall stable stage distribution of basic life history stages in a
#' historical matrix are not provided (the \code{\link{stablestage3.lefkoMat}()}
#' historical estimates these on the basis of stage description information
#' provided in the \code{lefkoMat} object used as input in that function). This
#' provided in the handle large and sparse matrices, and so can be used with
#' large historical matrices, IPMs, age x stage matrices, as well as smaller
#' ahistorical matrices.
#'
#' @name stablestage3.list
#'
#' @param mats A list of population projection matrices, all in either class
#' \code{matrix} or class \code{dgCMatrix}.
#' @param stochastic A logical value indicating whether to use deterministic
#' (\code{FALSE}) or stochastic (\code{TRUE}) analysis. Defaults to
#' \code{FALSE}.
#' @param times An integer variable indicating number of occasions to project if
#' using stochastic analysis. Defaults to 10000.
#' @param tweights An optional numeric vector or matrix denoting the
#' probabilities of choosing each matrix in a stochastic projection. If a matrix
#' is input, then a first-order Markovian environment is assumed, in which the
#' probability of choosing a specific annual matrix depends on which annual
#' matrix is currently chosen. If a vector is input, then the choice of annual
#' matrix is assumed to be independent of the current matrix. Defaults to equal
#' weighting among matrices.
#' @param seed A number to use as a random number seed in stochastic projection.
#' @param force_sparse A text string indicating whether to use sparse matrix
#' encoding (\code{"yes"}) when supplied with standard matrices. Defaults to
#' \code{"auto"}, in which case sparse matrix encoding is used with square
#' matrices with at least 50 rows and no more than 50\% of elements with values
#' greater than zero.
#' @param ... Other parameters.
#'
#' @return This function returns a list of vector data frames characterizing the
#' stable stage distributions for stages of each population projection matrix.
#'
#' @section Notes:
#' Speed can sometimes be increased by shifting from automatic sparse matrix
#' determination to forced dense or sparse matrix projection. This will most
#' likely occur when matrices have between 30 and 300 rows and columns.
#' Defaults work best when matrices are very small and dense, or very large and
#' sparse.
#'
#'
#' @examples
#' data(lathyrus)
#'
#' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0)
#' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm")
#' repvector <- c(0, 0, 0, 0, 0, 1, 0)
#' obsvector <- c(0, 1, 1, 1, 1, 1, 0)
#' matvector <- c(0, 0, 1, 1, 1, 1, 1)
#' immvector <- c(1, 1, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5)
#'
#' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   immstatus = immvector, indataset = indataset, binhalfwidth = binvec,
#'   propstatus = propvector)
#'
#' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
#'   patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
#'   juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88",
#'   nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea",
#'   censorcol = "Missing1988", censorkeep = NA, censor = TRUE)
#'
#' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"),
#'   stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"),
#'   stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"),
#'   eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"),
#'   eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"),
#'   eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"),
#'   givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA),
#'   type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1),
#'   stageframe = lathframe, historical = TRUE)
#'
#' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all",
#'   stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3,
#'   yearcol = "year2", indivcol = "individ")
#'
#' ehrlen3mean <- lmean(ehrlen3)
#' stablestage3(ehrlen3mean$A) #' #' @export stablestage3.list <- function(mats, stochastic = FALSE, times = 10000, tweights = NA, seed = NA, force_sparse = "auto", ...) { sparsemethod <- 0 sparse_initial <- FALSE w_list <- Xlist <- NULL if (is(mats[[1]], "dgCMatrix")) sparse_initial <- TRUE if (!sparse_initial) { if (is.logical(force_sparse)) { if (force_sparse) { sparsemethod <- 1 } else sparsemethod <- 0 } else if (is.element(tolower(force_sparse), c("y", "yes", "yea", "yeah", "t", "true", "ja", "tak"))) { sparsemethod <- 1 } else if (is.element(tolower(force_sparse), c("n", "no", "non", "nah", "f", "false", "nein", "nie"))) { sparsemethod <- 0 } else { elements_total <- length(mats[[1]]) dense_elements <- length(which(mats[[1]] != 0)) if ((dense_elements / elements_total) < 0.5 & elements_total > 2500) { sparsemethod <- 1 } else sparsemethod <- 0 } } list_length <- length(mats) massive_panic <- apply(as.matrix(c(1:list_length)), 1, function(X) { if (is.matrix(mats[[X]])) { outchunk <- 1 } else if (is(mats[[X]], "dgCMatrix")) { outchunk <- 2 } else outchunk <- 3 }) if (!all(massive_panic == 1) & !all(massive_panic == 2)) { stop("Matrix list must be entirely in standard matrix format, or entirely in dgCMatrix format.", call. = FALSE) } if (!stochastic) { w_list <- lapply(mats, function(X) { if (is.matrix(X)) { output <- .ss3matrix(X, sparsemethod) } else if (is(X, "dgCMatrix")) { output <- .ss3matrix_sp(X) } else { stop("Unrecognized input in object mats.", FALSE) } }) } else { if (!is.na(seed[1])) { set.seed(seed[1]) } used_slots <- c(1:list_length) if (!any(is.na(tweights))) { tw_error_1 <- "Option tweights must be NA, a numeric vector equal to the" tw_error_2 <- "number of matrices supplied." tw_error <- paste(tw_error_1, tw_error_2) if (is.matrix(tweights)) { theprophecy <- rep(0, times) used_weights <- tweights[, 1] for (i in c(1:times)) { used_weights <- used_weights / sum(used_weights) chosen_timepath <- sample(used_slots, 1, replace = TRUE, prob = used_weights) theprophecy[i] <- chosen_timepath[1] - 1 used_weights <- tweights[, chosen_timepath[1]] } } else if (is.numeric(tweights)) { if (length(tweights) != length(used_slots)) { if (length(tweights) == list_length) { used_weights <- tweights[used_slots] / sum(tweights[used_slots]) } else { stop(tw_error, call. = FALSE) } } else { used_weights <- tweights / sum(tweights) } theprophecy <- sample(used_slots, times, replace = TRUE, prob = used_weights) - 1 } else { stop(tw_error, call. = FALSE) } } else { used_weights <- rep(1, length(used_slots)) used_weights <- used_weights / sum(used_weights) theprophecy <- sample(used_slots, times, replace = TRUE, prob = used_weights) - 1 } starter <- if (!sparse_initial) { .ss3matrix(mats[[1]], sparsemethod) } else .ss3matrix_sp(mats[[1]]) theseventhmatrix <- if (!sparse_initial) { .proj3(starter, mats, theprophecy, 1, 0, 0, TRUE, sparsemethod) } else .proj3sp(starter, mats, theprophecy, 1, 0, 0) ssonly <- theseventhmatrix[((dim(mats[[1]])[1]) + 1):(2 *(dim(mats[[1]])[1])),] if (times > 2000) { Xlist <- ssonly[,(times - 999):(times)] } else if (times > 500) { Xlist <- ssonly[,(times-199):(times)] } else { Xlist <- ssonly } w_list <- apply(Xlist, 1, mean) } return(w_list) } #' Estimate Reproductive Value #' #' \code{repvalue3()} is a generic function that estimates returns the #' reproductive values of stages in a population projection matrix or a set of #' matrices. The specifics of estimation vary with the class of input object. #' This function is made to handle very large and sparse matrices supplied as #' \code{lefkoMat} objects or as individual matrices, and can be used with large #' historical matrices, IPMs, age x stage matrices, as well as ahistorical #' matrices. #' #' @name repvalue3 #' #' @param mats A lefkoMat object, a population projection matrix, or a list of #' population projection matrices for which the reproductive value vector is #' desired. #' @param ... Other parameters. #' #' @return The value returned depends on the class of the \code{mats} argument. #' See related functions for details. #' #' @seealso \code{\link{repvalue3.lefkoMat}()} #' @seealso \code{\link{repvalue3.matrix}()} #' @seealso \code{\link{repvalue3.dgCMatrix}()} #' @seealso \code{\link{repvalue3.list}()} #' #' @examples #' # Lathyrus deterministic example #' data(lathyrus) #' #' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0) #' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm") #' repvector <- c(0, 0, 0, 0, 0, 1, 0) #' obsvector <- c(0, 1, 1, 1, 1, 1, 0) #' matvector <- c(0, 0, 1, 1, 1, 1, 1) #' immvector <- c(1, 1, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5) #' #' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' immstatus = immvector, indataset = indataset, binhalfwidth = binvec, #' propstatus = propvector) #' #' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988, #' patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9, #' juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88", #' fecacol = "Intactseed88", deadacol = "Dead1988", #' nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea", #' censorcol = "Missing1988", censorkeep = NA, censor = TRUE) #' #' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"), #' stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"), #' stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"), #' eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"), #' eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"), #' eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"), #' givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA), #' multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA), #' type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1), #' stageframe = lathframe, historical = TRUE) #' #' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all", #' stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3, #' yearcol = "year2", indivcol = "individ") #' #' ehrlen3mean <- lmean(ehrlen3) #' repvalue3(ehrlen3mean) #' #' # Cypripedium stochastic example #' data(cypdata) #' #' sizevector <- c(0, 0, 0, 0, 0, 0, 1, 2.5, 4.5, 8, 17.5) #' stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "XSm", "Sm", "Md", "Lg", #' "XLg") #' repvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1) #' obsvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1) #' matvector <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1) #' immvector <- c(0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 0, 0, 0, 0, 0.5, 0.5, 1, 1, 2.5, 7) #' #' cypframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' propstatus = propvector, immstatus = immvector, indataset = indataset, #' binhalfwidth = binvec) #' #' cypraw_v1 <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004, #' patchidcol = "patch", individcol = "plantid", blocksize = 4, #' sizeacol = "Inf2.04", sizebcol = "Inf.04", sizeccol = "Veg.04", #' repstracol = "Inf.04", repstrbcol = "Inf2.04", fecacol = "Pod.04", #' stageassign = cypframe_raw, stagesize = "sizeadded", NAas0 = TRUE, #' NRasRep = TRUE) #' #' # Here we use supplemental() to provide overwrite and reproductive info #' cypsupp2r <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "D", #' "XSm", "Sm", "SD", "P1"), #' stage2 = c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "rep", #' "rep"), #' eststage3 = c(NA, NA, NA, NA, NA, "D", "XSm", "Sm", NA, NA), #' eststage2 = c(NA, NA, NA, NA, NA, "XSm", "XSm", "XSm", NA, NA), #' givenrate = c(0.10, 0.20, 0.20, 0.20, 0.25, NA, NA, NA, NA, NA), #' multiplier = c(NA, NA, NA, NA, NA, NA, NA, NA, 0.5, 0.5), #' type =c(1, 1, 1, 1, 1, 1, 1, 1, 3, 3), #' stageframe = cypframe_raw, historical = FALSE) #' #' 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") #' #' repvalue3(cypmatrix2r, stochastic = TRUE) #' #' @export repvalue3 <- function(mats, ...) UseMethod("repvalue3") #' Estimate Reproductive Value Vectors of Matrices in a lefkoMat Object #' #' \code{repvalue3.lefkoMat()} returns the reproductive values for stages in a #' set of population projection matrices provided as a \code{lefkoMat} object. #' This function can handle large and sparse matrices, and so can be used with #' large historical matrices, IPMs, age x stage matrices, as well as ahistorical #' matrices. #' #' @name repvalue3.lefkoMat #' #' @param mats An object of class \code{lefkoMat} object. #' @param stochastic A logical value indicating whether to use deterministic #' (\code{FALSE}) or stochastic (\code{TRUE}) analysis. Defaults to #' \code{FALSE}. #' @param times An integer variable indicating number of occasions to project if #' using stochastic analysis. Defaults to 10000. #' @param tweights An optional numeric vector or matrix denoting the #' probabilities of choosing each matrix in a stochastic projection. If a matrix #' is input, then a first-order Markovian environment is assumed, in which the #' probability of choosing a specific annual matrix depends on which annual #' matrix is currently chosen. If a vector is input, then the choice of annual #' matrix is assumed to be independent of the current matrix. Defaults to equal #' weighting among matrices. #' @param seed A number to use as a random number seed. #' @param force_sparse A text string indicating whether to use sparse matrix #' encoding (\code{"yes"}) when supplied with standard matrices. Defaults to #' \code{"auto"}, in which case sparse matrix encoding is used with square #' matrices with at least 50 rows and no more than 50\% of elements with values #' greater than zero. #' @param ... Other parameters. #' #' @return This function returns the asymptotic reproductive value vectors if #' deterministic analysis is chosen, and long-run mean reproductive value #' vectors if stochastic analysis is chosen. #' #' The output depends on whether the \code{lefkoMat} object used as input is #' ahistorical or historical, and whether the analysis is deterministic or #' stochastic. If deterministic and ahistorical, then a single data frame is #' output, which includes the number of the matrix within the \code{A} element #' of the input \code{lefkoMat} object, followed by the stage id (numeric and #' assigned through \code{\link{sf_create}()}), the stage name, and the #' estimated proportion of the reproductive value vector (\code{rep_value}). If #' stochastic and ahistorical, then a single data frame is output starting with #' the number of the population-patch (\code{matrix_set}), a string #' concatenating the names of the population and the patch (\code{poppatch}), #' the assigned stage id number (\code{stage_id}), and the stage name #' (\code{stage}), and the long-run mean reproductive value vector #' (\code{rep_value}). #' #' If a historical matrix is used as input, then two data frames are output #' into a list object. The \code{hist} element describes the historical #' stage-pair reproductive values, while the \code{ahist} element describes the #' stage reproductive values. If deterministic, then \code{hist} contains a data #' frame including the matrix number (\code{matrix}), the numeric stage #' designations for stages in occasions \emph{t} and \emph{t}-1, #' (\code{stage_id_2} and \code{stage_id_1}, respectively), followed by the #' respective stage names (\code{stage_2} and \code{stage_1}), and ending with #' the estimated reproductive values (\code{rep_value}). The associated #' \code{ahist} element is as before. If stochastic, then the \code{hist} #' element contains a single data frame with the number of the population-patch #' (\code{matrix_set}), a string concatenating the names of the population and #' the patch (\code{poppatch}), the assigned stage id numbers in times \emph{t} #' and \emph{t}-1 (\code{stage_id_2} and \code{stage_id_2}, respectively), and #' the associated stage names (\code{stage_2} and \code{stage_1}, respectively), #' and the long-run mean reproductive values (\code{rep_value}). The associated #' \code{ahist} element is as before in the ahistorical, stochastic case. #' #' In addition to the data frames noted above, stochastic analysis will result #' in the additional output of a list of matrices containing the actual #' projected reproductive value vectors across all projected occasions, in the #' order of population-patch combinations in the \code{lefkoMat} input. #' #' @section Notes: #' In stochastic analysis, the projected mean reproductive value vector is the #' arithmetic mean across the final projected 1000 occasions if the simulation #' is at least 2000 projected occasions long. If between 500 and 2000 projected #' occasions long, then only the final 200 are used, and if fewer than 500 #' occasions are used, then all are used. Note that because reproductive values #' in stochastic simulations can change greatly in the initial portion of the #' run, we encourage a minimum 2000 projected occasions per simulation, with #' 10000 preferred. #' #' Speed can sometimes be increased by shifting from automatic sparse matrix #' determination to forced dense or sparse matrix projection. This will most #' likely occur when matrices have several hundred rows and columns. Defaults #' work best when matrices are very small and dense, or very large and sparse. #' #' @seealso \code{\link{repvalue3}()} #' @seealso \code{\link{repvalue3.matrix}()} #' @seealso \code{\link{repvalue3.dgCMatrix}()} #' @seealso \code{\link{repvalue3.list}()} #' #' @examples #' data(cypdata) #' #' sizevector <- c(0, 0, 0, 0, 0, 0, 1, 2.5, 4.5, 8, 17.5) #' stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "XSm", "Sm", "Md", "Lg", #' "XLg") #' repvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1) #' obsvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1) #' matvector <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1) #' immvector <- c(0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 0, 0, 0, 0, 0.5, 0.5, 1, 1, 2.5, 7) #' #' cypframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' propstatus = propvector, immstatus = immvector, indataset = indataset, #' binhalfwidth = binvec) #' #' cypraw_v1 <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004, #' patchidcol = "patch", individcol = "plantid", blocksize = 4, #' sizeacol = "Inf2.04", sizebcol = "Inf.04", sizeccol = "Veg.04", #' repstracol = "Inf.04", repstrbcol = "Inf2.04", fecacol = "Pod.04", #' stageassign = cypframe_raw, stagesize = "sizeadded", NAas0 = TRUE, #' NRasRep = TRUE) #' #' # Here we use supplemental() to provide overwrite and reproductive info #' cypsupp2r <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "D", #' "XSm", "Sm", "SD", "P1"), #' stage2 = c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "rep", #' "rep"), #' eststage3 = c(NA, NA, NA, NA, NA, "D", "XSm", "Sm", NA, NA), #' eststage2 = c(NA, NA, NA, NA, NA, "XSm", "XSm", "XSm", NA, NA), #' givenrate = c(0.10, 0.20, 0.20, 0.20, 0.25, NA, NA, NA, NA, NA), #' multiplier = c(NA, NA, NA, NA, NA, NA, NA, NA, 0.5, 0.5), #' type =c(1, 1, 1, 1, 1, 1, 1, 1, 3, 3), #' stageframe = cypframe_raw, historical = FALSE) #' #' 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") #' #' repvalue3(cypmatrix2r, stochastic = TRUE) #' #' @export repvalue3.lefkoMat <- function(mats, stochastic = FALSE, times = 10000, tweights = NA, seed = NA, force_sparse = "auto", ...) { matrix_set <- poppatch <- NULL sparsemethod <- sparse_auto <- sparse_input <- FALSE if (is(mats$A[[1]], "dgCMatrix")) sparse_input <- TRUE

if (!sparse_input) {
if (is.logical(force_sparse)) {
if (force_sparse) {
sparsemethod <- TRUE
} else sparsemethod <- FALSE
} else if (is.element(tolower(force_sparse), c("y", "yes", "yea", "yeah", "t", "true", "ja", "tak"))) {
sparsemethod <- TRUE
} else if (is.element(tolower(force_sparse), c("n", "no", "non", "nah", "f", "false", "nein", "nie"))) {
sparsemethod <- FALSE
} else {
if (is.element(tolower(force_sparse), c("au", "aut", "auto"))) sparse_auto <- TRUE

elements_total <- length(mats$A[[1]]) dense_elements <- length(which(mats$A[[1]] != 0))

if ((dense_elements / elements_total) < 0.5 & elements_total > 2500) {
sparsemethod <- 1
} else sparsemethod <- 0
}
}

if (!stochastic) {
baldrick <- if (is.matrix(mats$A)) { almost_final <- if (!sparse_input) { .rv3matrix(mats$A, sparsemethod)
} else .rv3matrix_sp(mats$A) } else if (is.list(mats$A)) {
final <- unlist(lapply(mats$A, function(X) { almost_final <- if (!sparse_input) { .rv3matrix(X, sparsemethod) } else .rv3matrix_sp(X) return(almost_final/almost_final[which(almost_final == min(almost_final[which(almost_final > 0)])[1])[1]]) })) } else { stop("Input not recognized.") } } else { if (!is.na(seed)) { set.seed(seed) } mats$labels$poppatch <- paste(mats$labels$pop, mats$labels$patch) used_poppatches <- as.list(unique(mats$labels$poppatch)) # Full stage distribution and reproductive value vector series for all occasions, as list # Stage distributions in top half of matrix; reproductive value is at bottom princegeorge <- lapply(used_poppatches, function(X) { used_slots <- which(mats$labels$poppatch == X) if (length(used_slots) < 2) { warning("Only 1 annual matrix found for some pop-patch combinations. Stochastic analysis requires multiple annual matrices per pop-patch combination.", call. = FALSE) } if (!any(is.na(tweights))) { tw_error_1 <- "Option tweights must be NA, a numeric vector equal to the" tw_error_2 <- "number of years or matrices supplied, or a square matrix" tw_error_3 <- "with dimensions equal to the number of matrices per patch." tw_error <- paste(tw_error_1, tw_error_2, tw_error_3) if (is.matrix(tweights)) { theprophecy <- rep(0, times) used_weights <- tweights[, 1] for (i in c(1:times)) { used_weights <- used_weights / sum(used_weights) chosen_timepath <- sample(used_slots, 1, replace = TRUE, prob = used_weights) theprophecy[i] <- chosen_timepath[1] - 1 used_weights <- tweights[, chosen_timepath[1]] } } else if (is.numeric(tweights)) { if (length(tweights) != length(used_slots)) { if (length(tweights) == length(mats$A)) {
used_weights <- tweights[used_slots] / sum(tweights[used_slots])
} else {
stop(tw_error, call. = FALSE)
}
} else {
used_weights <- tweights / sum(tweights)
}

theprophecy <- sample(used_slots, times, replace = TRUE,
prob = used_weights) - 1
} else {
stop(tw_error, call. = FALSE)
}
} else {
used_weights <- rep(1, length(used_slots))
used_weights <- used_weights / sum(used_weights)

theprophecy <- sample(used_slots, times, replace = TRUE,
prob = used_weights) - 1
}

if (!sparse_input) {
starter <- .ss3matrix(mats$A[[used_slots[1]]], sparsemethod) } else starter <- .ss3matrix_sp(mats$A[[used_slots[1]]])

theseventhmatrix <- if (!sparse_input) {
.proj3(starter, mats$A, theprophecy, 1, 0, 0, sparse_auto, sparsemethod) } else .proj3sp(starter, mats$A, theprophecy, 1, 0, 0)

almostall <- theseventhmatrix[((dim(mats$A[[1]])[1]) + 1):(3 * (dim(mats$A[[1]])[1])),]

return(almostall)
})

# Create mean distributions
baldrick <- unlist(
lapply(princegeorge, function(X) {
if (times > 2000) {
usedX1 <- X[1:(dim(mats$A[[1]])[1]), (times - 998):(times+1)] usedX2 <- X[((dim(mats$A[[1]])[1]) + 1):(2*(dim(mats$A[[1]])[1])),1:1000] } else if (times > 500) { usedX1 <- X[1:(dim(mats$A[[1]])[1]), (times - 198):(times+1)]
usedX2 <- X[((dim(mats$A[[1]])[1]) + 1):(2*(dim(mats$A[[1]])[1])), 1:200]
} else {
usedX1 <- X[1:(dim(mats$A[[1]])[1]),] usedX2 <- X[((dim(mats$A[[1]])[1]) + 1):(2*(dim(mats$A[[1]])[1])),] } meanX1 <- apply(usedX1, 1, mean) meanX2 <- apply(usedX2, 1, mean) meanX2 <- zapsmall(meanX2) meanX2 <- meanX2 / meanX2[(which(meanX2 > 0)[1])] meanX <- c(meanX1, meanX2) return(meanX) }) ) princegeorge <- lapply(princegeorge, function(X) { return(X[((dim(mats$A[[1]])[1]) + 1):(2*(dim(mats$A[[1]])[1])),]) }) } if (is.list(mats$A)) {
if (!stochastic) {
multiplier <- length(mats$A) } else { multiplier <- length(used_poppatches) } } else multiplier <- 1 if (all(is.na(mats$hstages))) {
labels_orig <- mats$ahstages[,1:2] mat_dims <- dim(mats$A[[1]])[1]
if (dim(labels_orig)[1] == mat_dims) {
labels <- labels_orig
} else {
newmult <- mat_dims / dim(labels_orig)[1]
if (mat_dims %% dim(labels_orig)[1] != 0) {
stop("Matrices do not appear to be ahistorical, historical, or age x stage.
Cannot proceed. Please make sure that matrix dimensions match stage descriptions.",
call. = FALSE)
}

if (!all(is.na(mats$agestages))) { check_min <- min(mats$agestages$age) } age_bit <- c(apply(as.matrix(c((0 + check_min):(newmult + check_min - 1))), 1, rep, dim(labels_orig)[1])) core_labels <- cbind.data.frame(age_bit, do.call("rbind.data.frame", replicate(newmult, labels_orig, simplify = FALSE))) core_labels$agestage_id <- c(1:length(core_labels$stage_id)) core_labels$agestage <- apply(as.matrix(core_labels$agestage_id), 1, function(X) { paste(core_labels$age_bit[X], core_labels$stage[X]) }) labels <- core_labels names(labels) <- c("age", "stage_id", "stage", "agestage_id", "agestage") } if (!stochastic) { modlabels <- cbind.data.frame(as.matrix(rep(c(1:multiplier), each = dim(labels)[1])), do.call("rbind.data.frame", apply(as.matrix(c(1:multiplier)), 1, function(X) return(labels)))) output <- cbind.data.frame(modlabels, baldrick) if (is.element("age", names(labels))) { names(output) <- c("matrix", "age", "stage_id", "stage", "agestage_id", "agestage", "rep_value") } else { names(output) <- c("matrix", "stage_id", "stage", "rep_value") } } else { modlabels <- cbind.data.frame(as.matrix(rep(unlist(used_poppatches), each = dim(labels)[1])), do.call("rbind.data.frame", apply(as.matrix(c(1:multiplier)), 1, function(X) return(labels)))) unipoppatches <- unique(modlabels[,1]) modnums <- apply(as.matrix(modlabels[,1]), 1, function(X) { return(which(unipoppatches == X)) }) output <- cbind.data.frame(modnums, modlabels, baldrick[(mat_dims+1):(2*mat_dims)]) if (is.element("age", names(labels))) { names(output) <- c("matrix_set", "poppatch", "age", "stage_id", "stage", "agestage_id", "agestage", "rep_value") } else { names(output) <- c("matrix_set", "poppatch", "stage_id", "stage", "rep_value") } } rownames(output) <- c(1:dim(output)[1]) } else { # This section translates historical results to ahistorical and then cleans everything up if (!stochastic) { ss3 <- stablestage3.lefkoMat(mats, force_sparse = force_sparse) #stablestage3.lefkoMat rahist <- ss3$ahist
rhist <-ss3$hist rhist$ss3sum <- apply(as.matrix(c(1:dim(rhist)[1])), 1, function(X) {
rahist$ss_prop[intersect(which(rahist$stage_id == rhist$stage_id_2[X]), which(rahist$matrix == rhist$matrix[X]))] }) rhist$sscorr <- rhist$ss_prop / rhist$ss3sum
rhist$sscorr[which(is.na(rhist$sscorr))] <- 0
rhist$rep_value <- baldrick rhist$rv3raw <- rhist$sscorr * rhist$rep_value

outputh <- rhist[,c("matrix", "stage_id_2", "stage_id_1", "stage_2", "stage_1", "rep_value")]

ahlabels <- mats$ahstages[,c("stage_id", "stage")] rv2 <- Re(c(apply(as.matrix(c(1:multiplier)), 1, function(X) { rightset <- subset(rhist, matrix == X) apply(as.matrix(ahlabels[,1]), 1, function(Y) { sum(rightset$rv3raw[which(rightset$stage_id_2 == Y)]) }) }))) outputah <- cbind.data.frame(rep(c(1:multiplier), each = length(ahlabels[,1])), rep(ahlabels[,1], multiplier), rep(ahlabels[,2], multiplier), rv2) names(outputah) <- c("matrix", "stage_id", "stage", "rep_value_unc") outputah$rep_value <- apply(as.matrix(c(1:dim(outputah)[1])), 1, function(X) {
matsub <- subset(outputah, matrix == outputah$matrix[X]) entrystage <- min(which(abs(matsub$rep_value_unc) > 0))
return(outputah$rep_value_unc[X] / matsub$rep_value_unc[entrystage])
})

outputah <- outputah[,c("matrix", "stage_id", "stage", "rep_value")]
rownames(outputah) <- c(1:dim(outputah)[1])

} else {

ss_unlisted <- apply(as.matrix(c(1:multiplier)), 1, function(X) {
return(baldrick[((2 * (X - 1) * (dim(mats$A[[1]])[1])) + 1): (((2 * (X - 1)) + 1) * (dim(mats$A[[1]])[1]))])
})
rv_unlisted <- apply(as.matrix(c(1:multiplier)), 1, function(X) {
return(baldrick[((((2 * (X - 1)) + 1) * (dim(mats$A[[1]])[1])) + 1): (X * 2 * (dim(mats$A[[1]])[1]))])
})

ss_sums <- apply(as.matrix(c(1:multiplier)), 1, function(X) {
sum_vecs <- apply(as.matrix(mats$hstages$stage_id_2), 1, function(Y) {
sum(ss_unlisted[(which(mats$hstages$stage_id_2 == Y)) + ((X - 1) * dim(mats$hstages)[1])]) }) return(sum_vecs) }) ahistsize <- length(mats$ahstages$stage_id) histsize <- length(mats$hstages$stage_id_2) poppatchcol <- rep(unlist(used_poppatches), each = ahistsize) unipoppatches <- unique(poppatchcol) poppatchnum <- apply(as.matrix(poppatchcol), 1, function(X) { return(which(unipoppatches == X)) }) rahist <- cbind.data.frame(matrix_set = poppatchnum, poppatch = poppatchcol, stage_id = rep(mats$ahstages$stage_id, times = multiplier), stage = rep(mats$ahstages$stage, times = multiplier)) poppatchcol <- rep(unlist(used_poppatches), each = histsize) poppatchnum <- apply(as.matrix(poppatchcol), 1, function(X) { return(which(unipoppatches == X)) }) rhist <- cbind.data.frame(matrix_set = poppatchnum, poppatch = poppatchcol, stage_id_2 = rep(mats$hstages$stage_id_2, times = multiplier), stage_id_1 = rep(mats$hstages$stage_id_1, times = multiplier), stage_2 = rep(mats$hstages$stage_2, times = multiplier), stage_1 = rep(mats$hstages$stage_1, times = multiplier)) rhist$ss_prop <- as.vector(ss_unlisted)
rhist$ss3sum <- as.vector(ss_sums) rhist$sscorr <- rhist$ss_prop / rhist$ss3sum
rhist$sscorr[which(is.na(rhist$sscorr))] <- 0
rhist$rep_value <- as.vector(rv_unlisted) rhist$rv3raw <- rhist$sscorr * rhist$rep_value

outputh <- rhist[,c("matrix_set", "poppatch", "stage_id_2", "stage_id_1",
"stage_2", "stage_1", "rep_value")]

ahlabels <- mats$ahstages[,c("stage_id", "stage")] rv2 <- Re(c(apply(as.matrix(unlist(used_poppatches)), 1, function(X) { rightset <- subset(rhist, poppatch == X) apply(as.matrix(ahlabels[,1]), 1, function(Y) { sum(rightset$rv3raw[which(rightset$stage_id_2 == Y)]) }) }))) poppatchcol <- rep(unlist(used_poppatches), each = ahistsize) outputah <- cbind.data.frame(rep(c(1:multiplier), each = length(ahlabels[,1])), poppatchcol, rep(ahlabels[,1], multiplier), rep(ahlabels[,2], multiplier), rv2) names(outputah) <- c("matrix_set", "poppatch", "stage_id", "stage", "rep_value_unc") outputah$rep_value <- apply(as.matrix(c(1:dim(outputah)[1])), 1, function(X) {
matsub <- subset(outputah, poppatch == outputah$poppatch[X]) if (!all(matsub$rep_value_unc == 0)) {
entrystage <- min(which(abs(matsub$rep_value_unc) > 0)) return(outputah$rep_value_unc[X] / matsub$rep_value_unc[entrystage]) } else return(0) }) outputah <- outputah[,c("matrix_set", "poppatch", "stage_id", "stage", "rep_value")] rownames(outputah) <- c(1:dim(outputah)[1]) } if (!stochastic) { output <-list(hist = outputh, ahist = outputah) } else { output <-list(hist = outputh, ahist = outputah, projections = princegeorge) } } return(output) } #' Estimate Reproductive Value Vector for a Single Population Projection Matrix #' #' \code{repvalue3.matrix()} returns the reproductive values for stages in a #' population projection matrix. The function makes no assumptions about whether #' the matrix is ahistorical and simply provides standard reproductive values #' corresponding to each row, meaning that the overall reproductive values of #' basic life history stages in a historical matrix are not provided (the #' \code{\link{repvalue3.lefkoMat}()} function estimates these on the basis of #' stage description information provided in the \code{lefkoMat} object used as #' input in that function). This function can handle large and sparse matrices, #' and so can be used with large historical matrices, IPMs, age x stage #' matrices, as well as smaller ahistorical matrices. #' #' @name repvalue3.matrix #' #' @param mats A population projection matrix. #' @param force_sparse A text string indicating whether to use sparse matrix #' encoding (\code{"yes"}) when supplied with standard matrices. Defaults to #' \code{"auto"}, in which case sparse matrix encoding is used with square #' matrices with at least 50 rows and no more than 50\% of elements with values #' greater than zero. #' @param ... Other parameters. #' #' @return This function returns a vector data frame characterizing the #' reproductive values for stages of a population projection matrix. This is #' given as the left eigenvector associated with largest real part of the #' dominant eigenvalue, divided by the first non-zero element of the left #' eigenvector. #' #' @section Notes: #' Speed can sometimes be increased by shifting from automatic sparse matrix #' determination to forced dense or sparse matrix projection. This will most #' likely occur when matrices have between 30 and 300 rows and columns. #' Defaults work best when matrices are very small and dense, or very large and #' sparse. #' #' @seealso \code{\link{repvalue3}()} #' @seealso \code{\link{repvalue3.lefkoMat}()} #' @seealso \code{\link{repvalue3.dgCMatrix}()} #' @seealso \code{\link{repvalue3.list}()} #' #' @examples #' data(lathyrus) #' #' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0) #' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm") #' repvector <- c(0, 0, 0, 0, 0, 1, 0) #' obsvector <- c(0, 1, 1, 1, 1, 1, 0) #' matvector <- c(0, 0, 1, 1, 1, 1, 1) #' immvector <- c(1, 1, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5) #' #' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' immstatus = immvector, indataset = indataset, binhalfwidth = binvec, #' propstatus = propvector) #' #' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988, #' patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9, #' juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88", #' fecacol = "Intactseed88", deadacol = "Dead1988", #' nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea", #' censorcol = "Missing1988", censorkeep = NA, censor = TRUE) #' #' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"), #' stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"), #' stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"), #' eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"), #' eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"), #' eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"), #' givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA), #' multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA), #' type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1), #' stageframe = lathframe, historical = TRUE) #' #' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all", #' stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3, #' yearcol = "year2", indivcol = "individ") #' #' ehrlen3mean <- lmean(ehrlen3) #' repvalue3(ehrlen3mean$A[[1]])
#'
#' @export
repvalue3.matrix <- function(mats, force_sparse = "auto", ...)
{
sparsemethod <- 0

if (is.logical(force_sparse)) {
if (force_sparse) {
sparsemethod <- 1
} else sparsemethod <- 0
} else if (is.element(tolower(force_sparse), c("y", "yes", "yea", "yeah", "t", "true", "ja", "tak"))) {
sparsemethod <- 1
} else if (is.element(tolower(force_sparse), c("n", "no", "non", "nah", "f", "false", "nein", "nie"))) {
sparsemethod <- 0
} else {
elements_total <- length(mats)
dense_elements <- length(which(mats != 0))

if ((dense_elements / elements_total) < 0.5 & elements_total > 2500) {
sparsemethod <- 1
} else sparsemethod <- 0
}

v <- .rv3matrix(mats, sparsemethod)

return(v)
}

#' Estimate Reproductive Value Vector for a Single Population Projection Matrix
#'
#' \code{repvalue3.dgCMatrix()} returns the reproductive values for stages in a
#' sparse population projection matrix. The function makes no assumptions about
#' whether the matrix is ahistorical and simply provides standard reproductive
#' values corresponding to each row, meaning that the overall reproductive
#' values of basic life history stages in a historical matrix are not provided
#' (the \code{\link{repvalue3.lefkoMat}()} function estimates these on the basis
#' of stage description information provided in the \code{lefkoMat} object used
#' as input in that function).
#'
#' @name repvalue3.dgCMatrix
#'
#' @param mats A population projection matrix.
#' @param ... Other parameters.
#'
#' @return This function returns a vector data frame characterizing the
#' reproductive values for stages of a population projection matrix. This is
#' given as the left eigenvector associated with largest real part of the
#' dominant eigenvalue, divided by the first non-zero element of the left
#' eigenvector.
#'
#' @section Notes:
#' Speed can sometimes be increased by shifting from automatic sparse matrix
#' determination to forced dense or sparse matrix projection. This will most
#' likely occur when matrices have several hundred rows and columns. Defaults
#' work best when matrices are very small and dense, or very large and sparse.
#'
#'
#' @examples
#' data(lathyrus)
#'
#' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0)
#' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm")
#' repvector <- c(0, 0, 0, 0, 0, 1, 0)
#' obsvector <- c(0, 1, 1, 1, 1, 1, 0)
#' matvector <- c(0, 0, 1, 1, 1, 1, 1)
#' immvector <- c(1, 1, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5)
#'
#' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   immstatus = immvector, indataset = indataset, binhalfwidth = binvec,
#'   propstatus = propvector)
#'
#' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
#'   patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
#'   juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88",
#'   nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea",
#'   censorcol = "Missing1988", censorkeep = NA, censor = TRUE)
#'
#' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"),
#'   stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"),
#'   stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"),
#'   eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"),
#'   eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"),
#'   eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"),
#'   givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA),
#'   type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1),
#'   stageframe = lathframe, historical = TRUE)
#'
#' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all",
#'   stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3,
#'   yearcol = "year2", indivcol = "individ", sparse_output = TRUE)
#'
#' repvalue3(ehrlen3$A[[1]]) #' #' @export repvalue3.dgCMatrix <- function(mats, ...) { v <- .rv3matrix_sp(mats) return(v) } #' Estimate Reproductive Value Vector for a List of Projection Matrices #' #' \code{repvalue3.list()} returns the reproductive values for stages in #' population projection matrices arranged in a general list. The function makes #' no assumptions about whether the matrix is ahistorical and simply provides #' standard reproductive values corresponding to each row, meaning that the #' overall reproductive values of basic life history stages in a historical #' matrix are not provided (the \code{\link{repvalue3.lefkoMat}()} function #' estimates these on the basis of stage description information provided in the #' \code{lefkoMat} object used as input in that function). This function can #' handle large and sparse matrices, and so can be used with large historical #' matrices, IPMs, age x stage matrices, as well as smaller ahistorical #' matrices. #' #' @name repvalue3.list #' #' @param mats A list of population projection matrices, all in either class #' \code{matrix} or class \code{dgCMatrix}. #' @param stochastic A logical value indicating whether to use deterministic #' (\code{FALSE}) or stochastic (\code{TRUE}) analysis. Defaults to #' \code{FALSE}. #' @param times An integer variable indicating number of occasions to project if #' using stochastic analysis. Defaults to 10000. #' @param tweights An optional numeric vector or matrix denoting the #' probabilities of choosing each matrix in a stochastic projection. If a matrix #' is input, then a first-order Markovian environment is assumed, in which the #' probability of choosing a specific annual matrix depends on which annual #' matrix is currently chosen. If a vector is input, then the choice of annual #' matrix is assumed to be independent of the current matrix. Defaults to equal #' weighting among matrices. #' @param seed A number to use as a random number seed in stochastic projection. #' @param force_sparse A text string indicating whether to use sparse matrix #' encoding (\code{"yes"}) when supplied with standard matrices. Defaults to #' \code{"auto"}, in which case sparse matrix encoding is used with square #' matrices with at least 50 rows and no more than 50\% of elements with values #' greater than zero. #' @param ... Other parameters. #' #' @return This function returns a list of vector data frames characterizing the #' reproductive values for stages of each population projection matrix. This is #' given as the left eigenvector associated with largest real part of the #' dominant eigenvalue, divided by the first non-zero element of the left #' eigenvector. #' #' @section Notes: #' Speed can sometimes be increased by shifting from automatic sparse matrix #' determination to forced dense or sparse matrix projection. This will most #' likely occur when matrices have several hundred rows and columns. Defaults #' work best when matrices are very small and dense, or very large and sparse. #' #' @seealso \code{\link{repvalue3}()} #' @seealso \code{\link{repvalue3.lefkoMat}()} #' @seealso \code{\link{repvalue3.dgCMatrix}()} #' @seealso \code{\link{repvalue3.matrix}()} #' #' @examples #' data(lathyrus) #' #' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0) #' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm") #' repvector <- c(0, 0, 0, 0, 0, 1, 0) #' obsvector <- c(0, 1, 1, 1, 1, 1, 0) #' matvector <- c(0, 0, 1, 1, 1, 1, 1) #' immvector <- c(1, 1, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5) #' #' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' immstatus = immvector, indataset = indataset, binhalfwidth = binvec, #' propstatus = propvector) #' #' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988, #' patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9, #' juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88", #' fecacol = "Intactseed88", deadacol = "Dead1988", #' nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea", #' censorcol = "Missing1988", censorkeep = NA, censor = TRUE) #' #' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"), #' stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"), #' stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"), #' eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"), #' eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"), #' eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"), #' givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA), #' multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA), #' type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1), #' stageframe = lathframe, historical = TRUE) #' #' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all", #' stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3, #' yearcol = "year2", indivcol = "individ") #' #' repvalue3(ehrlen3$A)
#'
#' @export
repvalue3.list <- function(mats, stochastic = FALSE, times = 10000,
tweights = NA, seed = NA, force_sparse = "auto", ...) {

sparsemethod <- 0
sparse_initial <- FALSE
v_list <- Xlist <- NULL

if (is(mats[[1]], "dgCMatrix")) sparse_initial <- TRUE

if (!sparse_initial) {
if (is.logical(force_sparse)) {
if (force_sparse) {
sparsemethod <- 1
} else sparsemethod <- 0
} else if (is.element(tolower(force_sparse), c("y", "yes", "yea", "yeah", "t", "true", "ja", "tak"))) {
sparsemethod <- 1
} else if (is.element(tolower(force_sparse), c("n", "no", "non", "nah", "f", "false", "nein", "nie"))) {
sparsemethod <- 0
} else {
elements_total <- length(mats[[1]])
dense_elements <- length(which(mats[[1]] != 0))

if ((dense_elements / elements_total) < 0.5 & elements_total > 2500) {
sparsemethod <- 1
} else sparsemethod <- 0
}
}

list_length <- length(mats)

massive_panic <- apply(as.matrix(c(1:list_length)), 1, function(X) {
if (is.matrix(mats[[X]])) {
outchunk <- 1
} else if (is(mats[[X]], "dgCMatrix")) {
outchunk <- 2
} else outchunk <- 3
})

if (!all(massive_panic == 1) & !all(massive_panic == 2)) {
stop("Matrix list must be entirely in standard matrix format, or entirely in dgCMatrix format.",
call. = FALSE)
}

if (!stochastic) {
v_list <- lapply(mats, function(X) {
if (is.matrix(X)) {
output <- .rv3matrix(X, sparsemethod)
} else if (is(X, "dgCMatrix")) {
output <- .rv3matrix_sp(X)
} else {
stop("Unrecognized input in object mats.", FALSE)
}
})
} else {
if (!is.na(seed[1])) {
set.seed(seed[1])
}

used_slots <- c(1:list_length)
if (!any(is.na(tweights))) {
tw_error_1 <- "Option tweights must be NA, a numeric vector equal to the"
tw_error_2 <- "number of matrices supplied."
tw_error <- paste(tw_error_1, tw_error_2)

if (is.matrix(tweights)) {
theprophecy <- rep(0, times)
used_weights <- tweights[, 1]

for (i in c(1:times)) {
used_weights <- used_weights / sum(used_weights)
chosen_timepath <- sample(used_slots, 1, replace = TRUE, prob = used_weights)
theprophecy[i] <- chosen_timepath[1] - 1
used_weights <- tweights[, chosen_timepath[1]]
}
} else if (is.numeric(tweights)) {
if (length(tweights) != length(used_slots)) {
if (length(tweights) == list_length) {
used_weights <- tweights[used_slots] / sum(tweights[used_slots])
} else {
stop(tw_error, call. = FALSE)
}
} else {
used_weights <- tweights / sum(tweights)
}

theprophecy <- sample(used_slots, times, replace = TRUE,
prob = used_weights) - 1
} else {
stop(tw_error, call. = FALSE)
}
} else {
used_weights <- rep(1, length(used_slots))
used_weights <- used_weights / sum(used_weights)

theprophecy <- sample(used_slots, times, replace = TRUE,
prob = used_weights) - 1
}

starter <- if (!sparse_initial) {
.rv3matrix(mats[[1]], sparsemethod)
} else .rv3matrix_sp(mats[[1]])

theseventhmatrix <- if (!sparse_initial) {
.proj3(starter, mats, theprophecy, 1, 0, 0, TRUE, sparsemethod)
} else .proj3sp(starter, mats, theprophecy, 1, 0, 0)
almostall <- theseventhmatrix[((dim(mats[[1]])[1]) + 1):(3 * (dim(mats[[1]])[1])),]

if (times > 2000) {
usedX1 <- almostall[1:(dim(mats[[1]])[1]), (times - 998):(times+1)]
usedX2 <- almostall[((dim(mats[[1]])[1]) + 1):(2*(dim(mats[[1]])[1])),1:1000]
} else if (times > 500) {
usedX1 <- almostall[1:(dim(mats[[1]])[1]), (times - 198):(times+1)]
usedX2 <- almostall[((dim(mats[[1]])[1]) + 1):(2*(dim(mats[[1]])[1])), 1:200]
} else {
usedX1 <- almostall[1:(dim(mats[[1]])[1]),]
usedX2 <- almostall[((dim(mats[[1]])[1]) + 1):(2*(dim(mats[[1]])[1])),]
}
meanX1 <- apply(usedX1, 1, mean)
meanX2 <- apply(usedX2, 1, mean)
meanX2 <- zapsmall(meanX2)
meanX2 <- meanX2 / meanX2[(which(meanX2 > 0)[1])]
v_list <- c(meanX1, meanX2)
}

return(v_list)
}

#' Estimate Sensitivity of Population Growth Rate to Matrix Elements
#'
#' \code{sensitivity3()} is a generic function that returns the sensitivity of
#' the population growth rate to the elements of the matrices in a matrix
#' population model. Currently, this function estimates both deterministic and
#' stochastic sensitivities, where the growth rate is \eqn{\lambda} in the
#' former case and the log of the stochastic \eqn{\lambda} in the latter case.
#' This function is made to handle very large and sparse matrices supplied as
#' \code{lefkoMat} objects, as lists of matrices, and as individual matrices.
#'
#' @name sensitivity3
#'
#' @param mats A lefkoMat object, or population projection matrix, for which
#' the stable stage distribution is desired.
#' @param ... Other parameters
#'
#' @return The value returned depends on the class of the \code{mats} argument.
#'
#'
#' @examples
#' data(lathyrus)
#'
#' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0)
#' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm")
#' repvector <- c(0, 0, 0, 0, 0, 1, 0)
#' obsvector <- c(0, 1, 1, 1, 1, 1, 0)
#' matvector <- c(0, 0, 1, 1, 1, 1, 1)
#' immvector <- c(1, 1, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5)
#'
#' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   immstatus = immvector, indataset = indataset, binhalfwidth = binvec,
#'   propstatus = propvector)
#'
#' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
#'   patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
#'   juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88",
#'   nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea",
#'   censorcol = "Missing1988", censorkeep = NA, censor = TRUE)
#'
#' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"),
#'   stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"),
#'   stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"),
#'   eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"),
#'   eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"),
#'   eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"),
#'   givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA),
#'   type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1),
#'   stageframe = lathframe, historical = TRUE)
#'
#' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all",
#'   stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3,
#'   yearcol = "year2", indivcol = "individ")
#'
#' sensitivity3(ehrlen3)
#'
#' @export
sensitivity3 <- function(mats, ...) UseMethod("sensitivity3")

#' Estimate Sensitivity of Population Growth Rate of a lefkoMat Object
#'
#' \code{sensitivity3.lefkoMat()} returns the sensitivities of population growth
#' rate to elements of all \code{$A} matrices in an object of class #' \code{lefkoMat}. If deterministic, then \eqn{\lambda} is taken as the #' population growth rate. If stochastic, then the log of stochastic #' \eqn{\lambda}, or the log stochastic growth rate, is taken as the population #' growth rate. This function can handle large and sparse matrices, and so can #' be used with large historical matrices, IPMs, age x stage matrices, as well #' as smaller ahistorical matrices. #' #' @name sensitivity3.lefkoMat #' #' @param mats An object of class \code{lefkoMat}. #' @param stochastic A logical value determining whether to conduct a #' deterministic (FALSE) or stochastic (TRUE) sensitivity analysis. Defaults to #' FALSE. #' @param times The number of occasions to project forward in stochastic #' simulation. Defaults to \code{10000}. #' @param tweights An optional numeric vector or matrix denoting the #' probabilities of choosing each matrix in a stochastic projection. If a matrix #' is input, then a first-order Markovian environment is assumed, in which the #' probability of choosing a specific annual matrix depends on which annual #' matrix is currently chosen. If a vector is input, then the choice of annual #' matrix is assumed to be independent of the current matrix. Defaults to equal #' weighting among matrices. #' @param seed A number to use as a random number seed in stochastic projection. #' @param sparse A text string indicating whether to use sparse matrix encoding #' (\code{"yes"}) or dense matrix encoding (\code{"no"}). Defaults to #' \code{"auto"}, in which case sparse matrix encoding is used with square #' matrices with at least 50 rows and no more than 50\% of elements with values #' greater than zero. #' @param append_mats A logical value indicating whether to include the original #' A, U, and F matrices in the output \code{lefkoSens} object. #' @param ... Other parameters. #' #' @return This function returns an object of class \code{lefkoSens}, which is a #' list of 8 elements. The first, \code{h_sensmats}, is a list of historical #' sensitivity matrices (\code{NULL} if an ahMPM is used as input). The second, #' \code{ah_elasmats}, is a list of either ahistorical sensitivity matrices if #' an ahMPM is used as input, or, if an hMPM is used as input, then the result #' is a list of ahistorical matrices based on the equivalent historical #' dependencies assumed in the input historical matrices. The third element, #' \code{hstages}, is a data frame showing historical stage pairs (\code{NULL} #' if an ahMPM used as input). The fourth element, \code{agestages}, show the #' order of age-stage combinations, if age-by-stage MPMs have been supplied. The #' fifth element, \code{ahstages}, is a data frame showing the order of #' ahistorical stages. The last 3 elements are the A, U, and F portions of the #' input. #' #' @section Notes: #' All sensitivity matrix outputs from this function are in standard matrix #' format. #' #' Deterministic sensitivities are estimated as eqn. 9.14 in Caswell (2001, #' Matrix Population Models). Stochastic sensitivities are estimated as eqn. #' 14.97 in Caswell (2001). Note that stochastic sensitivities are of the log of #' the stochastic \eqn{\lambda}. #' #' Speed can sometimes be increased by shifting from automatic sparse matrix #' determination to forced dense or sparse matrix projection. This will most #' likely occur when matrices have between 30 and 300 rows and columns. #' Defaults work best when matrices are very small and dense, or very large and #' sparse. #' #' The \code{time_weights} and \code{steps} arguments are now deprecated. #' Instead, please use the \code{tweights} and \code{times} arguments. #' #' @seealso \code{\link{sensitivity3}()} #' @seealso \code{\link{sensitivity3.matrix}()} #' @seealso \code{\link{sensitivity3.dgCMatrix}()} #' @seealso \code{\link{sensitivity3.list}()} #' #' @examples #' data(lathyrus) #' #' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0) #' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm") #' repvector <- c(0, 0, 0, 0, 0, 1, 0) #' obsvector <- c(0, 1, 1, 1, 1, 1, 0) #' matvector <- c(0, 0, 1, 1, 1, 1, 1) #' immvector <- c(1, 1, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5) #' #' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' immstatus = immvector, indataset = indataset, binhalfwidth = binvec, #' propstatus = propvector) #' #' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988, #' patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9, #' juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88", #' fecacol = "Intactseed88", deadacol = "Dead1988", #' nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea", #' censorcol = "Missing1988", censorkeep = NA, censor = TRUE) #' #' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"), #' stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"), #' stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"), #' eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"), #' eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"), #' eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"), #' givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA), #' multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA), #' type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1), #' stageframe = lathframe, historical = TRUE) #' #' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all", #' stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3, #' yearcol = "year2", indivcol = "individ") #' #' sensitivity3(ehrlen3, stochastic = TRUE) #' #' @export sensitivity3.lefkoMat <- function(mats, stochastic = FALSE, times = 10000, tweights = NA, seed = NA, sparse = "auto", append_mats = FALSE, ...) { sparsemethod <- 0 sparse_input <- FALSE possible_args <- c("mats", "stochastic", "times", "tweights", "historical", "seed", "force_sparse", "append_mats", "time_weights", "steps", "sparse") further_args <- list(...) further_args_names <- names(further_args) if (length(setdiff(further_args_names, possible_args)) > 0) { stop("Some arguments not recognized.", call. = FALSE) } if ("time_weights" %in% further_args_names) { stop("Argument time_weights is deprecated. Please use argument tweights instead", call. = FALSE) } if ("steps" %in% further_args_names) { stop("Argument steps is deprecated. Please use argument times instead", call. = FALSE) } if (is(mats$A[[1]], "dgCMatrix")) sparse_input <- TRUE

if (is.logical(sparse)) {
if (sparse) {
sparsemethod <- 1
} else sparsemethod <- 0
} else if (is.element(tolower(sparse), c("y", "yes", "yea", "yeah", "t", "true", "ja", "tak"))) {
sparsemethod <- 1
} else if (is.element(tolower(sparse), c("n", "no", "non", "nah", "f", "false", "nein", "nie"))) {
sparsemethod <- 0
} else {
elements_total <- length(mats$A[[1]]) dense_elements <- if (!sparse_input ) { length(which(mats$A[[1]] != 0))
} else {
sum(diff(mats$A[[1]]@p)) } if ((dense_elements / elements_total) < 0.5 & elements_total > 2500) { sparsemethod <- 1 } else sparsemethod <- 0 } if (stochastic & length(mats$A) < 2) {
stop("Stochastic sensitivity estimation cannot be completed with fewer than 2 annual matrices.",
call. = FALSE)
}

if (!is.element("agestages", names(mats))) {
mats$agestages <- NA } if (!stochastic) { # Deterministic sensitivity analysis message("Running deterministic analysis...") baldrick <- if (all(is.na(mats$hstages))) {
if (!sparse_input) {
lapply(mats$A, .sens3matrix, sparsemethod) } else { lapply(mats$A, .sens3matrix_spinp)
}
} else {
if (!sparse_input) {
lapply(mats$A, .sens3hlefko, mats$ahstages, mats$hstages) } else { lapply(mats$A, .sens3hlefko_sp, mats$ahstages, mats$hstages)
}
}

hlabels <- mats$hstages ahlabels <- mats$ahstages
agelabels <- mats$agestages new_labels <- mats$labels

if (all(is.na(mats$hstages))) { output <- list(h_sensmats = NULL, ah_sensmats = baldrick, hstages = NULL, ahstages = ahlabels, agestages = agelabels, labels = new_labels) } else { he_list <- lapply(baldrick, function(X) {X$h_smat})
ahe_list <- lapply(baldrick, function(X) {X$ah_smat}) output <- list(h_sensmats = he_list, ah_sensmats = ahe_list, hstages = hlabels, ahstages = ahlabels, agestages = agelabels, labels = new_labels) } } else { # Stochastic sensitivity analysis if (!is.na(seed[1])) { set.seed(seed[1]) } if(!any(is.na(tweights))) { message("Running stochastic analysis with tweights option...") returned_cubes <- .stoch_senselas(mats, times = times, historical = FALSE, style = 1, sparsemethod, lefkoProj = TRUE, tweights = tweights) } else { message("Running stochastic analysis...") returned_cubes <- .stoch_senselas(mats, times = times, historical = FALSE, style = 1, sparsemethod, lefkoProj = TRUE) } old_labels <- mats$labels
old_labels$concat <- paste(old_labels$pop, old_labels$patch) labs_unique <- unique(old_labels$concat)
lab_indices <- apply(as.matrix(labs_unique), 1, function(X) {
found_guys <- which(old_labels$concat == X) min_guy <- min(found_guys) return(min_guy) }) new_labels <- mats$labels[lab_indices, c(1, 2)]

if (!all(is.na(mats$hstages))) { output <- list(h_sensmats = returned_cubes[[1]], ah_sensmats = returned_cubes[[2]], hstages = mats$hstages, agestages = mats$agestages, ahstages = mats$ahstages, labels = new_labels)
} else {
output <- list(h_sensmats = NULL, ah_sensmats = returned_cubes[[1]],
hstages = mats$hstages, agestages = mats$agestages,
ahstages = mats$ahstages, labels = new_labels) } } if (append_mats) { output$A <- mats$A output$U <- mats$U output$F <- mats$F } class(output) <- "lefkoSens" return(output) } #' Estimate Sensitivity of Population Growth Rate of a Single Matrix #' #' \code{sensitivity3.matrix()} returns the sensitivities of \eqn{\lambda} to #' elements of a single matrix. Because this handles only one matrix, the #' sensitivities are inherently deterministic and based on the dominant eigen #' value as the best metric of the population growth rate. This function can #' handle large and sparse matrices, and so can be used with large historical #' matrices, IPMs, age x stage matrices, as well as smaller ahistorical #' matrices. #' #' @name sensitivity3.matrix #' #' @param mats An object of class \code{matrix}. #' @param sparse A text string indicating whether to use sparse matrix encoding #' (\code{"yes"}) or dense matrix encoding (\code{"no"}). Defaults to #' \code{"auto"}, in which case sparse matrix encoding is used with square #' matrices with at least 50 rows and no more than 50\% of elements with values #' greater than zero. #' @param ... Other parameters. #' #' @return This function returns a single deterministic sensitivity matrix. #' #' @section Notes: #' All sensitivity matrix outputs from this function are in standard matrix #' format. #' #' Speed can sometimes be increased by shifting from automatic sparse matrix #' determination to forced dense or sparse matrix projection. This will most #' likely occur when matrices have between 30 and 300 rows and columns. #' Defaults work best when matrices are very small and dense, or very large and #' sparse. #' #' @seealso \code{\link{sensitivity3}()} #' @seealso \code{\link{sensitivity3.lefkoMat}()} #' @seealso \code{\link{sensitivity3.dgCMatrix}()} #' @seealso \code{\link{sensitivity3.list}()} #' #' @examples #' data(lathyrus) #' #' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0) #' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm") #' repvector <- c(0, 0, 0, 0, 0, 1, 0) #' obsvector <- c(0, 1, 1, 1, 1, 1, 0) #' matvector <- c(0, 0, 1, 1, 1, 1, 1) #' immvector <- c(1, 1, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5) #' #' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' immstatus = immvector, indataset = indataset, binhalfwidth = binvec, #' propstatus = propvector) #' #' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988, #' patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9, #' juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88", #' fecacol = "Intactseed88", deadacol = "Dead1988", #' nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea", #' censorcol = "Missing1988", censorkeep = NA, censor = TRUE) #' #' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"), #' stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"), #' stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"), #' eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"), #' eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"), #' eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"), #' givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA), #' multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA), #' type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1), #' stageframe = lathframe, historical = TRUE) #' #' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all", #' stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3, #' yearcol = "year2", indivcol = "individ") #' #' ehrlen3mean <- lmean(ehrlen3) #' sensitivity3(ehrlen3mean$A[[1]])
#'
#' @export
sensitivity3.matrix <- function(mats, sparse = "auto", ...)
{
sparsemethod <- 0

possible_args <- c("mats", "sparse")
further_args <- list(...)
further_args_names <- names(further_args)
if (length(setdiff(further_args_names, possible_args)) > 0) {
stop("Some arguments not recognized.", call. = FALSE)
}

if (is.logical(sparse)) {
if (sparse) {
sparsemethod <- 1
} else sparsemethod <- 0
} else if (is.element(tolower(sparse), c("y", "yes", "yea", "yeah", "t", "true", "ja", "tak"))) {
sparsemethod <- 1
} else if (is.element(tolower(sparse), c("n", "no", "non", "nah", "f", "false", "nein", "nie"))) {
sparsemethod <- 0
} else {
elements_total <- length(mats)
dense_elements <- length(which(mats != 0))

if ((dense_elements / elements_total) < 0.5 & elements_total > 2500) {
sparsemethod <- 1
} else sparsemethod <- 0
}

wcorr <- .sens3matrix(mats, FALSE)

if (sparsemethod == 1) wcorr <- as(wcorr, "dgCMatrix")

return(wcorr)
}

#' Estimate Sensitivity of Population Growth Rate of a Single Matrix
#'
#' \code{sensitivity3.dgCMatrix()} returns the sensitivities of \eqn{\lambda} to
#' elements of a single, sparse matrix. Because this handles only one matrix,
#' sensitivities are inherently deterministic and based on the dominant eigen
#' value as the best metric of the population growth rate.
#'
#' @name sensitivity3.dgCMatrix
#'
#' @param mats An object of class \code{dgCMatrix}.
#' @param sparse A text string indicating whether to use sparse matrix encoding
#' (\code{"yes"}) or dense matrix encoding (\code{"no"}). Defaults to
#' \code{"auto"}, in which case sparse matrix encoding is used with square
#' matrices with at least 50 rows and no more than 50\% of elements with values
#' greater than zero.
#' @param ... Other parameters.
#'
#' @return This function returns a single deterministic sensitivity matrix.
#'
#' @section Notes:
#' All sensitivity matrix outputs from this function are in standard matrix
#' format.
#'
#'
#' @examples
#' data(lathyrus)
#'
#' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0)
#' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm")
#' repvector <- c(0, 0, 0, 0, 0, 1, 0)
#' obsvector <- c(0, 1, 1, 1, 1, 1, 0)
#' matvector <- c(0, 0, 1, 1, 1, 1, 1)
#' immvector <- c(1, 1, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5)
#'
#' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   immstatus = immvector, indataset = indataset, binhalfwidth = binvec,
#'   propstatus = propvector)
#'
#' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
#'   patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
#'   juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88",
#'   nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea",
#'   censorcol = "Missing1988", censorkeep = NA, censor = TRUE)
#'
#' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"),
#'   stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"),
#'   stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"),
#'   eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"),
#'   eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"),
#'   eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"),
#'   givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA),
#'   type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1),
#'   stageframe = lathframe, historical = TRUE)
#'
#' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all",
#'   stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3,
#'   yearcol = "year2", indivcol = "individ", sparse_output = TRUE)
#'
#' sensitivity3(ehrlen3$A[[1]]) #' #' @export sensitivity3.dgCMatrix <- function(mats, sparse = "auto", ...) { possible_args <- c("mats", "sparse") further_args <- list(...) further_args_names <- names(further_args) if (length(setdiff(further_args_names, possible_args)) > 0) { stop("Some arguments not recognized.", call. = FALSE) } if (is.logical(sparse)) { if (sparse) { sparsemethod <- 1 } else sparsemethod <- 0 } else if (is.element(tolower(sparse), c("y", "yes", "yea", "yeah", "t", "true", "ja", "tak"))) { sparsemethod <- 1 } else if (is.element(tolower(sparse), c("n", "no", "non", "nah", "f", "false", "nein", "nie"))) { sparsemethod <- 0 } else { elements_total <- length(mats) dense_elements <- sum(diff(mats@p)) if ((dense_elements / elements_total) < 0.5 & elements_total > 2500) { sparsemethod <- 1 } else sparsemethod <- 0 } wcorr <- .sens3matrix_spinp(mats) if (sparsemethod == 1) { wcorr <- as(wcorr, "dgCMatrix") } return(wcorr) } #' Estimate Sensitivity of Population Growth Rate of a List of Matrices #' #' \code{sensitivity3.list()} returns the sensitivities of population growth #' rate to elements of matrices supplied in a list. The sensitivity analysis can #' be deterministic or stochastic, but if the latter then at least two A #' matrices must be included in the list. This function can handle large and #' sparse matrices, and so can be used with large historical matrices, IPMs, #' age x stage matrices, as well as smaller ahistorical matrices. #' #' @name sensitivity3.list #' #' @param mats An object of class \code{matrix}. #' @param stochastic A logical value determining whether to conduct a #' deterministic (FALSE) or stochastic (TRUE) sensitivity analysis. Defaults to #' FALSE. #' @param times The number of occasions to project forward in stochastic #' simulation. Defaults to 10,000. #' @param tweights An optional numeric vector or matrix denoting the #' probabilities of choosing each matrix in a stochastic projection. If a matrix #' is input, then a first-order Markovian environment is assumed, in which the #' probability of choosing a specific annual matrix depends on which annual #' matrix is currently chosen. If a vector is input, then the choice of annual #' matrix is assumed to be independent of the current matrix. Defaults to equal #' weighting among matrices. #' @param historical A logical value indicating whether matrices are historical. #' Defaults to \code{FALSE}. #' @param seed A number to use as a random number seed in stochastic projection. #' @param sparse A text string indicating whether to use sparse matrix encoding #' (\code{"yes"}) or dense matrix encoding (\code{"no"}). Defaults to #' \code{"auto"}, in which case sparse matrix encoding is used with square #' matrices with at least 50 rows and no more than 50\% of elements with values #' greater than zero. #' @param append_mats A logical value indicating whether to include the original #' matrices input as object \code{mats} in the output \code{lefkoSense} object. #' Defaults to FALSE. #' @param ... Other parameters. #' #' @return This function returns an object of class \code{lefkoSens}, which is a #' list of 8 elements. The first, \code{h_sensmats}, is a list of historical #' sensitivity matrices (\code{NULL} if an ahMPM is used as input). The second, #' \code{ah_elasmats}, is a list of ahistorical sensitivity matrices if an ahMPM #' is used as input (\code{NULL} if an hMPM is used as input). The third #' element, \code{hstages}, the fourth element, \code{agestages}, and the fifth #' element, \code{ahstages}, are \code{NULL}. The last 3 elements include the #' original A matrices supplied (as the \code{A} element), followed by #' \code{NULL}s for the U and F elements. #' #' @section Notes: #' All sensitivity matrix outputs from this function are in standard matrix #' format. #' #' Deterministic sensitivities are estimated as eqn. 9.14 in Caswell (2001, #' Matrix Population Models). Stochastic sensitivities are estimated as eqn. #' 14.97 in Caswell (2001). Note that stochastic sensitivities are with regard #' to the log of the stochastic \eqn{\lambda}. #' #' Currently, this function does not estimate equivalent ahistorical stochastic #' sensitivities for input historical matrices, due to the lack of guidance #' input on the order of stages (guidance is provided within \code{lefkoMat} #' objects). #' #' Speed can sometimes be increased by shifting from automatic sparse matrix #' determination to forced dense or sparse matrix projection. This will most #' likely occur when matrices have between 30 and 300 rows and columns. #' Defaults work best when matrices are very small and dense, or very large and #' sparse. #' #' The \code{time_weights} and \code{steps} arguments are now deprecated. #' Instead, please use the \code{tweights} and \code{times} arguments. #' #' @seealso \code{\link{sensitivity3}()} #' @seealso \code{\link{sensitivity3.lefkoMat}()} #' @seealso \code{\link{sensitivity3.matrix}()} #' @seealso \code{\link{sensitivity3.dgCMatrix}()} #' #' @examples #' # Lathyrus example #' data(lathyrus) #' #' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0) #' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm") #' repvector <- c(0, 0, 0, 0, 0, 1, 0) #' obsvector <- c(0, 1, 1, 1, 1, 1, 0) #' matvector <- c(0, 0, 1, 1, 1, 1, 1) #' immvector <- c(1, 1, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5) #' #' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' immstatus = immvector, indataset = indataset, binhalfwidth = binvec, #' propstatus = propvector) #' #' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988, #' patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9, #' juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88", #' fecacol = "Intactseed88", deadacol = "Dead1988", #' nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea", #' censorcol = "Missing1988", censorkeep = NA, censor = TRUE) #' #' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"), #' stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"), #' stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"), #' eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"), #' eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"), #' eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"), #' givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA), #' multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA), #' type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1), #' stageframe = lathframe, historical = TRUE) #' #' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all", #' stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3, #' yearcol = "year2", indivcol = "individ") #' #' sensitivity3(ehrlen3$A)
#'
#' # Cypripedium example
#' data(cypdata)
#'
#' sizevector <- c(0, 0, 0, 0, 0, 0, 1, 2.5, 4.5, 8, 17.5)
#' stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "XSm", "Sm", "Md", "Lg",
#'   "XLg")
#' repvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
#' obsvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
#' matvector <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)
#' immvector <- c(0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 0, 0, 0, 0, 0.5, 0.5, 1, 1, 2.5, 7)
#'
#' cypframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   propstatus = propvector, immstatus = immvector, indataset = indataset,
#'   binhalfwidth = binvec)
#'
#' cypraw_v1 <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004,
#'   patchidcol = "patch", individcol = "plantid", blocksize = 4,
#'   sizeacol = "Inf2.04", sizebcol = "Inf.04", sizeccol = "Veg.04",
#'   repstracol = "Inf.04", repstrbcol = "Inf2.04", fecacol = "Pod.04",
#'   stageassign = cypframe_raw, stagesize = "sizeadded", NAas0 = TRUE,
#'   NRasRep = TRUE)
#'
#' cypsupp2r <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "D",
#'     "XSm", "Sm", "SD", "P1"),
#'   stage2 = c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "rep",
#'     "rep"),
#'   eststage3 = c(NA, NA, NA, NA, NA, "D", "XSm", "Sm", NA, NA),
#'   eststage2 = c(NA, NA, NA, NA, NA, "XSm", "XSm", "XSm", NA, NA),
#'   givenrate = c(0.10, 0.20, 0.20, 0.20, 0.25, NA, NA, NA, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, NA, NA, NA, NA, 0.5, 0.5),
#'   type =c(1, 1, 1, 1, 1, 1, 1, 1, 3, 3),
#'   stageframe = cypframe_raw, historical = FALSE)
#'
#' cypmatrix2r <- rlefko2(data = cypraw_v1, stageframe = cypframe_raw,
#'   year = "all", patch = "all", stages = c("stage3", "stage2", "stage1"),
#'   yearcol = "year2", patchcol = "patchid", indivcol = "individ")
#'
#' sensitivity3(cypmatrix2r$A) #' #' @export sensitivity3.list <- function(mats, stochastic = FALSE, times = 10000, tweights = NA, historical = FALSE, seed = NA, sparse = "auto", append_mats = FALSE, ...) { sparsemethod <- 0 sparse_input <- FALSE possible_args <- c("mats", "stochastic", "times", "tweights", "historical", "seed", "force_sparse", "append_mats", "time_weights", "steps", "sparse") further_args <- list(...) further_args_names <- names(further_args) if (length(setdiff(further_args_names, possible_args)) > 0) { stop("Some arguments not recognized.", call. = FALSE) } if ("time_weights" %in% further_args_names) { stop("Argument time_weights is deprecated. Please use argument tweights instead", call. = FALSE) } if ("steps" %in% further_args_names) { stop("Argument steps is deprecated. Please use argument times instead", call. = FALSE) } if (is(mats[[1]], "dgCMatrix")) sparse_input <- TRUE if (is.logical(sparse)) { if (sparse) { sparsemethod <- 1 } else sparsemethod <- 0 } else if (is.element(tolower(sparse), c("y", "yes", "yea", "yeah", "t", "true", "ja", "tak"))) { sparsemethod <- 1 } else if (is.element(tolower(sparse), c("n", "no", "non", "nah", "f", "false", "nein", "nie"))) { sparsemethod <- 0 } else { elements_total <- length(mats[[1]]) dense_elements <- if (!sparse_input ) { length(which(mats[[1]] != 0)) } else { sum(diff(mats[[1]]@p)) } if ((dense_elements / elements_total) < 0.5 & elements_total > 2500) { sparsemethod <- 1 } else sparsemethod <- 0 } if(length(setdiff(unlist(lapply(mats, class)), c("matrix", "array", "dgCMatrix"))) > 0) { stop("Input list must be composed only of numeric matrices.", call. = FALSE) } if (stochastic & length(mats) < 2) { stop("Stochastic sensitivity estimation cannot be completed with fewer than 2 annual matrices.", call. = FALSE) } if (!stochastic) { # Deterministic sensitivity analysis message("Running deterministic analysis...") baldrick <- if (!sparse_input) { lapply(mats, .sens3matrix, sparsemethod) } else { lapply(mats, .sens3matrix_spinp) } if (historical) { output <- list(h_sensmats = baldrick, ah_sensmats = NULL, hstages = NULL, agestages = NULL, ahstages = NULL, U = NULL, F = NULL) } else { output <- list(h_sensmats = NULL, ah_sensmats = baldrick, hstages = NULL, agestages = NULL, ahstages = NULL, U = NULL, F = NULL) } } else { # Stochastic sensitivity analysis if (!is.na(seed[1])) { set.seed(seed[1]) } if(!any(is.na(tweights))) { message("Running stochastic analysis with tweights option...") returned_cube <- .stoch_senselas(mats, times = times, historical = historical, style = 1, sparsemethod, lefkoProj = FALSE, tweights = tweights)[[1]] } else { message("Running stochastic analysis...") returned_cube <- .stoch_senselas(mats, times = times, historical = historical, style = 1, sparsemethod, lefkoProj = FALSE)[[1]] } if (historical) { output <- list(h_sensmats = list(returned_cube[[1]]), ah_sensmats = NULL, hstages = NULL, agestages = NULL, ahstages = NULL, U = NULL, F = NULL) } else { output <- list(h_sensmats = NULL, ah_sensmats = list(returned_cube[[1]]), hstages = NULL, agestages = NULL, ahstages = NULL, U = NULL, F = NULL) } } if (append_mats) { output$A <- mats
}

class(output) <- "lefkoSens"

return(output)
}

#' Estimate Elasticity of Population Growth Rate to Matrix Elements
#'
#' \code{elasticity3()} is a generic function that returns the elasticity of
#' the population growth rate to the elements of the matrices in a matrix
#' population model. Currently, this function estimates both deterministic and
#' stochastic elasticities, where the growth rate is \eqn{\lambda} in the former
#' case and the log of the stochastic \eqn{\lambda} in the latter case. This
#' function is made to handle very large and sparse matrices supplied as
#' \code{lefkoMat} objects, as lists of matrices, and as individual matrices.
#'
#' @name elasticity3
#'
#' @param mats A lefkoMat object, a population projection matrix, or a list of
#' population projection matrices for which the stable stage distribution is
#' desired.
#' @param ... Other parameters.
#'
#' @return The value returned depends on the class of the \code{mats} argument.
#'
#'
#' @examples
#' # Lathyrus example
#' data(lathyrus)
#'
#' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0)
#' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm")
#' repvector <- c(0, 0, 0, 0, 0, 1, 0)
#' obsvector <- c(0, 1, 1, 1, 1, 1, 0)
#' matvector <- c(0, 0, 1, 1, 1, 1, 1)
#' immvector <- c(1, 1, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5)
#'
#' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   immstatus = immvector, indataset = indataset, binhalfwidth = binvec,
#'   propstatus = propvector)
#'
#' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
#'   patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
#'   juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88",
#'   nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea",
#'   censorcol = "Missing1988", censorkeep = NA, censor = TRUE)
#'
#' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"),
#'   stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"),
#'   stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"),
#'   eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"),
#'   eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"),
#'   eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"),
#'   givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA),
#'   type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1),
#'   stageframe = lathframe, historical = TRUE)
#'
#' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all",
#'   stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3,
#'   yearcol = "year2", indivcol = "individ")
#'
#' ehrlen3mean <- lmean(ehrlen3)
#' elasticity3(ehrlen3mean)
#'
#' # Cypripedium example
#' data(cypdata)
#'
#' sizevector <- c(0, 0, 0, 0, 0, 0, 1, 2.5, 4.5, 8, 17.5)
#' stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "XSm", "Sm", "Md", "Lg",
#'   "XLg")
#' repvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
#' obsvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1)
#' matvector <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)
#' immvector <- c(0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 0, 0, 0, 0, 0.5, 0.5, 1, 1, 2.5, 7)
#'
#' cypframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   propstatus = propvector, immstatus = immvector, indataset = indataset,
#'   binhalfwidth = binvec)
#'
#' cypraw_v1 <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004,
#'   patchidcol = "patch", individcol = "plantid", blocksize = 4,
#'   sizeacol = "Inf2.04", sizebcol = "Inf.04", sizeccol = "Veg.04",
#'   repstracol = "Inf.04", repstrbcol = "Inf2.04", fecacol = "Pod.04",
#'   stageassign = cypframe_raw, stagesize = "sizeadded", NAas0 = TRUE,
#'   NRasRep = TRUE)
#'
#' cypsupp2r <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "D",
#'     "XSm", "Sm", "SD", "P1"),
#'   stage2 = c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "rep",
#'     "rep"),
#'   eststage3 = c(NA, NA, NA, NA, NA, "D", "XSm", "Sm", NA, NA),
#'   eststage2 = c(NA, NA, NA, NA, NA, "XSm", "XSm", "XSm", NA, NA),
#'   givenrate = c(0.10, 0.20, 0.20, 0.20, 0.25, NA, NA, NA, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, NA, NA, NA, NA, 0.5, 0.5),
#'   type = c(1, 1, 1, 1, 1, 1, 1, 1, 3, 3),
#'   stageframe = cypframe_raw, historical = FALSE)
#'
#' cypmatrix2r <- rlefko2(data = cypraw_v1, stageframe = cypframe_raw,
#'   year = "all", patch = "all", stages = c("stage3", "stage2", "stage1"),
#'   yearcol = "year2", patchcol = "patchid", indivcol = "individ")
#'
#' elasticity3(cypmatrix2r)
#'
#' @export
elasticity3 <- function(mats, ...) UseMethod("elasticity3")

#' Estimate Elasticity of Population Growth Rate of a lefkoMat Object
#'
#' \code{elasticity3.lefkoMat()} returns the elasticities of population growth
#' rate to elements of all \code{$A} matrices in an object of class #' \code{lefkoMat}. If deterministic, then \eqn{\lambda} is taken as the #' population growth rate. If stochastic, then stochastic \eqn{\lambda}, or #' the stochastic growth rate, is taken as the population growth rate. This #' function can handle large and sparse matrices, and so can be used with large #' historical matrices, IPMs, age x stage matrices, as well as smaller #' ahistorical matrices. #' #' @name elasticity3.lefkoMat #' #' @param mats An object of class \code{lefkoMat}. #' @param stochastic A logical value determining whether to conduct a #' deterministic (FALSE) or stochastic (TRUE) elasticity analysis. Defaults to #' FALSE. #' @param times The number of occasions to project forward in stochastic #' simulation. Defaults to 10,000. #' @param tweights An optional numeric vector or matrix denoting the #' probabilities of choosing each matrix in a stochastic projection. If a matrix #' is input, then a first-order Markovian environment is assumed, in which the #' probability of choosing a specific annual matrix depends on which annual #' matrix is currently chosen. If a vector is input, then the choice of annual #' matrix is assumed to be independent of the current matrix. Defaults to equal #' weighting among matrices. #' @param seed A number to use as a random number seed in stochastic projection. #' @param sparse A text string indicating whether to use sparse matrix encoding #' (\code{"yes"}) or dense matrix encoding (\code{"no"}). Defaults to #' \code{"auto"}, in which case sparse matrix encoding is used with square #' matrices with at least 50 rows and no more than 50\% of elements with values #' greater than zero. #' @param append_mats A logical value indicating whether to include the original #' A, U, and F matrices in the output \code{lefkoElas} object. #' @param ... Other parameters. #' #' @return This function returns an object of class \code{lefkoElas}, which is a #' list with 8 elements. The first, \code{h_elasmats}, is a list of historical #' elasticity matrices (\code{NULL} if an ahMPM is used as input). The second, #' \code{ah_elasmats}, is a list of either ahistorical elasticity matrices if an #' ahMPM is used as input, or, if an hMPM is used as input, then the result is a #' list of elasticity matrices in which historical elasticities have been summed #' by the stage in occasions \emph{t} and \emph{t}+1 to produce #' historically-corrected elasticity matrices, which are equivalent in dimension #' to ahistorical elasticity matrices but reflect the effects of stage in #' occasion \emph{t}-1. The third element, \code{hstages}, is a data frame #' showing historical stage pairs (NULL if ahMPM used as input). The fourth #' element, \code{agestages}, shows age-stage combinations in the order used in #' age-by-stage MPMs, if suppled. The fifth element, \code{ahstages}, is a data #' frame showing the order of ahistorical stages. The last 3 elements are the A, #' U, and F portions of the input. #' #' @section Notes: #' Deterministic elasticities are estimated as eqn. 9.72 in Caswell (2001, #' Matrix Population Models). Stochastic elasticities are estimated as eqn. #' 14.99 in Caswell (2001). Note that stochastic elasticities are of the #' stochastic \eqn{\lambda}, while stochastic sensitivities are with regard to #' the log of the stochastic \eqn{\lambda}. #' #' Speed can sometimes be increased by shifting from automatic sparse matrix #' determination to forced dense or sparse matrix projection. This will most #' likely occur when matrices have between 30 and 300 rows and columns. #' Defaults work best when matrices are very small and dense, or very large and #' sparse. #' #' The \code{time_weights}, \code{steps}, and \code{force_sparse} arguments are #' now deprecated. Instead, please use the \code{tweights}, \code{times}, and #' \code{sparse} arguments. #' #' @seealso \code{\link{elasticity3}()} #' @seealso \code{\link{elasticity3.dgCMatrix}()} #' @seealso \code{\link{elasticity3.matrix}()} #' @seealso \code{\link{elasticity3.list}()} #' @seealso \code{\link{summary.lefkoElas}()} #' #' @examples #' # Lathyrus example #' data(lathyrus) #' #' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0) #' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm") #' repvector <- c(0, 0, 0, 0, 0, 1, 0) #' obsvector <- c(0, 1, 1, 1, 1, 1, 0) #' matvector <- c(0, 0, 1, 1, 1, 1, 1) #' immvector <- c(1, 1, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5) #' #' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' immstatus = immvector, indataset = indataset, binhalfwidth = binvec, #' propstatus = propvector) #' #' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988, #' patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9, #' juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88", #' fecacol = "Intactseed88", deadacol = "Dead1988", #' nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea", #' censorcol = "Missing1988", censorkeep = NA, censor = TRUE) #' #' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"), #' stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"), #' stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"), #' eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"), #' eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"), #' eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"), #' givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA), #' multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA), #' type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1), #' stageframe = lathframe, historical = TRUE) #' #' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all", #' stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3, #' yearcol = "year2", indivcol = "individ") #' #' elasticity3(ehrlen3, stochastic = TRUE) #' #' # Cypripedium example #' data(cypdata) #' #' sizevector <- c(0, 0, 0, 0, 0, 0, 1, 2.5, 4.5, 8, 17.5) #' stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "XSm", "Sm", "Md", "Lg", #' "XLg") #' repvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1) #' obsvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1) #' matvector <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1) #' immvector <- c(0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 0, 0, 0, 0, 0.5, 0.5, 1, 1, 2.5, 7) #' #' cypframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' propstatus = propvector, immstatus = immvector, indataset = indataset, #' binhalfwidth = binvec) #' #' cypraw_v1 <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004, #' patchidcol = "patch", individcol = "plantid", blocksize = 4, #' sizeacol = "Inf2.04", sizebcol = "Inf.04", sizeccol = "Veg.04", #' repstracol = "Inf.04", repstrbcol = "Inf2.04", fecacol = "Pod.04", #' stageassign = cypframe_raw, stagesize = "sizeadded", NAas0 = TRUE, #' NRasRep = TRUE) #' #' cypsupp2r <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "D", #' "XSm", "Sm", "SD", "P1"), #' stage2 = c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "rep", #' "rep"), #' eststage3 = c(NA, NA, NA, NA, NA, "D", "XSm", "Sm", NA, NA), #' eststage2 = c(NA, NA, NA, NA, NA, "XSm", "XSm", "XSm", NA, NA), #' givenrate = c(0.10, 0.20, 0.20, 0.20, 0.25, NA, NA, NA, NA, NA), #' multiplier = c(NA, NA, NA, NA, NA, NA, NA, NA, 0.5, 0.5), #' type =c(1, 1, 1, 1, 1, 1, 1, 1, 3, 3), #' stageframe = cypframe_raw, historical = FALSE) #' #' 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") #' #' elasticity3(cypmatrix2r) #' #' @export elasticity3.lefkoMat <- function(mats, stochastic = FALSE, times = 10000, tweights = NA, seed = NA, sparse = "auto", append_mats = FALSE, ...) { sparsemethod <- 0 sparse_input <- FALSE possible_args <- c("mats", "stochastic", "times", "tweights", "historical", "seed", "force_sparse", "sparse", "append_mats", "time_weights", "steps") further_args <- list(...) further_args_names <- names(further_args) if (length(setdiff(further_args_names, possible_args)) > 0) { stop("Some arguments not recognized.", call. = FALSE) } if ("time_weights" %in% further_args_names) { stop("Argument time_weights is deprecated. Please use argument tweights instead", call. = FALSE) } if ("steps" %in% further_args_names) { stop("Argument steps is deprecated. Please use argument times instead", call. = FALSE) } if ("force_sparse" %in% further_args_names) { stop("Argument force_sparse is deprecated. Please use argument sparse instead", call. = FALSE) } if (is(mats$A[[1]], "dgCMatrix")) sparse_input <- TRUE

if (is.logical(sparse)) {
if (sparse) {
sparsemethod <- 1
} else sparsemethod <- 0
} else if (is.element(tolower(sparse), c("y", "yes", "yea", "yeah", "t", "true", "ja", "tak"))) {
sparsemethod <- 1
} else if (is.element(tolower(sparse), c("n", "no", "non", "nah", "f", "false", "nein", "nie"))) {
sparsemethod <- 0
} else {
elements_total <- length(mats$A[[1]]) dense_elements <- if (!sparse_input ) { length(which(mats$A[[1]] != 0))
} else {
sum(diff(mats$A[[1]]@p)) } if ((dense_elements / elements_total) < 0.5 & elements_total > 2500) { sparsemethod <- 1 } else sparsemethod <- 0 } if (!is.element("agestages", names(mats))) { mats$agestages <- NA
}

if (!stochastic) {
# Deterministic elasticity analysis
message("Running deterministic analysis...")

baldrick <- if (is.matrix(mats$A)) { if (!sparse_input) { .elas3matrix(mats$A, sparsemethod)
} else {
.elas3sp_matrix(mats$A) } } else if (is.list(mats$A)) {

if (all(is.na(mats$hstages))) { if(!sparse_input) { lapply(mats$A, .elas3matrix, sparsemethod)
} else {
lapply(mats$A, .elas3sp_matrix) } } else { if (!sparse_input) { lapply(mats$A, .elas3hlefko, mats$ahstages, mats$hstages)
} else {
lapply(mats$A, .elas3sp_hlefko, mats$ahstages, mats$hstages) } } } else { stop("Input not recognized.") } if (is.list(mats$A)) {
multiplier <- length(mats$A) } else multiplier <- 1 hlabels <- mats$hstages
ahlabels <- mats$ahstages agelabels <- mats$agestages
new_labels <- mats$labels if (all(is.na(mats$hstages))) {

output <- list(h_elasmats = NULL, ah_elasmats = baldrick, hstages = NULL,
agestages = mats$agestages, ahstages = ahlabels, agestages = agelabels, labels = new_labels) } else { he_list <- lapply(baldrick, function(X) {X$h_emat})
ahe_list <- lapply(baldrick, function(X) {X$ah_emat}) hlabels <- mats$hstages
ahlabels <- mats$ahstages #Originally only the first two columns output <- list(h_elasmats = he_list, ah_elasmats = ahe_list, hstages = hlabels, agestages = mats$agestages, ahstages = ahlabels,
labels = new_labels)
}
} else {
# Stochastic elasticity analysis
if (!is.na(seed[1])) {
set.seed(seed[1])
}

if(!any(is.na(tweights))) {
message("Running stochastic analysis with tweights option...")

returned_cubes <- .stoch_senselas(mats, times = times, historical = FALSE,
style = 2, sparsemethod, lefkoProj = TRUE, tweights = tweights)
} else {
message("Running stochastic analysis...")

returned_cubes <- .stoch_senselas(mats, times = times, historical = FALSE,
style = 2, sparsemethod, lefkoProj = TRUE)
}

old_labels <- mats$labels old_labels$concat <- paste(old_labels$pop, old_labels$patch)
labs_unique <- unique(old_labels$concat) lab_indices <- apply(as.matrix(labs_unique), 1, function(X) { found_guys <- which(old_labels$concat == X)
min_guy <- min(found_guys)

return(min_guy)
})
new_labels <- mats$labels[lab_indices, c(1, 2)] if (!all(is.na(mats$hstages))) {
output <- list(h_elasmats = returned_cubes[[1]], ah_elasmats = returned_cubes[[2]],
hstages = mats$hstages, agestages = mats$agestages,
ahstages = mats$ahstages, labels = new_labels) } else { output <- list(h_elasmats = NULL, ah_elasmats = returned_cubes[[1]], hstages = mats$hstages, agestages = mats$agestages, ahstages = mats$ahstages, labels = new_labels)
}
}

if (append_mats) {
output$A <- mats$A
output$U <- mats$U
output$F <- mats$F
}

class(output) <- "lefkoElas"

return(output)
}

#' Estimate Elasticity of Population Growth Rate of a Single Matrix
#'
#' \code{elasticity3.matrix()} returns the elasticities of lambda to elements
#' of a single matrix.  Because this handles only one matrix, the elasticities
#' are inherently deterministic and based on the dominant eigen value as the
#' best metric of the population growth rate. This function can handle large and
#' sparse matrices, and so can be used with large historical matrices, IPMs,
#' age x stage matrices, as well as smaller ahistorical matrices.
#'
#' @name elasticity3.matrix
#'
#' @param mats An object of class \code{matrix}.
#' @param sparse A text string indicating whether to use sparse matrix encoding
#' (\code{"yes"}) or dense matrix encoding (\code{"no"}). Defaults to
#' \code{"auto"}, in which case sparse matrix encoding is used with square
#' matrices with at least 50 rows and no more than 50\% of elements with values
#' greater than zero.
#' @param ... Other parameters.
#'
#' @return This function returns a single elasticity matrix.
#'
#' @section Notes:
#' Speed can sometimes be increased by shifting from automatic sparse matrix
#' determination to forced dense or sparse matrix projection. This will most
#' likely occur when matrices have between 30 and 300 rows and columns.
#' Defaults work best when matrices are very small and dense, or very large and
#' sparse.
#'
#' The \code{force_sparse} argument is now deprecated. Please use \code{sparse}
#'
#'
#' @examples
#' # Lathyrus example
#' data(lathyrus)
#'
#' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0)
#' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm")
#' repvector <- c(0, 0, 0, 0, 0, 1, 0)
#' obsvector <- c(0, 1, 1, 1, 1, 1, 0)
#' matvector <- c(0, 0, 1, 1, 1, 1, 1)
#' immvector <- c(1, 1, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5)
#'
#' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   immstatus = immvector, indataset = indataset, binhalfwidth = binvec,
#'   propstatus = propvector)
#'
#' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
#'   patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
#'   juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88",
#'   nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea",
#'   censorcol = "Missing1988", censorkeep = NA, censor = TRUE)
#'
#' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"),
#'   stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"),
#'   stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"),
#'   eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"),
#'   eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"),
#'   eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"),
#'   givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA),
#'   type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1),
#'   stageframe = lathframe, historical = TRUE)
#'
#' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all",
#'   stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3,
#'   yearcol = "year2", indivcol = "individ")
#'
#' ehrlen3mean <- lmean(ehrlen3)
#' elasticity3(ehrlen3mean$A[[1]]) #' #' @export elasticity3.matrix <- function(mats, sparse = "auto", ...) { sparsemethod <- 0 possible_args <- c("mats", "force_sparse", "sparse") further_args <- list(...) further_args_names <- names(further_args) if (length(setdiff(further_args_names, possible_args)) > 0) { stop("Some arguments not recognized.", call. = FALSE) } if ("force_sparse" %in% further_args_names) { stop("Argument force_sparse is deprecated. Please use argument sparse instead", call. = FALSE) } if (is.logical(sparse[1])) { if (sparse[1]) { sparsemethod <- 1 } else sparsemethod <- 0 } else if (is.element(tolower(sparse), c("y", "yes", "yea", "yeah", "t", "true", "ja", "tak"))) { sparsemethod <- 1 } else if (is.element(tolower(sparse), c("n", "no", "non", "nah", "f", "false", "nein", "nie"))) { sparsemethod <- 0 } else { elements_total <- length(mats) dense_elements <- length(which(mats != 0)) if ((dense_elements / elements_total) < 0.5 & elements_total > 2500) { sparsemethod <- 1 } else sparsemethod <- 0 } wcorr <- .elas3matrix(mats, FALSE) if (sparsemethod == 1) wcorr <- as(wcorr, "dgCMatrix") return(wcorr) } #' Estimate Elasticity of Population Growth Rate of a Single Sparse Matrix #' #' \code{elasticity3.dgCMatrix()} returns the elasticities of lambda to elements #' of a single matrix. Because this handles only one matrix, the elasticities #' are inherently deterministic and based on the dominant eigen value as the #' best metric of the population growth rate. #' #' @name elasticity3.dgCMatrix #' #' @param mats An object of class \code{dgCMatrix}. #' @param sparse A text string indicating whether to use sparse matrix encoding #' (\code{"yes"}) or dense matrix encoding (\code{"no"}). Defaults to #' \code{"auto"}, in which case sparse matrix encoding is used with square #' matrices with at least 50 rows and no more than 50\% of elements with values #' greater than zero. #' @param ... Other parameters. #' #' @return This function returns a single elasticity matrix in \code{dgCMatrix} #' format. #' #' @seealso \code{\link{elasticity3}()} #' @seealso \code{\link{elasticity3.lefkoMat}()} #' @seealso \code{\link{elasticity3.list}()} #' @seealso \code{\link{elasticity3.matrix}()} #' @seealso \code{\link{summary.lefkoElas}()} #' #' @examples #' # Lathyrus example #' data(lathyrus) #' #' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0) #' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm") #' repvector <- c(0, 0, 0, 0, 0, 1, 0) #' obsvector <- c(0, 1, 1, 1, 1, 1, 0) #' matvector <- c(0, 0, 1, 1, 1, 1, 1) #' immvector <- c(1, 1, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5) #' #' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' immstatus = immvector, indataset = indataset, binhalfwidth = binvec, #' propstatus = propvector) #' #' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988, #' patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9, #' juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88", #' fecacol = "Intactseed88", deadacol = "Dead1988", #' nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea", #' censorcol = "Missing1988", censorkeep = NA, censor = TRUE) #' #' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"), #' stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"), #' stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"), #' eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"), #' eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"), #' eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"), #' givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA), #' multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA), #' type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1), #' stageframe = lathframe, historical = TRUE) #' #' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all", #' stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3, #' yearcol = "year2", indivcol = "individ") #' #' ehrlen3mean <- lmean(ehrlen3) #' elasticity3(ehrlen3mean$A[[1]])
#'
#' @export
elasticity3.dgCMatrix <- function(mats, sparse = "auto", ...)
{

possible_args <- c("mats", "sparse")
further_args <- list(...)
further_args_names <- names(further_args)
if (length(setdiff(further_args_names, possible_args)) > 0) {
stop("Some arguments not recognized.", call. = FALSE)
}

if (is.logical(sparse)) {
if (sparse) {
sparsemethod <- 1
} else sparsemethod <- 0
} else if (is.element(tolower(sparse), c("y", "yes", "yea", "yeah", "t", "true", "ja", "tak"))) {
sparsemethod <- 1
} else if (is.element(tolower(sparse), c("n", "no", "non", "nah", "f", "false", "nein", "nie"))) {
sparsemethod <- 0
} else {
elements_total <- length(mats)
dense_elements <- sum(diff(mats@p))

if ((dense_elements / elements_total) < 0.5 & elements_total > 2500) {
sparsemethod <- 1
} else sparsemethod <- 0
}

wcorr <- .elas3sp_matrix(mats)

if (sparsemethod == 0) wcorr <- as.matrix(wcorr)

return(wcorr)
}

#' Estimate Elasticity of Population Growth Rate of a List of Matrices
#'
#' \code{elasticity3.list()} returns the elasticities of lambda to elements
#' of a single matrix. This function can handle large and sparse matrices, and
#' so can be used with large historical matrices, IPMs, age x stage matrices,
#' as well as smaller ahistorical matrices.
#'
#' @name elasticity3.list
#'
#' @param mats A list of objects of class \code{matrix} or \code{dgCMatrix}.
#' @param stochastic A logical value determining whether to conduct a
#' deterministic (FALSE) or stochastic (TRUE) elasticity analysis. Defaults to
#' FALSE.
#' @param times The number of occasions to project forward in stochastic
#' simulation. Defaults to 10,000.
#' @param tweights An optional numeric vector or matrix denoting the
#' probabilities of choosing each matrix in a stochastic projection. If a matrix
#' is input, then a first-order Markovian environment is assumed, in which the
#' probability of choosing a specific annual matrix depends on which annual
#' matrix is currently chosen. If a vector is input, then the choice of annual
#' matrix is assumed to be independent of the current matrix. Defaults to equal
#' weighting among matrices.
#' @param historical A logical value denoting whether the input matrices are
#' historical. Defaults to FALSE.
#' @param seed A number to use as a random number seed in stochastic projection.
#' @param sparse A text string indicating whether to use sparse matrix encoding
#' (\code{"yes"}) or dense matrix encoding (\code{"no"}). Defaults to
#' \code{"auto"}, in which case sparse matrix encoding is used with square
#' matrices with at least 50 rows and no more than 50\% of elements with values
#' greater than zero.
#' @param append_mats A logical value indicating whether to include the original
#' matrices input as object \code{mats} in the output \code{lefkoElas} object.
#' @param ... Other parameters.
#'
#' @return This function returns an object of class \code{lefkoElas}, which is a
#' list with 8 elements. The first, \code{h_elasmats}, is a list of historical
#' elasticity matrices, though in the standard list case it returns a NULL
#' value. The second, \code{ah_elasmats}, is a list of ahistorical elasticity
#' matrices. The third element, \code{hstages}, the fourth element,
#' \code{agestages}, and the fifth element, \code{ahstages}, are set to NULL.
#' The last 3 elements are the original A matrices in element A, followed by
#' NULL values for the U and F elements.
#'
#' @section Notes:
#' Deterministic elasticities are estimated as eqn. 9.72 in Caswell (2001,
#' Matrix Population Models). Stochastic elasticities are estimated as eqn.
#' 14.99 in Caswell (2001). Note that stochastic elasticities are of stochastic
#' \eqn{\lambda}, while stochastic sensitivities are with regard to the log of
#' the stochastic \eqn{\lambda}.
#'
#' Speed can sometimes be increased by shifting from automatic sparse matrix
#' determination to forced dense or sparse matrix projection. This will most
#' likely occur when matrices have between 30 and 300 rows and columns.
#' Defaults work best when matrices are very small and dense, or very large and
#' sparse.
#'
#' The \code{time_weights}, \code{steps}, and \code{force_sparse} arguments are
#' \code{sparse} arguments.
#'
#'
#' @examples
#' # Lathyrus example
#' data(lathyrus)
#'
#' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0)
#' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm")
#' repvector <- c(0, 0, 0, 0, 0, 1, 0)
#' obsvector <- c(0, 1, 1, 1, 1, 1, 0)
#' matvector <- c(0, 0, 1, 1, 1, 1, 1)
#' immvector <- c(1, 1, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5)
#'
#' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   immstatus = immvector, indataset = indataset, binhalfwidth = binvec,
#'   propstatus = propvector)
#'
#' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
#'   patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
#'   juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88",
#'   nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea",
#'   censorcol = "Missing1988", censorkeep = NA, censor = TRUE)
#'
#' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"),
#'   stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"),
#'   stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"),
#'   eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"),
#'   eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"),
#'   eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"),
#'   givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA),
#'   type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1),
#'   stageframe = lathframe, historical = TRUE)
#'
#' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all",
#'   stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3,
#'   yearcol = "year2", indivcol = "individ")
#'
#' elasticity3(ehrlen3$A, stochastic = TRUE) #' #' # Cypripedium example #' data(cypdata) #' #' sizevector <- c(0, 0, 0, 0, 0, 0, 1, 2.5, 4.5, 8, 17.5) #' stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "XSm", "Sm", "Md", "Lg", #' "XLg") #' repvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1) #' obsvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1) #' matvector <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1) #' immvector <- c(0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 0, 0, 0, 0, 0.5, 0.5, 1, 1, 2.5, 7) #' #' cypframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' propstatus = propvector, immstatus = immvector, indataset = indataset, #' binhalfwidth = binvec) #' #' cypraw_v1 <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004, #' patchidcol = "patch", individcol = "plantid", blocksize = 4, #' sizeacol = "Inf2.04", sizebcol = "Inf.04", sizeccol = "Veg.04", #' repstracol = "Inf.04", repstrbcol = "Inf2.04", fecacol = "Pod.04", #' stageassign = cypframe_raw, stagesize = "sizeadded", NAas0 = TRUE, #' NRasRep = TRUE) #' #' cypsupp2r <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "D", #' "XSm", "Sm", "SD", "P1"), #' stage2 = c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "rep", #' "rep"), #' eststage3 = c(NA, NA, NA, NA, NA, "D", "XSm", "Sm", NA, NA), #' eststage2 = c(NA, NA, NA, NA, NA, "XSm", "XSm", "XSm", NA, NA), #' givenrate = c(0.10, 0.20, 0.20, 0.20, 0.25, NA, NA, NA, NA, NA), #' multiplier = c(NA, NA, NA, NA, NA, NA, NA, NA, 0.5, 0.5), #' type =c(1, 1, 1, 1, 1, 1, 1, 1, 3, 3), #' stageframe = cypframe_raw, historical = FALSE) #' #' 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") #' #' elasticity3(cypmatrix2r$A)
#'
#' @export
elasticity3.list <- function(mats, stochastic = FALSE, times = 10000,
tweights = NA, historical = FALSE, seed = NA, sparse = "auto",
append_mats = FALSE, ...) {

sparsemethod <- 0
sparse_input <- FALSE

possible_args <- c("mats", "stochastic", "times", "tweights", "historical",
"seed", "force_sparse", "append_mats", "time_weights", "steps")
further_args <- list(...)
further_args_names <- names(further_args)
if (length(setdiff(further_args_names, possible_args)) > 0) {
stop("Some arguments not recognized.", call. = FALSE)
}
if ("time_weights" %in% further_args_names) {
call. = FALSE)
}
if ("steps" %in% further_args_names) {
call. = FALSE)
}
if ("force_sparse" %in% further_args_names) {
call. = FALSE)
}

if (is(mats[[1]], "dgCMatrix")) sparse_input <- TRUE

if (is.logical(sparse)) {
if (sparse) {
sparsemethod <- 1
} else sparsemethod <- 0
} else if (is.element(tolower(sparse), c("y", "yes", "yea", "yeah", "t", "true", "ja", "tak"))) {
sparsemethod <- 1
} else if (is.element(tolower(sparse), c("n", "no", "non", "nah", "f", "false", "nein", "nie"))) {
sparsemethod <- 0
} else {
elements_total <- length(mats[[1]])
dense_elements <- if (!sparse_input ) {
length(which(mats[[1]] != 0))
} else {
sum(diff(mats[[1]]@p))
}

if ((dense_elements / elements_total) < 0.5 & elements_total > 2500) {
sparsemethod <- 1
} else sparsemethod <- 0
}

if(length(setdiff(unlist(lapply(mats, class)), c("matrix", "array", "dgCMatrix"))) > 0) {
stop("Input list must be composed only of numeric matrices.", call. = FALSE)
}
if (stochastic & length(mats) < 2) {
stop("Stochastic elasticity estimation cannot be completed with fewer than 2 annual matrices.", call. = FALSE)
}

if (!stochastic) {
# Deterministic elasticity analysis
message("Running deterministic analysis...")

baldrick <- if (!sparse_input) {
lapply(mats, .elas3matrix, sparsemethod)
} else {
lapply(mats, .elas3sp_matrix)
}

if (historical) {
output <- list(h_elasmats = baldrick, ah_elasmats = NULL, hstages = NULL,
ahstages = NULL)
} else {
output <- list(h_elasmats = NULL, ah_elasmats = baldrick, hstages = NULL,
ahstages = NULL)
}
} else {
# Stochastic elasticity analysis
if (!is.na(seed[1])) {
set.seed(seed[1])
}

if(!any(is.na(tweights))) {
message("Running stochastic analysis with tweights option...")

returned_cube <- .stoch_senselas(mats, times = times, historical = historical,
style = 2, sparsemethod, lefkoProj = FALSE, tweights = tweights)[[1]]
} else {
message("Running stochastic analysis...")

returned_cube <- .stoch_senselas(mats, times = times, historical = historical,
style = 2, sparsemethod, lefkoProj = FALSE)[[1]]
}

if (historical) {
output <- list(h_elasmats = list(returned_cube[[1]]), ah_elasmats = NULL,
hstages = NULL, ahstages = NULL)
} else {
output <- list(h_elasmats = NULL, ah_elasmats = list(returned_cube[[1]]),
hstages = NULL, ahstages = NULL)
}
}

if (append_mats) {
output$A <- mats } class(output) <- "lefkoElas" return(output) } #' Conduct a Life Table Response Experiment #' #' \code{ltre3()} returns a set of matrices of one-way LTRE (life table response #' experiment), stochastic LTRE (sLTRE) matrices, or small noise approximation #' LTRE (sna-LTRE) contributions. #' #' @name ltre3 #' #' @param mats An object of class \code{lefkoMat}. #' @param refmats A reference lefkoMat object, or matrix, for use as the #' control. Default is \code{NA}, which sets to the same object as \code{mats}. #' @param ref A numeric value indicating which matrix or matrices in #' \code{refmats} to use as the control. The numbers used must correspond to the #' number of the matrices in the \code{labels} element of the associated #' \code{lefkoMat} object. The default setting, \code{NA}, uses all entries in #' \code{refmats}. #' @param stochastic A logical value determining whether to conduct a #' deterministic (\code{FALSE}) or stochastic (\code{TRUE}) elasticity analysis. #' Defaults to \code{FALSE}. #' @param times The number of occasions to project forward in stochastic #' simulation. Defaults to \code{10000}. #' @param burnin The number of initial steps to ignore in stochastic projection #' when calculating stochastic elasticities. Must be smaller than \code{steps}. #' Defaults to \code{3000}. #' @param tweights An optional numeric vector or matrix denoting the #' probabilities of choosing each matrix in a stochastic projection. If a matrix #' is input, then a first-order Markovian environment is assumed, in which the #' probability of choosing a specific annual matrix depends on which annual #' matrix is currently chosen. If a vector is input, then the choice of annual #' matrix is assumed to be independent of the current matrix. Defaults to equal #' weighting among matrices. Note that SNA-LTRE analysis cannot take matrix #' input. #' @param sparse A text string indicating whether to use sparse matrix encoding #' (\code{"yes"}) or dense matrix encoding (\code{"no"}). Defaults to #' \code{"auto"}, in which case sparse matrix encoding is used with square #' matrices with at least 50 rows and no more than 50\% of elements with values #' greater than zero. #' @param seed Optional numeric value corresponding to the random seed for #' stochastic simulation. #' @param append_mats A logical value denoting whether to include the original #' \code{A}, \code{U}, and \code{F} matrices in the returned \code{lefkoLTRE} #' object. Defaults to \code{FALSE}. #' @param sna_ltre A logical value indicating whether to treat stochastic LTRE #' via the sna-LTRE approach from Davison et al. (2019) (\code{TRUE}), or the #' stochastic LTRE approximation from Davison et al. (2010) (\code{FALSE}). #' Defaults to \code{FALSE}. #' @param tol A numeric value indicating a lower positive limit to matrix #' element values when applied to stochastic and small noise approximation LTRE #' estimation protocols. Matrix element values lower than this will be treated #' as \code{0.0} values. Defaults to \code{1e-30}. #' @param ... Other parameters. #' #' @return This function returns an object of class \code{lefkoLTRE}. This #' includes a list of LTRE matrices as object \code{cont_mean} if a #' deterministic LTRE is called for, or a list of mean-value LTRE matrices as #' object \code{cont_mean} and a list of SD-value LTRE matrices as object #' \code{cont_sd} if a stochastic LTRE is called for. If a small-noise #' approximation LTRE (SNA-LTRE) is performed, then the output includes six #' objects: \code{cont_mean}, which provides the contributions of shifts in mean #' matrix elements; \code{cont_elas}, which provides the contributions of shifts #' in the elasticities of matrix elements; \code{cont_cv}, which provides the #' contributions of temporal variation in matrix elements; \code{cont_corr}, #' which provides the contributions of temporal correlations in matrix elements; #' \code{r_values_m}, which provides a vector of log deterministic lambda values #' for treatment populations; and \code{r_values_ref}, which provides the log #' deterministic lambda of the mean reference matrix.This is followed by the #' stageframe as object \code{ahstages}, the order of historical stages as #' object \code{hstages}, the age-by-stage order as object \code{agestages}, the #' order of matrices as object \code{labels}, and, if requested, the original A, #' U, and F matrices. #' #' @section Notes: #' Deterministic LTRE is one-way, fixed, and based on the sensitivities of the #' matrix midway between each input matrix and the reference matrix, per Caswell #' (2001, Matrix Population Models, Sinauer Associates, MA, USA). Stochastic #' LTRE is performed via two methods. The stochastic LTRE approximation is #' simulated per Davison et al. (2010) Journal of Ecology 98:255-267 #' (doi: 10.1111/j.1365-2745.2009.01611.x). The small noise approximation #' (sna-LTRE) is analyzed per Davison et al. (2019) Ecological Modelling 408: #' 108760 (doi: 10.1016/j.ecolmodel.2019.108760). #' #' All stochastic and small noise approximation LTREs conducted without #' reference matrices are performed as spatial tests of the population dynamics #' among patches. #' #' Default behavior for stochastic LTRE uses the full population provided in #' \code{mats} as the reference if no \code{refmats} and \code{ref} is provided. #' If no \code{refmats} is provided but \code{ref} is, then the matrices noted #' in \code{ref} are used as the reference matrix set. Year and patch order is #' utilized from object \code{mats}, but not from object \code{refmats}, in #' which each matrix is assumed to represent a different year from one #' population. This function cannot currently handle multiple populations within #' the same \code{mats} object (although such analysis is possible if these #' populations are designated as patches instead). #' #' If \code{sparse = "auto"}, the default, then sparse matrix encoding #' will be used if the size of the input matrices is at least 50 columns by 50 #' rows for deterministic and stochastic LTREs and 10 columns by 10 rows for #' small noise approximation LTREs, in all cases as long as 50\% of the elements #' in the first matrix are non-zero. #' #' Stochastic LTREs do not test for the impact of temporal change in vital #' rates. An MPM with a single population, a single patch, and only annual #' matrices will produce contributions of 0 to stochastic \eqn{\lambda}. #' #' Speed can sometimes be increased by shifting from automatic sparse matrix #' determination to forced dense or sparse matrix projection. This will most #' likely occur when matrices have between 10 and 300 rows and columns. #' Defaults work best when matrices are very small and dense, or very large and #' sparse. #' #' SNA-LTRE analysis cannot test the impact of first-order Markovian #' environments. However, different random weightings of annual matrices are #' allowed if given in vector format. #' #' The \code{time_weights}, \code{steps}, \code{force_sparse}, and \code{rseed} #' arguments are now deprecated. Instead, please use the \code{tweights}, #' \code{times}, \code{sparse}, and \code{seed} arguments. #' #' @seealso \code{\link{summary.lefkoLTRE}()} #' #' @examples #' data(cypdata) #' #' sizevector <- c(0, 0, 0, 0, 0, 0, 1, 2.5, 4.5, 8, 17.5) #' stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "XSm", "Sm", "Md", "Lg", #' "XLg") #' repvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1) #' obsvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1) #' matvector <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1) #' immvector <- c(0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 0, 0, 0, 0, 0.5, 0.5, 1, 1, 2.5, 7) #' #' cypframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' propstatus = propvector, immstatus = immvector, indataset = indataset, #' binhalfwidth = binvec) #' #' cypraw_v1 <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004, #' patchidcol = "patch", individcol = "plantid", blocksize = 4, #' sizeacol = "Inf2.04", sizebcol = "Inf.04", sizeccol = "Veg.04", #' repstracol = "Inf.04", repstrbcol = "Inf2.04", fecacol = "Pod.04", #' stageassign = cypframe_raw, stagesize = "sizeadded", NAas0 = TRUE, #' NRasRep = TRUE) #' #' cypsupp2r <- supplemental(stage3 = c("SD", "P1", "P2", "P3", "SL", "D", #' "XSm", "Sm", "SD", "P1"), #' stage2 = c("SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "rep", #' "rep"), #' eststage3 = c(NA, NA, NA, NA, NA, "D", "XSm", "Sm", NA, NA), #' eststage2 = c(NA, NA, NA, NA, NA, "XSm", "XSm", "XSm", NA, NA), #' givenrate = c(0.10, 0.20, 0.20, 0.20, 0.25, NA, NA, NA, NA, NA), #' multiplier = c(NA, NA, NA, NA, NA, NA, NA, NA, 0.5, 0.5), #' type =c(1, 1, 1, 1, 1, 1, 1, 1, 3, 3), #' stageframe = cypframe_raw, historical = FALSE) #' #' 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") #' #' ltre3(cypmatrix2r, sna_ltre = TRUE) #' #' @export ltre3 <- function(mats, refmats = NA, ref = NA, stochastic = FALSE, times = 10000, burnin = 3000, tweights = NA, sparse = "auto", seed = NA, append_mats = FALSE, sna_ltre = FALSE, tol = 1e-30, ...) { sparsemethod <- 0 sparse_input <- FALSE possible_args <- c("mats", "stochastic", "times", "tweights", "rseed", "seed", "force_sparse", "sparse", "append_mats", "time_weights", "steps", "refmats", "ref", "sna_ltre", "tol") further_args <- list(...) further_args_names <- names(further_args) if (length(setdiff(further_args_names, possible_args)) > 0) { stop("Some arguments not recognized.", call. = FALSE) } if ("time_weights" %in% further_args_names) { stop("Argument time_weights is deprecated. Please use argument tweights instead", call. = FALSE) } if ("steps" %in% further_args_names) { stop("Argument steps is deprecated. Please use argument times instead", call. = FALSE) } if ("force_sparse" %in% further_args_names) { stop("Argument force_sparse is deprecated. Please use argument sparse instead", call. = FALSE) } if ("rseed" %in% further_args_names) { stop("Argument rseed is deprecated. Please use argument seed instead", call. = FALSE) } if (!is(mats, "lefkoMat")) stop("Function ltre3() requires a lefkoMat object as input.", call. = FALSE) if (is(mats$A[[1]], "dgCMatrix")) sparse_input = TRUE

if (is.logical(sparse)) {
if (sparse) {
sparsemethod <- 1
} else sparsemethod <- 0
} else if (is.element(tolower(sparse), c("y", "yes", "yea", "yeah", "t", "true", "ja", "tak"))) {
sparsemethod <- 1
} else if (is.element(tolower(sparse), c("n", "no", "non", "nah", "f", "false", "nein", "nie"))) {
sparsemethod <- 0
} else {
elements_total <- length(mats$A[[1]]) dense_elements <- if (!sparse_input ) { length(which(mats$A[[1]] != 0))
} else {
sum(diff(mats$A[[1]]@p)) } if ((dense_elements / elements_total) < 0.5 & elements_total > 2500) { sparsemethod <- 1 } else sparsemethod <- 0 } if (!is.element("agestages", names(mats))) { mats$agestages <- NULL
}

if (all(is.na(refmats))) {
warning("Matrices input as mats will also be used in reference matrix calculation.",
call. = FALSE)
} else if (is(refmats, "lefkoMat")) {
refmats <- refmats$A } else if (is.list(refmats)) { if (!is.matrix(refmats[[1]]) & !is(refmats[[1]], "dgCMatrix")) { stop("Object refmats not recognized. Use only objects of class lefkoMat, list, matrix, or dgCMatrix.", call. = FALSE) } } else if (is.matrix(refmats) | is(refmats, "dgCMatrix")) { refmats <- list(refmats) } else { stop("Object refmats not recognized. Use only objects of class lefkoMat, list, or matrix.", call. = FALSE) } if (all(is.numeric(ref))) { if (length(ref) > 1) { if (!stochastic) { message("This function currently conducts deterministic LTREs against only single matrices, or mean matrices. Will use the mean.") } } } else if (all(is.na(ref))) { message("Using all refmats matrices in reference matrix calculation.") if (!all(is.na(refmats))) { ref <- c(1:length(refmats)) } else { ref <- c(1:length(mats$A))
}
}

if (sna_ltre & !stochastic) {
stochastic <- TRUE
}

if (!stochastic & length(ref) > 1) {
meanout <- 1
} else if (stochastic & length(ref) == 1) {
stop("Stochastic LTRE requires multiple matrices in object refmats.", call. = FALSE)
} else if (stochastic & length(mats$A) == 1) { stop("Stochastic LTRE requires multiple matrices in object mats.", call. = FALSE) } else { meanout <- 0; } if (all(is.na(refmats))) { if (any(ref > length(mats$A)) | any(ref < 1)) {
stop("Numbers of matrices provided in object ref must be integers between 1 and
the total number of matrices in object mats (or refmats, if given).",
call. = FALSE)
}

} else {
if (dim(mats$A[[1]])[1] != dim(refmats[[1]])[1] | dim(mats$A[[1]])[2] != dim(refmats[[1]])[2]) {
stop("Matrices used in objects mats and refmats must be of equal dimension.",
call. = FALSE)
}
}

if (!stochastic) {
# Deterministic LTRE analysis

if (all(is.na(refmats))) {
baldrick <- .ltre3matrix(mats$A, refnum = ref, mean = meanout, sparse = sparsemethod) } else { baldrick <- .ltre3matrix(mats$A, refnum = ref, refmats_ = refmats,
mean = meanout, sparse = sparsemethod)
}

baldrick$ahstages <- mats$ahstages
baldrick$agestages <- mats$agestages
baldrick$hstages <- mats$hstages
baldrick$labels <- mats$labels

} else {
# Stochastic and SNA LTRE analysis
if (burnin >= times) {
stop("Option burnin must be smaller than option steps.", call. = FALSE)
}

if (!is.na(seed[1])) {
set.seed(seed[1])
}

if (!sna_ltre) {
if (all(is.na(refmats))) {
if (all(is.na(tweights))) {
baldrick <- .sltre3matrix(mats$A, labels = mats$labels, refnum = ref,
steps = times, burnin = burnin, sparse = sparsemethod,
tol_used = tol)
} else {
baldrick <- .sltre3matrix(mats$A, labels = mats$labels, refnum = ref,
tweights_ = tweights, steps = times, burnin = burnin,
sparse = sparsemethod, tol_used = tol)
}
} else {
if (all(is.na(tweights))) {
baldrick <- .sltre3matrix(mats$A, labels = mats$labels, refnum = ref,
refmats_ = refmats, steps = times, burnin = burnin,
sparse = sparsemethod, tol_used = tol)
} else {
baldrick <- .sltre3matrix(mats$A, labels = mats$labels, refnum = ref,
refmats_ = refmats, tweights_ = tweights, steps = times,
burnin = burnin, sparse = sparsemethod, tol_used = tol)
}
}
} else {
if (is.matrix(tweights)) {
stop("SNA-LTRE analysis cannot proceed with matrix input in the tweights argument. Use vectors only.",
call. = FALSE)
}

if (all(is.na(refmats))) {
if (all(is.na(tweights))) {
baldrick <- .snaltre3matrix(mats$A, labels = mats$labels,
refnum = ref, sparse = sparsemethod, tol_used = tol)
} else {
baldrick <- .snaltre3matrix(mats$A, labels = mats$labels,
refnum = ref, tweights_ = tweights, sparse = sparsemethod,
tol_used = tol)
}
} else {
if (all(is.na(tweights))) {
baldrick <- .snaltre3matrix(mats$A, labels = mats$labels,
refnum = ref, refmats_ = refmats, sparse = sparsemethod,
tol_used = tol)
} else {
baldrick <- .snaltre3matrix(mats$A, labels = mats$labels,
refnum = ref, refmats_ = refmats, tweights_ = tweights,
sparse = sparsemethod, tol_used = tol)
}
}
}

alabels <- mats$labels[,c("pop", "patch")] alabels <- unique(alabels) baldrick$ahstages <- mats$ahstages baldrick$agestages <- mats$agestages baldrick$hstages <- mats$hstages baldrick$labels <- alabels
}

if (append_mats) {
baldrick$A <- mats$A
baldrick$U <- mats$U
baldrick$F <- mats$F
}

return(baldrick)
}

#' Summarize lefkoElas Objects
#'
#' Function \code{summary.lefkoElas()} summarizes \code{lefkoElas} objects.
#' Particularly, it breaks down elasticity values by the kind of ahistorical
#' and, if applicable, historical transition.
#'
#' @name summary.lefkoElas
#'
#' @param object A \code{lefkoElas} object.
#' @param ... Other parameters currently not utilized.
#'
#' @return A list composed of 2 data frames. The first, \code{hist}, is a data
#' frame showing the summed elasticities for all 16 kinds of historical
#' transition per matrix, with each column corresponding to each elasticity
#' matrix in order. The second, \code{ahist}, is a data frame showing the
#' summed elasticities for all 4 kinds of ahistorical transition per matrix,
#' with each column corresponding to each elasticity matrix in order.
#'
#' @examples
#' data(lathyrus)
#'
#' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0)
#' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm")
#' repvector <- c(0, 0, 0, 0, 0, 1, 0)
#' obsvector <- c(0, 1, 1, 1, 1, 1, 0)
#' matvector <- c(0, 0, 1, 1, 1, 1, 1)
#' immvector <- c(1, 1, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5)
#'
#' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   immstatus = immvector, indataset = indataset, binhalfwidth = binvec,
#'   propstatus = propvector)
#'
#' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
#'   patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
#'   juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88",
#'   nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea",
#'   censorcol = "Missing1988", censorkeep = NA, censor = TRUE)
#'
#' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"),
#'   stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"),
#'   stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"),
#'   eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"),
#'   eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"),
#'   eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"),
#'   givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA),
#'   type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1),
#'   stageframe = lathframe, historical = TRUE)
#'
#' lathsupp2 <- supplemental(stage3 = c("Sd", "Sdl", "Sd", "Sdl"),
#'   stage2 = c("Sd", "Sd", "rep", "rep"),
#'   givenrate = c(0.345, 0.054, NA, NA),
#'   multiplier = c(NA, NA, 0.345, 0.054),
#'   type = c(1, 1, 3, 3), stageframe = lathframe, historical = FALSE)
#'
#' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all",
#'   stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3,
#'   yearcol = "year2", indivcol = "individ")
#'
#' ehrlen2 <- rlefko2(data = lathvert, stageframe = lathframe, year = "all",
#'   stages = c("stage3", "stage2"), supplement = lathsupp2,
#'   yearcol = "year2", indivcol = "individ")
#'
#' ehrlen3elas <- elasticity3(ehrlen3)
#' ehrlen2elas <- elasticity3(ehrlen2)
#'
#' summary(ehrlen3elas)
#' summary(ehrlen2elas)
#'
#' @export
summary.lefkoElas <- function(object, ...) {

sparse_input <- FALSE
num_h_mats <- length(object$h_elasmats) num_ah_mats <- length(object$ah_elasmats)

used_emats <- if(num_h_mats == 0) object$ah_elasmats else object$h_elasmats

used_iterations <- if(num_h_mats > 0) num_h_mats else num_ah_mats
if (is(used_emats[[1]], "dgCMatrix")) sparse_input <- TRUE

if (num_h_mats == 0) {
if (!all(is.null(object$agestages))) { if (!all(is.na(object$agestages))) {
if (is.element("stage_id", names(object$agestages))) { new_ahstages_list <- apply(as.matrix(c(1:length(object$agestages$stage_id))), 1, function(X) { return(object$ahstages[which(object$ahstages$stage_id == object$agestages$stage_id[X]),])
})
new_ahstages <- do.call("rbind.data.frame", new_ahstages_list)
indices <- .bambi2(new_ahstages)
} else {
indices <- .bambi2(object$ahstages) } } else { indices <- .bambi2(object$ahstages)
}
} else {
if (all(is.null(object$ahstages))) { mat_size <- dim(object$ah_elasmats[[1]])[1]
new_sf <- sf_skeleton(mat_size, standard = FALSE)

indices <- .bambi2(new_sf)
} else {
indices <- .bambi2(object$ahstages) } } } else { indices <- .bambi3(object$ahstages, object$hstages) } for (i in c(1:used_iterations)) { if (!sparse_input) { trialguy <- .demolition3(used_emats[[i]], indices) } else { trialguy <- .demolition3sp(used_emats[[i]], indices) } if (i == 1) { if (num_h_mats > 0) hist <- trialguy$hist[,c(1,2)]
ahist <- trialguy$ahist[,c(1,2)] if (num_h_mats > 0) names(hist)[2] <- "matrix1" names(ahist)[2] <- "matrix1" } else { if (num_h_mats > 0) hist <- cbind.data.frame(hist, trialguy$hist[,2])
ahist <- cbind.data.frame(ahist, trialguy$ahist[,2]) if (num_h_mats > 0) names(hist)[(i+1)] <- paste0("matrix", i) names(ahist)[(i+1)] <- paste0("matrix", i) } } if (num_h_mats > 0) { output <- list(hist = hist, ahist = ahist) } else { output <- list(hist = NULL, ahist = ahist) } return (output) } #' Summarize lefkoLTRE Objects #' #' Function \code{summary.lefkoLTRE()} summarizes \code{lefkoLTRE} objects. #' Particularly, it breaks down LTRE contributions by the kind of ahistorical #' and, if applicable, historical transition. #' #' @name summary.lefkoLTRE #' #' @param object A \code{lefkoLTRE} object. #' @param ... Other parameters currently not utilized. #' #' @return A list of data frames. In all cases, the first data frame is one #' showing the positive, negative, and total contributions of elements in #' each LTRE contribution matrix. If not a SNA-LTRE, then there are an #' additional two (if deterministic) or four (if stochastic) data frames. If #' deterministic, then \code{hist_det} is a data frame showing the summed LTRE #' contributions for all 16 kinds of historical transition per matrix, with each #' column corresponding to each A matrix in order, followed by all summed #' positive and all summed negative contributions. Object \code{ahist_det} is a #' data frame showing the summed LTRE contributions for all four kinds of #' ahistorical transition per matrix, with order as before, followed by summed #' positive and summed negative contributions. If stochastic, then #' \code{hist_mean} and \code{hist_sd} are the summed LTRE contributions for the #' mean vital rates and variability in vital rates, respectively, according to #' all 16 historical transition types, followed by summed positive and negative #' contributions, and \code{ahist_mean} and \code{ahist_sd} are the equivalent #' ahistorical versions. The output for the SNA-LTRE also includes the #' logs of the deterministic lambda estimated through function \code{ltre3()}. #' #' @examples #' data(lathyrus) #' #' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0) #' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm") #' repvector <- c(0, 0, 0, 0, 0, 1, 0) #' obsvector <- c(0, 1, 1, 1, 1, 1, 0) #' matvector <- c(0, 0, 1, 1, 1, 1, 1) #' immvector <- c(1, 1, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5) #' #' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' immstatus = immvector, indataset = indataset, binhalfwidth = binvec, #' propstatus = propvector) #' #' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988, #' patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9, #' juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88", #' fecacol = "Intactseed88", deadacol = "Dead1988", #' nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea", #' censorcol = "Missing1988", censorkeep = NA, censor = TRUE) #' #' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl", "mat"), #' stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep", "Sdl"), #' stage1 = c("Sd", "rep", "Sd", "rep", "npr", "npr", "Sd"), #' eststage3 = c(NA, NA, NA, NA, NA, NA, "mat"), #' eststage2 = c(NA, NA, NA, NA, NA, NA, "Sdl"), #' eststage1 = c(NA, NA, NA, NA, NA, NA, "NotAlive"), #' givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA, NA), #' multiplier = c(NA, NA, NA, NA, 0.345, 0.054, NA), #' type = c(1, 1, 1, 1, 3, 3, 1), type_t12 = c(1, 2, 1, 2, 1, 1, 1), #' stageframe = lathframe, historical = TRUE) #' #' lathsupp2 <- supplemental(stage3 = c("Sd", "Sdl", "Sd", "Sdl"), #' stage2 = c("Sd", "Sd", "rep", "rep"), #' givenrate = c(0.345, 0.054, NA, NA), #' multiplier = c(NA, NA, 0.345, 0.054), #' type = c(1, 1, 3, 3), stageframe = lathframe, historical = FALSE) #' #' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, year = "all", #' stages = c("stage3", "stage2", "stage1"), supplement = lathsupp3, #' yearcol = "year2", indivcol = "individ") #' #' ehrlen2 <- rlefko2(data = lathvert, stageframe = lathframe, year = "all", #' stages = c("stage3", "stage2"), supplement = lathsupp2, #' yearcol = "year2", indivcol = "individ") #' #' ehrlen3ltre <- ltre3(ehrlen3) #' summary(ehrlen3ltre) #' #' @export summary.lefkoLTRE <- function(object, ...) { trialguy1 <- trialguy2 <- NULL hist1 <- hist2 <- ahist1 <- ahist2 <- NULL sparse_input <- FALSE if (is.element("cont_cv", names(object))) { ltretype <- 3 # sna-LTRE } else if (is.element("cont_sd", names(object))) { ltretype <- 2 # Stochastic LTRE } else if (is.element("cont_mean", names(object))) { ltretype <- 1 # Deterministic LTRE } else { stop("Input object is of unrecognized type. Use only lefkoLTRE obects with this function.", call. = FALSE) } numstages <- dim(object$cont_mean[[1]])[1]
if (is(object$cont_mean[[1]], "dgCMatrix")) sparse_input <- TRUE if(!all(is.na(object$hstages))) {
if(numstages == dim(object$hstages)[1]) { historical <- TRUE } else { historical <- FALSE } } else { historical <- FALSE } if (!historical) { if (!all(is.null(object$agestages))) {
if (!all(is.na(object$agestages))) { if (is.element("stage_id", names(object$agestages))) {
new_ahstages_list <- apply(as.matrix(c(1:length(object$agestages$stage_id))),
1, function(X) {
return(object$ahstages[which(object$ahstages$stage_id == object$agestages$stage_id[X]),]) }) new_ahstages <- do.call("rbind.data.frame", new_ahstages_list) indices <- .bambi2(new_ahstages) } else { indices <- .bambi2(object$ahstages)
}
} else {
indices <- .bambi2(object$ahstages) } } else { if (all(is.null(object$ahstages))) {
mat_size <- dim(object$ah_elasmats[[1]])[1] new_sf <- sf_skeleton(mat_size, standard = FALSE) indices <- .bambi2(new_sf) } else { indices <- .bambi2(object$ahstages)
}
}
} else {
indices <- .bambi3(object$ahstages, object$hstages)
}

used_iterations <- length(object$cont_mean) # General summary for all types of LTRE general_df <- .demolition4(object) # Additional summaries for LTRE and sLTRE for (i in c(1:used_iterations)) { if (!sparse_input) { trialguy1 <- .demolition3(object$cont_mean[[i]], indices)
} else {
trialguy1 <- .demolition3sp(object$cont_mean[[i]], indices) } if (ltretype == 2) { if (!sparse_input) { trialguy2 <- .demolition3(object$cont_sd[[i]], indices)
} else {
trialguy2 <- .demolition3sp(object$cont_sd[[i]], indices) } } if (i == 1) { hist1 <- trialguy1$hist
ahist1 <- trialguy1$ahist if (historical) names(hist1)[2] <- "matrix1" if (historical) names(hist1)[3] <- "matrix1_pos" if (historical) names(hist1)[4] <- "matrix1_neg" names(ahist1)[2] <- "matrix1" names(ahist1)[3] <- "matrix1_pos" names(ahist1)[4] <- "matrix1_neg" if (ltretype == 2) { hist2 <- trialguy2$hist
ahist2 <- trialguy2$ahist if (historical) names(hist2)[2] <- "matrix1" if (historical) names(hist2)[3] <- "matrix1_pos" if (historical) names(hist2)[4] <- "matrix1_neg" names(ahist2)[2] <- "matrix1" names(ahist2)[3] <- "matrix1_pos" names(ahist2)[4] <- "matrix1_neg" } } else { if (historical) hist1 <- cbind.data.frame(hist1, trialguy1$hist[,c(2:4)])
ahist1 <- cbind.data.frame(ahist1, trialguy1$ahist[,c(2:4)]) if (historical) names(hist1)[(((i-1)*3)+2)] <- paste0("matrix", i) if (historical) names(hist1)[(((i-1)*3)+3)] <- paste0("matrix", i, "_pos") if (historical) names(hist1)[(((i-1)*3)+4)] <- paste0("matrix", i, "_neg") names(ahist1)[(((i-1)*3)+2)] <- paste0("matrix", i) names(ahist1)[(((i-1)*3)+3)] <- paste0("matrix", i, "_pos") names(ahist1)[(((i-1)*3)+4)] <- paste0("matrix", i, "_neg") if (ltretype == 2) { if (historical) hist2 <- cbind.data.frame(hist2, trialguy2$hist[,c(2:4)])
ahist2 <- cbind.data.frame(ahist2, trialguy2$ahist[,c(2:4)]) if (historical) names(hist2)[(((i-1)*3)+2)] <- paste0("matrix", i) if (historical) names(hist2)[(((i-1)*3)+3)] <- paste0("matrix", i, "_pos") if (historical) names(hist2)[(((i-1)*3)+4)] <- paste0("matrix", i, "_neg") names(ahist2)[(((i-1)*3)+2)] <- paste0("matrix", i) names(ahist2)[(((i-1)*3)+3)] <- paste0("matrix", i, "_pos") names(ahist2)[(((i-1)*3)+4)] <- paste0("matrix", i, "_neg") } } } output <- if (ltretype == 1) { list(overall = general_df, hist_mean = hist1, ahist_mean = ahist1) } else if (ltretype == 2) { list(overall = general_df, hist_mean = hist1, hist_sd = hist2, ahist_mean = ahist1, ahist_sd = ahist2) } else { list(overall = general_df, hist_mean = hist1, ahist_mean = ahist1, r_values_m = object$r_values_m, r_value_ref = object$r_value_ref) } return (output) } #' Summarize lefkoProj Objects #' #' Function \code{summary.lefkoProj()} summarizes \code{lefkoProj} objects. #' Particularly, it breaks down the data frames provided in the #' \code{projection} element in ways meaningful for those running simulations. #' #' @name summary.lefkoProj #' #' @param object A \code{lefkoProj} object. #' @param threshold A threshold population size to be searched for in #' projections. Defaults to 1. #' @param inf_alive A logical value indicating whether to treat infinitely #' large population size as indicating that the population is still extant. #' If \code{FALSE}, then the population is considered extinct. Defaults to #' \code{TRUE}. #' @param milepost A numeric vector indicating at which points in the projection #' to assess detailed results. Can be input as integer values, in which case #' each number must be between 1 and the total number of occasions projected in #' each projection, or decimals between 0 and 1, which would then be translated #' into the corresponding projection steps of the total. Defaults to #' \code{c(0, 0.25, 0.50, 0.75, 1.00)}. #' @param ext_time A logical value indicating whether to output extinction times #' per population-patch. Defaults to \code{FALSE}. #' @param ... Other parameters currently not utilized. #' #' @return Apart from a statement of the results, this function outputs a list #' with the following elements: #' \item{milepost_sums}{A data frame showing the number of replicates at each #' of the milepost times that is above the threshold population/patch size.} #' \item{extinction_times}{A dataframe showing the numbers of replicates going #' extinct (\code{ext_reps}) and mean extinction time (\code{ext_time}) per #' population-patch. If \code{ext_time = FALSE}, then only outputs \code{NA}.} #' #' @section Notes: #' The \code{inf_alive} and \code{ext_time} options both assess whether #' replicates have reached a value of \code{NaN} or \code{Inf}. If #' \code{inf_alive = TRUE} or \code{ext_time = TRUE} and one of these values is #' found, then the replicate is counted in the \code{milepost_sums} object if #' the last numeric value in the replicate is above the \code{threshold} value, #' and is counted as extant and not extinct if the last numeric value in the #' replicate is above the extinction threshold of a single individual. #' #' Extinction time is calculated on the basis of whether the replicate ever #' falls below a single individual. A replicate with a positive population size #' below 0.0 that manages to rise above 1.0 individual is still considered to #' have gone extinct the first time it crossed below 1.0. #' #' If the input \code{lefkoProj} object is a mixture of two or more other #' \code{lefkoProj} objects, then mileposts will be given relative to the #' maximum number of time steps noted. #' #' @examples #' # Lathyrus example #' data(lathyrus) #' #' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0) #' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm") #' repvector <- c(0, 0, 0, 0, 0, 1, 0) #' obsvector <- c(0, 1, 1, 1, 1, 1, 0) #' matvector <- c(0, 0, 1, 1, 1, 1, 1) #' immvector <- c(1, 1, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5) #' #' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' immstatus = immvector, indataset = indataset, binhalfwidth = binvec, #' propstatus = propvector) #' #' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988, #' patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9, #' juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88", #' fecacol = "Intactseed88", deadacol = "Dead1988", #' nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea", #' censorcol = "Missing1988", censorkeep = NA, censor = TRUE) #' #' lathrepm <- matrix(0, 7, 7) #' lathrepm[1, 6] <- 0.345 #' lathrepm[2, 6] <- 0.054 #' #' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl"), #' stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep"), #' stage1 = c("Sd", "rep", "Sd", "rep", "all", "all"), #' givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA), #' multiplier = c(NA, NA, NA, NA, 0.345, 0.054), #' type = c(1, 1, 1, 1, 3, 3), type_t12 = c(1, 2, 1, 2, 1, 1), #' stageframe = lathframe, historical = TRUE) #' #' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe, #' year = c(1989, 1990), stages = c("stage3", "stage2", "stage1"), #' repmatrix = lathrepm, supplement = lathsupp3, yearcol = "year2", #' indivcol = "individ") #' #' lathproj <- projection3(ehrlen3, nreps = 5, stochastic = TRUE) #' summary(lathproj) #' #' # Cypripedium example #' data(cypdata) #' #' sizevector <- c(0, 0, 0, 0, 0, 0, 1, 2.5, 4.5, 8, 17.5) #' stagevector <- c("SD", "P1", "P2", "P3", "SL", "D", "XSm", "Sm", "Md", "Lg", #' "XLg") #' repvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1) #' obsvector <- c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1) #' matvector <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1) #' immvector <- c(0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0) #' propvector <- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #' indataset <- c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1) #' binvec <- c(0, 0, 0, 0, 0, 0.5, 0.5, 1, 1, 2.5, 7) #' #' cypframe_raw <- sf_create(sizes = sizevector, stagenames = stagevector, #' repstatus = repvector, obsstatus = obsvector, matstatus = matvector, #' propstatus = propvector, immstatus = immvector, indataset = indataset, #' binhalfwidth = binvec) #' #' cypraw_v1 <- verticalize3(data = cypdata, noyears = 6, firstyear = 2004, #' patchidcol = "patch", individcol = "plantid", blocksize = 4, #' sizeacol = "Inf2.04", sizebcol = "Inf.04", sizeccol = "Veg.04", #' repstracol = "Inf.04", repstrbcol = "Inf2.04", fecacol = "Pod.04", #' stageassign = cypframe_raw, stagesize = "sizeadded", NAas0 = TRUE, #' NRasRep = TRUE) #' #' cypsupp3r <- supplemental(stage3 = c("SD", "SD", "P1", "P1", "P2", "P3", "SL", #' "D", "XSm", "Sm", "D", "XSm", "Sm", "mat", "mat", "mat", "SD", "P1"), #' stage2 = c("SD", "SD", "SD", "SD", "P1", "P2", "P3", "SL", "SL", "SL", "SL", #' "SL", "SL", "D", "XSm", "Sm", "rep", "rep"), #' stage1 = c("SD", "rep", "SD", "rep", "SD", "P1", "P2", "P3", "P3", "P3", #' "SL", "SL", "SL", "SL", "SL", "SL", "mat", "mat"), #' eststage3 = c(NA, NA, NA, NA, NA, NA, NA, "D", "XSm", "Sm", "D", "XSm", "Sm", #' "mat", "mat", "mat", NA, NA), #' eststage2 = c(NA, NA, NA, NA, NA, NA, NA, "XSm", "XSm", "XSm", "XSm", "XSm", #' "XSm", "D", "XSm", "Sm", NA, NA), #' eststage1 = c(NA, NA, NA, NA, NA, NA, NA, "XSm", "XSm", "XSm", "XSm", "XSm", #' "XSm", "XSm", "XSm", "XSm", NA, NA), #' givenrate = c(0.1, 0.1, 0.2, 0.2, 0.2, 0.2, 0.25, 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, 0.5, 0.5), #' type = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3), #' type_t12 = c(1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), #' stageframe = cypframe_raw, historical = TRUE) #' #' cypmatrix3r <- rlefko3(data = cypraw_v1, stageframe = cypframe_raw, #' year = "all", patch = "all", stages = c("stage3", "stage2", "stage1"), #' size = c("size3added", "size2added", "size1added"), #' supplement = cypsupp3r, yearcol = "year2", #' patchcol = "patchid", indivcol = "individ") #' #' cypstoch <- projection3(cypmatrix3r, nreps = 5, stochastic = TRUE) #' summary(cypstoch, ext_time = TRUE) #' #' @export summary.lefkoProj <- function(object, threshold = 1, inf_alive = TRUE, milepost = c(0, 0.25, 0.50, 0.75, 1.00), ext_time = FALSE, ...) { num_reps_vec <- num_times_vec <- NULL appended <- FALSE max_times <- max_reps <- 1L ave_times <- ave_reps <- 1.0 poppatches <- length(object$labels[,1])

if (is.matrix(object$control)) { unique_indices <- unique(object$control[,1])
dup_counts <- apply(as.matrix(unique_indices), 1, function(X) {
found <- length(which(object$control[,1] == X)) return(found) }) appended <- TRUE max_reps <- max(dup_counts, na.rm = TRUE) ave_reps <- mean(dup_counts, na.rm = TRUE) max_times <- max(object$control[,3], na.rm = TRUE)
total_reps_across <- sum(object$control[,2], na.rm = TRUE) times_contribs <- apply(as.matrix(c(1:(dim(object$control)[1]))), 1, function(X) {
mean_div <- object$control[X, 2] / total_reps_across current_contrib <- object$control[X, 3] * mean_div

return(current_contrib)
})
ave_times <- sum(times_contribs)

num_reps_vec <- apply(as.matrix(c(1:(dim(object$labels)[1]))), 1, function(X) { found_reps <- sum(object$control[which(object$control[,1] == X), 2], na.rm = TRUE) }) num_times_vec <- apply(as.matrix(c(1:(dim(object$labels)[1]))), 1, function(X) {
found_times <- max(object$control[which(object$control[,1] == X), 3], na.rm = TRUE)
})
} else {
max_reps <- object$control[1] max_times <- ave_times <- object$control[2]

num_reps_vec <- rep(max_reps, poppatches)
num_times_vec <- rep(max_times, poppatches)
}

if (any(milepost < 0)) {
stop("Option milepost may not take negative values.", call. = FALSE)
}
if (any(milepost > max_times)) {
stop("Option milepost may not take values higher than the number of actual
number of projected occasions.", call. = FALSE)
}

if (inf_alive | ext_time) {
for (i in c(1:poppatches)) {
for (j in c(1:num_reps_vec[i])) {
for (k in c(1:(num_times_vec[i] + 1))) {
if ((is.nan(object$pop_size[[i]][j, k]) | is.infinite(object$pop_size[[i]][j, k])) & k > 1) {
object$pop_size[[i]][j, k] <- object$pop_size[[i]][j, (k - 1)]   #. max_found
}
}
}
}
}

if (ext_time) {
the_numbers <- apply(as.matrix(c(1:poppatches)), 1, function(X) {
freemasonry <- apply(as.matrix(c(1:num_reps_vec[X])), 1, function(Y) {
ext_points <- which(object$pop_size[[X]][Y,] < 1) if (length(ext_points) > 0) return (min(ext_points)) else return (NA) }) ext_varmints <- length(which(!is.na(freemasonry) & !is.nan(freemasonry))) if (ext_varmints > 0) { ext_time <- mean(freemasonry, na.rm = TRUE) } else { ext_time <- NA } return (c(ext_varmints, ext_time)) }) the_numbers <- t(the_numbers) the_numbers <- as.data.frame(the_numbers) colnames(the_numbers) <- c("ext_reps", "ext_time") if (dim(object$labels)[1] > 1) {
row_labels <- apply(object$labels, 1, function(X) { paste(X[1], X[2]) }) rownames(the_numbers) <- row_labels } } else { the_numbers <- NA } for (i in c(1:poppatches)) { if (any(milepost > num_times_vec[i])) { stop("Entered milepost values are outside the allowable range.", call. = FALSE) } } milepost_sums <- apply(as.matrix(c(1:poppatches)), 1, function (X) { used_milepost <- milepost if (all(milepost >= 0) & all(milepost <= 1)) { used_milepost <- floor(used_milepost * num_times_vec[X]) + 1 } else if (any(milepost == 0)) { used_milepost <- used_milepost + 1 } if (num_reps_vec[X] > 1) { phew <- apply(as.matrix(object$pop_size[[X]][,used_milepost]), 2, function(Y) {
above_vector <- which(as.vector(Y) >= threshold)

return(length(above_vector))
})
return(phew)
} else {
phew <- apply(as.matrix(object$pop_size[[X]][,used_milepost]), 1, function(Y) { above_vector <- which(as.vector(Y) >= threshold) return(length(above_vector)) }) return(phew) } }) if (is.matrix(milepost_sums)) { rownames(milepost_sums) <- milepost col_labels <- apply(object$labels, 1, function(X) {
paste(X[1], X[2])
})
colnames(milepost_sums) <- col_labels
} else {
rownames(milepost_sums) <- milepost
}

writeLines(paste0("\nThe input lefkoProj object covers ", poppatches,
" population-patches."), con = stdout())
if (appended) {
writeLines(paste0("It is an appended projection, including an average and maximum of ",
format(ave_times, digits = 5), " and ", max_times, " steps per "))
writeLines(paste0("replicate, and an average and maximum of ",
format(ave_reps, digits = 5), " and ", max_reps, " replicates."))
} else {
writeLines(paste0("It is a single projection including ", max_times,
" projected steps per replicate, and ", max_reps, " replicates, respectively."),
con = stdout())
}
writeLines(paste0("The number of replicates with population size above the threshold size of ",
threshold, " is as in"), con = stdout())
writeLines(paste0("the following matrix, with pop-patches given by column and milepost times given by row: \n"),
con = stdout())

output <- list(milepost_sums = milepost_sums, extinction_times = the_numbers)

return (output)
}

#' Plot Projection Simulations
#'
#' Function \code{plot.lefkoProj()} produces plots of \code{lefkoProj} objects.
#' Acts as a convenient wrapper for the \code{plot.default()} function.
#'
#' @name plot.lefkoProj
#'
#' @param x A \code{lefkoProj} object.
#' @param variable The focus variable of the plot to produce. Defaults to
#' \code{"popsize"}, which produces line plots of the \code{popsize} element in
#' object \code{x}.
#' @param style A string denoting ther kind of plot to produce. Currently
#' limited to \code{"timeseries"}, which shows \code{variable} against time on
#' the x axis. Other choices include \code{"statespace"}, which plots
#' \code{variable} at one time on the x axis against the same variable in the
#' next time on the y axis.
#' @param repl The replicate to plot. Defaults to \code{"all"}, in which case
#' all replicates are plotted.
#' @param patch The patch to plot, as labeled in the \code{labels} element in
#' object \code{x}. Defaults to \code{"pop"}, in which case only the final
#' population-level projection is plotted. Can also be set to \code{"all"}, in
#' which case projections for all patches and population in the \code{labels}
#' element are plotted.
#' @param auto_ylim A logical value indicating whether the maximum of the y axis
#' should be determined automatically. Defaults to \code{TRUE}, but reverts to
#' \code{FALSE} if any setting for \code{ylim} is given.
#' @param auto_col A logical value indicating whether to shift the color of
#' lines associated with each patch automatically. Defaults to \code{TRUE}, but
#' reverts to \code{FALSE} if any setting for \code{col} is given.
#' @param auto_lty A logical value indicating whether to shift the line type
#' associated with each replicate automatically. Defaults to \code{TRUE}, but
#' reverts to \code{FALSE} if any setting for \code{lty} is given.
#' @param auto_title A logical value indicating whether to add a title to each
#' plot. The plot is composed of the concatenated population and patch names.
#' Defaults to \code{FALSE}.
#' @param ... Other parameters used by functions \code{plot.default()} and
#' \code{lines()}.
#'
#' @return A plot of the results of a \code{\link{projection3}()} run.
#'
#' @section Notes:
#' Output plots are currently limited to time series and state space plots of
#' population size.
#'
#' The default settings will preferentially plot any projections marked as
#' \code{0} in the \code{patch} portion of the \code{labels} element of the
#' input MPM. This can produce confusing results if a mean MPM resulting from
#' the \code{lmean()} function is used as input and the \code{add_mean} setting
#' is set to the default, which is \code{TRUE}.
#'
#' @examples
#' data(lathyrus)
#'
#' sizevector <- c(0, 100, 13, 127, 3730, 3800, 0)
#' stagevector <- c("Sd", "Sdl", "VSm", "Sm", "VLa", "Flo", "Dorm")
#' repvector <- c(0, 0, 0, 0, 0, 1, 0)
#' obsvector <- c(0, 1, 1, 1, 1, 1, 0)
#' matvector <- c(0, 0, 1, 1, 1, 1, 1)
#' immvector <- c(1, 1, 0, 0, 0, 0, 0)
#' propvector <- c(1, 0, 0, 0, 0, 0, 0)
#' indataset <- c(0, 1, 1, 1, 1, 1, 1)
#' binvec <- c(0, 100, 11, 103, 3500, 3800, 0.5)
#'
#' lathframe <- sf_create(sizes = sizevector, stagenames = stagevector,
#'   repstatus = repvector, obsstatus = obsvector, matstatus = matvector,
#'   immstatus = immvector, indataset = indataset, binhalfwidth = binvec,
#'   propstatus = propvector)
#'
#' lathvert <- verticalize3(lathyrus, noyears = 4, firstyear = 1988,
#'   patchidcol = "SUBPLOT", individcol = "GENET", blocksize = 9,
#'   juvcol = "Seedling1988", sizeacol = "Volume88", repstracol = "FCODE88",
#'   nonobsacol = "Dormant1988", stageassign = lathframe, stagesize = "sizea",
#'   censorcol = "Missing1988", censorkeep = NA, censor = TRUE)
#'
#' lathrepm <- matrix(0, 7, 7)
#' lathrepm[1, 6] <- 0.345
#' lathrepm[2, 6] <- 0.054
#'
#' lathsupp3 <- supplemental(stage3 = c("Sd", "Sd", "Sdl", "Sdl", "Sd", "Sdl"),
#'   stage2 = c("Sd", "Sd", "Sd", "Sd", "rep", "rep"),
#'   stage1 = c("Sd", "rep", "Sd", "rep", "all", "all"),
#'   givenrate = c(0.345, 0.345, 0.054, 0.054, NA, NA),
#'   multiplier = c(NA, NA, NA, NA, 0.345, 0.054),
#'   type = c(1, 1, 1, 1, 3, 3), type_t12 = c(1, 2, 1, 2, 1, 1),
#'   stageframe = lathframe, historical = TRUE)
#'
#' ehrlen3 <- rlefko3(data = lathvert, stageframe = lathframe,
#'   year = c(1989, 1990), stages = c("stage3", "stage2", "stage1"),
#'   repmatrix = lathrepm, supplement = lathsupp3, yearcol = "year2",
#'   indivcol = "individ")
#'
#' lathproj <- projection3(ehrlen3, nreps = 5, stochastic = TRUE)
#' plot(lathproj)
#'
#' @export
plot.lefkoProj <- function(x, variable = "popsize", style = "time",
repl = "all", patch = "pop", auto_ylim = TRUE, auto_col = TRUE,
auto_lty = TRUE, auto_title = FALSE, ...) {

appended <- FALSE

further_args <- list(...)

if (length(further_args) == 0) further_args <- list()

if (!is.element("type", names(further_args))) {
further_args$type <- "l" } if (is.element("col", names(further_args))) { auto_col <- FALSE } if (is.element("lty", names(further_args))) { auto_lty <- FALSE } if (is.element("ylim", names(further_args))) { auto_ylim <- FALSE } if (is.element("main", names(further_args))) { auto_title <- FALSE } basal_args <- further_args actual_patches <- c(1:length(x$labels$patch)) if (is.matrix(x$control)) {
appended <- TRUE
actual_replicates <- apply(as.matrix(actual_patches), 1, function(X) {
num_reps <- sum(x$control[which(x$control[, 1] == X), 2])
})
} else {
actual_replicates <- rep(x$control[1], length(x$labels$patch)) } if (length(grep("pop", variable)) > 0 | length(grep("size", variable)) > 0) { variable <- "popsize" } if (length(grep("stat", style)) > 0) { style <- "statespace" if (!is.element("ylab", names(further_args))) { further_args$ylab <- "State in time t+1"
}
if (!is.element("xlab", names(further_args))) {
further_args$xlab <- "State in time t" } } else if (length(grep("tim", style)) > 0) { style <- "timeseries" if (!is.element("ylab", names(further_args))) { further_args$ylab <- "Population size"
}
if (!is.element("xlab", names(further_args))) {
further_args$xlab <- "Time" } } if (all(is.na(patch))) { patch <- actual_patches } else if (is.character(patch)) { if (any(grep("al", patch))) { patch <- actual_patches } else if (length(patch) == 1) { if (grep("po", patch)) { if (length(actual_patches) == 1) { patch <- actual_patches } else { patch <- which(x$labels$patch == 0) } } } else { patch <- as.numeric(patch) if (any(is.na(patch))) { stop("Setting patch not understood. Please do not combine text with numbers.", call. = FALSE) } if (!all(is.element(patch, actual_patches))) { stop("Setting patch not understood. Please check the labels element of the input lefkoProj object.", call. = FALSE) } } } else if (is.numeric(patch)) { if (any(!is.element(patch, actual_patches))) { stop("Setting patch not understood. Please check the labels element of your input lefkoProj object.", call. = FALSE) } } if (is.character(repl)) { if (any(grep("al", repl))) { repl <- actual_replicates } else { stop("Setting repl not understood.", call. = FALSE) } } if (variable != "popsize") { stop("Function plot.lefkoProj() does not currently handle plots of elements other than popsize.", call. = FALSE) } used_col <- 1 for (i in patch) { if (auto_title) { used_string <- paste("pop", x$labels[i, 1], "patch", x$labels[i, 2]) further_args$main <- used_string
}

used_lty <- 1

if (auto_ylim) {
further_args$ylim <- c(0, max(x$pop_size[[i]], na.rm = TRUE))
}
if (auto_col) {
further_args$col <- used_col basal_args$col <- used_col
}
if (auto_lty) {
further_args$lty <- used_lty basal_args$lty <- used_lty
}

if (style == "timeseries") {
c_xy <- xy.coords(x = c(1:length(x$pop_size[[i]][1,])), y = x$pop_size[[i]][1,])
further_args$x <- c_xy } else if (style == "statespace") { c_xy <- xy.coords(x = x$pop_size[[i]][1,1:(dim(x$pop_size[[i]])[2] - 1)], y = x$pop_size[[i]][1,2:dim(x$pop_size[[i]])[2]]) further_args$x <- c_xy
}

do.call("plot.default", further_args)

if (repl[i] > 1) {
for (j in c(2:actual_replicates[i])) {
if (style == "timeseries") {
basal_args$x <- c(1:length(x$pop_size[[i]][j,]))
basal_args$y <- x$pop_size[[i]][j,]

} else if (style == "statespace") {
basal_args$y <- x$pop_size[[i]][j,2:dim(x$pop_size[[i]])[2]] basal_args$x <- x$pop_size[[i]][j,1:(dim(x$pop_size[[i]])[2] - 1)]

}
do.call("lines", basal_args)
}
used_lty <- used_lty + 1;
}
used_col <- used_col + 1;
if (used_col > length(palette())) used_col <- 1;
}
}


## Try the lefko3 package in your browser

Any scripts or data that you put into this service are public.

lefko3 documentation built on May 29, 2024, 6:39 a.m.