From fec78622db551acddabb8d0cdab36422419a17a3 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Mon, 3 Aug 2026 14:10:14 +0200 Subject: [PATCH 01/36] Added print and plot methods for ipmr uncertainty objects --- R/generics.R | 120 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 120 insertions(+) diff --git a/R/generics.R b/R/generics.R index 92311ad..553885b 100644 --- a/R/generics.R +++ b/R/generics.R @@ -1018,6 +1018,49 @@ print.ipmr_vital_rate_funs <- function(x, ...) { } +#' Print method for ipmr uncertainty objects +#' +#' @param x An object of class `ipmr_uncertainty` +#' @param ... Unused +#' +#' @return \code{x} invisibly. +#' +#' @export +print.ipmr_uncertainty <- function(x, ...) { + + cat("\n") + cat("====================================\n") + cat(" ipmr_uncertainty object\n") + cat("====================================\n\n") + + ## ---- Basic info + cat("Number of parameter samples:", nrow(x$lambdas), "\n") + cat("Number of parameters:", nrow(x$params_uncert), "\n\n") + + ## ---- Lambda uncertainty + cat("Lambda variance:\n") + print(x$mod_uncert) + cat("\n") + + ## ---- Parameter-level summary + cat("Top contributing parameters (by variance contribution):\n") + top <- x$params_uncert[order(-x$params_uncert$variance_prop), ] + print(utils::head(top[, c("parameter", "vital_rate", "variance_prop")], 5)) + cat("\n") + + ## ---- Vital-rate summary + if (!is.null(x$vr_uncert)) { + cat("Vital-rate contributions:\n") + print(x$vr_uncert) + cat("\n") + } + + cat("Use `plot()` to visualize uncertainty decomposition.\n") + cat("====================================\n") + + invisible(x) +} + # Lambda------------ #' @title Compute the per-capita growth rate for an IPM object #' @rdname lambda @@ -1880,6 +1923,83 @@ plot.general_di_det_ipm <- function(x = NULL, y = NULL, } +#' Plot method for ipmr uncertainty objects +#' +#' @param x An object of class `ipmr_uncertainty` +#' @param type Type of plot: +#' \itemize{ +#' \item "param" - parameter-level 3-panel plot +#' \item "vr" - vital-rate aggregated contributions +#' \item "stacked" - stacked contributions vs lambda variance +#' \item "all" - returns all plots as a list +#' } +#' @param vr_colors A named character vector of colors used to represent +#' vital rate categories in plots. Names must match the values in the +#' \code{vital_rate} column of the uncertainty output. If \code{NULL}, +#' a default colorblind-friendly palette is used. +#' @param ... Additional arguments passed to plotting functions +#' +#' @return A \code{ggplot} object or a list of plots. +#' +#' @export +plot.ipmr_uncertainty <- function(x, + type = c("all", "param", "vr", "stacked"), + vr_colors = NULL, + ...) { + + type <- match.arg(type) + + if (!inherits(x, "ipmr_uncertainty")) { + stop("`x` must be an object of class 'ipmr_uncertainty'.") + } + + default_cols <- c( + "#D0873C", # soft orange + "#8FB339", # olive green + "#7A5195", # muted purple + "#88CCEE", # light blue + "#332288", # deep blue + "#44AA99", # teal + "#882255", # dark red + "#DDCC77", # sand + "#C94C4C", # muted red + "#6C757D" # grey + ) + + vr_levels <- unique(x$params_uncert$vital_rate) + + if (is.null(vr_colors)) { + vr_colors <- setNames( + rep(default_cols, length.out = length(vr_levels)), + vr_levels + ) + } else { + if (is.null(names(vr_colors))) { + stop("`vr_colors` must be a named vector with names matching vital rates.") + } + + missing_cols <- setdiff(vr_levels, names(vr_colors)) + if (length(missing_cols) > 0) { + stop("Missing colors for vital rates: ", + paste(missing_cols, collapse = ", ")) + } + } + + # Route to internal plotting helpers + plots <- list( + param = .plot_uncert_param(x, vr_colors = vr_colors, ...), + vr = .plot_uncert_vr(x, vr_colors = vr_colors, ...), + stacked = .plot_uncert_stacked(x, vr_colors = vr_colors, ...) + ) + + switch(type, + param = plots$param, + vr = plots$vr, + stacked = plots$stacked, + all = plots + ) +} + # right_ev ---------------- #' @rdname eigenvectors From 8a11e54a30dd94eb65c11422d2dc56c5eb056ec1 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Mon, 3 Aug 2026 14:16:11 +0200 Subject: [PATCH 02/36] Add helper functions for sensitivity and uncertainty --- R/utils-uncertainty.R | 481 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 481 insertions(+) create mode 100644 R/utils-uncertainty.R diff --git a/R/utils-uncertainty.R b/R/utils-uncertainty.R new file mode 100644 index 0000000..eb7c5db --- /dev/null +++ b/R/utils-uncertainty.R @@ -0,0 +1,481 @@ +#' Extract dominant eigenvalue +#' +#' Returns the dominant eigenvalue using the same ordering assumption +#' currently used in sensitivity calculations. +#' +#' @noRd +dominant_eigen <- function(mat) { + eig <- eigen(mat) + Re(eig$values[which.max(Mod(eig$values))]) +} + + +#' Validate inputs for sensitivity and uncertainty functions +#' +#' Returns usable objects to be passed to the function-specific validations. +#' +#' @noRd +validate_ipm_base <- function(ipm) { + + if (missing(ipm) || is.null(ipm)) { + stop("`ipm` must be provided and cannot be NULL.") + } + + if (!is.list(ipm)) { + stop("`ipm` must be a valid ipmr object.") + } + + cls <- strsplit(class(ipm)[1], "_")[[1]] + + if (length(cls) < 3) { + stop("Malformed IPM class structure.") + } + + di_dd <- cls[2] + det_stoch <- cls[3] + + if (is.na(di_dd) || di_dd != "di") { + stop("Only density-independent IPMs are currently supported.") + } + + if (is.na(det_stoch) || det_stoch != "det") { + stop("Only deterministic IPMs are currently supported.") + } + + pars_all <- tryCatch( + parameters(ipm), + error = function(e) { + stop("Could not extract parameters from IPM.") + } + ) + + if (!is.list(pars_all) || is.null(names(pars_all))) { + stop("`parameters(ipm)` must return a named list.") + } + + list( + ipm = ipm, + pars_all = pars_all + ) +} + + +#' Validate inputs to sensitivity function +#' +#' Returns usable objects for the sensitivity function. +#' +#' @noRd +validate_ipm_sensitivity <- function(ipm, pars, kernels, delta) { + + base <- validate_ipm_base(ipm) + + pars_all <- base$pars_all + ipm <- base$ipm + + if (is.null(names(pars_all))) { + stop("Parameters must be named.") + } + + if (is.null(pars)) { + pars <- names(pars_all) + } else { + if (!is.character(pars)) { + stop("`pars` must be a character vector of parameter names.") + } + + if (length(pars) == 0) { + stop("`pars` cannot be empty.") + } + + if (!all(pars %in% names(pars_all))) { + bad <- pars[!pars %in% names(pars_all)] + stop("Unknown parameter(s): ", paste(bad, collapse = ", ")) + } + + pars <- unique(pars) + } + + if (any(pars_all[pars] == 0)) { + warning("Zero-valued parameters detected; elasticity may be undefined.") + } + + ## ---- kernels + + if (missing(kernels) || is.null(kernels)) { + stop("`kernels` must be provided.") + } + + if (!is.character(kernels) || length(kernels) == 0) { + stop("`kernels` must be a non-empty character vector.") + } + + # Check for duplicates + dup_kerns <- unique(kernels[duplicated(kernels)]) + + if (length(dup_kerns) > 0) { + warning( + "Duplicate kernel name(s) detected in `kernels`: ", + paste(dup_kerns, collapse = ", "), + ". This may lead to incorrect kernel construction if not intentional." + ) + } + + available_kernels <- names(ipm$sub_kernels) + + if (!all(kernels %in% available_kernels)) { + bad <- setdiff(kernels, available_kernels) + stop( + "Unknown kernel name(s): ", + paste(bad, collapse = ", "), + ". Available kernels are: ", + paste(available_kernels, collapse = ", ") + ) + } + + n_expected <- length(ipm$sub_kernels) + + if (length(kernels) != n_expected) { + warning( + "`kernels` has length ", length(kernels), + " but IPM contains ", n_expected, " subkernels. ", + "Ensure kernels are provided in correct row-major order." + ) + } + + ## ---- delta + + if (!is.numeric(delta) || length(delta) != 1 || is.na(delta)) { + stop("`delta` must be a single numeric value.") + } + + if (delta <= 0) { + stop("`delta` must be > 0.") + } + + ## return cleaned objects for main function + list( + ipm = ipm, + pars = pars, + pars_all = pars_all, + kernels = kernels, + delta = delta + ) +} + + +#' Create a new ipmr uncertainty object +#' +#' Internal constructor for ipmr uncertainty outputs. +#' +#' @param ipm The original IPM +#' @param lambdas Data frame of lambda values +#' @param mod_uncert Numeric variance in lambda +#' @param params_uncert Data frame of parameter summaries +#' @param vr_uncert Data frame of uncertainty contributions summed by vital rate +#' +#' @noRd +new_ipmr_uncertainty <- function(ipm, lambdas, mod_uncert, params_uncert, vr_uncert) { + + structure( + list( + ipm = ipm, + lambdas = lambdas, + mod_uncert = mod_uncert, + params_uncert = params_uncert, + vr_uncert = vr_uncert + ), + class = "ipmr_uncertainty" + ) +} + +#' Validate inputs to uncertainty function +#' +#' +#' @noRd +validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, + delta, cores) { + + base <- validate_ipm_base(ipm) + + pars_all <- base$pars_all + ipm <- base$ipm + + if (is.null(pars)) { + pars <- names(pars_all) + } + + if (!is.character(pars)) { + stop("`pars` must be a character vector.") + } + + if (!all(pars %in% names(pars_all))) { + stop("Unknown parameter(s): ", + paste(pars[!pars %in% names(pars_all)], collapse = ", ")) + } + + ## ---- samples + if (missing(samples) || !is.data.frame(samples)) { + stop("`samples` must be a data.frame.") + } + + if (nrow(samples) < 1) { + stop("`samples` must have at least one row. We suggest using no fewer than 100 samples.") + } + + missing_in_samples <- setdiff(names(pars_all), colnames(samples)) + if (length(missing_in_samples) > 0) { + stop("Samples missing required parameter columns: ", + paste(missing_in_samples, collapse = ", ")) + } + + extra_in_samples <- setdiff(colnames(samples), names(pars_all)) + if (length(extra_in_samples) > 0) { + warning("Extra columns in `samples` not used: ", + paste(extra_in_samples, collapse = ", ")) + } + + ## ---- kernels + if (missing(kernels) || is.null(kernels)) { + stop("`kernels` must be provided.") + } + + if (!is.character(kernels) || length(kernels) == 0) { + stop("`kernels` must be a non-empty character vector.") + } + + # Check for duplicates + dup_kerns <- unique(kernels[duplicated(kernels)]) + + if (length(dup_kerns) > 0) { + warning( + "Duplicate kernel name(s) detected in `kernels`: ", + paste(dup_kerns, collapse = ", "), + ". This may lead to incorrect kernel construction if not intentional." + ) + } + + available_kernels <- names(ipm$sub_kernels) + + if (!all(kernels %in% available_kernels)) { + bad <- setdiff(kernels, available_kernels) + stop( + "Unknown kernel name(s): ", + paste(bad, collapse = ", "), + ". Available kernels are: ", + paste(available_kernels, collapse = ", ") + ) + } + + n_expected <- length(ipm$sub_kernels) + + if (length(kernels) != n_expected) { + warning( + "`kernels` has length ", length(kernels), + " but IPM contains ", n_expected, " subkernels. ", + "Ensure kernels are provided in correct row-major order." + ) + } + + ## ---- vr_table + if (!identical(vr_table, "template")) { + + if (!is.data.frame(vr_table)) { + stop("`vr_table` must be a data.frame or 'template'.") + } + + if (!all(c("parameter", "vital_rate") %in% colnames(vr_table))) { + stop("`vr_table` must contain columns: parameter, vital_rate") + } + + if (any(duplicated(vr_table$parameter))) { + stop("Each parameter must appear only once in `vr_table`.") + } + + if (!all(pars %in% vr_table$parameter)) { + stop("Missing parameters in `vr_table`.") + } + } + + ## ---- delta + if (!is.numeric(delta) || length(delta) != 1 || is.na(delta) || delta <= 0) { + stop("`delta` must be a single numeric value > 0.") + } + + + ## ---- cores + if (!is.numeric(cores) || length(cores) != 1 || cores < 1) { + stop("`cores` must be a single integer >= 1.") + } + + cores <- as.integer(cores) + + if (cores < 1) { + stop("`cores` must be >= 1.") + } + + max_cores <- parallel::detectCores(logical = TRUE) + if (cores > max_cores) { + stop( + "`cores` (", cores, ") exceeds available logical cores (", + max_cores, ")." + ) + } + + slurm <- Sys.getenv("SLURM_JOB_ID") != "" + pbs <- Sys.getenv("PBS_JOBID") != "" + + if ((slurm || pbs) && cores == max_cores) { + warning( + "Detected HPC scheduler environment (Slurm/PBS). ", + "Using all detected cores may exceed allocated resources." + ) + } + + list( + ipm = ipm, + pars = pars, + samples = samples, + kernels = kernels, + vr_table = vr_table, + delta = delta, + cores = as.integer(cores) + ) +} + +#' Plot three-panel uncertainty plot +#' +#' @noRd +.plot_uncert_param <- function(x, vr_colors, ...) { + + df <- x$params_uncert + + # Order parameters by contribution + df$parameter <- factor( + df$parameter, + levels = rev(df$parameter[order(df$vital_rate)]) + ) + + library(ggplot2) + library(patchwork) + + base_theme <- theme_minimal() + + theme( + plot.title = element_text(size = 10, face = "bold", hjust = 0.5), + plot.margin = margin(5.5, 10, 5.5, 5.5) + ) + + title_wrap <- function(x) str_wrap(x, width = 28) + + p1 <- ggplot(df, aes(x = cv, y = parameter, color = vital_rate)) + + geom_vline(xintercept = 0) + + geom_point(size = 2.5) + + geom_segment(aes(x = 0, xend = cv, yend = parameter), linewidth = 1) + + labs(x = "Coefficient of variation", + y = "Parameter", + title = title_wrap("Parameter uncertainty")) + + scale_color_manual(values = vr_colors, guide = "none") + + base_theme + + p2 <- ggplot(df, aes(x = elasticity, y = parameter, color = vital_rate)) + + geom_vline(xintercept = 0) + + geom_point(size = 2.5) + + geom_segment(aes(x = 0, xend = elasticity, yend = parameter), linewidth = 1) + + labs(x = "Elasticity", + y = NULL, + title = title_wrap("Model sensitivity to parameter")) + + scale_color_manual(values = vr_colors, guide = "none") + + base_theme + + p3 <- ggplot(df, aes(x = variance_prop, y = parameter, color = vital_rate)) + + geom_vline(xintercept = 0) + + geom_point(size = 2.5) + + geom_segment(aes(x = 0, xend = variance_prop, yend = parameter), linewidth = 1) + + labs(x = "Uncertainty contribution", + y = NULL, + title = title_wrap("Parameter contribution to uncertainty in λ"), + color = "Vital rate") + + scale_color_manual(values = vr_colors) + + base_theme + + p1 + p2 + p3 + plot_layout(ncol = 3, axes = 'collect_y') +} + +#' Plot uncertainty aggregated by vital rate +#' +#' @noRd +.plot_uncert_vr <- function(x, vr_colors, ...) { + + df <- x$vr_uncert + + df$vital_rate <- factor( + df$vital_rate, + levels = df$vital_rate[order(df$variance_sum)] + ) + + library(ggplot2) + + ggplot(df, aes(x = variance_sum, y = vital_rate, fill = vital_rate)) + + geom_col() + + labs( + x = "Total uncertainty contribution", + y = "Vital rate" + ) + + scale_fill_manual(values = vr_colors) + + theme_minimal() + + guides(fill = "none") +} + +#' Plot model uncertainty and vital rate uncertainty contributions +#' +#' @noRd +.plot_uncert_stacked <- function(x, vr_colors, ...) { + + df <- x$vr_uncert + + total_var <- x$mod_uncert + + library(ggplot2) + + p <- ggplot(df, aes(x = 1, y = variance_sum, fill = vital_rate)) + + geom_col(width = 0.6) + + + # Outline box for total lambda variance + geom_rect( + aes( + xmin = 0.7, + xmax = 1.3, + ymin = 0, + ymax = total_var, + color = "Total variance in λ" + ), + inherit.aes = FALSE, + fill = NA, + linewidth = 1 + ) + + + labs( + x = NULL, + y = "Variance in λ", + title = "Uncertainty decomposition by vital rate", + fill = "Vital rate" + ) + + + scale_fill_manual(values = vr_colors) + + scale_color_manual( + name = NULL, + values = c("Total variance in λ" = "black") + ) + + + guides( + fill = guide_legend(order = 1), + color = guide_legend(order = 2) + ) + + + theme_minimal() + + theme( + axis.text.x = element_blank(), + axis.ticks.x = element_blank() + ) + + p +} From 77c655cc25358ae978b99b7247117346608e7290 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Mon, 3 Aug 2026 14:18:16 +0200 Subject: [PATCH 03/36] Add sensitivity function --- R/sensitivity.R | 110 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 110 insertions(+) create mode 100644 R/sensitivity.R diff --git a/R/sensitivity.R b/R/sensitivity.R new file mode 100644 index 0000000..5106660 --- /dev/null +++ b/R/sensitivity.R @@ -0,0 +1,110 @@ +#' @title Sensitivity and elasticity analysis of population growth rate +#' +#' @description Computes numerical sensitivities and elasticities of the +#' dominant eigenvalue (\eqn{\lambda}) with respect to model parameters. +#' +#' @param ipm An \code{ipmr} IPM object. +#' @param pars A character vector of parameter names to perturb. Defaults to +#' all parameters in \code{ipm}. +#' @param kernels A character vector specifying the kernel structure in row +#' major order. +#' @param delta A numeric scalar specifying the perturbation size used for +#' numerical differentiation. Default is 1e-4. +#' +#'@details +#' This function uses a central difference approximation to estimate the +#' sensitivity of \eqn{\lambda} to each parameter in \code{pars}. +#' +#' Currently implemented for deterministic ("det"), density-independent +#' ("di") IPMs constructed with \code{ipmr}. Stochastic ("stoch") and +#' density-dependent ("dd") IPMs are not supported. +#' +#' @return A named list with two elements: +#' \itemize{ +#' \item \code{sensitivity}: Numeric vector of sensitivities of +#' \eqn{\lambda} to each parameter. +#' \item \code{elasticity}: Numeric vector of elasticities of +#' \eqn{\lambda} to each parameter. +#' } +#' +#' @export +sensitivity <- function(ipm, pars = NULL, kernels, delta = 1e-4) { + + ## ---- Input checks + + val <- validate_ipm_sensitivity(ipm, pars, kernels, delta) + + ipm <- val$ipm + pars <- val$pars + pars_all <- val$pars_all + kernels <- val$kernels + delta <- val$delta + + ## ---- Setup + + npar <- length(pars) + spar <- numeric(npar) + vpar <- pars_all + + kernel <- make_iter_kernel(ipm, mega_mat = kernels) + + if (!is.list(kernel) || is.null(kernel[[1]])) { + stop("Kernel construction failed.") + } + + eig <- tryCatch( + eigen(kernel[[1]]), + error = function(e) stop("Eigen decomposition failed for base kernel.") + ) + + lambda_all <- dominant_eigen(kernel[[1]]) + + ## ---- Loop over each parameter + + for(i in seq_len(npar)) { + + par_now <- pars[i] + + # ---- DOWN perturbation + vpar <- pars_all + vpar[[par_now]] <- vpar[[par_now]] - delta + + ipm_down <- ipm$proto_ipm + parameters(ipm_down) <- vpar + + ker_down <- ipm_down %>% + make_ipm(iterate = FALSE) %>% + make_iter_kernel(mega_mat = kernels) + + lambda_down <- dominant_eigen(ker_down[[1]]) + + # ---- UP perturbation + vpar <- pars_all + vpar[[par_now]] <- vpar[[par_now]] + delta + + ipm_up <- ipm$proto_ipm + parameters(ipm_up) <- vpar + + ker_up <- ipm_up %>% + make_ipm(iterate = FALSE) %>% + make_iter_kernel(mega_mat = kernels) + + lambda_up <- dominant_eigen(ker_up[[1]]) + + # ---- Central difference + spar[i] <- (lambda_up - lambda_down) / (2 * delta) + } + + ## ---- Output + + sens <- spar + elas <- spar * abs(as.numeric(pars_all[pars]) / lambda_all) + + names(sens) <- pars + names(elas) <- pars + + return(list( + sensitivity = sens, + elasticity = elas + )) +} From 04edd3eb0ebbc1dba9f50664b8b67d67e1e7606d Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Mon, 3 Aug 2026 14:21:36 +0200 Subject: [PATCH 04/36] Add uncertainty function --- R/uncertainty.R | 240 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 240 insertions(+) create mode 100644 R/uncertainty.R diff --git a/R/uncertainty.R b/R/uncertainty.R new file mode 100644 index 0000000..4663c86 --- /dev/null +++ b/R/uncertainty.R @@ -0,0 +1,240 @@ +#' @title Uncertainty analysis of IPMs +#' +#' @description Performs uncertainty analysis of IPMs. +#' +#' @param ipm An \code{ipmr} IPM object. +#' @param pars A character vector of parameter names on which to perform the +#' uncertainty analysis. Defaults to all parameters in \code{ipm}. +#' @param samples A data.frame containing sampled parameter values. +#' @param kernels A character vector specifying the kernel structure in row +#' major order. +#' @param vr_table A data.frame with columns 'parameter' and 'vital_rate' +#' mapping each parameter to a user-defined vital rate category. Every +#' parameter in \code{pars} must be present exactly once in this table. +#' @param delta A numeric scalar specifying the perturbation size used for +#' numerical differentiation when calculating sensitivities. Default is 1e-4. +#' @param cores An integer specifying the number of cores between which the +#' computations will be distributed. Defaults to 1. +#' +#'@details +#' Currently implemented for deterministic ("det"), density-independent +#' ("di") IPMs constructed with \code{ipmr}. Stochastic ("stoch") and +#' density-dependent ("dd") IPMs are not supported. +#' +#' Sampled parameter values may either be drawn from the joint distribution or +#' iteratively refit to bootstrapped data. The \code{samples} dataframe must +#' contain columns for all parameters needed to build the IPM, including those +#' which do not vary between samples. +#' +#' Because parameters may appear in multiple kernels or represent multiple +#' biological processes, \code{ipmr} does not attempt to infer vital rate +#' categories automatically. Users must explicitly define this categorization. +#' A template for \code{vr_table} can be generated by calling the function with +#' \code{vr_table = "template"}. +#' +#' Sensitivities within each sampled model are computed using a central +#' difference approximation (see \code{\link{sensitivity}}). The \code{delta} +#' argument controls the perturbation size used in this approximation. +#' +#' On systems with multiple cores, computations can be parallelized using the +#' \code{cores} argument. The function does not automatically infer scheduler +#' resource limits on high-performance computing (HPC) systems. +#' +#' Users running analyses on shared clusters (e.g. via Slurm, PBS, or similar +#' job schedulers) should ensure that the number of cores requested does not +#' exceed the allocation specified by the job submission script. Using all +#' available system cores may lead to oversubscription and reduced performance +#' of concurrent jobs. +#' +#' @return +#' +#' @examples +#' # Generate a template for defining vital rates +#' vr <- uncertainty(ipm, vr_table = "template") +#' +#' # Fill in the vital_rate column manually, then run: +#' # vr$vital_rate <- c("survival", "growth", "reproduction", ...) +#' +#' # Run uncertainty analysis +#' # out <- uncertainty(ipm, pars, samples, kernels, vr_table = vr) +#' +#' @export +uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, + delta = 1e-4, cores = 1) { + + if (isTRUE(vr_table == "template")) { + if (is.null(pars)) { + pars <- names(parameters(ipm)) + } + return(data.frame( + parameter = pars, + vital_rate = NA_character_ + )) + } + + ## ---- Input checks + val <- validate_ipm_uncertainty(ipm, pars, samples, kernels, + vr_table, delta, cores) + + ipm <- val$ipm + pars <- val$pars + samples <- val$samples + kernels <- val$kernels + vr_table <- val$vr_table + cores <- val$cores + + npar <- length(pars) + nsam <- nrow(samples) + + ## ---- Parallel setup + base_ipm <- ipm$proto_ipm + cluster <- parallel::makePSOCKcluster(cores) + on.exit(parallel::stopCluster(cluster), add = TRUE) + + parallel::clusterEvalQ(cluster, { + library(ipmr) + }) + + parallel::clusterExport( + cluster, + c("base_ipm", "samples", "pars", "kernels", "delta"), + envir = environment() + ) + + ## ---- Worker + main_process_opt <- function(i) { + + ipm_temp <- base_ipm + parameters(ipm_temp) <- samples[i, ] + + ipm_temp <- make_ipm(ipm_temp, iterate = FALSE) + + ker_base <- make_iter_kernel(ipm_temp, mega_mat = kernels) + mat_base <- ker_base[[1]] + + lambda_base <- Re(eigen(mat_base, only.values = TRUE)$values[1]) + + pars_all <- parameters(ipm_temp) + + sens <- numeric(length(pars)) + names(sens) <- pars + + ## ---- Sensitivity loop + for (j in seq_along(pars)) { + + par_now <- pars[j] + + ## DOWN perturbation + pars_down <- pars_all + pars_down[[par_now]] <- pars_down[[par_now]] - delta + + ipm_down <- base_ipm + parameters(ipm_down) <- pars_down + + ipm_down <- make_ipm(ipm_down, iterate = FALSE) + ker_down <- make_iter_kernel(ipm_down, mega_mat = kernels) + + lambda_down <- Re(eigen(ker_down[[1]], only.values = TRUE)$values[1]) + + ## UP perturbation + pars_up <- pars_all + pars_up[[par_now]] <- pars_up[[par_now]] + delta + + ipm_up <- base_ipm + parameters(ipm_up) <- pars_up + + ipm_up <- make_ipm(ipm_up, iterate = FALSE) + ker_up <- make_iter_kernel(ipm_up, mega_mat = kernels) + + lambda_up <- Re(eigen(ker_up[[1]], only.values = TRUE)$values[1]) + + sens[j] <- (lambda_up - lambda_down) / (2 * delta) + } + + elas <- sens * abs(as.numeric(pars_all[pars]) / lambda_base) + + list( + lambda = data.frame(lambda = lambda_base), + sens = as.data.frame(t(sens)), + elas = as.data.frame(t(elas)) + ) + } + + ## ---- Run parallel + main <- parallel::parLapply(cluster, seq_len(nsam), main_process_opt) + + ## ---- Combine + lam_sam <- dplyr::bind_rows(lapply(main, `[[`, "lambda")) + sens_sam <- dplyr::bind_rows(lapply(main, `[[`, "sens")) + elas_sam <- dplyr::bind_rows(lapply(main, `[[`, "elas")) + + var_lambda <- var(lam_sam$lambda) + + ## ---- Summary stats + stats_temp <- data.frame( + parameter = pars, + sensitivity = colMeans(sens_sam), + elasticity = colMeans(elas_sam), + variance = NA, + se = NA, + cv = NA + ) + + par_idx <- match(pars, colnames(samples)) + + for (i in seq_len(npar)) { + dat_temp <- samples[, par_idx[i]] + stats_temp[i, "variance"] <- var(dat_temp) + stats_temp[i, "se"] <- sd(dat_temp) / sqrt(length(dat_temp)) + stats_temp[i, "cv"] <- sd(dat_temp) / mean(dat_temp) * 100 + } + + ## ---- Covariance contributions + cov_temp <- cov(samples[, pars]) + + cov_long <- reshape2::melt(cov_temp) + colnames(cov_long) <- c("par1", "par2", "covariance") + + sens_lookup <- stats_temp[, c("parameter", "sensitivity")] + + cov_long <- cov_long %>% + left_join( + sens_lookup %>% rename(sens1 = sensitivity), + by = c("par1" = "parameter") + ) %>% + left_join( + sens_lookup %>% rename(sens2 = sensitivity), + by = c("par2" = "parameter") + ) + + cov_long$variance_part <- cov_long$covariance * cov_long$sens1 * cov_long$sens2 + + var_cont <- cov_long |> + dplyr::group_by(par1) |> + dplyr::summarize(variance_sum = sum(variance_part), .groups = "drop") + + var_cont$variance_prop <- var_cont$variance_sum / var_lambda + + stats_temp <- dplyr::left_join(stats_temp, var_cont, + by = dplyr::join_by(parameter == par1)) + + ## ---- Add vital rates + stats_temp <- dplyr::left_join(stats_temp, vr_table, by = "parameter") + + ## ---- Add summary by vital rates + vr_uncert <- stats_temp %>% + dplyr::group_by(vital_rate) %>% + dplyr::summarise( + variance_sum = sum(variance_sum, na.rm = TRUE), + .groups = "drop" + ) + + ## ---- Return + new_ipmr_uncertainty( + ipm = ipm, + lambdas = lam_sam, + mod_uncert = var_lambda, + params_uncert = stats_temp, + vr_uncert = vr_uncert + ) +} From d6efca990964396091547fcabdc56387a073fd7f Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Mon, 3 Aug 2026 15:18:34 +0200 Subject: [PATCH 05/36] Add descriptive warning when sensitivity analysis perturbation may change sign on parameters --- R/utils-uncertainty.R | 51 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/R/utils-uncertainty.R b/R/utils-uncertainty.R index eb7c5db..ce7566f 100644 --- a/R/utils-uncertainty.R +++ b/R/utils-uncertainty.R @@ -59,6 +59,36 @@ validate_ipm_base <- function(ipm) { ) } +#' Warn if finite-difference perturbation may change parameter sign +#' +#' Internal helper used by sensitivity and uncertainty validation. +#' +#' @param par_names Character vector of parameter names. +#' @param par_values Numeric vector of parameter values corresponding to +#' `par_names`. +#' @param delta Finite-difference perturbation size. +#' +#' @noRd +warn_delta_crossing <- function(par_names, par_values = NULL, delta) { + + if (!is.null(par_values)) { + par_names <- unique(par_names[abs(par_values) <= delta]) + } + + if (length(par_names) > 0) { + warning( + "The following parameter(s) are within `delta` of zero and may change ", + "sign during perturbation: ", + paste(par_names, collapse = ", "), + ". This may result in failed or invalid subkernel construction or failed ", + "sensitivity/uncertainty calculations." + ) + } + + invisible(NULL) +} + + #' Validate inputs to sensitivity function #' @@ -152,6 +182,14 @@ validate_ipm_sensitivity <- function(ipm, pars, kernels, delta) { stop("`delta` must be > 0.") } + ## ---- check if delta changes sign on any parameters + + warn_delta_crossing( + par_names = pars, + par_values = pars_all[pars], + delta = delta + ) + ## return cleaned objects for main function list( ipm = ipm, @@ -301,6 +339,19 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, stop("`delta` must be a single numeric value > 0.") } + ## ---- check if delta changes sign on any parameters + problem_pars <- pars[ + vapply( + pars, + function(p) any(abs(samples[[p]]) <= delta), + logical(1) + ) + ] + + warn_delta_crossing( + par_names = problem_pars, + delta = delta + ) ## ---- cores if (!is.numeric(cores) || length(cores) != 1 || cores < 1) { From 134781cf10b5d68098574409698c324f238fbb5b Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Mon, 3 Aug 2026 15:18:55 +0200 Subject: [PATCH 06/36] Add tests for uncertainty tools --- tests/testthat/test-plot-uncertainty.R | 322 ++++++++++++++++++++ tests/testthat/test-sensitivity.R | 332 ++++++++++++++++++++ tests/testthat/test-uncertainty.R | 382 ++++++++++++++++++++++++ tests/testthat/test-utils-uncertainty.R | 6 + 4 files changed, 1042 insertions(+) create mode 100644 tests/testthat/test-plot-uncertainty.R create mode 100644 tests/testthat/test-sensitivity.R create mode 100644 tests/testthat/test-uncertainty.R create mode 100644 tests/testthat/test-utils-uncertainty.R diff --git a/tests/testthat/test-plot-uncertainty.R b/tests/testthat/test-plot-uncertainty.R new file mode 100644 index 0000000..9d3c098 --- /dev/null +++ b/tests/testthat/test-plot-uncertainty.R @@ -0,0 +1,322 @@ +# Function to build different IPMs to play with + +make_test_ipm <- function(type = "simple_di_det") { + + inv_logit <- function(int, slope, x) { + 1 / (1 + exp(-(int + slope * x))) + } + + if (type == "simple_di_det") { + + data_list <- list( + s_int = 2.2, s_slope = 0.25, + g_int = 0.2, g_slope = 1.02, sd_g = 0.7, + f_r_int = 0.003, f_r_slope = 0.015, + f_s_int = 1.3, f_s_slope = 0.075, + mu_fd = 2, sd_fd = 0.3 + ) + + ipm <- init_ipm(sim_gen = "simple", di_dd = "di", det_stoch = "det") %>% + define_kernel( + "P", + formula = s * g, + family = "CC", + s = inv_logit(s_int, s_slope, dbh_1), + g = dnorm(dbh_2, mu_g, sd_g), + mu_g = g_int + g_slope * dbh_1, + data_list = data_list, + states = list(c("dbh")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "g") + ) %>% + define_kernel( + "F", + formula = f_r * f_s * f_d, + family = "CC", + f_r = inv_logit(f_r_int, f_r_slope, dbh_1), + f_s = exp(f_s_int + f_s_slope * dbh_1), + f_d = dnorm(dbh_2, mu_fd, sd_fd), + data_list = data_list, + states = list(c("dbh")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "f_d") + ) %>% + define_impl( + make_impl_args_list( + c("P", "F"), + int_rule = rep("midpoint", 2), + state_start = rep("dbh", 2), + state_end = rep("dbh", 2) + ) + ) %>% + define_domains(dbh = c(0, 50, 100)) %>% + define_pop_state(n_dbh = runif(100)) %>% + make_ipm( + usr_funs = list(inv_logit = inv_logit), + iterate = TRUE, + iterations = 50, + normalize_pop_size = FALSE + ) + + } else if (type == "simple_dd_det") { + + data_list <- list( + s_int = 2.2, s_slope = 0.25, dd_slope = -0.022, + g_int = 0.2, g_slope = 1.02, sd_g = 0.7, + f_r_int = 0.03, f_r_slope = 0.015, + f_s_int = 1.3, f_s_slope = 0.075, + mu_fd = 0.5, sd_fd = 0.2 + ) + + ipm <- init_ipm(sim_gen = "simple", di_dd = "dd", det_stoch = "det") %>% + define_kernel( + "P", + formula = s * g, + family = "CC", + s = plogis(s_int + s_slope * dbh_1 + dd_slope * sum(n_dbh_t)), + g = dnorm(dbh_2, mu_g, sd_g), + mu_g = g_int + g_slope * dbh_1, + data_list = data_list, + states = list(c("dbh")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "g") + ) %>% + define_kernel( + "F", + formula = f_r * f_s * f_d, + family = "CC", + f_r = plogis(f_r_int + f_r_slope * dbh_1), + f_s = exp(f_s_int + f_s_slope * dbh_1 + dd_slope * sum(n_dbh_t)), + f_d = dnorm(dbh_2, mu_fd, sd_fd), + data_list = data_list, + states = list(c("dbh")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "f_d") + ) %>% + define_impl( + make_impl_args_list( + c("P", "F"), + int_rule = rep("midpoint", 2), + state_start = rep("dbh", 2), + state_end = rep("dbh", 2) + ) + ) %>% + define_domains(dbh = c(0, 50, 100)) %>% + define_pop_state(n_dbh = runif(100)) %>% + make_ipm( + iterate = TRUE, + iterations = 50, + normalize_pop_size = FALSE + ) + + } else if (type == "simple_di_stoch") { + + inv_logit <- function(int, slope, x) { + 1 / (1 + exp(-(int + slope * x))) + } + + mvt_wrapper <- function(r_means, r_sigma, nms) { + out <- mvtnorm::rmvnorm(1, r_means, r_sigma) + out <- as.list(out) + names(out) <- nms + out + } + + data_list <- list( + s_slope = 0.2, + g_slope = 0.99, + g_sd = 0.2, + f_r_slope = 0.003, + f_s_slope = 0.01, + f_d_mu = 2, + f_d_sd = 0.75 + ) + + r_means <- c( + s_int_yr = 0.8, + g_int_yr = 0.1, + f_r_int_yr = 0.3, + f_s_int_yr = 0.01 + ) + + r_sigma <- diag(4) * 0.1 + + ipm <- init_ipm( + sim_gen = "simple", + di_dd = "di", + det_stoch = "stoch", + kern_param = "param" + ) %>% + define_kernel( + "P", + formula = s * g, + family = "CC", + g_mu = g_int_yr + g_slope * surf_area_1, + s = inv_logit(s_int_yr, s_slope, surf_area_1), + g = dnorm(surf_area_2, g_mu, g_sd), + data_list = data_list, + states = list(c("surf_area")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "g") + ) %>% + define_kernel( + "F", + formula = f_r * f_s * f_d, + family = "CC", + f_r = inv_logit(f_r_int_yr, f_r_slope, surf_area_1), + f_s = exp(f_s_int_yr + f_s_slope * surf_area_1), + f_d = dnorm(surf_area_2, f_d_mu, f_d_sd), + data_list = data_list, + states = list(c("surf_area")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "f_d") + ) %>% + define_impl( + make_impl_args_list( + c("P", "F"), + int_rule = rep("midpoint", 2), + state_start = rep("surf_area", 2), + state_end = rep("surf_area", 2) + ) + ) %>% + define_domains(surf_area = c(0, 10, 100)) %>% + define_env_state( + env_params = mvt_wrapper( + r_means, r_sigma, + names(r_means) + ), + data_list = list(r_means = r_means, r_sigma = r_sigma) + ) %>% + define_pop_state(n_surf_area = runif(100)) %>% + make_ipm( + usr_funs = list( + inv_logit = inv_logit, + mvt_wrapper = mvt_wrapper + ), + iterate = TRUE, + iterations = 50, + normalize_pop_size = TRUE + ) + + } else { + stop("Unknown type") + } + + return(ipm) +} + +test_that("plot.ipmr_uncertainty runs", { + + ipm <- make_test_ipm("simple_di_det") + + samples <- as.data.frame(matrix( + rep(unlist(parameters(ipm)), each = 5), + nrow = 5, + byrow = TRUE + )) + colnames(samples) <- names(parameters(ipm)) + + vr <- data.frame( + parameter = names(parameters(ipm)), + vital_rate = rep("test", length(parameters(ipm))) + ) + + res <- uncertainty(ipm, samples = samples, kernels = c("P", "F"), vr_table = vr) + + expect_s3_class(plot(res, type = "param"), "ggplot") + expect_s3_class(plot(res, type = "vr"), "ggplot") + expect_s3_class(plot(res, type = "stacked"), "ggplot") +}) + +test_that("plot type 'all' returns list", { + + ipm <- make_test_ipm("simple_di_det") + + samples <- as.data.frame(matrix( + rep(unlist(parameters(ipm)), each = 5), + nrow = 5, + byrow = TRUE + )) + colnames(samples) <- names(parameters(ipm)) + + vr <- data.frame( + parameter = names(parameters(ipm)), + vital_rate = rep("test", length(parameters(ipm))) + ) + + res <- uncertainty(ipm, samples = samples, kernels = c("P", "F"), vr_table = vr) + + out <- plot(res, type = "all") + + expect_type(out, "list") + expect_true(all(c("param", "vr", "stacked") %in% names(out))) +}) + +test_that("plot errors on unnamed vr_colors", { + + ipm <- make_test_ipm("simple_di_det") + + samples <- as.data.frame(matrix( + rep(unlist(parameters(ipm)), each = 5), + nrow = 5, + byrow = TRUE + )) + colnames(samples) <- names(parameters(ipm)) + + vr <- data.frame( + parameter = names(parameters(ipm)), + vital_rate = rep("test", length(parameters(ipm))) + ) + + res <- uncertainty(ipm, samples = samples, kernels = c("P", "F"), vr_table = vr) + + expect_error( + plot(res, vr_colors = c("red", "blue")), + "must be a named vector" + ) +}) + +test_that("plot errors on missing vr_colors entries", { + + ipm <- make_test_ipm("simple_di_det") + + samples <- as.data.frame(matrix( + rep(unlist(parameters(ipm)), each = 5), + nrow = 5, + byrow = TRUE + )) + colnames(samples) <- names(parameters(ipm)) + + vr <- data.frame( + parameter = names(parameters(ipm)), + vital_rate = rep("test", length(parameters(ipm))) + ) + + res <- uncertainty(ipm, samples = samples, kernels = c("P", "F"), vr_table = vr) + + expect_error( + plot(res, vr_colors = c(foo = "red")), + "Missing colors" + ) +}) + +test_that("print.ipmr_uncertainty runs", { + + ipm <- make_test_ipm("simple_di_det") + + samples <- as.data.frame(matrix( + rep(unlist(parameters(ipm)), each = 5), + nrow = 5, + byrow = TRUE + )) + colnames(samples) <- names(parameters(ipm)) + + vr <- data.frame( + parameter = names(parameters(ipm)), + vital_rate = rep("test", length(parameters(ipm))) + ) + + res <- uncertainty(ipm, samples = samples, kernels = c("P", "F"), vr_table = vr) + + expect_invisible(print(res)) +}) diff --git a/tests/testthat/test-sensitivity.R b/tests/testthat/test-sensitivity.R new file mode 100644 index 0000000..4d71609 --- /dev/null +++ b/tests/testthat/test-sensitivity.R @@ -0,0 +1,332 @@ +# Function to build different IPMs to play with + +make_test_ipm <- function(type = "simple_di_det") { + + inv_logit <- function(int, slope, x) { + 1 / (1 + exp(-(int + slope * x))) + } + + if (type == "simple_di_det") { + + data_list <- list( + s_int = 2.2, s_slope = 0.25, + g_int = 0.2, g_slope = 1.02, sd_g = 0.7, + f_r_int = 0.003, f_r_slope = 0.015, + f_s_int = 1.3, f_s_slope = 0.075, + mu_fd = 2, sd_fd = 0.3 + ) + + ipm <- init_ipm(sim_gen = "simple", di_dd = "di", det_stoch = "det") %>% + define_kernel( + "P", + formula = s * g, + family = "CC", + s = inv_logit(s_int, s_slope, dbh_1), + g = dnorm(dbh_2, mu_g, sd_g), + mu_g = g_int + g_slope * dbh_1, + data_list = data_list, + states = list(c("dbh")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "g") + ) %>% + define_kernel( + "F", + formula = f_r * f_s * f_d, + family = "CC", + f_r = inv_logit(f_r_int, f_r_slope, dbh_1), + f_s = exp(f_s_int + f_s_slope * dbh_1), + f_d = dnorm(dbh_2, mu_fd, sd_fd), + data_list = data_list, + states = list(c("dbh")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "f_d") + ) %>% + define_impl( + make_impl_args_list( + c("P", "F"), + int_rule = rep("midpoint", 2), + state_start = rep("dbh", 2), + state_end = rep("dbh", 2) + ) + ) %>% + define_domains(dbh = c(0, 50, 100)) %>% + define_pop_state(n_dbh = runif(100)) %>% + make_ipm( + usr_funs = list(inv_logit = inv_logit), + iterate = TRUE, + iterations = 50, + normalize_pop_size = FALSE + ) + + } else if (type == "simple_dd_det") { + + data_list <- list( + s_int = 2.2, s_slope = 0.25, dd_slope = -0.022, + g_int = 0.2, g_slope = 1.02, sd_g = 0.7, + f_r_int = 0.03, f_r_slope = 0.015, + f_s_int = 1.3, f_s_slope = 0.075, + mu_fd = 0.5, sd_fd = 0.2 + ) + + ipm <- init_ipm(sim_gen = "simple", di_dd = "dd", det_stoch = "det") %>% + define_kernel( + "P", + formula = s * g, + family = "CC", + s = plogis(s_int + s_slope * dbh_1 + dd_slope * sum(n_dbh_t)), + g = dnorm(dbh_2, mu_g, sd_g), + mu_g = g_int + g_slope * dbh_1, + data_list = data_list, + states = list(c("dbh")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "g") + ) %>% + define_kernel( + "F", + formula = f_r * f_s * f_d, + family = "CC", + f_r = plogis(f_r_int + f_r_slope * dbh_1), + f_s = exp(f_s_int + f_s_slope * dbh_1 + dd_slope * sum(n_dbh_t)), + f_d = dnorm(dbh_2, mu_fd, sd_fd), + data_list = data_list, + states = list(c("dbh")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "f_d") + ) %>% + define_impl( + make_impl_args_list( + c("P", "F"), + int_rule = rep("midpoint", 2), + state_start = rep("dbh", 2), + state_end = rep("dbh", 2) + ) + ) %>% + define_domains(dbh = c(0, 50, 100)) %>% + define_pop_state(n_dbh = runif(100)) %>% + make_ipm( + iterate = TRUE, + iterations = 50, + normalize_pop_size = FALSE + ) + + } else if (type == "simple_di_stoch") { + + inv_logit <- function(int, slope, x) { + 1 / (1 + exp(-(int + slope * x))) + } + + mvt_wrapper <- function(r_means, r_sigma, nms) { + out <- mvtnorm::rmvnorm(1, r_means, r_sigma) + out <- as.list(out) + names(out) <- nms + out + } + + data_list <- list( + s_slope = 0.2, + g_slope = 0.99, + g_sd = 0.2, + f_r_slope = 0.003, + f_s_slope = 0.01, + f_d_mu = 2, + f_d_sd = 0.75 + ) + + r_means <- c( + s_int_yr = 0.8, + g_int_yr = 0.1, + f_r_int_yr = 0.3, + f_s_int_yr = 0.01 + ) + + r_sigma <- diag(4) * 0.1 + + ipm <- init_ipm( + sim_gen = "simple", + di_dd = "di", + det_stoch = "stoch", + kern_param = "param" + ) %>% + define_kernel( + "P", + formula = s * g, + family = "CC", + g_mu = g_int_yr + g_slope * surf_area_1, + s = inv_logit(s_int_yr, s_slope, surf_area_1), + g = dnorm(surf_area_2, g_mu, g_sd), + data_list = data_list, + states = list(c("surf_area")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "g") + ) %>% + define_kernel( + "F", + formula = f_r * f_s * f_d, + family = "CC", + f_r = inv_logit(f_r_int_yr, f_r_slope, surf_area_1), + f_s = exp(f_s_int_yr + f_s_slope * surf_area_1), + f_d = dnorm(surf_area_2, f_d_mu, f_d_sd), + data_list = data_list, + states = list(c("surf_area")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "f_d") + ) %>% + define_impl( + make_impl_args_list( + c("P", "F"), + int_rule = rep("midpoint", 2), + state_start = rep("surf_area", 2), + state_end = rep("surf_area", 2) + ) + ) %>% + define_domains(surf_area = c(0, 10, 100)) %>% + define_env_state( + env_params = mvt_wrapper( + r_means, r_sigma, + names(r_means) + ), + data_list = list(r_means = r_means, r_sigma = r_sigma) + ) %>% + define_pop_state(n_surf_area = runif(100)) %>% + make_ipm( + usr_funs = list( + inv_logit = inv_logit, + mvt_wrapper = mvt_wrapper + ), + iterate = TRUE, + iterations = 50, + normalize_pop_size = TRUE + ) + + } else { + stop("Unknown type") + } + + return(ipm) +} + +ipm <- make_test_ipm("simple_di_det") +ipm_dd <- make_test_ipm("simple_dd_det") +ipm_stoch <- make_test_ipm("simple_di_stoch") + +test_that("sensitivity errors on missing ipm", { + expect_error( + sensitivity(pars = NULL, kernels = NULL, delta = 1e-4), + "`ipm` must be provided and cannot be NULL." + ) +}) + +test_that("sensitivity rejects non-list ipm", { + expect_error( + sensitivity(ipm = 5, pars = NULL, kernels = NULL), + "must be a valid ipmr object" + ) +}) + +test_that("sensitivity runs on simple_di_det", { + + ipm <- make_test_ipm("simple_di_det") + + res <- sensitivity(ipm, kernels = c("P", "F")) + + expect_type(res, "list") + expect_true(all(c("sensitivity", "elasticity") %in% names(res))) + + expect_true(is.numeric(res$sensitivity)) + expect_true(is.numeric(res$elasticity)) +}) + +test_that("sensitivity respects parameter subset", { + + ipm <- make_test_ipm("simple_di_det") + + res <- sensitivity( + ipm, + pars = c("s_int", "g_int"), + kernels = c("P", "F") + ) + + expect_equal(names(res$sensitivity), c("s_int", "g_int")) + expect_equal(length(res$sensitivity), 2) +}) + +test_that("sensitivity rejects density-dependent IPMs", { + + ipm <- make_test_ipm("simple_dd_det") + + expect_error( + sensitivity(ipm, kernels = c("P", "F")), + "This function currently supports density-independent IPMs only." + ) +}) + +test_that("sensitivity rejects stochastic IPMs", { + + ipm <- make_test_ipm("simple_di_stoch") + + expect_error( + sensitivity(ipm, kernels = c("P", "F")), + "This function currently supports deterministic IPMs only." + ) +}) + +test_that("sensitivity errors on unknown parameter", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + sensitivity(ipm, pars = "not_a_param", kernels = c("P", "F")), + "Unknown parameter" + ) +}) + +test_that("elasticity is consistent with sensitivity definition", { + + ipm <- make_test_ipm("simple_di_det") + + res <- sensitivity(ipm, kernels = c("P", "F")) + + pars_all <- parameters(ipm) + lambda <- { + ker <- make_iter_kernel(ipm, mega_mat = c("P", "F")) + Re(eigen(ker[[1]])$values[1]) + } + + expected_elas <- res$sensitivity * abs(as.numeric(pars_all[names(res$sensitivity)]) / lambda) + + expect_equal( + res$elasticity, + expected_elas, + tolerance = 1e-6 + ) +}) + +test_that("duplicate kernels trigger warning", { + + ipm <- make_test_ipm("simple_di_det") + + expect_warning( + sensitivity(ipm, kernels = c("P", "P")), + "Duplicate kernel name" + ) +}) + +test_that("kernel length mismatch triggers warning", { + + ipm <- make_test_ipm("simple_di_det") + + expect_warning( + sensitivity(ipm, kernels = c("P")), + "has length" + ) +}) + +test_that("unknown kernel errors clearly", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + sensitivity(ipm, kernels = c("P", "NOT_REAL")), + "Unknown kernel" + ) +}) diff --git a/tests/testthat/test-uncertainty.R b/tests/testthat/test-uncertainty.R new file mode 100644 index 0000000..2aa35ac --- /dev/null +++ b/tests/testthat/test-uncertainty.R @@ -0,0 +1,382 @@ +# Function to build different IPMs to play with -------------------------------- + +make_test_ipm <- function(type = "simple_di_det") { + + inv_logit <- function(int, slope, x) { + 1 / (1 + exp(-(int + slope * x))) + } + + if (type == "simple_di_det") { + + data_list <- list( + s_int = 2.2, s_slope = 0.25, + g_int = 0.2, g_slope = 1.02, sd_g = 0.7, + f_r_int = 0.003, f_r_slope = 0.015, + f_s_int = 1.3, f_s_slope = 0.075, + mu_fd = 2, sd_fd = 0.3 + ) + + ipm <- init_ipm(sim_gen = "simple", di_dd = "di", det_stoch = "det") %>% + define_kernel( + "P", + formula = s * g, + family = "CC", + s = inv_logit(s_int, s_slope, dbh_1), + g = dnorm(dbh_2, mu_g, sd_g), + mu_g = g_int + g_slope * dbh_1, + data_list = data_list, + states = list(c("dbh")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "g") + ) %>% + define_kernel( + "F", + formula = f_r * f_s * f_d, + family = "CC", + f_r = inv_logit(f_r_int, f_r_slope, dbh_1), + f_s = exp(f_s_int + f_s_slope * dbh_1), + f_d = dnorm(dbh_2, mu_fd, sd_fd), + data_list = data_list, + states = list(c("dbh")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "f_d") + ) %>% + define_impl( + make_impl_args_list( + c("P", "F"), + int_rule = rep("midpoint", 2), + state_start = rep("dbh", 2), + state_end = rep("dbh", 2) + ) + ) %>% + define_domains(dbh = c(0, 50, 100)) %>% + define_pop_state(n_dbh = runif(100)) %>% + make_ipm( + usr_funs = list(inv_logit = inv_logit), + iterate = TRUE, + iterations = 50, + normalize_pop_size = FALSE + ) + + } else if (type == "simple_dd_det") { + + data_list <- list( + s_int = 2.2, s_slope = 0.25, dd_slope = -0.022, + g_int = 0.2, g_slope = 1.02, sd_g = 0.7, + f_r_int = 0.03, f_r_slope = 0.015, + f_s_int = 1.3, f_s_slope = 0.075, + mu_fd = 0.5, sd_fd = 0.2 + ) + + ipm <- init_ipm(sim_gen = "simple", di_dd = "dd", det_stoch = "det") %>% + define_kernel( + "P", + formula = s * g, + family = "CC", + s = plogis(s_int + s_slope * dbh_1 + dd_slope * sum(n_dbh_t)), + g = dnorm(dbh_2, mu_g, sd_g), + mu_g = g_int + g_slope * dbh_1, + data_list = data_list, + states = list(c("dbh")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "g") + ) %>% + define_kernel( + "F", + formula = f_r * f_s * f_d, + family = "CC", + f_r = plogis(f_r_int + f_r_slope * dbh_1), + f_s = exp(f_s_int + f_s_slope * dbh_1 + dd_slope * sum(n_dbh_t)), + f_d = dnorm(dbh_2, mu_fd, sd_fd), + data_list = data_list, + states = list(c("dbh")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "f_d") + ) %>% + define_impl( + make_impl_args_list( + c("P", "F"), + int_rule = rep("midpoint", 2), + state_start = rep("dbh", 2), + state_end = rep("dbh", 2) + ) + ) %>% + define_domains(dbh = c(0, 50, 100)) %>% + define_pop_state(n_dbh = runif(100)) %>% + make_ipm( + iterate = TRUE, + iterations = 50, + normalize_pop_size = FALSE + ) + + } else if (type == "simple_di_stoch") { + + inv_logit <- function(int, slope, x) { + 1 / (1 + exp(-(int + slope * x))) + } + + mvt_wrapper <- function(r_means, r_sigma, nms) { + out <- mvtnorm::rmvnorm(1, r_means, r_sigma) + out <- as.list(out) + names(out) <- nms + out + } + + data_list <- list( + s_slope = 0.2, + g_slope = 0.99, + g_sd = 0.2, + f_r_slope = 0.003, + f_s_slope = 0.01, + f_d_mu = 2, + f_d_sd = 0.75 + ) + + r_means <- c( + s_int_yr = 0.8, + g_int_yr = 0.1, + f_r_int_yr = 0.3, + f_s_int_yr = 0.01 + ) + + r_sigma <- diag(4) * 0.1 + + ipm <- init_ipm( + sim_gen = "simple", + di_dd = "di", + det_stoch = "stoch", + kern_param = "param" + ) %>% + define_kernel( + "P", + formula = s * g, + family = "CC", + g_mu = g_int_yr + g_slope * surf_area_1, + s = inv_logit(s_int_yr, s_slope, surf_area_1), + g = dnorm(surf_area_2, g_mu, g_sd), + data_list = data_list, + states = list(c("surf_area")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "g") + ) %>% + define_kernel( + "F", + formula = f_r * f_s * f_d, + family = "CC", + f_r = inv_logit(f_r_int_yr, f_r_slope, surf_area_1), + f_s = exp(f_s_int_yr + f_s_slope * surf_area_1), + f_d = dnorm(surf_area_2, f_d_mu, f_d_sd), + data_list = data_list, + states = list(c("surf_area")), + evict_cor = TRUE, + evict_fun = truncated_distributions("norm", "f_d") + ) %>% + define_impl( + make_impl_args_list( + c("P", "F"), + int_rule = rep("midpoint", 2), + state_start = rep("surf_area", 2), + state_end = rep("surf_area", 2) + ) + ) %>% + define_domains(surf_area = c(0, 10, 100)) %>% + define_env_state( + env_params = mvt_wrapper( + r_means, r_sigma, + names(r_means) + ), + data_list = list(r_means = r_means, r_sigma = r_sigma) + ) %>% + define_pop_state(n_surf_area = runif(100)) %>% + make_ipm( + usr_funs = list( + inv_logit = inv_logit, + mvt_wrapper = mvt_wrapper + ), + iterate = TRUE, + iterations = 50, + normalize_pop_size = TRUE + ) + + } else { + stop("Unknown type") + } + + return(ipm) +} + +ipm <- make_test_ipm("simple_di_det") +ipm_dd <- make_test_ipm("simple_dd_det") +ipm_stoch <- make_test_ipm("simple_di_stoch") + +# Tests for structural correctness --------------------------------------------- + +test_that("uncertainty errors on missing IPM", { + + expect_error( + uncertainty(pars = NULL, samples = NULL, kernels = NULL, vr_table = NULL), + "`ipm` must be provided" + ) +}) + +test_that("uncertainty rejects non-list IPM", { + + expect_error( + uncertainty(ipm = 5, samples = data.frame(), kernels = NULL, vr_table = NULL), + "valid ipmr object" + ) +}) + +test_that("uncertainty rejects DD IPMs", { + + ipm <- make_test_ipm("simple_dd_det") + + expect_error( + uncertainty(ipm, samples = data.frame(), kernels = c("P","F"), vr_table = NULL), + "density-independent" + ) +}) + +test_that("uncertainty rejects stochastic IPMs", { + + ipm <- make_test_ipm("simple_di_stoch") + + expect_error( + uncertainty(ipm, samples = data.frame(), kernels = c("P","F"), vr_table = NULL), + "deterministic" + ) +}) + +test_that("uncertainty returns correct structure", { + + ipm <- make_test_ipm("simple_di_det") + + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) + + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" + + res <- uncertainty( + ipm, + samples = samples, + kernels = c("P","F"), + vr_table = vr + ) + + expect_s3_class(res, "ipmr_uncertainty") + + expect_true(all(c("lambdas", "mod_uncert", "params_uncert") %in% + names(res))) +}) + +test_that("uncertainty respects parameter subset", { + + ipm <- make_test_ipm("simple_di_det") + + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) + + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" + + res <- uncertainty( + ipm, + pars = c("s_int", "g_int"), + samples = samples, + kernels = c("P","F"), + vr_table = vr + ) + + expect_true(all(res$params_uncert$parameter %in% + c("s_int", "g_int"))) +}) + +# Tests for numerical consistency ---------------------------------------------- + +test_that("uncertainty computes lambda variance", { + + ipm <- make_test_ipm("simple_di_det") + + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) + + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" + + res <- uncertainty( + ipm, + samples = samples, + kernels = c("P","F"), + vr_table = vr + ) + + expect_true(is.numeric(res$mod_uncert)) + expect_true(res$mod_uncert >= 0) +}) + +test_that("uncertainty sensitivities match sensitivity()", { + + ipm <- make_test_ipm("simple_di_det") + + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) + + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" + + unc <- uncertainty( + ipm, + samples = samples, + kernels = c("P","F"), + vr_table = vr + ) + + sens <- sensitivity(ipm, kernels = c("P","F"))$sensitivity + + expect_equal( + unc$params_uncert$sensitivity, + sens[names(unc$params_uncert$parameter)], + tolerance = 1e-6 + ) +}) + +test_that("vr_table join works correctly", { + + ipm <- make_test_ipm("simple_di_det") + + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) + + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "growth" + + res <- uncertainty( + ipm, + samples = samples, + kernels = c("P","F"), + vr_table = vr + ) + + expect_true("vital_rate" %in% names(res$params_uncert)) +}) + +test_that("uncertainty passes kernel warnings", { + + ipm <- make_test_ipm("simple_di_det") + + samples <- as.data.frame(matrix( + rep(unlist(parameters(ipm)), each = 3), + nrow = 3, + byrow = TRUE + )) + colnames(samples) <- names(parameters(ipm)) + + vr <- data.frame( + parameter = names(parameters(ipm)), + vital_rate = "test" + ) + + expect_warning( + uncertainty(ipm, + samples = samples, + kernels = c("P", "P"), + vr_table = vr), + "Duplicate kernel name" + ) +}) diff --git a/tests/testthat/test-utils-uncertainty.R b/tests/testthat/test-utils-uncertainty.R new file mode 100644 index 0000000..939e75d --- /dev/null +++ b/tests/testthat/test-utils-uncertainty.R @@ -0,0 +1,6 @@ +test_that("dominant_eigen returns largest eigenvalue", { + m <- matrix(c(2, 0, + 0, 1), nrow = 2) + + expect_equal(dominant_eigen(m), 2) +}) From c16c3cd0ebb30f5d8fa85e9958fa981324d32035 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Mon, 3 Aug 2026 18:26:13 +0200 Subject: [PATCH 07/36] Add vignette introducing uncertainty analysis tools --- vignettes/uncertainty-analysis.Rmd | 202 +++++++++++++++++++++++++++++ 1 file changed, 202 insertions(+) create mode 100644 vignettes/uncertainty-analysis.Rmd diff --git a/vignettes/uncertainty-analysis.Rmd b/vignettes/uncertainty-analysis.Rmd new file mode 100644 index 0000000..4c87b47 --- /dev/null +++ b/vignettes/uncertainty-analysis.Rmd @@ -0,0 +1,202 @@ +--- +title: "Quantifying Parameter Uncertainty with ipmr" +output: + rmarkdown::html_vignette: + toc: yes + toc_depth: 3 +vignette: > + %\VignetteIndexEntry{Quantifying Parameter Uncertainty with ipmr} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +## Overview + +Integral projection models (IPMs) are commonly used to estimate population growth rates and quantify the demographic importance of different vital rates. However, parameter estimates used to construct IPMs are themselves uncertain, and this uncertainty propagates into model predictions. The `uncertainty()` function estimates how uncertainty in individual parameters contributes to uncertainty in the population growth rate (λ) using finite-difference sensitivities and the covariance structure of sampled parameter estimates. + +This vignette demonstrates how to: +- prepare sampled parameter values +- perform an uncertainty analysis +- interpret the output +- visualize parameter and vital rate contributions to uncertainty + +## Quick start guide {#quick} + +To perform an uncertainty analysis with `uncertainty()`, you will need four required inputs and three optional inputs. + +### `ipm` + +An `ipmr` model object returned by `make_ipm()`. Currently, `uncertainty()` supports **simple** and **general**, **d**ensity-**i**ndependent, **det**erministic IPMs. Stochastic and density-dependent IPMs are not currently supported. + +### `pars` + +A character vector giving the names of the parameters to examine in the uncertainty analysis. In most applications, this should include every parameter that varies among the sampled parameter sets. + +If `pars = NULL`, all parameters used by the IPM are included. Parameters that do not vary across samples will simply have zero estimated uncertainty contribution. + +### `samples` + +A dataframe containing sampled parameter values. Each row should represent one complete parameter set (one sampled model), and each column should contain the sampled values of a single model parameter. Column names must exactly match the parameter names used to construct the `ipmr` object. + +Sampled parameter values may be obtained by bootstrap resampling demographic data and refitting the vital rate models, or by drawing from the joint posterior distribution of a Bayesian model. Guidance for preparing these samples is provided in Section @ref(bootstrap). + +All parameters required to construct the IPM must be included in `samples`, even if some parameters are held constant across all sampled models. We recommend using at least 1,000 sampled parameter sets to obtain stable estimates of parameter uncertainty. + +### `kernels` + +A character vector specifying the subkernel structure of the `ipmr` object. The kernel names must match those used to construct the `ipmr` object and must be supplied in row-major order. + +### `vr_table` (optional) + +An optional dataframe mapping each parameter to the vital rate to which it corresponds. This table is used to summarize parameter uncertainty contributions by vital rate and to support downstream visualization with the `plot()` function. If omitted, uncertainty contributions are only reported at the parameter level. Calling `uncertainty(ipm, vr_table = "template")` generates a template dataframe that can be manually filled in for all parameters in the IPM. + +### `delta` (optional) + +A numeric scalar specifying the the perturbation size used for numerical differentiation when calculating parameter sensitivities. The default is `delta = 1e-4`. For some model constructions, selecting a delta value greater than the absolute value of the smallest parameter value may result in invalid subkernels. Care should be taken to select a sufficiently small value for delta such that adding or subtracting this value from all sampled parameter values does not change the sign on the parameter. + +### `cores` (optional) + +The number of processor cores to use during the analysis. The default (`cores = 1`) performs the analysis sequentially. For larger IPMs or many sampled parameter sets, parallel execution can substantially reduce computation time. + +We recommend using + +```{r} +parallelly::availableCores(omit = 1) +``` + +to determine the number of available cores while leaving one core free for other system processes. + +Once these objects have been prepared, the uncertainty analysis can be performed with a single function call: + +```{r} +unc <- uncertainty( + ipm = ipm, + samples = samples, + pars = pars, + kernels = kernels, + vr_table = vr_table, + cores = parallelly::availableCores(omit = 1) +) +``` + +The resulting object contains the sampled population growth rates, estimates of total model uncertainty, and summaries of uncertainty contributions at both the parameter and (optionally) vital rate levels. + +## Understanding the output {#output} + +`uncertainty()` returns a list containing the original IPM, the sampled population growth rates, and summaries of uncertainty at both the parameter and vital rate levels. + +### `ipm` + +The original IPM object supplied to `uncertainty()`. This is returned unchanged for convenience and reproducibility. + +### `lambdas` + +A numeric vector containing the population growth rates (λ) calculated from each of the sampled parameter sets. The variance of this vector represents the estimated overall model uncertainty. + +### `mod_uncert` + +A single numeric value giving the estimated model uncertainty, calculated as the variance of the sampled λ values (`var(lambdas))`). + +### `params_uncert` + +A dataframe summarizing parameter-level uncertainty. Each row corresponds to a single model parameter and includes: +- the mean sensitivity and elasticity across all sampled models +- the variance, standard error, and coefficient of variation of the sampled parameter values +- the estimated contribution of the parameter to the uncertainty in λ +- the proportion of total model uncertainty explained by the parameter +- the vital rate the parameter corresponds to (if provided) + +The parameter uncertainty contributions are calculated by combining the parameter covariance matrix with the average sensitivities estimated across the sampled models. + +### `vr_uncert` + +If a `vr_table` is supplied, `vr_uncert` summarizes parameter uncertainty contributions by vital rate. Parameter contributions are summed within each vital rate, allowing users to identify which demographic processes contribute most to uncertainty in the population growth rate. If `vr_table` is omitted, this element contains `NA`. + +## Obtaining sampled parameter values {#bootstrap} + +The `samples` argument should contain a representative sample from the joint distribution of the model parameters. Because `uncertainty()` combines parameter sensitivities with the covariance structure of the sampled parameter values, it is important that correlations among parameters are preserved. + +Sampled parameter values may be obtained using either Bayesian or frequentist methods. + +### Bayesian models + +When vital rate models are fit within a Bayesian framework, sampled parameter values can be obtained directly by drawing from the joint posterior distribution. Each posterior draw represents one complete parameter set and can be supplied directly to `uncertainty()` after arranging the draws into a dataframe whose column names match the parameter names used in the IPM. + +### Bootstrap resampling + +For models fit using maximum likelihood or other frequentist methods, we recommend generating sampled parameter values using bootstrap resampling. Briefly, bootstrap resampling repeatedly samples the original demographic data with replacement, refits each vital rate model to the resampled data, and stores the resulting parameter estimates. Repeating this procedure produces a distribution of parameter values while preserving the covariance structure among parameters estimated from the same resampled dataset. + +Although the exact implementation will depend on the study design and modelling framework, the following recommendations will help produce sampled parameter values that are compatible with `uncertainty()`. + +#### Keep constants fixed across bootstrap samples + +Some quantities used to construct an IPM, such as the lower and upper integration limits (`L` and `U`), are determinated by the observed range of the original dataset rather than estimated by a statistical model. These quantities should generally be held constant across bootstrap samples. + +Recalculating such values for every bootstrap replicate may produce artificial variability arising from differences in the sampled observations rather than genuine uncertainty in the demographic process. + +#### Use a consistent model structure + +If model selection is performed before constructing the IPM, it should be performed only once using the original dataset. The selected model structure should then be refit to every bootstrap sample. + +Allowing the model structure to change among bootstrap replicates introduces model selection uncertainty in addition to parameter uncertainty. Because `uncertainty()` is designed to quantify parameter uncertainty, all sampled parameter sets should be generated from the same underlying model structure. + +#### Consider the sampling design + +Bootstrap resampling should reflect the biological and statistical structure of the original study whenever possible/ Depending on the study design, it may be appropriate to stratify sampling by factors such as year, site, or experimental block so that each bootstrap replicate preserves important features of the data. However, it is important to consider the tradeoff between stratification structure and the resulting sample size of each resampled stratified unit. + +The appropriate resampling strategy depends on the size and structure of the dataset, and no single approach is appropriate for all applications. Users should select a strategy that is consistent with the assumptions of their demographic models and the biological processes being studied. + +#### Take caution when independently resampling multiple datasets + +Integral projection models frequently combine information from multiple demographic datasets. For example, survival, growth, and reproduction may be estimated from a long-term monitoring study, while fecundity or recruitment may be estimated from a separate, one-off experiment. + +When these datasets are bootstrapped independently, parameter covariances between demographic processes cannot be accurately estimated because each dataset is resampled separately. Consequently, the covariance terms used in the uncertainty analysis may not reflect biologically-relevant covariance between these terms, or the true covariance may be much higher than estimated, and the resulting uncertainty contributions must be interpreted with caution. + +This limitation is not unique to `uncertainty()` but reflects the information available in the sampled parameter values. Whenever possible, uncertainty analyses should therefore be based on parameter samples generated from a single integrated dataset or from a modelling framework that estimates the full joint parameter distribution. + + +## Plotting the results {#viz} + +The output of `uncertainty()` can be passed directly to `plot()` to visualize the estimated contributions of parameters or vital rates to model uncertainty using `ggplot2`. + +By default, `plot()` returns a list containing all available visualizations. Individual plots may also be specified using the `type` argument. `type = "vr"` and `type = "stacked"` are only available if `vr_table` is specified. + +### Parameter-level uncertainty + +```{r, eval = F, echo = T} +plot(uncert_out, type = "param") +``` + +The parameter-level plot summarizes the uncertainty analysis for every parameter specified in the analysis. Three panels are displayed: +- **Parameter uncertainty**, represented as the coefficient of variation (CV) of the sampled parameter values +- **Model sensitivity**, shown as the mean elasticity of λ to each parameter across all sampled models +- **Parameter contribution to uncertainty**, showing the estimated proportion of model uncertainty attributable to each parameter + +Parameters are grouped and colored according to the vital rate to which they correspond, when supplied, making it easy to distinguish whether uncertainty arises from variation within a vital rate or differences among vital rates. + +### Vital rate uncertainty + +```{r, eval = F. echo = T} +plot(uncert_out, type = "vr") +``` + +When a `vr_table` is supplied to `uncertainty()`, parameter uncertainty contributions are summed within each vital rate. This plot provides a higher-level summary by showing the total contribution of each demographic process to uncertainty in λ. + +This visualization is useful for identifying which processes of the life cycle contribute most to uncertainty in model projections. + +### Model uncertainty decomposition + +```{r, eval = F, echo = T} +plot(uncert_out, type = "stacked") +``` + +The stacked uncertainty plot summarizes the uncertainty contributions of all vital rates within a single bar. Colored segments represent the uncertainty attributed to each vital rate, while the outlined box indicates the total variance in λ estimated from the sampled models. + +This visualization provides a quick assessment of how much of the total model uncertainty is explained by the parameter uncertainty captured in the analysis. + +### Customizing colors + +By default, the plotting functions use a colorblind-friendly palette to distinguish vital rates. Custom colors can be supplied using the `vr_colors` argument as a named character vector, where the names correspond to the values in the `vital_rate` column of the uncertainty output. + + +Note: when exporting figures as PDF files using `ggsave()`, use `device = cairo_pdf` to preserve the λ character in the figure text. From 8b43dea4002281e8c26025c39036fa49f5127196 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 10:47:48 +0200 Subject: [PATCH 08/36] Update documentation for uncertainty tools --- DESCRIPTION | 2 +- NAMESPACE | 4 ++ man/plot.ipmr_uncertainty.Rd | 32 ++++++++++++++ man/print.ipmr_uncertainty.Rd | 19 ++++++++ man/sensitivity.Rd | 41 ++++++++++++++++++ man/uncertainty.Rd | 81 +++++++++++++++++++++++++++++++++++ 6 files changed, 178 insertions(+), 1 deletion(-) create mode 100644 man/plot.ipmr_uncertainty.Rd create mode 100644 man/print.ipmr_uncertainty.Rd create mode 100644 man/sensitivity.Rd create mode 100644 man/uncertainty.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 7b99f72..a02c260 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -71,10 +71,10 @@ Imports: utils, Rcpp Depends: R (>= 3.4) -RoxygenNote: 7.2.2 Config/testthat/parallel: true Config/testthat/edition: 3 LinkingTo: Rcpp URL: https://padrinoDB.github.io/ipmr/, https://github.com/padrinoDB/ipmr BugReports: https://github.com/padrinoDB/ipmr/issues +Config/roxygen2/version: 8.0.0 diff --git a/NAMESPACE b/NAMESPACE index ba2d3c0..03b91fd 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -54,6 +54,7 @@ S3method(parameters,default) S3method(parameters,proto_ipm) S3method(plot,general_di_det_ipm) S3method(plot,ipmr_matrix) +S3method(plot,ipmr_uncertainty) S3method(plot,simple_di_det_ipm) S3method(plot,simple_di_stoch_kern_ipm) S3method(plot,simple_di_stoch_param_ipm) @@ -71,6 +72,7 @@ S3method(print,ipmr_matrix) S3method(print,ipmr_mp_mesh) S3method(print,ipmr_parameters) S3method(print,ipmr_pop_state) +S3method(print,ipmr_uncertainty) S3method(print,ipmr_vital_rate_exprs) S3method(print,ipmr_vital_rate_funs) S3method(print,proto_ipm) @@ -136,9 +138,11 @@ export(quo_name) export(quos) export(right_ev) export(right_mult) +export(sensitivity) export(sym) export(syms) export(truncated_distributions) +export(uncertainty) export(use_vr_model) export(vital_rate_exprs) export(vital_rate_funs) diff --git a/man/plot.ipmr_uncertainty.Rd b/man/plot.ipmr_uncertainty.Rd new file mode 100644 index 0000000..8d08946 --- /dev/null +++ b/man/plot.ipmr_uncertainty.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/generics.R +\name{plot.ipmr_uncertainty} +\alias{plot.ipmr_uncertainty} +\title{Plot method for ipmr uncertainty objects} +\usage{ +\method{plot}{ipmr_uncertainty}(x, type = c("all", "param", "vr", "stacked"), vr_colors = NULL, ...) +} +\arguments{ +\item{x}{An object of class `ipmr_uncertainty`} + +\item{type}{Type of plot: +\itemize{ + \item "param" - parameter-level 3-panel plot + \item "vr" - vital-rate aggregated contributions + \item "stacked" - stacked contributions vs lambda variance + \item "all" - returns all plots as a list +}} + +\item{vr_colors}{A named character vector of colors used to represent +vital rate categories in plots. Names must match the values in the +\code{vital_rate} column of the uncertainty output. If \code{NULL}, +a default colorblind-friendly palette is used.} + +\item{...}{Additional arguments passed to plotting functions} +} +\value{ +A \code{ggplot} object or a list of plots. +} +\description{ +Plot method for ipmr uncertainty objects +} diff --git a/man/print.ipmr_uncertainty.Rd b/man/print.ipmr_uncertainty.Rd new file mode 100644 index 0000000..6687bc7 --- /dev/null +++ b/man/print.ipmr_uncertainty.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/generics.R +\name{print.ipmr_uncertainty} +\alias{print.ipmr_uncertainty} +\title{Print method for ipmr uncertainty objects} +\usage{ +\method{print}{ipmr_uncertainty}(x, ...) +} +\arguments{ +\item{x}{An object of class `ipmr_uncertainty`} + +\item{...}{Unused} +} +\value{ +\code{x} invisibly. +} +\description{ +Print method for ipmr uncertainty objects +} diff --git a/man/sensitivity.Rd b/man/sensitivity.Rd new file mode 100644 index 0000000..e2501a4 --- /dev/null +++ b/man/sensitivity.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/sensitivity.R +\name{sensitivity} +\alias{sensitivity} +\title{Sensitivity and elasticity analysis of population growth rate} +\usage{ +sensitivity(ipm, pars = NULL, kernels, delta = 1e-04) +} +\arguments{ +\item{ipm}{An \code{ipmr} IPM object.} + +\item{pars}{A character vector of parameter names to perturb. Defaults to +all parameters in \code{ipm}.} + +\item{kernels}{A character vector specifying the kernel structure in row +major order.} + +\item{delta}{A numeric scalar specifying the perturbation size used for +numerical differentiation. Default is 1e-4.} +} +\value{ +A named list with two elements: +\itemize{ + \item \code{sensitivity}: Numeric vector of sensitivities of + \eqn{\lambda} to each parameter. + \item \code{elasticity}: Numeric vector of elasticities of + \eqn{\lambda} to each parameter. +} +} +\description{ +Computes numerical sensitivities and elasticities of the +dominant eigenvalue (\eqn{\lambda}) with respect to model parameters. +} +\details{ +This function uses a central difference approximation to estimate the +sensitivity of \eqn{\lambda} to each parameter in \code{pars}. + +Currently implemented for deterministic ("det"), density-independent +("di") IPMs constructed with \code{ipmr}. Stochastic ("stoch") and +density-dependent ("dd") IPMs are not supported. +} diff --git a/man/uncertainty.Rd b/man/uncertainty.Rd new file mode 100644 index 0000000..d298392 --- /dev/null +++ b/man/uncertainty.Rd @@ -0,0 +1,81 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/uncertainty.R +\name{uncertainty} +\alias{uncertainty} +\title{Uncertainty analysis of IPMs} +\usage{ +uncertainty( + ipm, + pars = NULL, + samples, + kernels, + vr_table, + delta = 1e-04, + cores = 1 +) +} +\arguments{ +\item{ipm}{An \code{ipmr} IPM object.} + +\item{pars}{A character vector of parameter names on which to perform the +uncertainty analysis. Defaults to all parameters in \code{ipm}.} + +\item{samples}{A data.frame containing sampled parameter values.} + +\item{kernels}{A character vector specifying the kernel structure in row +major order.} + +\item{vr_table}{A data.frame with columns 'parameter' and 'vital_rate' +mapping each parameter to a user-defined vital rate category. Every +parameter in \code{pars} must be present exactly once in this table.} + +\item{delta}{A numeric scalar specifying the perturbation size used for +numerical differentiation when calculating sensitivities. Default is 1e-4.} + +\item{cores}{An integer specifying the number of cores between which the +computations will be distributed. Defaults to 1.} +} +\description{ +Performs uncertainty analysis of IPMs. +} +\details{ +Currently implemented for deterministic ("det"), density-independent +("di") IPMs constructed with \code{ipmr}. Stochastic ("stoch") and +density-dependent ("dd") IPMs are not supported. + +Sampled parameter values may either be drawn from the joint distribution or +iteratively refit to bootstrapped data. The \code{samples} dataframe must +contain columns for all parameters needed to build the IPM, including those +which do not vary between samples. + +Because parameters may appear in multiple kernels or represent multiple +biological processes, \code{ipmr} does not attempt to infer vital rate +categories automatically. Users must explicitly define this categorization. +A template for \code{vr_table} can be generated by calling the function with +\code{vr_table = "template"}. + +Sensitivities within each sampled model are computed using a central +difference approximation (see \code{\link{sensitivity}}). The \code{delta} +argument controls the perturbation size used in this approximation. + +On systems with multiple cores, computations can be parallelized using the +\code{cores} argument. The function does not automatically infer scheduler +resource limits on high-performance computing (HPC) systems. + +Users running analyses on shared clusters (e.g. via Slurm, PBS, or similar +job schedulers) should ensure that the number of cores requested does not +exceed the allocation specified by the job submission script. Using all +available system cores may lead to oversubscription and reduced performance +of concurrent jobs. +} +\examples{ +# Generate a template for defining vital rates +vr <- uncertainty(ipm, vr_table = "template") + +# Fill in the vital_rate column manually, then run: +# vr$vital_rate <- c("survival", "growth", "reproduction", ...) + +# Run uncertainty analysis +# out <- uncertainty(ipm, pars, samples, kernels, vr_table = vr) + +} From 76482d91ab4e7ee353b037f7e05978581205091e Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 11:00:22 +0200 Subject: [PATCH 09/36] Update chunk options --- vignettes/uncertainty-analysis.Rmd | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/vignettes/uncertainty-analysis.Rmd b/vignettes/uncertainty-analysis.Rmd index 4c87b47..de05d27 100644 --- a/vignettes/uncertainty-analysis.Rmd +++ b/vignettes/uncertainty-analysis.Rmd @@ -60,7 +60,7 @@ The number of processor cores to use during the analysis. The default (`cores = We recommend using -```{r} +```{r, eval = F, echo = T} parallelly::availableCores(omit = 1) ``` @@ -68,7 +68,7 @@ to determine the number of available cores while leaving one core free for other Once these objects have been prepared, the uncertainty analysis can be performed with a single function call: -```{r} +```{r, eval = F, echo = T} unc <- uncertainty( ipm = ipm, samples = samples, @@ -176,7 +176,7 @@ Parameters are grouped and colored according to the vital rate to which they cor ### Vital rate uncertainty -```{r, eval = F. echo = T} +```{r, eval = F, echo = T} plot(uncert_out, type = "vr") ``` From 2844cac72615044951ea23369688a6871c29b4bf Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 11:22:11 +0200 Subject: [PATCH 10/36] Replace lambda character with ASCII code --- R/utils-uncertainty.R | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/R/utils-uncertainty.R b/R/utils-uncertainty.R index ce7566f..6772635 100644 --- a/R/utils-uncertainty.R +++ b/R/utils-uncertainty.R @@ -443,7 +443,7 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, geom_segment(aes(x = 0, xend = variance_prop, yend = parameter), linewidth = 1) + labs(x = "Uncertainty contribution", y = NULL, - title = title_wrap("Parameter contribution to uncertainty in λ"), + title = title_wrap("Parameter contribution to uncertainty in \u03BB"), color = "Vital rate") + scale_color_manual(values = vr_colors) + base_theme @@ -497,7 +497,7 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, xmax = 1.3, ymin = 0, ymax = total_var, - color = "Total variance in λ" + color = "Total variance in \u03BB" ), inherit.aes = FALSE, fill = NA, @@ -506,7 +506,7 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, labs( x = NULL, - y = "Variance in λ", + y = "Variance in \u03BB", title = "Uncertainty decomposition by vital rate", fill = "Vital rate" ) + @@ -514,7 +514,7 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, scale_fill_manual(values = vr_colors) + scale_color_manual( name = NULL, - values = c("Total variance in λ" = "black") + values = c("Total variance in \u03BB" = "black") ) + guides( From 87e5d9fcb0b4d4445f3e868996970fc1b052952e Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 11:22:37 +0200 Subject: [PATCH 11/36] Fix spelling error --- vignettes/uncertainty-analysis.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/uncertainty-analysis.Rmd b/vignettes/uncertainty-analysis.Rmd index de05d27..62c4950 100644 --- a/vignettes/uncertainty-analysis.Rmd +++ b/vignettes/uncertainty-analysis.Rmd @@ -130,7 +130,7 @@ Although the exact implementation will depend on the study design and modelling #### Keep constants fixed across bootstrap samples -Some quantities used to construct an IPM, such as the lower and upper integration limits (`L` and `U`), are determinated by the observed range of the original dataset rather than estimated by a statistical model. These quantities should generally be held constant across bootstrap samples. +Some quantities used to construct an IPM, such as the lower and upper integration limits (`L` and `U`), are determined by the observed range of the original dataset rather than estimated by a statistical model. These quantities should generally be held constant across bootstrap samples. Recalculating such values for every bootstrap replicate may produce artificial variability arising from differences in the sampled observations rather than genuine uncertainty in the demographic process. From 3ffc40cb5ab739748f2297ea9703479fff0d1d15 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 11:23:45 +0200 Subject: [PATCH 12/36] Don't run examples --- R/uncertainty.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/uncertainty.R b/R/uncertainty.R index 4663c86..1a7696c 100644 --- a/R/uncertainty.R +++ b/R/uncertainty.R @@ -49,6 +49,7 @@ #' @return #' #' @examples +#' \dontrun{ #' # Generate a template for defining vital rates #' vr <- uncertainty(ipm, vr_table = "template") #' @@ -57,6 +58,7 @@ #' #' # Run uncertainty analysis #' # out <- uncertainty(ipm, pars, samples, kernels, vr_table = vr) +#' } #' #' @export uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, From fb398ab4367d74cd9d8ea51c50bf1f6b5270a492 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 11:41:10 +0200 Subject: [PATCH 13/36] Fix imports --- R/uncertainty.R | 18 ++--- R/utils-uncertainty.R | 148 +++++++++++++++++++----------------------- 2 files changed, 75 insertions(+), 91 deletions(-) diff --git a/R/uncertainty.R b/R/uncertainty.R index 1a7696c..19bfe98 100644 --- a/R/uncertainty.R +++ b/R/uncertainty.R @@ -170,7 +170,7 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, sens_sam <- dplyr::bind_rows(lapply(main, `[[`, "sens")) elas_sam <- dplyr::bind_rows(lapply(main, `[[`, "elas")) - var_lambda <- var(lam_sam$lambda) + var_lambda <- stats::var(lam_sam$lambda) ## ---- Summary stats stats_temp <- data.frame( @@ -186,13 +186,13 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, for (i in seq_len(npar)) { dat_temp <- samples[, par_idx[i]] - stats_temp[i, "variance"] <- var(dat_temp) - stats_temp[i, "se"] <- sd(dat_temp) / sqrt(length(dat_temp)) - stats_temp[i, "cv"] <- sd(dat_temp) / mean(dat_temp) * 100 + stats_temp[i, "variance"] <- stats::var(dat_temp) + stats_temp[i, "se"] <- stats::sd(dat_temp) / sqrt(length(dat_temp)) + stats_temp[i, "cv"] <- stats::sd(dat_temp) / mean(dat_temp) * 100 } ## ---- Covariance contributions - cov_temp <- cov(samples[, pars]) + cov_temp <- stats::cov(samples[, pars]) cov_long <- reshape2::melt(cov_temp) colnames(cov_long) <- c("par1", "par2", "covariance") @@ -200,12 +200,12 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, sens_lookup <- stats_temp[, c("parameter", "sensitivity")] cov_long <- cov_long %>% - left_join( - sens_lookup %>% rename(sens1 = sensitivity), + dplyr::left_join( + sens_lookup %>% dplyr::rename(sens1 = sensitivity), by = c("par1" = "parameter") ) %>% - left_join( - sens_lookup %>% rename(sens2 = sensitivity), + dplyr::left_join( + sens_lookup %>% dplyr::rename(sens2 = sensitivity), by = c("par2" = "parameter") ) diff --git a/R/utils-uncertainty.R b/R/utils-uncertainty.R index 6772635..0aa795d 100644 --- a/R/utils-uncertainty.R +++ b/R/utils-uncertainty.R @@ -406,49 +406,52 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, levels = rev(df$parameter[order(df$vital_rate)]) ) - library(ggplot2) - library(patchwork) - - base_theme <- theme_minimal() + - theme( - plot.title = element_text(size = 10, face = "bold", hjust = 0.5), - plot.margin = margin(5.5, 10, 5.5, 5.5) + base_theme <- ggplot2::theme_minimal() + + ggplot2::theme( + plot.title = ggplot2::element_text(size = 10, face = "bold", hjust = 0.5), + plot.margin = ggplot2::margin(5.5, 10, 5.5, 5.5) ) title_wrap <- function(x) str_wrap(x, width = 28) - p1 <- ggplot(df, aes(x = cv, y = parameter, color = vital_rate)) + - geom_vline(xintercept = 0) + - geom_point(size = 2.5) + - geom_segment(aes(x = 0, xend = cv, yend = parameter), linewidth = 1) + - labs(x = "Coefficient of variation", - y = "Parameter", - title = title_wrap("Parameter uncertainty")) + - scale_color_manual(values = vr_colors, guide = "none") + + p1 <- ggplot2::ggplot(df, ggplot2::aes(x = cv, y = parameter, + color = vital_rate)) + + ggplot2::geom_vline(xintercept = 0) + + ggplot2::geom_point(size = 2.5) + + ggplot2::geom_segment(ggplot2::aes(x = 0, xend = cv, yend = parameter), + linewidth = 1) + + ggplot2::labs(x = "Coefficient of variation", + y = "Parameter", + title = ggplot2::title_wrap("Parameter uncertainty")) + + ggplot2::scale_color_manual(values = vr_colors, guide = "none") + base_theme - p2 <- ggplot(df, aes(x = elasticity, y = parameter, color = vital_rate)) + - geom_vline(xintercept = 0) + - geom_point(size = 2.5) + - geom_segment(aes(x = 0, xend = elasticity, yend = parameter), linewidth = 1) + - labs(x = "Elasticity", - y = NULL, - title = title_wrap("Model sensitivity to parameter")) + - scale_color_manual(values = vr_colors, guide = "none") + + p2 <- ggplot2::ggplot(df, ggplot2::aes(x = elasticity, y = parameter, + color = vital_rate)) + + ggplot2::geom_vline(xintercept = 0) + + ggplot2::geom_point(size = 2.5) + + ggplot2::geom_segment(ggplot2::aes(x = 0, xend = elasticity, + yend = parameter), linewidth = 1) + + ggplot2::labs(x = "Elasticity", + y = NULL, + title = title_wrap("Model sensitivity to parameter")) + + ggplot2::scale_color_manual(values = vr_colors, guide = "none") + base_theme - p3 <- ggplot(df, aes(x = variance_prop, y = parameter, color = vital_rate)) + - geom_vline(xintercept = 0) + - geom_point(size = 2.5) + - geom_segment(aes(x = 0, xend = variance_prop, yend = parameter), linewidth = 1) + - labs(x = "Uncertainty contribution", - y = NULL, - title = title_wrap("Parameter contribution to uncertainty in \u03BB"), - color = "Vital rate") + - scale_color_manual(values = vr_colors) + + p3 <- ggplot2::ggplot(df, ggplot2::aes(x = variance_prop, y = parameter, + color = vital_rate)) + + ggplot2::geom_vline(xintercept = 0) + + ggplot2::geom_point(size = 2.5) + + ggplot2::geom_segment(ggplot2::aes(x = 0, xend = variance_prop, + yend = parameter), linewidth = 1) + + ggplot2::labs(x = "Uncertainty contribution", + y = NULL, + title = title_wrap("Parameter contribution to uncertainty in \u03BB"), + color = "Vital rate") + + ggplot2::scale_color_manual(values = vr_colors) + base_theme - p1 + p2 + p3 + plot_layout(ncol = 3, axes = 'collect_y') + p1 + p2 + p3 + patchwork::plot_layout(ncol = 3, axes = 'collect_y') } #' Plot uncertainty aggregated by vital rate @@ -463,17 +466,14 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, levels = df$vital_rate[order(df$variance_sum)] ) - library(ggplot2) - - ggplot(df, aes(x = variance_sum, y = vital_rate, fill = vital_rate)) + - geom_col() + - labs( - x = "Total uncertainty contribution", - y = "Vital rate" - ) + - scale_fill_manual(values = vr_colors) + - theme_minimal() + - guides(fill = "none") + ggplot2::ggplot(df, ggplot2::aes(x = variance_sum, y = vital_rate, + fill = vital_rate)) + + ggplot2::geom_col() + + ggplot2::labs( x = "Total uncertainty contribution", + y = "Vital rate" ) + + ggplot2::scale_fill_manual(values = vr_colors) + + ggplot2::theme_minimal() + + ggplot2::guides(fill = "none") } #' Plot model uncertainty and vital rate uncertainty contributions @@ -485,48 +485,32 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, total_var <- x$mod_uncert - library(ggplot2) - - p <- ggplot(df, aes(x = 1, y = variance_sum, fill = vital_rate)) + - geom_col(width = 0.6) + + p <- ggplot2::ggplot(df, ggplot2::aes(x = 1, y = variance_sum, + fill = vital_rate)) + + ggplot2::geom_col(width = 0.6) + # Outline box for total lambda variance - geom_rect( - aes( - xmin = 0.7, - xmax = 1.3, - ymin = 0, - ymax = total_var, - color = "Total variance in \u03BB" - ), + ggplot2::geom_rect( + ggplot2::aes( xmin = 0.7, + xmax = 1.3, + ymin = 0, + ymax = total_var, + color = "Total variance in \u03BB"), inherit.aes = FALSE, fill = NA, - linewidth = 1 - ) + - - labs( - x = NULL, - y = "Variance in \u03BB", - title = "Uncertainty decomposition by vital rate", - fill = "Vital rate" - ) + - - scale_fill_manual(values = vr_colors) + - scale_color_manual( - name = NULL, - values = c("Total variance in \u03BB" = "black") - ) + - - guides( - fill = guide_legend(order = 1), - color = guide_legend(order = 2) - ) + - - theme_minimal() + - theme( - axis.text.x = element_blank(), - axis.ticks.x = element_blank() - ) + linewidth = 1 ) + + ggplot2::labs( x = NULL, + y = "Variance in \u03BB", + title = "Uncertainty decomposition by vital rate", + fill = "Vital rate" ) + + ggplot2::scale_fill_manual(values = vr_colors) + + ggplot2::scale_color_manual( name = NULL, + values = c("Total variance in \u03BB" = "black") ) + + ggplot2::guides( fill = ggplot2::guide_legend(order = 1), + color = guide_legend(order = 2) ) + + ggplot2::theme_minimal() + + ggplot2::theme( axis.text.x = ggplot2::element_blank(), + axis.ticks.x = ggplot2::element_blank() ) p } From b9943fe904469964a46f2d7b8f0511a2fb862201 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 11:45:26 +0200 Subject: [PATCH 14/36] Define global variables --- R/utils-uncertainty.R | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/R/utils-uncertainty.R b/R/utils-uncertainty.R index 0aa795d..fdb1caf 100644 --- a/R/utils-uncertainty.R +++ b/R/utils-uncertainty.R @@ -414,36 +414,36 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, title_wrap <- function(x) str_wrap(x, width = 28) - p1 <- ggplot2::ggplot(df, ggplot2::aes(x = cv, y = parameter, - color = vital_rate)) + + p1 <- ggplot2::ggplot(df, ggplot2::aes(x = .data$cv, y = .data$parameter, + color = .data$vital_rate)) + ggplot2::geom_vline(xintercept = 0) + ggplot2::geom_point(size = 2.5) + - ggplot2::geom_segment(ggplot2::aes(x = 0, xend = cv, yend = parameter), - linewidth = 1) + + ggplot2::geom_segment(ggplot2::aes(x = 0, xend = .data$cv, + yend = .data$parameter), linewidth = 1) + ggplot2::labs(x = "Coefficient of variation", y = "Parameter", title = ggplot2::title_wrap("Parameter uncertainty")) + ggplot2::scale_color_manual(values = vr_colors, guide = "none") + base_theme - p2 <- ggplot2::ggplot(df, ggplot2::aes(x = elasticity, y = parameter, - color = vital_rate)) + + p2 <- ggplot2::ggplot(df, ggplot2::aes(x = .data$elasticity, y = .data$parameter, + color = .data$vital_rate)) + ggplot2::geom_vline(xintercept = 0) + ggplot2::geom_point(size = 2.5) + - ggplot2::geom_segment(ggplot2::aes(x = 0, xend = elasticity, - yend = parameter), linewidth = 1) + + ggplot2::geom_segment(ggplot2::aes(x = 0, xend = .data$elasticity, + yend = .data$parameter), linewidth = 1) + ggplot2::labs(x = "Elasticity", y = NULL, title = title_wrap("Model sensitivity to parameter")) + ggplot2::scale_color_manual(values = vr_colors, guide = "none") + base_theme - p3 <- ggplot2::ggplot(df, ggplot2::aes(x = variance_prop, y = parameter, - color = vital_rate)) + + p3 <- ggplot2::ggplot(df, ggplot2::aes(x = .data$variance_prop, y = .data$parameter, + color = .data$vital_rate)) + ggplot2::geom_vline(xintercept = 0) + ggplot2::geom_point(size = 2.5) + - ggplot2::geom_segment(ggplot2::aes(x = 0, xend = variance_prop, - yend = parameter), linewidth = 1) + + ggplot2::geom_segment(ggplot2::aes(x = 0, xend = .data$variance_prop, + yend = .data$parameter), linewidth = 1) + ggplot2::labs(x = "Uncertainty contribution", y = NULL, title = title_wrap("Parameter contribution to uncertainty in \u03BB"), @@ -466,8 +466,8 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, levels = df$vital_rate[order(df$variance_sum)] ) - ggplot2::ggplot(df, ggplot2::aes(x = variance_sum, y = vital_rate, - fill = vital_rate)) + + ggplot2::ggplot(df, ggplot2::aes(x = .data$variance_sum, y = .data$vital_rate, + fill = .data$vital_rate)) + ggplot2::geom_col() + ggplot2::labs( x = "Total uncertainty contribution", y = "Vital rate" ) + @@ -486,8 +486,8 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, total_var <- x$mod_uncert - p <- ggplot2::ggplot(df, ggplot2::aes(x = 1, y = variance_sum, - fill = vital_rate)) + + p <- ggplot2::ggplot(df, ggplot2::aes(x = 1, y = .data$variance_sum, + fill = .data$vital_rate)) + ggplot2::geom_col(width = 0.6) + # Outline box for total lambda variance ggplot2::geom_rect( From d954fd3375e1598452cce0daf93864c0d8c0b764 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 11:46:26 +0200 Subject: [PATCH 15/36] Fix imports --- R/utils-uncertainty.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/utils-uncertainty.R b/R/utils-uncertainty.R index fdb1caf..f98413a 100644 --- a/R/utils-uncertainty.R +++ b/R/utils-uncertainty.R @@ -412,7 +412,7 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, plot.margin = ggplot2::margin(5.5, 10, 5.5, 5.5) ) - title_wrap <- function(x) str_wrap(x, width = 28) + title_wrap <- function(x) stringr::str_wrap(x, width = 28) p1 <- ggplot2::ggplot(df, ggplot2::aes(x = .data$cv, y = .data$parameter, color = .data$vital_rate)) + From b88af0969d353a03f39bd385d755af72bd58c949 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 12:01:17 +0200 Subject: [PATCH 16/36] Don't run examples --- R/uncertainty.R | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/R/uncertainty.R b/R/uncertainty.R index 19bfe98..2433376 100644 --- a/R/uncertainty.R +++ b/R/uncertainty.R @@ -49,16 +49,14 @@ #' @return #' #' @examples -#' \dontrun{ #' # Generate a template for defining vital rates -#' vr <- uncertainty(ipm, vr_table = "template") +#' # vr <- uncertainty(ipm, vr_table = "template") #' #' # Fill in the vital_rate column manually, then run: #' # vr$vital_rate <- c("survival", "growth", "reproduction", ...) #' #' # Run uncertainty analysis #' # out <- uncertainty(ipm, pars, samples, kernels, vr_table = vr) -#' } #' #' @export uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, From 8f73d97c760df6afe9787697250238da6ad2e9c2 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 12:01:34 +0200 Subject: [PATCH 17/36] Remove lambda symbol from plot legend --- R/utils-uncertainty.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/utils-uncertainty.R b/R/utils-uncertainty.R index f98413a..c771f25 100644 --- a/R/utils-uncertainty.R +++ b/R/utils-uncertainty.R @@ -495,7 +495,7 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, xmax = 1.3, ymin = 0, ymax = total_var, - color = "Total variance in \u03BB"), + color = "Total variance in lambda"), inherit.aes = FALSE, fill = NA, linewidth = 1 ) + @@ -505,7 +505,7 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, fill = "Vital rate" ) + ggplot2::scale_fill_manual(values = vr_colors) + ggplot2::scale_color_manual( name = NULL, - values = c("Total variance in \u03BB" = "black") ) + + values = c("Total variance in lambda" = "black") ) + ggplot2::guides( fill = ggplot2::guide_legend(order = 1), color = guide_legend(order = 2) ) + ggplot2::theme_minimal() + From 8d0593784768cfa034fbea64d86a99363bf80b97 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 12:03:09 +0200 Subject: [PATCH 18/36] Fix imports --- R/utils-uncertainty.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/utils-uncertainty.R b/R/utils-uncertainty.R index c771f25..8ea1515 100644 --- a/R/utils-uncertainty.R +++ b/R/utils-uncertainty.R @@ -507,7 +507,7 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, ggplot2::scale_color_manual( name = NULL, values = c("Total variance in lambda" = "black") ) + ggplot2::guides( fill = ggplot2::guide_legend(order = 1), - color = guide_legend(order = 2) ) + + color = ggplot2::guide_legend(order = 2) ) + ggplot2::theme_minimal() + ggplot2::theme( axis.text.x = ggplot2::element_blank(), axis.ticks.x = ggplot2::element_blank() ) From 034cf8ddfdf936d6d324c5ac921ddedfd83cce15 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 13:08:32 +0200 Subject: [PATCH 19/36] Don't run example --- man/uncertainty.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/man/uncertainty.Rd b/man/uncertainty.Rd index d298392..8a2ebb7 100644 --- a/man/uncertainty.Rd +++ b/man/uncertainty.Rd @@ -70,7 +70,7 @@ of concurrent jobs. } \examples{ # Generate a template for defining vital rates -vr <- uncertainty(ipm, vr_table = "template") +# vr <- uncertainty(ipm, vr_table = "template") # Fill in the vital_rate column manually, then run: # vr$vital_rate <- c("survival", "growth", "reproduction", ...) From 831dc10ce8f7b7ec9d8e5e568600e9990cc4a585 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 13:08:57 +0200 Subject: [PATCH 20/36] Update imports --- DESCRIPTION | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index a02c260..a63a1f5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -69,7 +69,12 @@ Imports: rlang (>= 0.3.0), stats, utils, - Rcpp + Rcpp, + dplyr, + ggplot2, + patchwork, + reshape2, + stringr Depends: R (>= 3.4) Config/testthat/parallel: true Config/testthat/edition: 3 From b42dcdf7209ec24317643cbc8af2f8cfaf0f5b83 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 13:17:20 +0200 Subject: [PATCH 21/36] Define global variables for uncertainty analysis --- R/globals.R | 9 +++++++++ 1 file changed, 9 insertions(+) create mode 100644 R/globals.R diff --git a/R/globals.R b/R/globals.R new file mode 100644 index 0000000..9299a54 --- /dev/null +++ b/R/globals.R @@ -0,0 +1,9 @@ +utils::globalVariables( + c( + "par1", + "parameter", + "variance_part", + "variance_sum", + "vital_rate" + ) +) From 09ce26c9124261cfdb2733c718e6923e376d08bb Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 13:41:11 +0200 Subject: [PATCH 22/36] Fix imports --- R/uncertainty.R | 14 +++++++++----- R/utils-uncertainty.R | 10 +++++++++- 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/R/uncertainty.R b/R/uncertainty.R index 2433376..fa1c55c 100644 --- a/R/uncertainty.R +++ b/R/uncertainty.R @@ -58,6 +58,10 @@ #' # Run uncertainty analysis #' # out <- uncertainty(ipm, pars, samples, kernels, vr_table = vr) #' +#' @importFrom dplyr bind_rows left_join rename group_by summarize join_by +#' @importFrom stats var sd cov +#' @importFrom reshape2 melt +#' #' @export uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, delta = 1e-4, cores = 1) { @@ -210,8 +214,8 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, cov_long$variance_part <- cov_long$covariance * cov_long$sens1 * cov_long$sens2 var_cont <- cov_long |> - dplyr::group_by(par1) |> - dplyr::summarize(variance_sum = sum(variance_part), .groups = "drop") + dplyr::group_by(.data$par1) |> + dplyr::summarize(variance_sum = sum(.data$variance_part), .groups = "drop") var_cont$variance_prop <- var_cont$variance_sum / var_lambda @@ -223,9 +227,9 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, ## ---- Add summary by vital rates vr_uncert <- stats_temp %>% - dplyr::group_by(vital_rate) %>% - dplyr::summarise( - variance_sum = sum(variance_sum, na.rm = TRUE), + dplyr::group_by(.data$vital_rate) %>% + dplyr::summarize( + variance_sum = sum(.data$variance_sum, na.rm = TRUE), .groups = "drop" ) diff --git a/R/utils-uncertainty.R b/R/utils-uncertainty.R index 8ea1515..fd9fd65 100644 --- a/R/utils-uncertainty.R +++ b/R/utils-uncertainty.R @@ -396,6 +396,10 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, #' Plot three-panel uncertainty plot #' #' @noRd +#' @importFrom patchwork plot_layout +#' @importFrom ggplot2 theme_minimal theme element_text margin aes ggplot +#' geom_vline geom_point geom_segment labs scale_color_manual +#' @importFrom stringr str_wrap .plot_uncert_param <- function(x, vr_colors, ...) { df <- x$params_uncert @@ -422,7 +426,7 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, yend = .data$parameter), linewidth = 1) + ggplot2::labs(x = "Coefficient of variation", y = "Parameter", - title = ggplot2::title_wrap("Parameter uncertainty")) + + title = title_wrap("Parameter uncertainty")) + ggplot2::scale_color_manual(values = vr_colors, guide = "none") + base_theme @@ -457,6 +461,8 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, #' Plot uncertainty aggregated by vital rate #' #' @noRd +#' @importFrom ggplot2 ggplot aes geom_col labs scale_fill_manual theme_minimal +#' guides .plot_uncert_vr <- function(x, vr_colors, ...) { df <- x$vr_uncert @@ -479,6 +485,8 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, #' Plot model uncertainty and vital rate uncertainty contributions #' #' @noRd +#' @importFrom ggplot2 ggplot aes geom_col geom_rect labs scale_fill_manual +#' scale_color_manual guides guide_legend theme theme_minimal element_blank .plot_uncert_stacked <- function(x, vr_colors, ...) { df <- x$vr_uncert From 6d6a2a82a5796e2236bed658186fb9050404578a Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 13:43:56 +0200 Subject: [PATCH 23/36] Fix imports --- NAMESPACE | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/NAMESPACE b/NAMESPACE index 03b91fd..5e9d330 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -148,6 +148,29 @@ export(vital_rate_exprs) export(vital_rate_funs) importFrom(Rcpp,evalCpp) importFrom(Rcpp,sourceCpp) +importFrom(dplyr,bind_rows) +importFrom(dplyr,group_by) +importFrom(dplyr,join_by) +importFrom(dplyr,left_join) +importFrom(dplyr,rename) +importFrom(dplyr,summarize) +importFrom(ggplot2,aes) +importFrom(ggplot2,element_blank) +importFrom(ggplot2,element_text) +importFrom(ggplot2,geom_col) +importFrom(ggplot2,geom_point) +importFrom(ggplot2,geom_rect) +importFrom(ggplot2,geom_segment) +importFrom(ggplot2,geom_vline) +importFrom(ggplot2,ggplot) +importFrom(ggplot2,guide_legend) +importFrom(ggplot2,guides) +importFrom(ggplot2,labs) +importFrom(ggplot2,margin) +importFrom(ggplot2,scale_color_manual) +importFrom(ggplot2,scale_fill_manual) +importFrom(ggplot2,theme) +importFrom(ggplot2,theme_minimal) importFrom(grDevices,grey) importFrom(grDevices,rainbow) importFrom(graphics,abline) @@ -158,11 +181,13 @@ importFrom(graphics,layout) importFrom(graphics,par) importFrom(magrittr,"%>%") importFrom(methods,hasArg) +importFrom(patchwork,plot_layout) importFrom(purrr,flatten) importFrom(purrr,map) importFrom(purrr,map2) importFrom(purrr,map_chr) importFrom(purrr,map_dbl) +importFrom(reshape2,melt) importFrom(rlang,":=") importFrom(rlang,.data) importFrom(rlang,as_environment) @@ -196,10 +221,14 @@ importFrom(rlang,quos) importFrom(rlang,sym) importFrom(rlang,syms) importFrom(stats,bw.nrd0) +importFrom(stats,cov) importFrom(stats,dnorm) importFrom(stats,na.omit) importFrom(stats,runif) +importFrom(stats,sd) importFrom(stats,setNames) +importFrom(stats,var) +importFrom(stringr,str_wrap) importFrom(utils,globalVariables) importFrom(utils,packageDescription) useDynLib(ipmr, .registration=TRUE) From c284fa7d34eb1dd97e717fb11fd811ad0f8107ae Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 13:54:22 +0200 Subject: [PATCH 24/36] Fix call to parameter values --- R/utils-uncertainty.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/utils-uncertainty.R b/R/utils-uncertainty.R index fd9fd65..0bbe8c3 100644 --- a/R/utils-uncertainty.R +++ b/R/utils-uncertainty.R @@ -186,7 +186,7 @@ validate_ipm_sensitivity <- function(ipm, pars, kernels, delta) { warn_delta_crossing( par_names = pars, - par_values = pars_all[pars], + par_values = as.numeric(pars_all[pars]), delta = delta ) From 0baafc22a3c54296904a41dcd78574d487964f1b Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 14:22:30 +0200 Subject: [PATCH 25/36] Fix expected error messages --- tests/testthat/test-sensitivity.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-sensitivity.R b/tests/testthat/test-sensitivity.R index 4d71609..0ecc455 100644 --- a/tests/testthat/test-sensitivity.R +++ b/tests/testthat/test-sensitivity.R @@ -256,7 +256,7 @@ test_that("sensitivity rejects density-dependent IPMs", { expect_error( sensitivity(ipm, kernels = c("P", "F")), - "This function currently supports density-independent IPMs only." + "Only density-independent IPMs are currently supported." ) }) @@ -266,7 +266,7 @@ test_that("sensitivity rejects stochastic IPMs", { expect_error( sensitivity(ipm, kernels = c("P", "F")), - "This function currently supports deterministic IPMs only." + "Only deterministic IPMs are currently supported." ) }) From d86c429ab3c93020108ce448d518aceb33f062ad Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 14:22:43 +0200 Subject: [PATCH 26/36] Update join syntax --- R/uncertainty.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/uncertainty.R b/R/uncertainty.R index fa1c55c..3f55f55 100644 --- a/R/uncertainty.R +++ b/R/uncertainty.R @@ -220,7 +220,7 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, var_cont$variance_prop <- var_cont$variance_sum / var_lambda stats_temp <- dplyr::left_join(stats_temp, var_cont, - by = dplyr::join_by(parameter == par1)) + by = c("parameter" == "par1")) ## ---- Add vital rates stats_temp <- dplyr::left_join(stats_temp, vr_table, by = "parameter") From b3ced50810a3147d9421210242818d0eee9ce6e3 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 14:39:31 +0200 Subject: [PATCH 27/36] Fix join syntax --- R/uncertainty.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/uncertainty.R b/R/uncertainty.R index 3f55f55..5bcaf12 100644 --- a/R/uncertainty.R +++ b/R/uncertainty.R @@ -220,7 +220,7 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, var_cont$variance_prop <- var_cont$variance_sum / var_lambda stats_temp <- dplyr::left_join(stats_temp, var_cont, - by = c("parameter" == "par1")) + by = c("parameter" = "par1")) ## ---- Add vital rates stats_temp <- dplyr::left_join(stats_temp, vr_table, by = "parameter") From 75b5bf3f92041b27bcc52b327331b5a64b901696 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 14:39:53 +0200 Subject: [PATCH 28/36] Fix syntax in expect equal --- tests/testthat/test-uncertainty.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-uncertainty.R b/tests/testthat/test-uncertainty.R index 2aa35ac..4ea9fe9 100644 --- a/tests/testthat/test-uncertainty.R +++ b/tests/testthat/test-uncertainty.R @@ -332,7 +332,7 @@ test_that("uncertainty sensitivities match sensitivity()", { expect_equal( unc$params_uncert$sensitivity, - sens[names(unc$params_uncert$parameter)], + sens[unc$params_uncert$parameter], tolerance = 1e-6 ) }) From f656d8ec1c8282200835ba2fca782f77956345e9 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Tue, 4 Aug 2026 14:59:31 +0200 Subject: [PATCH 29/36] Fix expect equal syntax --- tests/testthat/test-uncertainty.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-uncertainty.R b/tests/testthat/test-uncertainty.R index 4ea9fe9..4645d2a 100644 --- a/tests/testthat/test-uncertainty.R +++ b/tests/testthat/test-uncertainty.R @@ -331,8 +331,8 @@ test_that("uncertainty sensitivities match sensitivity()", { sens <- sensitivity(ipm, kernels = c("P","F"))$sensitivity expect_equal( - unc$params_uncert$sensitivity, - sens[unc$params_uncert$parameter], + unname(unc$params_uncert$sensitivity), + unname(sens[unc$params_uncert$parameter]), tolerance = 1e-6 ) }) From facf4c0bab8646256f83d7f5ee6422d34c0b0968 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Wed, 5 Aug 2026 19:19:54 +0200 Subject: [PATCH 30/36] Add functionality to accommodate multiple methods of finite-difference sensitivity calculations --- R/sensitivity.R | 164 ++++++++++++---- R/uncertainty.R | 137 +++++++++---- R/utils-uncertainty.R | 302 ++++++++++++++++++++++++++--- tests/testthat/test-sensitivity.R | 192 ++++++++++++++++++ tests/testthat/test-uncertainty.R | 194 ++++++++++++++++++ vignettes/uncertainty-analysis.Rmd | 2 +- 6 files changed, 879 insertions(+), 112 deletions(-) diff --git a/R/sensitivity.R b/R/sensitivity.R index 5106660..b168b1b 100644 --- a/R/sensitivity.R +++ b/R/sensitivity.R @@ -1,7 +1,8 @@ #' @title Sensitivity and elasticity analysis of population growth rate #' #' @description Computes numerical sensitivities and elasticities of the -#' dominant eigenvalue (\eqn{\lambda}) with respect to model parameters. +#' dominant eigenvalue (\eqn{\lambda}) with respect to model parameters using +#' finite differences. #' #' @param ipm An \code{ipmr} IPM object. #' @param pars A character vector of parameter names to perturb. Defaults to @@ -10,41 +11,105 @@ #' major order. #' @param delta A numeric scalar specifying the perturbation size used for #' numerical differentiation. Default is 1e-4. +#' @param bounds An optional named list specifying lower and/or upper bounds +#' for model parameters when calculating numerical sensitivities. Each element +#' should be a numeric vector of length two giving the lower and upper bounds, +#' respectively (e.g., \code{list(surv = c(0, 1), recruit = c(0, Inf))}). +#' Parameters omitted from this list are assumed to be unbounded and are +#' perturbed using the central finite difference. For parameters with +#' specified bounds, the function automatically uses a forward or backward +#' finite difference whenever perturbation by \code{delta} would cross a +#' boundary; otherwise, the central finite difference is used. #' #'@details -#' This function uses a central difference approximation to estimate the -#' sensitivity of \eqn{\lambda} to each parameter in \code{pars}. +#' This function uses perturbation to estimate the sensitivity of \eqn{\lambda} +#' to each parameter in \code{pars}. #' #' Currently implemented for deterministic ("det"), density-independent #' ("di") IPMs constructed with \code{ipmr}. Stochastic ("stoch") and #' density-dependent ("dd") IPMs are not supported. #' -#' @return A named list with two elements: +#' ## Finite differences and bounded parameters +#' +#' Sensitivities are estimated numerically by perturbing each parameter by +#' \code{delta}. By default, a central finite difference is used for all +#' parameters. However, some demographic parameters are biologically constrained +#' (e.g., probabilities bounded between 0 and 1, or fecundity parameters bounded +#' below by 0). If perturbing a parameter by \code{delta} would cross one of +#' these user-specified bounds, the function automatically computes the +#' appropriate one-sided finite difference (forward or backward) rather than +#' evaluating the model outside the feasible parameter space. +#' +#' One-sided finite differences approximate the same local derivative as the +#' central finite difference, but are used only when a symmetric perturbation +#' would leave the biologically feasible parameter space. Consequently, +#' sensitivities calculated using one-sided and central finite differences are +#' directly comparable and may be used together in downstream uncertainty +#' analyses. +#' +#' Supplying \code{bounds} is recommended whenever model parameters are subject +#' to biological constraints and numerical perturbation could otherwise produce +#' invalid parameter values or biologically impossible kernels. +#' +#' @return A named list with three elements: #' \itemize{ -#' \item \code{sensitivity}: Numeric vector of sensitivities of -#' \eqn{\lambda} to each parameter. -#' \item \code{elasticity}: Numeric vector of elasticities of -#' \eqn{\lambda} to each parameter. +#' \item \code{sensitivity}: A named numeric vector giving the numerical +#' sensitivities of the dominant eigenvalue (\eqn{\lambda}) to each +#' parameter. +#' \item \code{elasticity}: A named numeric vector giving the elasticities +#' of the dominant eigenvalue (\eqn{\lambda}) to each parameter. +#' \item \code{difference_method}: A named character vector indicating the +#' finite-difference scheme used to calculate the sensitivity for each +#' parameter. This component is provided to facilitate reproducibility and +#' to identify parameters whose sensitivities were evaluated using one-sided +#' finite differences because of user-specified parameter bounds. Values are +#' \code{"central"}, \code{"forward"}, or \code{"backward"}. #' } #' +#' @examples +#' # Probability bounded between 0 and 1 +#' # bounds <- list( +#' # surv = c(0, 1) +#' # ) +#' +#' # Recruitment bounded below by zero +#' # bounds <- list( +#' # recruit = c(0, Inf) +#' # ) +#' +#' # Multiple constrained parameters +#' # bounds <- list( +#' # surv = c(0, 1), +#' # recruit = c(0, Inf), +#' # fec = c(0, Inf) +#' # ) +#' +#' # Compute sensitivities using defined parameter boundaries +#' # sens <- sensitivity( +#' # ipm, +#' # kernels = c("P", "F"), +#' # bounds = bounds ) +#' #' @export -sensitivity <- function(ipm, pars = NULL, kernels, delta = 1e-4) { +sensitivity <- function(ipm, pars = NULL, kernels, delta = 1e-4, + bounds = NULL) { ## ---- Input checks - val <- validate_ipm_sensitivity(ipm, pars, kernels, delta) + val <- validate_ipm_sensitivity(ipm, pars, kernels, delta, bounds) ipm <- val$ipm pars <- val$pars pars_all <- val$pars_all kernels <- val$kernels delta <- val$delta + bounds <- val$bounds ## ---- Setup - npar <- length(pars) - spar <- numeric(npar) - vpar <- pars_all + npar <- length(pars) + spar <- numeric(npar) + diff_method <- character(npar) kernel <- make_iter_kernel(ipm, mega_mat = kernels) @@ -52,47 +117,64 @@ sensitivity <- function(ipm, pars = NULL, kernels, delta = 1e-4) { stop("Kernel construction failed.") } - eig <- tryCatch( - eigen(kernel[[1]]), - error = function(e) stop("Eigen decomposition failed for base kernel.") - ) - lambda_all <- dominant_eigen(kernel[[1]]) + ## ---- Helper to perturb a single parameter and calculate lambda + + perturb_lambda <- function(par, value) { + vpar <- pars_all + vpar[[par]] <- value + ipm_tmp <- ipm$proto_ipm + parameters(ipm_tmp) <- vpar + ker <- ipm_tmp %>% + make_ipm(iterate = FALSE) %>% + make_iter_kernel(mega_mat = kernels) + dominant_eigen(ker[[1]]) + } + ## ---- Loop over each parameter - for(i in seq_len(npar)) { + for (i in seq_len(npar)) { par_now <- pars[i] + theta <- pars_all[[par_now]] - # ---- DOWN perturbation - vpar <- pars_all - vpar[[par_now]] <- vpar[[par_now]] - delta + ## ---- Define bounds (default none) - ipm_down <- ipm$proto_ipm - parameters(ipm_down) <- vpar + if (!is.null(bounds) && par_now %in% names(bounds)) { + lower <- bounds[[par_now]][1] + upper <- bounds[[par_now]][2] + } else { + lower <- -Inf + upper <- Inf + } - ker_down <- ipm_down %>% - make_ipm(iterate = FALSE) %>% - make_iter_kernel(mega_mat = kernels) + use_forward <- theta - delta < lower + use_backward <- theta + delta > upper - lambda_down <- dominant_eigen(ker_down[[1]]) + ## ---- Forward difference - # ---- UP perturbation - vpar <- pars_all - vpar[[par_now]] <- vpar[[par_now]] + delta + if (use_forward) { + lambda_up <- perturb_lambda(par_now, theta + delta) - ipm_up <- ipm$proto_ipm - parameters(ipm_up) <- vpar + spar[i] <- (lambda_up - lambda_all)/delta + diff_method[i] <- "forward" - ker_up <- ipm_up %>% - make_ipm(iterate = FALSE) %>% - make_iter_kernel(mega_mat = kernels) + ## ---- Backward difference + } else if (use_backward) { + lambda_down <- perturb_lambda(par_now, theta - delta) + + spar[i] <- (lambda_all - lambda_down)/delta + diff_method[i] <- "backward" - lambda_up <- dominant_eigen(ker_up[[1]]) + ## ---- Central difference + } else { + lambda_down <- perturb_lambda(par_now, theta - delta) + lambda_up <- perturb_lambda(par_now, theta + delta) - # ---- Central difference - spar[i] <- (lambda_up - lambda_down) / (2 * delta) + spar[i] <- (lambda_up - lambda_down)/(2*delta) + diff_method[i] <- "central" + } } ## ---- Output @@ -102,9 +184,11 @@ sensitivity <- function(ipm, pars = NULL, kernels, delta = 1e-4) { names(sens) <- pars names(elas) <- pars + names(diff_method) <- pars return(list( sensitivity = sens, - elasticity = elas + elasticity = elas, + difference_method = diff_method )) } diff --git a/R/uncertainty.R b/R/uncertainty.R index 5bcaf12..66cda32 100644 --- a/R/uncertainty.R +++ b/R/uncertainty.R @@ -13,6 +13,15 @@ #' parameter in \code{pars} must be present exactly once in this table. #' @param delta A numeric scalar specifying the perturbation size used for #' numerical differentiation when calculating sensitivities. Default is 1e-4. +#' @param bounds An optional named list specifying lower and/or upper bounds +#' for model parameters when calculating numerical sensitivities. Each element +#' should be a numeric vector of length two giving the lower and upper bounds, +#' respectively (e.g., \code{list(surv = c(0, 1), recruit = c(0, Inf))}). +#' Parameters omitted from this list are assumed to be unbounded and are +#' perturbed using the central finite difference. For parameters with +#' specified bounds, the function automatically uses a forward or backward +#' finite difference whenever perturbation by \code{delta} would cross a +#' boundary; otherwise, the central finite difference is used. #' @param cores An integer specifying the number of cores between which the #' computations will be distributed. Defaults to 1. #' @@ -32,9 +41,12 @@ #' A template for \code{vr_table} can be generated by calling the function with #' \code{vr_table = "template"}. #' -#' Sensitivities within each sampled model are computed using a central -#' difference approximation (see \code{\link{sensitivity}}). The \code{delta} -#' argument controls the perturbation size used in this approximation. +#' The finite-difference scheme used in \code{\link{sensitivity}} is used when +#' estimating parameter sensitivities for the uncertainty decomposition. +#' Consequently, supplying biological parameter bounds ensures that +#' sensitivities are evaluated only within the feasible parameter space +#' throughout the uncertainty analysis. The \code{delta} argument controls the +#' perturbation size used in the sensitivity approximation. #' #' On systems with multiple cores, computations can be parallelized using the #' \code{cores} argument. The function does not automatically infer scheduler @@ -47,6 +59,19 @@ #' of concurrent jobs. #' #' @return +#' An object of class \code{"ipmr_uncertainty"} containing: +#' \itemize{ +#' \item \code{ipm}: The complete IPM object as specified in \code{ipm}. +#' \item \code{lambdas}: Population growth rates (\eqn{\lambda}) calculated +#' for each sampled parameter set. +#' \item \code{mod_uncert}: The variance in \eqn{\lambda} across all sampled +#' parameter sets. +#' \item \code{params_uncert}: A data frame summarizing parameter-level +#' sensitivities, elasticities, finite-difference methods, sampling +#' variance, coefficients of variation, and contributions to uncertainty. +#' \item \code{vr_uncert}: A data frame summarizing variance contributions +#' aggregated by user-defined vital rate, if given. +#' } #' #' @examples #' # Generate a template for defining vital rates @@ -64,7 +89,7 @@ #' #' @export uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, - delta = 1e-4, cores = 1) { + delta = 1e-4, bounds = NULL, cores = 1) { if (isTRUE(vr_table == "template")) { if (is.null(pars)) { @@ -78,13 +103,15 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, ## ---- Input checks val <- validate_ipm_uncertainty(ipm, pars, samples, kernels, - vr_table, delta, cores) + vr_table, delta, bounds, cores) ipm <- val$ipm pars <- val$pars samples <- val$samples kernels <- val$kernels vr_table <- val$vr_table + delta <- val$delta + bounds <- val$bounds cores <- val$cores npar <- length(pars) @@ -101,7 +128,7 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, parallel::clusterExport( cluster, - c("base_ipm", "samples", "pars", "kernels", "delta"), + c("base_ipm", "samples", "pars", "kernels", "delta", "bounds"), envir = environment() ) @@ -120,50 +147,71 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, pars_all <- parameters(ipm_temp) + calc_lambda <- function(par, value) { + pars_tmp <- pars_all + pars_tmp[[par]] <- value + ipm_tmp <- base_ipm + parameters(ipm_tmp) <- pars_tmp + ipm_tmp <- make_ipm(ipm_tmp, iterate = FALSE) + ker <- make_iter_kernel(ipm_tmp, + mega_mat = kernels) + Re(eigen(ker[[1]], only.values = TRUE)$values[1]) + } + sens <- numeric(length(pars)) names(sens) <- pars + diff_method <- character(length(pars)) ## ---- Sensitivity loop for (j in seq_along(pars)) { - par_now <- pars[j] - - ## DOWN perturbation - pars_down <- pars_all - pars_down[[par_now]] <- pars_down[[par_now]] - delta - - ipm_down <- base_ipm - parameters(ipm_down) <- pars_down - - ipm_down <- make_ipm(ipm_down, iterate = FALSE) - ker_down <- make_iter_kernel(ipm_down, mega_mat = kernels) - - lambda_down <- Re(eigen(ker_down[[1]], only.values = TRUE)$values[1]) - - ## UP perturbation - pars_up <- pars_all - pars_up[[par_now]] <- pars_up[[par_now]] + delta - - ipm_up <- base_ipm - parameters(ipm_up) <- pars_up - - ipm_up <- make_ipm(ipm_up, iterate = FALSE) - ker_up <- make_iter_kernel(ipm_up, mega_mat = kernels) - - lambda_up <- Re(eigen(ker_up[[1]], only.values = TRUE)$values[1]) - - sens[j] <- (lambda_up - lambda_down) / (2 * delta) + theta <- pars_all[[par_now]] + if (!is.null(bounds) && + par_now %in% names(bounds)) { + lower <- bounds[[par_now]][1] + upper <- bounds[[par_now]][2] + } else { + lower <- -Inf + upper <- Inf + } + use_forward <- theta - delta < lower + use_backward <- theta + delta > upper + if (use_forward) { + lambda_up <- calc_lambda(par_now, theta + delta) + sens[j] <- (lambda_up - lambda_base) / delta + diff_method[j] <- "forward" + } else if (use_backward) { + lambda_down <- calc_lambda(par_now, theta - delta) + sens[j] <- (lambda_base - lambda_down) / delta + diff_method[j] <- "backward" + } else { + lambda_down <- calc_lambda(par_now, theta - delta) + lambda_up <- calc_lambda(par_now, theta + delta) + sens[j] <- (lambda_up - lambda_down) / (2 * delta) + diff_method[j] <- "central" + } } elas <- sens * abs(as.numeric(pars_all[pars]) / lambda_base) list( - lambda = data.frame(lambda = lambda_base), - sens = as.data.frame(t(sens)), - elas = as.data.frame(t(elas)) + lambda = data.frame(lambda = lambda_base), + sens = as.data.frame(t(sens)), + elas = as.data.frame(t(elas)), + diff_method = as.data.frame(t(diff_method), + stringsAsFactors = FALSE) ) } + summarize_difference_method <- function(x) { + tab <- sort(table(x), decreasing = TRUE) + if (length(tab) == 1) { + return(names(tab)) + } + pct <- format( round(100 * tab / sum(tab), 1),trim = TRUE, nsmall = 0 ) + paste0( "mixed (", paste0(pct, "% ", names(tab), collapse = ", "), ")" ) + } + ## ---- Run parallel main <- parallel::parLapply(cluster, seq_len(nsam), main_process_opt) @@ -171,19 +219,24 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, lam_sam <- dplyr::bind_rows(lapply(main, `[[`, "lambda")) sens_sam <- dplyr::bind_rows(lapply(main, `[[`, "sens")) elas_sam <- dplyr::bind_rows(lapply(main, `[[`, "elas")) + diff_sam <- dplyr::bind_rows(lapply(main, `[[`, "diff_method")) var_lambda <- stats::var(lam_sam$lambda) ## ---- Summary stats stats_temp <- data.frame( - parameter = pars, - sensitivity = colMeans(sens_sam), - elasticity = colMeans(elas_sam), - variance = NA, - se = NA, - cv = NA + parameter = pars, + sensitivity = colMeans(sens_sam), + elasticity = colMeans(elas_sam), + difference_method = NA, + variance = NA, + se = NA, + cv = NA ) + stats_temp$difference_method <- vapply(diff_sam, summarize_difference_method, + character(1)) + par_idx <- match(pars, colnames(samples)) for (i in seq_len(npar)) { diff --git a/R/utils-uncertainty.R b/R/utils-uncertainty.R index 0bbe8c3..ea44b80 100644 --- a/R/utils-uncertainty.R +++ b/R/utils-uncertainty.R @@ -59,7 +59,7 @@ validate_ipm_base <- function(ipm) { ) } -#' Warn if finite-difference perturbation may change parameter sign +#' Warn if finite-difference perturbation may cross parameter boundaries #' #' Internal helper used by sensitivity and uncertainty validation. #' @@ -67,25 +67,97 @@ validate_ipm_base <- function(ipm) { #' @param par_values Numeric vector of parameter values corresponding to #' `par_names`. #' @param delta Finite-difference perturbation size. +#' @param bounds Named list whose elements are numeric vectors of length two +#' defining parameter boundary constraints. #' #' @noRd -warn_delta_crossing <- function(par_names, par_values = NULL, delta) { +warn_boundary_crossing <- function(par_names, + par_values, + delta, + bounds = NULL) { + + ## ---- Standardize parameter values to a named list + + if (is.data.frame(par_values)) { + par_values <- as.list(par_values) + } else if (is.numeric(par_values)) { + par_values <- as.list(par_values) + } + + warn_pars <- character(0) + + ## ---- No bounds supplied (backwards-compatible behaviour) + + if (is.null(bounds)) { + + warn_pars <- par_names[ + vapply( + par_names, + function(p) any(abs(par_values[[p]]) <= delta), + logical(1) + ) + ] + + if (length(warn_pars) > 0) { + + warning( + "The following parameter(s) are within `delta` of zero and ", + "may change sign during perturbation: ", + paste(warn_pars, collapse = ", "), + ".\n\n", + "If these parameters are biologically constrained (e.g. ", + "probabilities or rates bounded below by zero), consider ", + "supplying the `bounds` argument so that one-sided finite ", + "differences can be used when necessary. Proceeding may ", + "result in failed or invalid subkernel construction or ", + "failed sensitivity/uncertainty calculations." + ) + + } + + return(invisible(NULL)) + + } + + ## ---- Bounds supplied + + for (par in names(bounds)) { + + # Skip parameters that are not being perturbed + if (!(par %in% par_names)) + next + + theta <- par_values[[par]] + + lower <- bounds[[par]][1] + upper <- bounds[[par]][2] + + crosses_lower <- any(theta - delta < lower) + crosses_upper <- any(theta + delta > upper) + + if (crosses_lower || crosses_upper) { + warn_pars <- c(warn_pars, par) + } - if (!is.null(par_values)) { - par_names <- unique(par_names[abs(par_values) <= delta]) } - if (length(par_names) > 0) { + warn_pars <- unique(warn_pars) + + if (length(warn_pars) > 0) { + warning( - "The following parameter(s) are within `delta` of zero and may change ", - "sign during perturbation: ", - paste(par_names, collapse = ", "), - ". This may result in failed or invalid subkernel construction or failed ", - "sensitivity/uncertainty calculations." + "Perturbation by `delta` crosses the specified parameter ", + "bounds for: ", + paste(warn_pars, collapse = ", "), + ".\n\n", + "One-sided finite differences should be used for these ", + "parameter(s) instead of the central difference." ) + } invisible(NULL) + } @@ -95,7 +167,7 @@ warn_delta_crossing <- function(par_names, par_values = NULL, delta) { #' Returns usable objects for the sensitivity function. #' #' @noRd -validate_ipm_sensitivity <- function(ipm, pars, kernels, delta) { +validate_ipm_sensitivity <- function(ipm, pars, kernels, delta, bounds) { base <- validate_ipm_base(ipm) @@ -182,21 +254,110 @@ validate_ipm_sensitivity <- function(ipm, pars, kernels, delta) { stop("`delta` must be > 0.") } - ## ---- check if delta changes sign on any parameters + ## ---- bounds + + if (!is.null(bounds)) { + if (!is.list(bounds)) { + stop("`bounds` must be NULL or a named list.") + } + if (is.null(names(bounds)) || any(names(bounds) == "")) { + stop("`bounds` must be a named list whose names correspond to parameter names.") + } + } + + bad_names <- setdiff(names(bounds), names(pars_all)) + if (length(bad_names) > 0) { + stop( + "Unknown parameter(s) supplied in `bounds`: ", + paste(bad_names, collapse = ", ") + ) + } + + bad_length <- vapply(bounds, length, integer(1)) != 2 + if (any(bad_length)) { + stop( + "Each element of `bounds` must be a numeric vector of length 2." + ) + } + + bad_numeric <- !vapply(bounds, is.numeric, logical(1)) + if (any(bad_numeric)) { + stop( + "Each element of `bounds` must be numeric." + ) + } + + has_na <- vapply( + bounds, + function(x) any(is.na(x)), + logical(1) ) + if (any(has_na)) { + stop( + "Bounds may not contain NA values. Use -Inf or Inf for unbounded parameters." + ) + } + + bad_order <- vapply( + bounds, + function(x) x[1] > x[2], + logical(1) ) + if (any(bad_order)) { + stop( + "Each lower bound must be less than or equal to its corresponding upper bound." + ) + } + + outside <- vapply( + names(bounds), + function(p) { + theta <- pars_all[[p]] + theta < bounds[[p]][1] || + theta > bounds[[p]][2] + }, + logical(1) + ) - warn_delta_crossing( + if (any(outside)) { + stop( + "Current parameter value(s) fall outside the supplied bounds: ", + paste(names(bounds)[outside], collapse = ", ") + ) + } + + too_large <- vapply( + names(bounds), + function(p) { + theta <- pars_all[[p]] + theta - delta < bounds[[p]][1] && + theta + delta > bounds[[p]][2] + }, + logical(1) + ) + + if (any(too_large)) { + warning( + "For the following parameter(s), `delta` exceeds the feasible ", + "parameter range: ", + paste(names(bounds)[too_large], collapse = ", "), + ". Numerical derivatives may be inaccurate." + ) + } + + warn_boundary_crossing( par_names = pars, par_values = as.numeric(pars_all[pars]), - delta = delta + delta = delta, + bounds = bounds ) ## return cleaned objects for main function list( - ipm = ipm, - pars = pars, + ipm = ipm, + pars = pars, pars_all = pars_all, - kernels = kernels, - delta = delta + kernels = kernels, + delta = delta, + bounds = bounds ) } @@ -231,7 +392,7 @@ new_ipmr_uncertainty <- function(ipm, lambdas, mod_uncert, params_uncert, vr_unc #' #' @noRd validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, - delta, cores) { + delta, bounds, cores) { base <- validate_ipm_base(ipm) @@ -339,18 +500,100 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, stop("`delta` must be a single numeric value > 0.") } - ## ---- check if delta changes sign on any parameters - problem_pars <- pars[ - vapply( - pars, - function(p) any(abs(samples[[p]]) <= delta), - logical(1) + ## ---- bounds + + if (!is.null(bounds)) { + if (!is.list(bounds)) { + stop("`bounds` must be NULL or a named list.") + } + if (is.null(names(bounds)) || any(names(bounds) == "")) { + stop("`bounds` must be a named list whose names correspond to parameter names.") + } + } + + bad_names <- setdiff(names(bounds), names(pars_all)) + if (length(bad_names) > 0) { + stop( + "Unknown parameter(s) supplied in `bounds`: ", + paste(bad_names, collapse = ", ") + ) + } + + bad_length <- vapply(bounds, length, integer(1)) != 2 + if (any(bad_length)) { + stop( + "Each element of `bounds` must be a numeric vector of length 2." ) - ] + } - warn_delta_crossing( - par_names = problem_pars, - delta = delta + bad_numeric <- !vapply(bounds, is.numeric, logical(1)) + if (any(bad_numeric)) { + stop( + "Each element of `bounds` must be numeric." + ) + } + + has_na <- vapply( + bounds, + function(x) any(is.na(x)), + logical(1) ) + if (any(has_na)) { + stop( + "Bounds may not contain NA values. Use -Inf or Inf for unbounded parameters." + ) + } + + bad_order <- vapply( + bounds, + function(x) x[1] > x[2], + logical(1) ) + if (any(bad_order)) { + stop( + "Each lower bound must be less than or equal to its corresponding upper bound." + ) + } + + outside <- vapply( + names(bounds), + function(p) { + theta <- pars_all[[p]] + theta < bounds[[p]][1] || + theta > bounds[[p]][2] + }, + logical(1) + ) + + if (any(outside)) { + stop( + "Current parameter value(s) fall outside the supplied bounds: ", + paste(names(bounds)[outside], collapse = ", ") + ) + } + + too_large <- vapply( + names(bounds), + function(p) { + theta <- pars_all[[p]] + theta - delta < bounds[[p]][1] && + theta + delta > bounds[[p]][2] + }, + logical(1) + ) + + if (any(too_large)) { + warning( + "For the following parameter(s), `delta` exceeds the feasible ", + "parameter range: ", + paste(names(bounds)[too_large], collapse = ", "), + ". Numerical derivatives may be inaccurate." + ) + } + + warn_boundary_crossing( + par_names = pars, + par_values = samples, + delta = delta, + bounds = bounds ) ## ---- cores @@ -389,6 +632,7 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, kernels = kernels, vr_table = vr_table, delta = delta, + bounds = bounds, cores = as.integer(cores) ) } diff --git a/tests/testthat/test-sensitivity.R b/tests/testthat/test-sensitivity.R index 0ecc455..3ba46fa 100644 --- a/tests/testthat/test-sensitivity.R +++ b/tests/testthat/test-sensitivity.R @@ -280,6 +280,198 @@ test_that("sensitivity errors on unknown parameter", { ) }) +test_that("sensitivity rejects non-list bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + sensitivity( + ipm, + kernels = c("P", "F"), + bounds = c(0, 1) + ), + "`bounds` must be NULL or a named list" + ) + +}) + +test_that("sensitivity rejects unnamed bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + sensitivity( + ipm, + kernels = c("P", "F"), + bounds = list(c(0, 1)) + ), + "must be a named list" + ) + +}) + +test_that("sensitivity rejects unknown parameter names in bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + sensitivity( + ipm, + kernels = c("P", "F"), + bounds = list( + fake_parameter = c(0, 1) + ) + ), + "Unknown parameter" + ) + +}) + +test_that("sensitivity rejects bounds of incorrect length 1", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + sensitivity( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c(0) + ) + ), + "length 2" + ) + +}) + +test_that("sensitivity rejects bounds of incorrect length 3", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + sensitivity( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c(0, 1, 2) + ) + ), + "length 2" + ) + +}) + +test_that("sensitivity rejects non-numeric bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + sensitivity( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c("a", "b") + ) + ), + "must be numeric" + ) + +}) + +test_that("sensitivity rejects missing values in bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + sensitivity( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c(NA, 1) + ) + ), + "may not contain NA" + ) + +}) + +test_that("sensitivity rejects reversed bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + sensitivity( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c(1, 0) + ) + ), + "lower bound" + ) + +}) + +test_that("sensitivity rejects parameter values outside supplied bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + sensitivity( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c(10, 20) + ) + ), + "outside the supplied bounds" + ) + +}) + +test_that("sensitivity accepts infinite bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_no_error( + + sensitivity( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c(-Inf, Inf), + f_s_int = c(0, Inf) + ) + ) + + ) + +}) + +test_that("sensitivity warns when delta exceeds feasible region", { + + ipm <- make_test_ipm("simple_di_det") + + pars <- parameters(ipm) + + expect_warning( + + sensitivity( + ipm, + kernels = c("P", "F"), + delta = 100, + bounds = list( + s_int = c(pars$s_int - 1, + pars$s_int + 1) + ) + ), + + "delta" + + ) + +}) + test_that("elasticity is consistent with sensitivity definition", { ipm <- make_test_ipm("simple_di_det") diff --git a/tests/testthat/test-uncertainty.R b/tests/testthat/test-uncertainty.R index 4645d2a..c27439a 100644 --- a/tests/testthat/test-uncertainty.R +++ b/tests/testthat/test-uncertainty.R @@ -290,6 +290,200 @@ test_that("uncertainty respects parameter subset", { c("s_int", "g_int"))) }) +test_that("uncertainty rejects non-list bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + uncertainty( + ipm, + kernels = c("P", "F"), + bounds = c(0, 1) + ), + "`bounds` must be NULL or a named list" + ) + +}) + +test_that("uncertainty rejects unnamed bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + uncertainty( + ipm, + kernels = c("P", "F"), + bounds = list(c(0, 1)) + ), + "must be a named list" + ) + +}) + +test_that("uncertainty rejects unknown parameter names in bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + uncertainty( + ipm, + kernels = c("P", "F"), + bounds = list( + fake_parameter = c(0, 1) + ) + ), + "Unknown parameter" + ) + +}) + +test_that("uncertainty rejects bounds of incorrect length 1", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + uncertainty( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c(0) + ) + ), + "length 2" + ) + +}) + +test_that("uncertainty rejects bounds of incorrect length 3", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + uncertainty( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c(0, 1, 2) + ) + ), + "length 2" + ) + +}) + +test_that("uncertainty rejects non-numeric bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + uncertainty( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c("a", "b") + ) + ), + "must be numeric" + ) + +}) + +test_that("uncertainty rejects missing values in bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + uncertainty( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c(NA, 1) + ) + ), + "may not contain NA" + ) + +}) + +test_that("uncertainty rejects reversed bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + uncertainty( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c(1, 0) + ) + ), + "lower bound" + ) + +}) + +test_that("uncertainty rejects parameter values outside supplied bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_error( + uncertainty( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c(10, 20) + ) + ), + "outside the supplied bounds" + ) + +}) + +test_that("uncertainty accepts infinite bounds", { + + ipm <- make_test_ipm("simple_di_det") + + expect_no_error( + + uncertainty( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c(-Inf, Inf), + f_s_int = c(0, Inf) + ) + ) + + ) + +}) + +test_that("uncertainty warns when delta exceeds feasible region", { + + ipm <- make_test_ipm("simple_di_det") + + pars <- parameters(ipm) + + expect_warning( + + uncertainty( + ipm, + kernels = c("P", "F"), + delta = 100, + bounds = list( + s_int = c(pars$s_int - 1, + pars$s_int + 1) + ) + ), + + "delta" + + ) + +}) + + + # Tests for numerical consistency ---------------------------------------------- test_that("uncertainty computes lambda variance", { diff --git a/vignettes/uncertainty-analysis.Rmd b/vignettes/uncertainty-analysis.Rmd index 62c4950..90ef9f7 100644 --- a/vignettes/uncertainty-analysis.Rmd +++ b/vignettes/uncertainty-analysis.Rmd @@ -100,7 +100,7 @@ A single numeric value giving the estimated model uncertainty, calculated as the ### `params_uncert` A dataframe summarizing parameter-level uncertainty. Each row corresponds to a single model parameter and includes: -- the mean sensitivity and elasticity across all sampled models +- the mean sensitivity and elasticity across all sampled models, and the finite-difference method used to calculate the sensitivities - the variance, standard error, and coefficient of variation of the sampled parameter values - the estimated contribution of the parameter to the uncertainty in λ - the proportion of total model uncertainty explained by the parameter From 4ad4d4d1da90f11836ec24a905f45bddf02c3ae3 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Wed, 5 Aug 2026 19:36:42 +0200 Subject: [PATCH 31/36] Add vr_table argument --- tests/testthat/test-uncertainty.R | 53 +++++++++++++++++++++++++------ 1 file changed, 43 insertions(+), 10 deletions(-) diff --git a/tests/testthat/test-uncertainty.R b/tests/testthat/test-uncertainty.R index c27439a..6cca0ef 100644 --- a/tests/testthat/test-uncertainty.R +++ b/tests/testthat/test-uncertainty.R @@ -293,11 +293,14 @@ test_that("uncertainty respects parameter subset", { test_that("uncertainty rejects non-list bounds", { ipm <- make_test_ipm("simple_di_det") + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" expect_error( uncertainty( ipm, kernels = c("P", "F"), + vr_table = vr, bounds = c(0, 1) ), "`bounds` must be NULL or a named list" @@ -308,12 +311,15 @@ test_that("uncertainty rejects non-list bounds", { test_that("uncertainty rejects unnamed bounds", { ipm <- make_test_ipm("simple_di_det") + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" expect_error( uncertainty( ipm, kernels = c("P", "F"), - bounds = list(c(0, 1)) + bounds = list(c(0, 1)), + vr_table = vr ), "must be a named list" ) @@ -323,6 +329,8 @@ test_that("uncertainty rejects unnamed bounds", { test_that("uncertainty rejects unknown parameter names in bounds", { ipm <- make_test_ipm("simple_di_det") + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" expect_error( uncertainty( @@ -330,7 +338,8 @@ test_that("uncertainty rejects unknown parameter names in bounds", { kernels = c("P", "F"), bounds = list( fake_parameter = c(0, 1) - ) + ), + vr_table = vr ), "Unknown parameter" ) @@ -340,6 +349,8 @@ test_that("uncertainty rejects unknown parameter names in bounds", { test_that("uncertainty rejects bounds of incorrect length 1", { ipm <- make_test_ipm("simple_di_det") + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" expect_error( uncertainty( @@ -347,7 +358,8 @@ test_that("uncertainty rejects bounds of incorrect length 1", { kernels = c("P", "F"), bounds = list( s_int = c(0) - ) + ), + vr_table = vr ), "length 2" ) @@ -357,6 +369,8 @@ test_that("uncertainty rejects bounds of incorrect length 1", { test_that("uncertainty rejects bounds of incorrect length 3", { ipm <- make_test_ipm("simple_di_det") + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" expect_error( uncertainty( @@ -364,7 +378,8 @@ test_that("uncertainty rejects bounds of incorrect length 3", { kernels = c("P", "F"), bounds = list( s_int = c(0, 1, 2) - ) + ), + vr_table = vr ), "length 2" ) @@ -374,6 +389,8 @@ test_that("uncertainty rejects bounds of incorrect length 3", { test_that("uncertainty rejects non-numeric bounds", { ipm <- make_test_ipm("simple_di_det") + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" expect_error( uncertainty( @@ -381,7 +398,8 @@ test_that("uncertainty rejects non-numeric bounds", { kernels = c("P", "F"), bounds = list( s_int = c("a", "b") - ) + ), + vr_table = vr ), "must be numeric" ) @@ -391,6 +409,8 @@ test_that("uncertainty rejects non-numeric bounds", { test_that("uncertainty rejects missing values in bounds", { ipm <- make_test_ipm("simple_di_det") + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" expect_error( uncertainty( @@ -398,7 +418,8 @@ test_that("uncertainty rejects missing values in bounds", { kernels = c("P", "F"), bounds = list( s_int = c(NA, 1) - ) + ), + vr_table = vr ), "may not contain NA" ) @@ -408,6 +429,8 @@ test_that("uncertainty rejects missing values in bounds", { test_that("uncertainty rejects reversed bounds", { ipm <- make_test_ipm("simple_di_det") + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" expect_error( uncertainty( @@ -415,7 +438,8 @@ test_that("uncertainty rejects reversed bounds", { kernels = c("P", "F"), bounds = list( s_int = c(1, 0) - ) + ), + vr_table = vr ), "lower bound" ) @@ -425,6 +449,8 @@ test_that("uncertainty rejects reversed bounds", { test_that("uncertainty rejects parameter values outside supplied bounds", { ipm <- make_test_ipm("simple_di_det") + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" expect_error( uncertainty( @@ -432,7 +458,8 @@ test_that("uncertainty rejects parameter values outside supplied bounds", { kernels = c("P", "F"), bounds = list( s_int = c(10, 20) - ) + ), + vr_table = vr ), "outside the supplied bounds" ) @@ -442,6 +469,8 @@ test_that("uncertainty rejects parameter values outside supplied bounds", { test_that("uncertainty accepts infinite bounds", { ipm <- make_test_ipm("simple_di_det") + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" expect_no_error( @@ -451,7 +480,8 @@ test_that("uncertainty accepts infinite bounds", { bounds = list( s_int = c(-Inf, Inf), f_s_int = c(0, Inf) - ) + ), + vr_table = vr ) ) @@ -461,6 +491,8 @@ test_that("uncertainty accepts infinite bounds", { test_that("uncertainty warns when delta exceeds feasible region", { ipm <- make_test_ipm("simple_di_det") + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" pars <- parameters(ipm) @@ -473,7 +505,8 @@ test_that("uncertainty warns when delta exceeds feasible region", { bounds = list( s_int = c(pars$s_int - 1, pars$s_int + 1) - ) + ), + vr_table = vr ), "delta" From dcecccb40dbcb3c3ab44fe3754c81edbf4637f66 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Wed, 5 Aug 2026 19:53:58 +0200 Subject: [PATCH 32/36] Add samples argument --- tests/testthat/test-uncertainty.R | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/tests/testthat/test-uncertainty.R b/tests/testthat/test-uncertainty.R index 6cca0ef..290603b 100644 --- a/tests/testthat/test-uncertainty.R +++ b/tests/testthat/test-uncertainty.R @@ -295,10 +295,12 @@ test_that("uncertainty rejects non-list bounds", { ipm <- make_test_ipm("simple_di_det") vr <- uncertainty(ipm, vr_table = "template") vr$vital_rate <- "test" + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) expect_error( uncertainty( ipm, + samples = samples, kernels = c("P", "F"), vr_table = vr, bounds = c(0, 1) @@ -313,10 +315,12 @@ test_that("uncertainty rejects unnamed bounds", { ipm <- make_test_ipm("simple_di_det") vr <- uncertainty(ipm, vr_table = "template") vr$vital_rate <- "test" + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) expect_error( uncertainty( ipm, + samples = samples, kernels = c("P", "F"), bounds = list(c(0, 1)), vr_table = vr @@ -331,10 +335,12 @@ test_that("uncertainty rejects unknown parameter names in bounds", { ipm <- make_test_ipm("simple_di_det") vr <- uncertainty(ipm, vr_table = "template") vr$vital_rate <- "test" + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) expect_error( uncertainty( ipm, + samples = samples, kernels = c("P", "F"), bounds = list( fake_parameter = c(0, 1) @@ -351,10 +357,12 @@ test_that("uncertainty rejects bounds of incorrect length 1", { ipm <- make_test_ipm("simple_di_det") vr <- uncertainty(ipm, vr_table = "template") vr$vital_rate <- "test" + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) expect_error( uncertainty( ipm, + samples = samples, kernels = c("P", "F"), bounds = list( s_int = c(0) @@ -371,10 +379,12 @@ test_that("uncertainty rejects bounds of incorrect length 3", { ipm <- make_test_ipm("simple_di_det") vr <- uncertainty(ipm, vr_table = "template") vr$vital_rate <- "test" + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) expect_error( uncertainty( ipm, + samples = samples, kernels = c("P", "F"), bounds = list( s_int = c(0, 1, 2) @@ -391,10 +401,12 @@ test_that("uncertainty rejects non-numeric bounds", { ipm <- make_test_ipm("simple_di_det") vr <- uncertainty(ipm, vr_table = "template") vr$vital_rate <- "test" + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) expect_error( uncertainty( ipm, + samples = samples, kernels = c("P", "F"), bounds = list( s_int = c("a", "b") @@ -411,10 +423,12 @@ test_that("uncertainty rejects missing values in bounds", { ipm <- make_test_ipm("simple_di_det") vr <- uncertainty(ipm, vr_table = "template") vr$vital_rate <- "test" + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) expect_error( uncertainty( ipm, + samples = samples, kernels = c("P", "F"), bounds = list( s_int = c(NA, 1) @@ -431,10 +445,12 @@ test_that("uncertainty rejects reversed bounds", { ipm <- make_test_ipm("simple_di_det") vr <- uncertainty(ipm, vr_table = "template") vr$vital_rate <- "test" + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) expect_error( uncertainty( ipm, + samples = samples, kernels = c("P", "F"), bounds = list( s_int = c(1, 0) @@ -451,10 +467,12 @@ test_that("uncertainty rejects parameter values outside supplied bounds", { ipm <- make_test_ipm("simple_di_det") vr <- uncertainty(ipm, vr_table = "template") vr$vital_rate <- "test" + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) expect_error( uncertainty( ipm, + samples = samples, kernels = c("P", "F"), bounds = list( s_int = c(10, 20) @@ -471,11 +489,13 @@ test_that("uncertainty accepts infinite bounds", { ipm <- make_test_ipm("simple_di_det") vr <- uncertainty(ipm, vr_table = "template") vr$vital_rate <- "test" + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) expect_no_error( uncertainty( ipm, + samples = samples, kernels = c("P", "F"), bounds = list( s_int = c(-Inf, Inf), @@ -493,6 +513,7 @@ test_that("uncertainty warns when delta exceeds feasible region", { ipm <- make_test_ipm("simple_di_det") vr <- uncertainty(ipm, vr_table = "template") vr$vital_rate <- "test" + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50))) pars <- parameters(ipm) @@ -500,6 +521,7 @@ test_that("uncertainty warns when delta exceeds feasible region", { uncertainty( ipm, + samples = samples, kernels = c("P", "F"), delta = 100, bounds = list( From e34a09b952e1debfa35013f792e902d4f19b15c1 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Wed, 5 Aug 2026 19:54:51 +0200 Subject: [PATCH 33/36] Add functionality for multiple finite-difference methods in sensitivity calculations --- man/sensitivity.Rd | 84 ++++++++++++++++++++++++++++++++++++++++------ man/uncertainty.Rd | 35 +++++++++++++++++-- 2 files changed, 106 insertions(+), 13 deletions(-) diff --git a/man/sensitivity.Rd b/man/sensitivity.Rd index e2501a4..4ba1953 100644 --- a/man/sensitivity.Rd +++ b/man/sensitivity.Rd @@ -4,7 +4,7 @@ \alias{sensitivity} \title{Sensitivity and elasticity analysis of population growth rate} \usage{ -sensitivity(ipm, pars = NULL, kernels, delta = 1e-04) +sensitivity(ipm, pars = NULL, kernels, delta = 1e-04, bounds = NULL) } \arguments{ \item{ipm}{An \code{ipmr} IPM object.} @@ -16,26 +16,90 @@ all parameters in \code{ipm}.} major order.} \item{delta}{A numeric scalar specifying the perturbation size used for -numerical differentiation. Default is 1e-4.} +numerical differentiation. Default is 1e-4. +@param bounds An optional named list specifying lower and/or upper bounds + for model parameters when calculating numerical sensitivities. Each element + should be a numeric vector of length two giving the lower and upper bounds, + respectively (e.g., \code{list(surv = c(0, 1), recruit = c(0, Inf))}). + Parameters omitted from this list are assumed to be unbounded and are + perturbed using the central finite difference. For parameters with + specified bounds, the function automatically uses a forward or backward + finite difference whenever perturbation by \code{delta} would cross a + boundary; otherwise, the central finite difference is used.} } \value{ -A named list with two elements: +A named list with three elements: \itemize{ - \item \code{sensitivity}: Numeric vector of sensitivities of - \eqn{\lambda} to each parameter. - \item \code{elasticity}: Numeric vector of elasticities of - \eqn{\lambda} to each parameter. + \item \code{sensitivity}: A named numeric vector giving the numerical + sensitivities of the dominant eigenvalue (\eqn{\lambda}) to each + parameter. + \item \code{elasticity}: A named numeric vector giving the elasticities + of the dominant eigenvalue (\eqn{\lambda}) to each parameter. + \item \code{difference_method}: A named character vector indicating the + finite-difference scheme used to calculate the sensitivity for each + parameter. This component is provided to facilitate reproducibility and + to identify parameters whose sensitivities were evaluated using one-sided + finite differences because of user-specified parameter bounds. Values are + \code{"central"}, \code{"forward"}, or \code{"backward"}. } } \description{ Computes numerical sensitivities and elasticities of the -dominant eigenvalue (\eqn{\lambda}) with respect to model parameters. +dominant eigenvalue (\eqn{\lambda}) with respect to model parameters using +finite differences. } \details{ -This function uses a central difference approximation to estimate the -sensitivity of \eqn{\lambda} to each parameter in \code{pars}. +This function uses perturbation to estimate the sensitivity of \eqn{\lambda} +to each parameter in \code{pars}. Currently implemented for deterministic ("det"), density-independent ("di") IPMs constructed with \code{ipmr}. Stochastic ("stoch") and density-dependent ("dd") IPMs are not supported. + +## Finite differences and bounded parameters + +Sensitivities are estimated numerically by perturbing each parameter by +\code{delta}. By default, a central finite difference is used for all +parameters. However, some demographic parameters are biologically constrained +(e.g., probabilities bounded between 0 and 1, or fecundity parameters bounded +below by 0). If perturbing a parameter by \code{delta} would cross one of +these user-specified bounds, the function automatically computes the +appropriate one-sided finite difference (forward or backward) rather than +evaluating the model outside the feasible parameter space. + +One-sided finite differences approximate the same local derivative as the +central finite difference, but are used only when a symmetric perturbation +would leave the biologically feasible parameter space. Consequently, +sensitivities calculated using one-sided and central finite differences are +directly comparable and may be used together in downstream uncertainty +analyses. + +Supplying \code{bounds} is recommended whenever model parameters are subject +to biological constraints and numerical perturbation could otherwise produce +invalid parameter values or biologically impossible kernels. +} +\examples{ +# Probability bounded between 0 and 1 +# bounds <- list( +# surv = c(0, 1) +# ) + +# Recruitment bounded below by zero +# bounds <- list( +# recruit = c(0, Inf) +# ) + +# Multiple constrained parameters +# bounds <- list( +# surv = c(0, 1), +# recruit = c(0, Inf), +# fec = c(0, Inf) +# ) + +# Compute sensitivities using defined parameter boundaries +# sens <- sensitivity( +# ipm, +# kernels = c("P", "F"), +# bounds = bounds ) + } diff --git a/man/uncertainty.Rd b/man/uncertainty.Rd index 8a2ebb7..2f5e9d3 100644 --- a/man/uncertainty.Rd +++ b/man/uncertainty.Rd @@ -11,6 +11,7 @@ uncertainty( kernels, vr_table, delta = 1e-04, + bounds = NULL, cores = 1 ) } @@ -32,9 +33,34 @@ parameter in \code{pars} must be present exactly once in this table.} \item{delta}{A numeric scalar specifying the perturbation size used for numerical differentiation when calculating sensitivities. Default is 1e-4.} +\item{bounds}{An optional named list specifying lower and/or upper bounds +for model parameters when calculating numerical sensitivities. Each element +should be a numeric vector of length two giving the lower and upper bounds, +respectively (e.g., \code{list(surv = c(0, 1), recruit = c(0, Inf))}). +Parameters omitted from this list are assumed to be unbounded and are +perturbed using the central finite difference. For parameters with +specified bounds, the function automatically uses a forward or backward +finite difference whenever perturbation by \code{delta} would cross a +boundary; otherwise, the central finite difference is used.} + \item{cores}{An integer specifying the number of cores between which the computations will be distributed. Defaults to 1.} } +\value{ +An object of class \code{"ipmr_uncertainty"} containing: +\itemize{ + \item \code{ipm}: The complete IPM object as specified in \code{ipm}. + \item \code{lambdas}: Population growth rates (\eqn{\lambda}) calculated + for each sampled parameter set. + \item \code{mod_uncert}: The variance in \eqn{\lambda} across all sampled + parameter sets. + \item \code{params_uncert}: A data frame summarizing parameter-level + sensitivities, elasticities, finite-difference methods, sampling + variance, coefficients of variation, and contributions to uncertainty. + \item \code{vr_uncert}: A data frame summarizing variance contributions + aggregated by user-defined vital rate, if given. +} +} \description{ Performs uncertainty analysis of IPMs. } @@ -54,9 +80,12 @@ categories automatically. Users must explicitly define this categorization. A template for \code{vr_table} can be generated by calling the function with \code{vr_table = "template"}. -Sensitivities within each sampled model are computed using a central -difference approximation (see \code{\link{sensitivity}}). The \code{delta} -argument controls the perturbation size used in this approximation. +The finite-difference scheme used in \code{\link{sensitivity}} is used when +estimating parameter sensitivities for the uncertainty decomposition. +Consequently, supplying biological parameter bounds ensures that +sensitivities are evaluated only within the feasible parameter space +throughout the uncertainty analysis. The \code{delta} argument controls the +perturbation size used in the sensitivity approximation. On systems with multiple cores, computations can be parallelized using the \code{cores} argument. The function does not automatically infer scheduler From f02b91eae4389ddfbe96f424cef35fb30f260ea0 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Fri, 7 Aug 2026 15:44:23 +0200 Subject: [PATCH 34/36] Prepare new version with multiple methods for sensitivity calculation --- .Rbuildignore | 2 + .gitignore | 2 + R/uncertainty.R | 1 + R/utils-uncertainty.R | 90 +++++++++- tests/testthat/test-sensitivity.R | 118 +++++++++++- tests/testthat/test-uncertainty.R | 290 +++++++++++++++++++++++++++++- 6 files changed, 494 insertions(+), 9 deletions(-) diff --git a/.Rbuildignore b/.Rbuildignore index a2c77ad..d11bb5a 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -28,3 +28,5 @@ ^tests/vignettes$ ^tests/vignettes/*$ ^collaborations/*$ +^doc$ +^Meta$ diff --git a/.gitignore b/.gitignore index 79cb8b1..7d8d095 100644 --- a/.gitignore +++ b/.gitignore @@ -31,3 +31,5 @@ src-i386/* paper/Drafts/submitted_revision* paper/Drafts/*.pdf collaborations/* +/doc/ +/Meta/ diff --git a/R/uncertainty.R b/R/uncertainty.R index 66cda32..7922095 100644 --- a/R/uncertainty.R +++ b/R/uncertainty.R @@ -101,6 +101,7 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, )) } + ## ---- Input checks val <- validate_ipm_uncertainty(ipm, pars, samples, kernels, vr_table, delta, bounds, cores) diff --git a/R/utils-uncertainty.R b/R/utils-uncertainty.R index ea44b80..38354bd 100644 --- a/R/utils-uncertainty.R +++ b/R/utils-uncertainty.R @@ -324,6 +324,42 @@ validate_ipm_sensitivity <- function(ipm, pars, kernels, delta, bounds) { ) } + for (par in names(bounds)) { + + theta <- pars_all[[par]] + lower <- bounds[[par]][1] + upper <- bounds[[par]][2] + + ## Parameter itself must lie inside bounds + if (theta < lower || theta > upper) { + stop( + "Parameter '", par, + "' has value ", theta, + ", which lies outside the supplied bounds [", + lower, ", ", upper, "].", + call. = FALSE + ) + } + + if (is.finite(lower) && + is.finite(upper) && + delta > (upper - lower)) { + + stop( + sprintf( + "Parameter '%s' has bounds [%g, %g], giving a feasible interval of %g, which is smaller than delta (%g). Choose a smaller value of delta.", + par, + lower, + upper, + upper - lower, + delta + ), + call. = FALSE + ) + } + + } + too_large <- vapply( names(bounds), function(p) { @@ -345,7 +381,7 @@ validate_ipm_sensitivity <- function(ipm, pars, kernels, delta, bounds) { warn_boundary_crossing( par_names = pars, - par_values = as.numeric(pars_all[pars]), + par_values = pars_all[pars], delta = delta, bounds = bounds ) @@ -570,6 +606,42 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, ) } + for (par in names(bounds)) { + + theta <- pars_all[[par]] + lower <- bounds[[par]][1] + upper <- bounds[[par]][2] + + ## Parameter itself must lie inside bounds + if (theta < lower || theta > upper) { + stop( + "Parameter '", par, + "' has value ", theta, + ", which lies outside the supplied bounds [", + lower, ", ", upper, "].", + call. = FALSE + ) + } + + if (is.finite(lower) && + is.finite(upper) && + delta > (upper - lower)) { + + stop( + sprintf( + "Parameter '%s' has bounds [%g, %g], giving a feasible interval of %g, which is smaller than delta (%g). Choose a smaller value of delta.", + par, + lower, + upper, + upper - lower, + delta + ), + call. = FALSE + ) + } + + } + too_large <- vapply( names(bounds), function(p) { @@ -589,6 +661,22 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, ) } + outside_samples <- vapply( + names(bounds), + function(p) { + vals <- samples[[p]] + any(vals < bounds[[p]][1] | vals > bounds[[p]][2]) + }, + logical(1) + ) + + if (any(outside_samples)) { + stop( + "Sampled value(s) fall outside the supplied bounds for: ", + paste(names(bounds)[outside_samples], collapse = ", ") + ) + } + warn_boundary_crossing( par_names = pars, par_values = samples, diff --git a/tests/testthat/test-sensitivity.R b/tests/testthat/test-sensitivity.R index 3ba46fa..b5ae728 100644 --- a/tests/testthat/test-sensitivity.R +++ b/tests/testthat/test-sensitivity.R @@ -209,6 +209,8 @@ ipm <- make_test_ipm("simple_di_det") ipm_dd <- make_test_ipm("simple_dd_det") ipm_stoch <- make_test_ipm("simple_di_stoch") +# Tests for structural correctness --------------------------------------------- + test_that("sensitivity errors on missing ipm", { expect_error( sensitivity(pars = NULL, kernels = NULL, delta = 1e-4), @@ -448,30 +450,107 @@ test_that("sensitivity accepts infinite bounds", { }) -test_that("sensitivity warns when delta exceeds feasible region", { +test_that("sensitivity errors when delta exceeds feasible interval", { ipm <- make_test_ipm("simple_di_det") pars <- parameters(ipm) - expect_warning( + expect_error( sensitivity( ipm, kernels = c("P", "F"), - delta = 100, + delta = 1, bounds = list( - s_int = c(pars$s_int - 1, - pars$s_int + 1) + s_int = c( + pars$s_int - 0.25, + pars$s_int + 0.25 + ) ) ), - "delta" + "Choose a smaller value of delta" ) }) +test_that("sensitivity returns difference method", { + + ipm <- make_test_ipm("simple_di_det") + + res <- sensitivity( + ipm, + kernels = c("P", "F") + ) + + expect_true("difference_method" %in% names(res)) +}) + +test_that("sensitivity records forward differences for bounded parameters", { + + ipm <- make_test_ipm("simple_di_det") + + pars <- parameters(ipm) + + res <- sensitivity( + ipm, + kernels = c("P", "F"), + bounds = list(s_int = c(pars$s_int, Inf)) + ) + + expect_warning( + sensitivity( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c(pars$s_int, Inf) + ) + ), + "parameter bounds" + ) + + expect_equal( + unname(res$difference_method["s_int"]), + "forward" + ) + +}) + +test_that("sensitivity records backward differences for bounded parameters", { + + ipm <- make_test_ipm("simple_di_det") + + pars <- parameters(ipm) + + res <- sensitivity( + ipm, + kernels = c("P", "F"), + bounds = list(s_int = c(-Inf, pars$s_int)) + ) + + expect_warning( + sensitivity( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c(pars$s_int, Inf) + ) + ), + "parameter bounds" + ) + + expect_equal( + unname(res$difference_method["s_int"]), + "backward" + ) + +}) + + +# Tests for numerical consistency ---------------------------------------------- + test_that("elasticity is consistent with sensitivity definition", { ipm <- make_test_ipm("simple_di_det") @@ -522,3 +601,30 @@ test_that("unknown kernel errors clearly", { "Unknown kernel" ) }) + +test_that("bounded sensitivity matches unbounded sensitivity when parameter is away from bounds", { + + ipm <- make_test_ipm("simple_di_det") + + pars <- parameters(ipm) + + res1 <- sensitivity( + ipm, + kernels = c("P", "F") + ) + + res2 <- sensitivity( + ipm, + kernels = c("P", "F"), + bounds = list( + s_int = c(0, Inf) + ) + ) + + expect_equal( + res1$sensitivity, + res2$sensitivity, + tolerance = 1e-6 + ) + +}) diff --git a/tests/testthat/test-uncertainty.R b/tests/testthat/test-uncertainty.R index 290603b..de1143a 100644 --- a/tests/testthat/test-uncertainty.R +++ b/tests/testthat/test-uncertainty.R @@ -209,6 +209,17 @@ ipm <- make_test_ipm("simple_di_det") ipm_dd <- make_test_ipm("simple_dd_det") ipm_stoch <- make_test_ipm("simple_di_stoch") +samples <- as.data.frame( + lapply(parameters(ipm), function(x) rep(x, 2)) +) + +vr <- data.frame( + parameter = names(parameters(ipm)), + vital_rate = "test" +) + +ker <- c("P", "F") + # Tests for structural correctness --------------------------------------------- test_that("uncertainty errors on missing IPM", { @@ -517,13 +528,13 @@ test_that("uncertainty warns when delta exceeds feasible region", { pars <- parameters(ipm) - expect_warning( + expect_error( uncertainty( ipm, samples = samples, kernels = c("P", "F"), - delta = 100, + delta = 5, bounds = list( s_int = c(pars$s_int - 1, pars$s_int + 1) @@ -537,6 +548,244 @@ test_that("uncertainty warns when delta exceeds feasible region", { }) +test_that("uncertainty returns difference method", { + + ipm <- make_test_ipm("simple_di_det") + + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 10))) + + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" + + res <- uncertainty( + ipm, + samples = samples, + kernels = c("P", "F"), + vr_table = vr + ) + + expect_true("difference_method" %in% names(res$params_uncert)) +}) + +test_that("uncertainty records forward differences for bounded parameters", { + + ipm <- make_test_ipm("simple_di_det") + + pars <- parameters(ipm) + + samples <- as.data.frame( + lapply(pars, function(x) rep(x, 10)) + ) + + samples$s_int[] <- pars$s_int + + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" + + res <- uncertainty( + ipm, + samples = samples, + kernels = c("P", "F"), + vr_table = vr, + bounds = list(s_int = c(pars$s_int, Inf)) + ) + + expect_equal( + res$params_uncert$difference_method[ + res$params_uncert$parameter == "s_int" + ], + "forward" + ) + +}) + +test_that("uncertainty records backward differences for bounded parameters", { + + ipm <- make_test_ipm("simple_di_det") + + pars <- parameters(ipm) + + samples <- as.data.frame( + lapply(pars, function(x) rep(x, 10)) + ) + + samples$s_int[] <- pars$s_int + + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" + + res <- uncertainty( + ipm, + samples = samples, + kernels = c("P", "F"), + vr_table = vr, + bounds = list(s_int = c(-Inf, pars$s_int)) + ) + + expect_equal( + res$params_uncert$difference_method[ + res$params_uncert$parameter == "s_int" + ], + "backward" + ) + +}) + +test_that("uncertainty records mixed difference methods", { + + ipm <- make_test_ipm("simple_di_det") + + pars <- parameters(ipm) + + samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 10))) + + ## Half of the samples sit on the lower bound, half do not + samples$s_int <- c(rep(pars$s_int, 5), rep(pars$s_int + 0.5, 5)) + + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" + + res <- uncertainty( + ipm, + samples = samples, + kernels = c("P", "F"), + vr_table = vr, + bounds = list( + s_int = c(pars$s_int, Inf) + ) + ) + + expect_match( + res$params_uncert$difference_method[ + res$params_uncert$parameter == "s_int" + ], + "^mixed" + ) + +}) + +test_that("uncertainty errors when delta exceeds feasible interval", { + + ipm <- make_test_ipm("simple_di_det") + + pars <- parameters(ipm) + + samples <- as.data.frame( + lapply(pars, function(x) rep(x, 5)) + ) + + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" + + expect_error( + + uncertainty( + ipm, + samples = samples, + kernels = c("P", "F"), + vr_table = vr, + delta = 1, + bounds = list( + s_int = c( + pars$s_int - 0.25, + pars$s_int + 0.25 + ) + ) + ), + + "Choose a smaller value of delta" + + ) + +}) + +test_that("uncertainty errors when sampled values cross parameter bounds", { + + ipm <- make_test_ipm("simple_di_det") + + pars <- parameters(ipm) + + samples <- as.data.frame( + lapply(pars, function(x) rep(x, 10)) + ) + + samples$s_int <- c( + rep(pars$s_int, 5), + rep(pars$s_int - 0.05, 5) + ) + + vr <- data.frame( + parameter = names(parameters(ipm)), + vital_rate = "test" ) + + expect_error( + + uncertainty( + ipm, + samples = samples, + kernels = c("P","F"), + vr_table = vr, + delta = 0.1, + bounds = list( + s_int = c(pars$s_int, Inf) + ) + ), + + "supplied bounds" + ) + +}) + +test_that("uncertainty rejects sampled values outside bounds", { + + pars <- parameters(ipm) + + samples <- as.data.frame( + lapply(pars, function(x) rep(x, 10)) + ) + + samples$s_int[1] <- pars$s_int - 0.5 + + expect_error( + + uncertainty( + ipm, + samples = samples, + kernels = c("P", "F"), + vr_table = vr, + bounds = list( + s_int = c(pars$s_int, Inf) + ) + ), + + "Sampled value\\(s\\) fall outside the supplied bounds" + + ) + +}) + +test_that("uncertainty template mode returns expected table", { + + ipm <- make_test_ipm("simple_di_det") + + vr <- uncertainty(ipm, vr_table = "template") + + expect_s3_class(vr, "data.frame") + + expect_equal( + names(vr), + c("parameter", "vital_rate") + ) + + expect_equal( + vr$parameter, + names(parameters(ipm)) + ) + + expect_true(all(is.na(vr$vital_rate))) + +}) + # Tests for numerical consistency ---------------------------------------------- @@ -586,6 +835,43 @@ test_that("uncertainty sensitivities match sensitivity()", { ) }) +test_that("bounded uncertainty sensitivities match bounded sensitivity()", { + + ipm <- make_test_ipm("simple_di_det") + + pars <- parameters(ipm) + + samples <- as.data.frame(lapply(pars, function(x) rep(x, 25))) + + vr <- uncertainty(ipm, vr_table = "template") + vr$vital_rate <- "test" + + bounds <- list( + s_int = c(pars$s_int, Inf) + ) + + unc <- uncertainty( + ipm, + samples = samples, + kernels = c("P", "F"), + vr_table = vr, + bounds = bounds + ) + + sens <- sensitivity( + ipm, + kernels = c("P", "F"), + bounds = bounds + ) + + expect_equal( + unname(unc$params_uncert$sensitivity), + unname(sens$sensitivity[unc$params_uncert$parameter]), + tolerance = 1e-6 + ) + +}) + test_that("vr_table join works correctly", { ipm <- make_test_ipm("simple_di_det") From 10fb95e7ec1bf119b29509e11f013872a6f4e913 Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Mon, 10 Aug 2026 15:18:49 +0200 Subject: [PATCH 35/36] Change kernels argument to mega_mat and update documentation to account for changes. --- DESCRIPTION | 1 + R/sensitivity.R | 43 +++--- R/uncertainty.R | 39 ++++-- R/utils-uncertainty.R | 121 ++++++---------- man/plot.ipmr_uncertainty.Rd | 10 +- man/print.ipmr_uncertainty.Rd | 2 +- man/sensitivity.Rd | 62 +++++---- man/uncertainty.Rd | 49 ++++--- tests/testthat/test-plot-uncertainty.R | 10 +- tests/testthat/test-sensitivity.R | 82 ++++------- tests/testthat/test-uncertainty.R | 85 ++++-------- vignettes/uncertainty-analysis.Rmd | 183 +++++++++++++++++++++++-- 12 files changed, 403 insertions(+), 284 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index a63a1f5..117fcd4 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -83,3 +83,4 @@ LinkingTo: URL: https://padrinoDB.github.io/ipmr/, https://github.com/padrinoDB/ipmr BugReports: https://github.com/padrinoDB/ipmr/issues Config/roxygen2/version: 8.0.0 +Config/roxygen2/markdown: TRUE diff --git a/R/sensitivity.R b/R/sensitivity.R index b168b1b..69c1060 100644 --- a/R/sensitivity.R +++ b/R/sensitivity.R @@ -7,19 +7,26 @@ #' @param ipm An \code{ipmr} IPM object. #' @param pars A character vector of parameter names to perturb. Defaults to #' all parameters in \code{ipm}. -#' @param kernels A character vector specifying the kernel structure in row -#' major order. +#' @param mega_mat A character vector or string specifying the structure of the +#' iteration matrix passed to \code{\link{make_iter_kernel}}. For simple IPMs, +#' \code{mega_mat} may be omitted or set to \code{NULL}. For general IPMs, +#' \code{mega_mat} must be supplied and specify the arrangement of matrix +#' blocks in row-major order. Individual blocks may contain sums of subkernels +#' when multiple kernels contribute to the same matrix block. For example, +#' \code{"c(P + F, seedbank,...)"} specifies an iteration matrix in which \code{P} +#' and \code{F} are summed in the same block. See \code{\link{make_iter_kernel}} +#' for additional details on specifying \code{mega_mat}. #' @param delta A numeric scalar specifying the perturbation size used for #' numerical differentiation. Default is 1e-4. -#' @param bounds An optional named list specifying lower and/or upper bounds -#' for model parameters when calculating numerical sensitivities. Each element -#' should be a numeric vector of length two giving the lower and upper bounds, -#' respectively (e.g., \code{list(surv = c(0, 1), recruit = c(0, Inf))}). -#' Parameters omitted from this list are assumed to be unbounded and are -#' perturbed using the central finite difference. For parameters with -#' specified bounds, the function automatically uses a forward or backward -#' finite difference whenever perturbation by \code{delta} would cross a -#' boundary; otherwise, the central finite difference is used. +#' @param bounds An optional named list specifying lower and/or upper bounds +#' for model parameters when calculating numerical sensitivities. Each element +#' should be a numeric vector of length two giving the lower and upper bounds, +#' respectively (e.g., \code{list(surv = c(0, 1), recruit = c(0, Inf))}). +#' Parameters omitted from this list are assumed to be unbounded and are +#' perturbed using the central finite difference. For parameters with +#' specified bounds, the function automatically uses a forward or backward +#' finite difference whenever perturbation by \code{delta} would cross a +#' boundary; otherwise, the central finite difference is used. #' #'@details #' This function uses perturbation to estimate the sensitivity of \eqn{\lambda} @@ -29,6 +36,10 @@ #' ("di") IPMs constructed with \code{ipmr}. Stochastic ("stoch") and #' density-dependent ("dd") IPMs are not supported. #' +#' For general IPMs, \code{mega_mat} should describe the complete block +#' structure of the iteration matrix. The syntax accepted by \code{make_iter_kernel} +#' is passed through without modification. +#' #' ## Finite differences and bounded parameters #' #' Sensitivities are estimated numerically by perturbing each parameter by @@ -91,17 +102,17 @@ #' # bounds = bounds ) #' #' @export -sensitivity <- function(ipm, pars = NULL, kernels, delta = 1e-4, +sensitivity <- function(ipm, pars = NULL, mega_mat = NULL, delta = 1e-4, bounds = NULL) { ## ---- Input checks - val <- validate_ipm_sensitivity(ipm, pars, kernels, delta, bounds) + val <- validate_ipm_sensitivity(ipm, pars, mega_mat, delta, bounds) ipm <- val$ipm pars <- val$pars pars_all <- val$pars_all - kernels <- val$kernels + mega_mat <- val$mega_mat delta <- val$delta bounds <- val$bounds @@ -111,7 +122,7 @@ sensitivity <- function(ipm, pars = NULL, kernels, delta = 1e-4, spar <- numeric(npar) diff_method <- character(npar) - kernel <- make_iter_kernel(ipm, mega_mat = kernels) + kernel <- make_iter_kernel(ipm, mega_mat = mega_mat) if (!is.list(kernel) || is.null(kernel[[1]])) { stop("Kernel construction failed.") @@ -128,7 +139,7 @@ sensitivity <- function(ipm, pars = NULL, kernels, delta = 1e-4, parameters(ipm_tmp) <- vpar ker <- ipm_tmp %>% make_ipm(iterate = FALSE) %>% - make_iter_kernel(mega_mat = kernels) + make_iter_kernel(mega_mat = mega_mat) dominant_eigen(ker[[1]]) } diff --git a/R/uncertainty.R b/R/uncertainty.R index 7922095..43f7ede 100644 --- a/R/uncertainty.R +++ b/R/uncertainty.R @@ -5,12 +5,19 @@ #' @param ipm An \code{ipmr} IPM object. #' @param pars A character vector of parameter names on which to perform the #' uncertainty analysis. Defaults to all parameters in \code{ipm}. -#' @param samples A data.frame containing sampled parameter values. -#' @param kernels A character vector specifying the kernel structure in row -#' major order. -#' @param vr_table A data.frame with columns 'parameter' and 'vital_rate' -#' mapping each parameter to a user-defined vital rate category. Every -#' parameter in \code{pars} must be present exactly once in this table. +#' @param samples A dataframe containing sampled parameter values. +#' @param mega_mat A character vector or string specifying the structure of the +#' iteration matrix passed to \code{\link{make_iter_kernel}}. For simple IPMs, +#' \code{mega_mat} may be omitted or set to \code{NULL}. For general IPMs, +#' \code{mega_mat} must be supplied and specify the arrangement of matrix +#' blocks in row-major order. Individual blocks may contain sums of subkernels +#' when multiple kernels contribute to the same matrix block. For example, +#' \code{"c(P + F, seedbank,...)"} specifies an iteration matrix in which \code{P} +#' and \code{F} are summed in the same block. See \code{\link{make_iter_kernel}} +#' for additional details on specifying \code{mega_mat}. +#' @param vr_table An optional dataframe with columns 'parameter' and +#' 'vital_rate' mapping each parameter to a user-defined vital rate category. +#' Every parameter in \code{pars} must be present exactly once in this table. #' @param delta A numeric scalar specifying the perturbation size used for #' numerical differentiation when calculating sensitivities. Default is 1e-4. #' @param bounds An optional named list specifying lower and/or upper bounds @@ -81,14 +88,20 @@ #' # vr$vital_rate <- c("survival", "growth", "reproduction", ...) #' #' # Run uncertainty analysis -#' # out <- uncertainty(ipm, pars, samples, kernels, vr_table = vr) +#' # mega_mat <- c("P", "F") +#' # out <- uncertainty(ipm, pars, samples, mega_mat, vr_table = vr) +#' +#' # For multiple kernels contributing to the same matrix block: +#' # mega_mat <- "c(P + repr, SB1_germ, SB2_germ, +#' # enter_SB1, SB1_SB1, SB2_SB1, +#' # enter_SB2, SB1_SB2, SB2_SB2)" #' #' @importFrom dplyr bind_rows left_join rename group_by summarize join_by #' @importFrom stats var sd cov #' @importFrom reshape2 melt #' #' @export -uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, +uncertainty <- function(ipm, pars = NULL, samples, mega_mat = NULL, vr_table, delta = 1e-4, bounds = NULL, cores = 1) { if (isTRUE(vr_table == "template")) { @@ -103,13 +116,13 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, ## ---- Input checks - val <- validate_ipm_uncertainty(ipm, pars, samples, kernels, + val <- validate_ipm_uncertainty(ipm, pars, samples, mega_mat, vr_table, delta, bounds, cores) ipm <- val$ipm pars <- val$pars samples <- val$samples - kernels <- val$kernels + mega_mat <- val$mega_mat vr_table <- val$vr_table delta <- val$delta bounds <- val$bounds @@ -129,7 +142,7 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, parallel::clusterExport( cluster, - c("base_ipm", "samples", "pars", "kernels", "delta", "bounds"), + c("base_ipm", "samples", "pars", "mega_mat", "delta", "bounds"), envir = environment() ) @@ -141,7 +154,7 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, ipm_temp <- make_ipm(ipm_temp, iterate = FALSE) - ker_base <- make_iter_kernel(ipm_temp, mega_mat = kernels) + ker_base <- make_iter_kernel(ipm_temp, mega_mat = mega_mat) mat_base <- ker_base[[1]] lambda_base <- Re(eigen(mat_base, only.values = TRUE)$values[1]) @@ -155,7 +168,7 @@ uncertainty <- function(ipm, pars = NULL, samples, kernels, vr_table, parameters(ipm_tmp) <- pars_tmp ipm_tmp <- make_ipm(ipm_tmp, iterate = FALSE) ker <- make_iter_kernel(ipm_tmp, - mega_mat = kernels) + mega_mat = mega_mat) Re(eigen(ker[[1]], only.values = TRUE)$values[1]) } diff --git a/R/utils-uncertainty.R b/R/utils-uncertainty.R index 38354bd..bc0eaff 100644 --- a/R/utils-uncertainty.R +++ b/R/utils-uncertainty.R @@ -167,7 +167,7 @@ warn_boundary_crossing <- function(par_names, #' Returns usable objects for the sensitivity function. #' #' @noRd -validate_ipm_sensitivity <- function(ipm, pars, kernels, delta, bounds) { +validate_ipm_sensitivity <- function(ipm, pars, mega_mat, delta, bounds) { base <- validate_ipm_base(ipm) @@ -201,47 +201,27 @@ validate_ipm_sensitivity <- function(ipm, pars, kernels, delta, bounds) { warning("Zero-valued parameters detected; elasticity may be undefined.") } - ## ---- kernels + ## ---- mega_mat + is_simple <- grepl("simple", class(ipm)[1]) - if (missing(kernels) || is.null(kernels)) { - stop("`kernels` must be provided.") - } - - if (!is.character(kernels) || length(kernels) == 0) { - stop("`kernels` must be a non-empty character vector.") - } + if (is_simple) { - # Check for duplicates - dup_kerns <- unique(kernels[duplicated(kernels)]) - - if (length(dup_kerns) > 0) { - warning( - "Duplicate kernel name(s) detected in `kernels`: ", - paste(dup_kerns, collapse = ", "), - ". This may lead to incorrect kernel construction if not intentional." - ) - } - - available_kernels <- names(ipm$sub_kernels) + # mega_mat may be NULL for simple IPMs + if (!is.null(mega_mat) && + (!is.character(mega_mat) || length(mega_mat) == 0)) { + stop("`mega_mat` must be a non-empty character vector or string, or NULL.") + } - if (!all(kernels %in% available_kernels)) { - bad <- setdiff(kernels, available_kernels) - stop( - "Unknown kernel name(s): ", - paste(bad, collapse = ", "), - ". Available kernels are: ", - paste(available_kernels, collapse = ", ") - ) - } + } else { - n_expected <- length(ipm$sub_kernels) + # General IPMs require an explicitly specified mega_mat + if (is.null(mega_mat)) { + stop("`mega_mat` must be provided for general IPMs.") + } - if (length(kernels) != n_expected) { - warning( - "`kernels` has length ", length(kernels), - " but IPM contains ", n_expected, " subkernels. ", - "Ensure kernels are provided in correct row-major order." - ) + if (!is.character(mega_mat) || length(mega_mat) == 0) { + stop("`mega_mat` must be a non-empty character vector or string.") + } } ## ---- delta @@ -388,12 +368,12 @@ validate_ipm_sensitivity <- function(ipm, pars, kernels, delta, bounds) { ## return cleaned objects for main function list( - ipm = ipm, - pars = pars, - pars_all = pars_all, - kernels = kernels, - delta = delta, - bounds = bounds + ipm = ipm, + pars = pars, + pars_all = pars_all, + mega_mat = mega_mat, + delta = delta, + bounds = bounds ) } @@ -427,7 +407,7 @@ new_ipmr_uncertainty <- function(ipm, lambdas, mod_uncert, params_uncert, vr_unc #' #' #' @noRd -validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, +validate_ipm_uncertainty <- function(ipm, pars, samples, mega_mat, vr_table, delta, bounds, cores) { base <- validate_ipm_base(ipm) @@ -469,46 +449,27 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, paste(extra_in_samples, collapse = ", ")) } - ## ---- kernels - if (missing(kernels) || is.null(kernels)) { - stop("`kernels` must be provided.") - } - - if (!is.character(kernels) || length(kernels) == 0) { - stop("`kernels` must be a non-empty character vector.") - } + ## ---- mega_mat + is_simple <- grepl("simple", class(ipm)[1]) - # Check for duplicates - dup_kerns <- unique(kernels[duplicated(kernels)]) - - if (length(dup_kerns) > 0) { - warning( - "Duplicate kernel name(s) detected in `kernels`: ", - paste(dup_kerns, collapse = ", "), - ". This may lead to incorrect kernel construction if not intentional." - ) - } + if (is_simple) { - available_kernels <- names(ipm$sub_kernels) + # mega_mat may be NULL for simple IPMs + if (!is.null(mega_mat) && + (!is.character(mega_mat) || length(mega_mat) == 0)) { + stop("`mega_mat` must be a non-empty character vector or string, or NULL.") + } - if (!all(kernels %in% available_kernels)) { - bad <- setdiff(kernels, available_kernels) - stop( - "Unknown kernel name(s): ", - paste(bad, collapse = ", "), - ". Available kernels are: ", - paste(available_kernels, collapse = ", ") - ) - } + } else { - n_expected <- length(ipm$sub_kernels) + # General IPMs require an explicitly specified mega_mat + if (is.null(mega_mat)) { + stop("`mega_mat` must be provided for general IPMs.") + } - if (length(kernels) != n_expected) { - warning( - "`kernels` has length ", length(kernels), - " but IPM contains ", n_expected, " subkernels. ", - "Ensure kernels are provided in correct row-major order." - ) + if (!is.character(mega_mat) || length(mega_mat) == 0) { + stop("`mega_mat` must be a non-empty character vector or string.") + } } ## ---- vr_table @@ -717,7 +678,7 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, kernels, vr_table, ipm = ipm, pars = pars, samples = samples, - kernels = kernels, + mega_mat = mega_mat, vr_table = vr_table, delta = delta, bounds = bounds, diff --git a/man/plot.ipmr_uncertainty.Rd b/man/plot.ipmr_uncertainty.Rd index 8d08946..f2f59d0 100644 --- a/man/plot.ipmr_uncertainty.Rd +++ b/man/plot.ipmr_uncertainty.Rd @@ -7,14 +7,14 @@ \method{plot}{ipmr_uncertainty}(x, type = c("all", "param", "vr", "stacked"), vr_colors = NULL, ...) } \arguments{ -\item{x}{An object of class `ipmr_uncertainty`} +\item{x}{An object of class \code{ipmr_uncertainty}} \item{type}{Type of plot: \itemize{ - \item "param" - parameter-level 3-panel plot - \item "vr" - vital-rate aggregated contributions - \item "stacked" - stacked contributions vs lambda variance - \item "all" - returns all plots as a list +\item "param" - parameter-level 3-panel plot +\item "vr" - vital-rate aggregated contributions +\item "stacked" - stacked contributions vs lambda variance +\item "all" - returns all plots as a list }} \item{vr_colors}{A named character vector of colors used to represent diff --git a/man/print.ipmr_uncertainty.Rd b/man/print.ipmr_uncertainty.Rd index 6687bc7..0b0ef02 100644 --- a/man/print.ipmr_uncertainty.Rd +++ b/man/print.ipmr_uncertainty.Rd @@ -7,7 +7,7 @@ \method{print}{ipmr_uncertainty}(x, ...) } \arguments{ -\item{x}{An object of class `ipmr_uncertainty`} +\item{x}{An object of class \code{ipmr_uncertainty}} \item{...}{Unused} } diff --git a/man/sensitivity.Rd b/man/sensitivity.Rd index 4ba1953..b15f362 100644 --- a/man/sensitivity.Rd +++ b/man/sensitivity.Rd @@ -4,7 +4,7 @@ \alias{sensitivity} \title{Sensitivity and elasticity analysis of population growth rate} \usage{ -sensitivity(ipm, pars = NULL, kernels, delta = 1e-04, bounds = NULL) +sensitivity(ipm, pars = NULL, mega_mat = NULL, delta = 1e-04, bounds = NULL) } \arguments{ \item{ipm}{An \code{ipmr} IPM object.} @@ -12,35 +12,43 @@ sensitivity(ipm, pars = NULL, kernels, delta = 1e-04, bounds = NULL) \item{pars}{A character vector of parameter names to perturb. Defaults to all parameters in \code{ipm}.} -\item{kernels}{A character vector specifying the kernel structure in row -major order.} +\item{mega_mat}{A character vector or string specifying the structure of the +iteration matrix passed to \code{\link{make_iter_kernel}}. For simple IPMs, +\code{mega_mat} may be omitted or set to \code{NULL}. For general IPMs, +\code{mega_mat} must be supplied and specify the arrangement of matrix +blocks in row-major order. Individual blocks may contain sums of subkernels +when multiple kernels contribute to the same matrix block. For example, +\code{"c(P + F, seedbank,...)"} specifies an iteration matrix in which \code{P} +and \code{F} are summed in the same block. See \code{\link{make_iter_kernel}} +for additional details on specifying \code{mega_mat}.} \item{delta}{A numeric scalar specifying the perturbation size used for -numerical differentiation. Default is 1e-4. -@param bounds An optional named list specifying lower and/or upper bounds - for model parameters when calculating numerical sensitivities. Each element - should be a numeric vector of length two giving the lower and upper bounds, - respectively (e.g., \code{list(surv = c(0, 1), recruit = c(0, Inf))}). - Parameters omitted from this list are assumed to be unbounded and are - perturbed using the central finite difference. For parameters with - specified bounds, the function automatically uses a forward or backward - finite difference whenever perturbation by \code{delta} would cross a - boundary; otherwise, the central finite difference is used.} +numerical differentiation. Default is 1e-4.} + +\item{bounds}{An optional named list specifying lower and/or upper bounds +for model parameters when calculating numerical sensitivities. Each element +should be a numeric vector of length two giving the lower and upper bounds, +respectively (e.g., \code{list(surv = c(0, 1), recruit = c(0, Inf))}). +Parameters omitted from this list are assumed to be unbounded and are +perturbed using the central finite difference. For parameters with +specified bounds, the function automatically uses a forward or backward +finite difference whenever perturbation by \code{delta} would cross a +boundary; otherwise, the central finite difference is used.} } \value{ A named list with three elements: \itemize{ - \item \code{sensitivity}: A named numeric vector giving the numerical - sensitivities of the dominant eigenvalue (\eqn{\lambda}) to each - parameter. - \item \code{elasticity}: A named numeric vector giving the elasticities - of the dominant eigenvalue (\eqn{\lambda}) to each parameter. - \item \code{difference_method}: A named character vector indicating the - finite-difference scheme used to calculate the sensitivity for each - parameter. This component is provided to facilitate reproducibility and - to identify parameters whose sensitivities were evaluated using one-sided - finite differences because of user-specified parameter bounds. Values are - \code{"central"}, \code{"forward"}, or \code{"backward"}. +\item \code{sensitivity}: A named numeric vector giving the numerical +sensitivities of the dominant eigenvalue (\eqn{\lambda}) to each +parameter. +\item \code{elasticity}: A named numeric vector giving the elasticities +of the dominant eigenvalue (\eqn{\lambda}) to each parameter. +\item \code{difference_method}: A named character vector indicating the +finite-difference scheme used to calculate the sensitivity for each +parameter. This component is provided to facilitate reproducibility and +to identify parameters whose sensitivities were evaluated using one-sided +finite differences because of user-specified parameter bounds. Values are +\code{"central"}, \code{"forward"}, or \code{"backward"}. } } \description{ @@ -56,7 +64,10 @@ Currently implemented for deterministic ("det"), density-independent ("di") IPMs constructed with \code{ipmr}. Stochastic ("stoch") and density-dependent ("dd") IPMs are not supported. -## Finite differences and bounded parameters +For general IPMs, \code{mega_mat} should describe the complete block +structure of the iteration matrix. The syntax accepted by \code{make_iter_kernel} +is passed through without modification. +\subsection{Finite differences and bounded parameters}{ Sensitivities are estimated numerically by perturbing each parameter by \code{delta}. By default, a central finite difference is used for all @@ -78,6 +89,7 @@ Supplying \code{bounds} is recommended whenever model parameters are subject to biological constraints and numerical perturbation could otherwise produce invalid parameter values or biologically impossible kernels. } +} \examples{ # Probability bounded between 0 and 1 # bounds <- list( diff --git a/man/uncertainty.Rd b/man/uncertainty.Rd index 2f5e9d3..845c02a 100644 --- a/man/uncertainty.Rd +++ b/man/uncertainty.Rd @@ -8,7 +8,7 @@ uncertainty( ipm, pars = NULL, samples, - kernels, + mega_mat = NULL, vr_table, delta = 1e-04, bounds = NULL, @@ -21,14 +21,21 @@ uncertainty( \item{pars}{A character vector of parameter names on which to perform the uncertainty analysis. Defaults to all parameters in \code{ipm}.} -\item{samples}{A data.frame containing sampled parameter values.} +\item{samples}{A dataframe containing sampled parameter values.} -\item{kernels}{A character vector specifying the kernel structure in row -major order.} +\item{mega_mat}{A character vector or string specifying the structure of the +iteration matrix passed to \code{\link{make_iter_kernel}}. For simple IPMs, +\code{mega_mat} may be omitted or set to \code{NULL}. For general IPMs, +\code{mega_mat} must be supplied and specify the arrangement of matrix +blocks in row-major order. Individual blocks may contain sums of subkernels +when multiple kernels contribute to the same matrix block. For example, +\code{"c(P + F, seedbank,...)"} specifies an iteration matrix in which \code{P} +and \code{F} are summed in the same block. See \code{\link{make_iter_kernel}} +for additional details on specifying \code{mega_mat}.} -\item{vr_table}{A data.frame with columns 'parameter' and 'vital_rate' -mapping each parameter to a user-defined vital rate category. Every -parameter in \code{pars} must be present exactly once in this table.} +\item{vr_table}{An optional dataframe with columns 'parameter' and +'vital_rate' mapping each parameter to a user-defined vital rate category. +Every parameter in \code{pars} must be present exactly once in this table.} \item{delta}{A numeric scalar specifying the perturbation size used for numerical differentiation when calculating sensitivities. Default is 1e-4.} @@ -49,16 +56,16 @@ computations will be distributed. Defaults to 1.} \value{ An object of class \code{"ipmr_uncertainty"} containing: \itemize{ - \item \code{ipm}: The complete IPM object as specified in \code{ipm}. - \item \code{lambdas}: Population growth rates (\eqn{\lambda}) calculated - for each sampled parameter set. - \item \code{mod_uncert}: The variance in \eqn{\lambda} across all sampled - parameter sets. - \item \code{params_uncert}: A data frame summarizing parameter-level - sensitivities, elasticities, finite-difference methods, sampling - variance, coefficients of variation, and contributions to uncertainty. - \item \code{vr_uncert}: A data frame summarizing variance contributions - aggregated by user-defined vital rate, if given. +\item \code{ipm}: The complete IPM object as specified in \code{ipm}. +\item \code{lambdas}: Population growth rates (\eqn{\lambda}) calculated +for each sampled parameter set. +\item \code{mod_uncert}: The variance in \eqn{\lambda} across all sampled +parameter sets. +\item \code{params_uncert}: A data frame summarizing parameter-level +sensitivities, elasticities, finite-difference methods, sampling +variance, coefficients of variation, and contributions to uncertainty. +\item \code{vr_uncert}: A data frame summarizing variance contributions +aggregated by user-defined vital rate, if given. } } \description{ @@ -105,6 +112,12 @@ of concurrent jobs. # vr$vital_rate <- c("survival", "growth", "reproduction", ...) # Run uncertainty analysis -# out <- uncertainty(ipm, pars, samples, kernels, vr_table = vr) +# mega_mat <- c("P", "F") +# out <- uncertainty(ipm, pars, samples, mega_mat, vr_table = vr) + +# For multiple kernels contributing to the same matrix block: +# mega_mat <- "c(P + repr, SB1_germ, SB2_germ, +# enter_SB1, SB1_SB1, SB2_SB1, +# enter_SB2, SB1_SB2, SB2_SB2)" } diff --git a/tests/testthat/test-plot-uncertainty.R b/tests/testthat/test-plot-uncertainty.R index 9d3c098..6a28d55 100644 --- a/tests/testthat/test-plot-uncertainty.R +++ b/tests/testthat/test-plot-uncertainty.R @@ -221,7 +221,7 @@ test_that("plot.ipmr_uncertainty runs", { vital_rate = rep("test", length(parameters(ipm))) ) - res <- uncertainty(ipm, samples = samples, kernels = c("P", "F"), vr_table = vr) + res <- uncertainty(ipm, samples = samples, mega_mat = c("P", "F"), vr_table = vr) expect_s3_class(plot(res, type = "param"), "ggplot") expect_s3_class(plot(res, type = "vr"), "ggplot") @@ -244,7 +244,7 @@ test_that("plot type 'all' returns list", { vital_rate = rep("test", length(parameters(ipm))) ) - res <- uncertainty(ipm, samples = samples, kernels = c("P", "F"), vr_table = vr) + res <- uncertainty(ipm, samples = samples, mega_mat = c("P", "F"), vr_table = vr) out <- plot(res, type = "all") @@ -268,7 +268,7 @@ test_that("plot errors on unnamed vr_colors", { vital_rate = rep("test", length(parameters(ipm))) ) - res <- uncertainty(ipm, samples = samples, kernels = c("P", "F"), vr_table = vr) + res <- uncertainty(ipm, samples = samples, mega_mat = c("P", "F"), vr_table = vr) expect_error( plot(res, vr_colors = c("red", "blue")), @@ -292,7 +292,7 @@ test_that("plot errors on missing vr_colors entries", { vital_rate = rep("test", length(parameters(ipm))) ) - res <- uncertainty(ipm, samples = samples, kernels = c("P", "F"), vr_table = vr) + res <- uncertainty(ipm, samples = samples, mega_mat = c("P", "F"), vr_table = vr) expect_error( plot(res, vr_colors = c(foo = "red")), @@ -316,7 +316,7 @@ test_that("print.ipmr_uncertainty runs", { vital_rate = rep("test", length(parameters(ipm))) ) - res <- uncertainty(ipm, samples = samples, kernels = c("P", "F"), vr_table = vr) + res <- uncertainty(ipm, samples = samples, mega_mat = c("P", "F"), vr_table = vr) expect_invisible(print(res)) }) diff --git a/tests/testthat/test-sensitivity.R b/tests/testthat/test-sensitivity.R index b5ae728..a8f4c4a 100644 --- a/tests/testthat/test-sensitivity.R +++ b/tests/testthat/test-sensitivity.R @@ -213,14 +213,14 @@ ipm_stoch <- make_test_ipm("simple_di_stoch") test_that("sensitivity errors on missing ipm", { expect_error( - sensitivity(pars = NULL, kernels = NULL, delta = 1e-4), + sensitivity(pars = NULL, mega_mat = NULL, delta = 1e-4), "`ipm` must be provided and cannot be NULL." ) }) test_that("sensitivity rejects non-list ipm", { expect_error( - sensitivity(ipm = 5, pars = NULL, kernels = NULL), + sensitivity(ipm = 5, pars = NULL, mega_mat = NULL), "must be a valid ipmr object" ) }) @@ -229,7 +229,7 @@ test_that("sensitivity runs on simple_di_det", { ipm <- make_test_ipm("simple_di_det") - res <- sensitivity(ipm, kernels = c("P", "F")) + res <- sensitivity(ipm, mega_mat = c("P", "F")) expect_type(res, "list") expect_true(all(c("sensitivity", "elasticity") %in% names(res))) @@ -245,7 +245,7 @@ test_that("sensitivity respects parameter subset", { res <- sensitivity( ipm, pars = c("s_int", "g_int"), - kernels = c("P", "F") + mega_mat = c("P", "F") ) expect_equal(names(res$sensitivity), c("s_int", "g_int")) @@ -257,7 +257,7 @@ test_that("sensitivity rejects density-dependent IPMs", { ipm <- make_test_ipm("simple_dd_det") expect_error( - sensitivity(ipm, kernels = c("P", "F")), + sensitivity(ipm, mega_mat = c("P", "F")), "Only density-independent IPMs are currently supported." ) }) @@ -267,7 +267,7 @@ test_that("sensitivity rejects stochastic IPMs", { ipm <- make_test_ipm("simple_di_stoch") expect_error( - sensitivity(ipm, kernels = c("P", "F")), + sensitivity(ipm, mega_mat = c("P", "F")), "Only deterministic IPMs are currently supported." ) }) @@ -277,7 +277,7 @@ test_that("sensitivity errors on unknown parameter", { ipm <- make_test_ipm("simple_di_det") expect_error( - sensitivity(ipm, pars = "not_a_param", kernels = c("P", "F")), + sensitivity(ipm, pars = "not_a_param", mega_mat = c("P", "F")), "Unknown parameter" ) }) @@ -289,7 +289,7 @@ test_that("sensitivity rejects non-list bounds", { expect_error( sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = c(0, 1) ), "`bounds` must be NULL or a named list" @@ -304,7 +304,7 @@ test_that("sensitivity rejects unnamed bounds", { expect_error( sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list(c(0, 1)) ), "must be a named list" @@ -319,7 +319,7 @@ test_that("sensitivity rejects unknown parameter names in bounds", { expect_error( sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( fake_parameter = c(0, 1) ) @@ -336,7 +336,7 @@ test_that("sensitivity rejects bounds of incorrect length 1", { expect_error( sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c(0) ) @@ -353,7 +353,7 @@ test_that("sensitivity rejects bounds of incorrect length 3", { expect_error( sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c(0, 1, 2) ) @@ -370,7 +370,7 @@ test_that("sensitivity rejects non-numeric bounds", { expect_error( sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c("a", "b") ) @@ -387,7 +387,7 @@ test_that("sensitivity rejects missing values in bounds", { expect_error( sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c(NA, 1) ) @@ -404,7 +404,7 @@ test_that("sensitivity rejects reversed bounds", { expect_error( sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c(1, 0) ) @@ -421,7 +421,7 @@ test_that("sensitivity rejects parameter values outside supplied bounds", { expect_error( sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c(10, 20) ) @@ -439,7 +439,7 @@ test_that("sensitivity accepts infinite bounds", { sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c(-Inf, Inf), f_s_int = c(0, Inf) @@ -460,7 +460,7 @@ test_that("sensitivity errors when delta exceeds feasible interval", { sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), delta = 1, bounds = list( s_int = c( @@ -482,7 +482,7 @@ test_that("sensitivity returns difference method", { res <- sensitivity( ipm, - kernels = c("P", "F") + mega_mat = c("P", "F") ) expect_true("difference_method" %in% names(res)) @@ -496,14 +496,14 @@ test_that("sensitivity records forward differences for bounded parameters", { res <- sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list(s_int = c(pars$s_int, Inf)) ) expect_warning( sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c(pars$s_int, Inf) ) @@ -526,14 +526,14 @@ test_that("sensitivity records backward differences for bounded parameters", { res <- sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list(s_int = c(-Inf, pars$s_int)) ) expect_warning( sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c(pars$s_int, Inf) ) @@ -555,7 +555,7 @@ test_that("elasticity is consistent with sensitivity definition", { ipm <- make_test_ipm("simple_di_det") - res <- sensitivity(ipm, kernels = c("P", "F")) + res <- sensitivity(ipm, mega_mat = c("P", "F")) pars_all <- parameters(ipm) lambda <- { @@ -572,36 +572,6 @@ test_that("elasticity is consistent with sensitivity definition", { ) }) -test_that("duplicate kernels trigger warning", { - - ipm <- make_test_ipm("simple_di_det") - - expect_warning( - sensitivity(ipm, kernels = c("P", "P")), - "Duplicate kernel name" - ) -}) - -test_that("kernel length mismatch triggers warning", { - - ipm <- make_test_ipm("simple_di_det") - - expect_warning( - sensitivity(ipm, kernels = c("P")), - "has length" - ) -}) - -test_that("unknown kernel errors clearly", { - - ipm <- make_test_ipm("simple_di_det") - - expect_error( - sensitivity(ipm, kernels = c("P", "NOT_REAL")), - "Unknown kernel" - ) -}) - test_that("bounded sensitivity matches unbounded sensitivity when parameter is away from bounds", { ipm <- make_test_ipm("simple_di_det") @@ -610,12 +580,12 @@ test_that("bounded sensitivity matches unbounded sensitivity when parameter is a res1 <- sensitivity( ipm, - kernels = c("P", "F") + mega_mat = c("P", "F") ) res2 <- sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c(0, Inf) ) diff --git a/tests/testthat/test-uncertainty.R b/tests/testthat/test-uncertainty.R index de1143a..5139d27 100644 --- a/tests/testthat/test-uncertainty.R +++ b/tests/testthat/test-uncertainty.R @@ -225,7 +225,7 @@ ker <- c("P", "F") test_that("uncertainty errors on missing IPM", { expect_error( - uncertainty(pars = NULL, samples = NULL, kernels = NULL, vr_table = NULL), + uncertainty(pars = NULL, samples = NULL, mega_mat = NULL, vr_table = NULL), "`ipm` must be provided" ) }) @@ -233,7 +233,7 @@ test_that("uncertainty errors on missing IPM", { test_that("uncertainty rejects non-list IPM", { expect_error( - uncertainty(ipm = 5, samples = data.frame(), kernels = NULL, vr_table = NULL), + uncertainty(ipm = 5, samples = data.frame(), mega_mat = NULL, vr_table = NULL), "valid ipmr object" ) }) @@ -243,7 +243,7 @@ test_that("uncertainty rejects DD IPMs", { ipm <- make_test_ipm("simple_dd_det") expect_error( - uncertainty(ipm, samples = data.frame(), kernels = c("P","F"), vr_table = NULL), + uncertainty(ipm, samples = data.frame(), mega_mat = c("P","F"), vr_table = NULL), "density-independent" ) }) @@ -253,7 +253,7 @@ test_that("uncertainty rejects stochastic IPMs", { ipm <- make_test_ipm("simple_di_stoch") expect_error( - uncertainty(ipm, samples = data.frame(), kernels = c("P","F"), vr_table = NULL), + uncertainty(ipm, samples = data.frame(), mega_mat = c("P","F"), vr_table = NULL), "deterministic" ) }) @@ -270,7 +270,7 @@ test_that("uncertainty returns correct structure", { res <- uncertainty( ipm, samples = samples, - kernels = c("P","F"), + mega_mat = c("P","F"), vr_table = vr ) @@ -293,7 +293,7 @@ test_that("uncertainty respects parameter subset", { ipm, pars = c("s_int", "g_int"), samples = samples, - kernels = c("P","F"), + mega_mat = c("P","F"), vr_table = vr ) @@ -312,7 +312,7 @@ test_that("uncertainty rejects non-list bounds", { uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), vr_table = vr, bounds = c(0, 1) ), @@ -332,7 +332,7 @@ test_that("uncertainty rejects unnamed bounds", { uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list(c(0, 1)), vr_table = vr ), @@ -352,7 +352,7 @@ test_that("uncertainty rejects unknown parameter names in bounds", { uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( fake_parameter = c(0, 1) ), @@ -374,7 +374,7 @@ test_that("uncertainty rejects bounds of incorrect length 1", { uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c(0) ), @@ -396,7 +396,7 @@ test_that("uncertainty rejects bounds of incorrect length 3", { uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c(0, 1, 2) ), @@ -418,7 +418,7 @@ test_that("uncertainty rejects non-numeric bounds", { uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c("a", "b") ), @@ -440,7 +440,7 @@ test_that("uncertainty rejects missing values in bounds", { uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c(NA, 1) ), @@ -462,7 +462,7 @@ test_that("uncertainty rejects reversed bounds", { uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c(1, 0) ), @@ -484,7 +484,7 @@ test_that("uncertainty rejects parameter values outside supplied bounds", { uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c(10, 20) ), @@ -507,7 +507,7 @@ test_that("uncertainty accepts infinite bounds", { uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = list( s_int = c(-Inf, Inf), f_s_int = c(0, Inf) @@ -533,7 +533,7 @@ test_that("uncertainty warns when delta exceeds feasible region", { uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), delta = 5, bounds = list( s_int = c(pars$s_int - 1, @@ -560,7 +560,7 @@ test_that("uncertainty returns difference method", { res <- uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), vr_table = vr ) @@ -585,7 +585,7 @@ test_that("uncertainty records forward differences for bounded parameters", { res <- uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), vr_table = vr, bounds = list(s_int = c(pars$s_int, Inf)) ) @@ -617,7 +617,7 @@ test_that("uncertainty records backward differences for bounded parameters", { res <- uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), vr_table = vr, bounds = list(s_int = c(-Inf, pars$s_int)) ) @@ -648,7 +648,7 @@ test_that("uncertainty records mixed difference methods", { res <- uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), vr_table = vr, bounds = list( s_int = c(pars$s_int, Inf) @@ -682,7 +682,7 @@ test_that("uncertainty errors when delta exceeds feasible interval", { uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), vr_table = vr, delta = 1, bounds = list( @@ -723,7 +723,7 @@ test_that("uncertainty errors when sampled values cross parameter bounds", { uncertainty( ipm, samples = samples, - kernels = c("P","F"), + mega_mat = c("P","F"), vr_table = vr, delta = 0.1, bounds = list( @@ -751,7 +751,7 @@ test_that("uncertainty rejects sampled values outside bounds", { uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), vr_table = vr, bounds = list( s_int = c(pars$s_int, Inf) @@ -802,7 +802,7 @@ test_that("uncertainty computes lambda variance", { res <- uncertainty( ipm, samples = samples, - kernels = c("P","F"), + mega_mat = c("P","F"), vr_table = vr ) @@ -822,11 +822,11 @@ test_that("uncertainty sensitivities match sensitivity()", { unc <- uncertainty( ipm, samples = samples, - kernels = c("P","F"), + mega_mat = c("P","F"), vr_table = vr ) - sens <- sensitivity(ipm, kernels = c("P","F"))$sensitivity + sens <- sensitivity(ipm, mega_mat = c("P","F"))$sensitivity expect_equal( unname(unc$params_uncert$sensitivity), @@ -853,14 +853,14 @@ test_that("bounded uncertainty sensitivities match bounded sensitivity()", { unc <- uncertainty( ipm, samples = samples, - kernels = c("P", "F"), + mega_mat = c("P", "F"), vr_table = vr, bounds = bounds ) sens <- sensitivity( ipm, - kernels = c("P", "F"), + mega_mat = c("P", "F"), bounds = bounds ) @@ -884,34 +884,9 @@ test_that("vr_table join works correctly", { res <- uncertainty( ipm, samples = samples, - kernels = c("P","F"), + mega_mat = c("P","F"), vr_table = vr ) expect_true("vital_rate" %in% names(res$params_uncert)) }) - -test_that("uncertainty passes kernel warnings", { - - ipm <- make_test_ipm("simple_di_det") - - samples <- as.data.frame(matrix( - rep(unlist(parameters(ipm)), each = 3), - nrow = 3, - byrow = TRUE - )) - colnames(samples) <- names(parameters(ipm)) - - vr <- data.frame( - parameter = names(parameters(ipm)), - vital_rate = "test" - ) - - expect_warning( - uncertainty(ipm, - samples = samples, - kernels = c("P", "P"), - vr_table = vr), - "Duplicate kernel name" - ) -}) diff --git a/vignettes/uncertainty-analysis.Rmd b/vignettes/uncertainty-analysis.Rmd index 90ef9f7..1e5f411 100644 --- a/vignettes/uncertainty-analysis.Rmd +++ b/vignettes/uncertainty-analysis.Rmd @@ -12,17 +12,156 @@ vignette: > ## Overview -Integral projection models (IPMs) are commonly used to estimate population growth rates and quantify the demographic importance of different vital rates. However, parameter estimates used to construct IPMs are themselves uncertain, and this uncertainty propagates into model predictions. The `uncertainty()` function estimates how uncertainty in individual parameters contributes to uncertainty in the population growth rate (λ) using finite-difference sensitivities and the covariance structure of sampled parameter estimates. +Integral projection models (IPMs) are commonly used to estimate population growth rates and quantify the demographic importance of different vital rates. However, parameter estimates used to construct IPMs are themselves uncertain, and this uncertainty propagates into model predictions. The `sensitivity()` function estimates how population growth rate (λ) responds to changes in individual model parameters, while `uncertainty()` uses these sensitivities together with the covariance structure of sampled parameter estimates to estimate how parameter uncertainty contributes to uncertainty in λ. This vignette demonstrates how to: +- perform a sensitivity analysis - prepare sampled parameter values - perform an uncertainty analysis - interpret the output -- visualize parameter and vital rate contributions to uncertainty +- visualize parameter and vital rate contributions to uncertainty + +## Sensitivity analysis {#sensitivity} + +The `sensitivity()` function calculates the sensitivity and elasticity of λ to changes in individual model parameters using finite differences. Sensitivity describes the absolute change in λ associated with a change in a parameter, whereas elasticity expresses this relationship on a proportional scale. + +### Performing the analysis with `sensitivity()` + +A sensitivity analysis requires an `ipmr` model object and a specification of the IPM matrix through the `mega_mat` argument. Parameters can optionally be specified using the pars argument. + +```{r, eval = F, echo = T} +sens <- sensitivity( + ipm = ipm, + pars = pars, + mega_mat = mega_mat +) +``` + +If `pars = NULL`, sensitivities are calculated for all parameters in the IPM. + +#### The `mega_mat` argument {#mega-mat} + +The `mega_mat` argument specifies how the subkernels of the IPM should be assembled into the iteration matrix used to calculate λ. It is passed directly to `make_iter_kernel()`, which constructs the complete iteration matrix from the subkernels stores in the `ipmr` object. + +For **simple IPMs**, the iteration kernel is constructed automatically by summing the model's subkernels. In this case, `mega_mat` does not need to specified explicitly. + +For **general IPMS**, `mega_mat` specifies the complete block structure of the iteration matrix. Matrix blocks must be supplied in **row-major order**, meaning that entries are specified from left to right across the first row, followed by the second row, and so on. + +For example, a two-state-class IPM with the block structure + +[ +\mathbf{K} = +\begin{pmatrix} +\mathbf{P}{11} & \mathbf{F}{12} \ +\mathbf{F}{21} & \mathbf{P}{22} +\end{pmatrix} +] + +would be specified as: + +```{r, eval = F, echo = T} +mega_mat <- c( + "P_11", "F_12", + "F_21", "P_22" +) +``` + +The four entries correspond to the four blocks of the matrix in row-major order: + +`P_11 F_12 +F_21 P_22` + +The entries in `mega_mat` must correspond to subkernels that exist in the ipmr object. The names should therefore match the kernel names used when constructing the IPM. + +##### Combining multiple subkernels within a block + +A single matrix block can contain the sum of multiple subkernels. This is useful when more than one demographic process contributes to the same transition between state classes. + +For example, suppose the upper-left block receives contributions from both `P` and `F`, while the lower-right block contains `P_22`. The iteration matrix could be written as + +[ +\mathbf{K} = +\begin{pmatrix} +\mathbf{P} + \mathbf{F} & \mathbf{B} \ +\mathbf{C} & \mathbf{P}_{22} +\end{pmatrix}. +] + +The corresponding mega_mat specification is: + +```{r, eval = F, echo = T} +mega_mat <- c( + "P + F", "B", + "C", "P_22" +) +``` + +Alternatively, the entire specification can be supplied as a character string: + +`mega_mat <- "c(P + F, B, C, P_22)"` + +Both forms use the same syntax internally. The string form can be convenient when the matrix specification is generated programmatically or stored as a single object. + +In most cases, users do not need to construct these matrices manually. The important requirement is that the `mega_mat` specification accurately represents the biological transitions encoded by the IPM and uses the same kernel names as the original model. + +##### Other specifications accepted by `make_iter_kernel()` + +As the `mega_mat` argument is passed only to `make_iter_kernel()`, other specifications are accepted as described in the documentation for \code{\link{make_iter_kernel}}. These include parameter set indices, the identity matrix `I`, and the zero matrix `0`. + +#### Numerical differentiation + +Sensitivities are calculated using finite differences. By default, `sensitivity()` perturbs each parameter by `delta = 1e-4` and uses a central finite difference: + +[ +\frac{\lambda(\theta + \delta) - \lambda(\theta - \delta)} +{2\delta}. +] + +The perturbation size should be small enough to provide a useful numerical approximation while remaining large enough to avoid numerical precision problems. + +For parameters with biological bounds, such as probabilities constrained to the interval [0, 1], subtracting `delta` may produce an invalid parameter value or kernel construction. The bounds argument can be used to specify these constraints: + +```{r, eval = F, echo = T} +sens <- sensitivity( + ipm = ipm, + mega_mat = c("P", "F"), + bounds = list( + s_int = c(-Inf, Inf), #default behavior; does not need to be specified + survival = c(0, 1) + ) +) +``` + +Bounds are specified as a named list, with each element containing the lower and upper bound for one parameter. Parameters not included in `bounds` are treated as unbounded. + +When a parameter is sufficiently close to a specified boundary that a central perturbation would cross that boundary, `sensitivity()` automatically uses a one-sided finite difference. For example, if a parameter is exactly at its lower bound, a forward difference is used: + +[ +\frac{\lambda(\theta + \delta) - \lambda(\theta)} +{\delta}. +] + +Similarly, a backward difference is used when perturbation toward the upper side would cross an upper bound. + +The finite-difference method used for each parameter is returned in the `difference_method` element of the result. + +#### Interpreting sensitivity output + +The result of `sensitivity()` contains sensitivity and elasticity estimates for each parameter: + +`sens$sensitivity +sens$elasticity +sens$difference_method` + +Sensitivity values describe the absolute effect of a parameter on λ, whereas elasticity values describe the proportional effect. Elasticity is therefore particularly useful when comparing parameters measured on different scales. + +The `difference_method` output indicates whether each parameter was evaluated using a central, forward, or backward finite difference. + +Sensitivity analysis can be useful as a standalone analysis for identifying parameters to which population growth is most responsive. It also provides the parameter sensitivities used by `uncertainty()` to decompose uncertainty in λ. ## Quick start guide {#quick} -To perform an uncertainty analysis with `uncertainty()`, you will need four required inputs and three optional inputs. +To perform an uncertainty analysis with `uncertainty()`, you will need four required inputs and four optional inputs. ### `ipm` @@ -42,17 +181,41 @@ Sampled parameter values may be obtained by bootstrap resampling demographic dat All parameters required to construct the IPM must be included in `samples`, even if some parameters are held constant across all sampled models. We recommend using at least 1,000 sampled parameter sets to obtain stable estimates of parameter uncertainty. -### `kernels` +### `mega_mat` + +A character vector or string specifying the subkernel structure of the `ipmr` object. For general IPMs, the `mega_mat` argument describes how subkernels are assembled into the complete iteration matrix. The subkernel names must match those used to construct the `ipmr` object and must be supplied in row-major order. -A character vector specifying the subkernel structure of the `ipmr` object. The kernel names must match those used to construct the `ipmr` object and must be supplied in row-major order. +The same `mega_mat` specification used for `sensitivity()` should be supplied to `uncertainty()` when both analyses are performed on the same IPM. See the [discussion of `mega_mat`]{#mega-mat} in the sensitivity analysis section for examples. ### `vr_table` (optional) -An optional dataframe mapping each parameter to the vital rate to which it corresponds. This table is used to summarize parameter uncertainty contributions by vital rate and to support downstream visualization with the `plot()` function. If omitted, uncertainty contributions are only reported at the parameter level. Calling `uncertainty(ipm, vr_table = "template")` generates a template dataframe that can be manually filled in for all parameters in the IPM. +An optional dataframe mapping each parameter to the vital rate to which it corresponds. This table is used to summarize parameter uncertainty contributions by vital rate and to support downstream visualization with the `plot()` function. If omitted, uncertainty contributions are only reported at the parameter level. + +Calling `uncertainty(ipm, vr_table = "template")` generates a template dataframe that can be manually filled in for all parameters in the IPM. ### `delta` (optional) -A numeric scalar specifying the the perturbation size used for numerical differentiation when calculating parameter sensitivities. The default is `delta = 1e-4`. For some model constructions, selecting a delta value greater than the absolute value of the smallest parameter value may result in invalid subkernels. Care should be taken to select a sufficiently small value for delta such that adding or subtracting this value from all sampled parameter values does not change the sign on the parameter. +A numeric scalar specifying the perturbation size used for numerical differentiation when calculating parameter sensitivities. The default is `delta = 1e-4`. + +For some model constructions, selecting a `delta` value greater than the absolute value of the smallest parameter value may result in invalid subkernels. Care should be taken to select a sufficiently small value for `delta` such that perturbation of the sampled parameter values remains within the biologically feasible parameter space. + +When parameter bounds are known, they can be supplied through the `bounds` argument. `uncertainty()` then uses the same finite-difference rules as `sensitivity()`, including forward or backward differences when a perturbation would cross a specified boundary. + +### `bounds` (optional) + +An optional named list specifying lower and upper bounds for model parameters. Each element should be a numeric vector of length two giving the lower and upper bounds, respectively. + +For example: +```{r, eval = F, echo = T} +bounds <- list( + survival = c(0, 1), + recruitment = c(0, Inf) +) +``` + +Parameters omitted from `bounds` are treated as unbounded. When sampled parameter values occur within `delta` of a boundary, `uncertainty()` automatically selects an appropriate one-sided finite difference rather than perturbing the parameter outside its feasible range. + +The supplied bounds must apply to every sampled value of the corresponding parameter. Sampled values outside the specified bounds are rejected because the resulting sensitivity calculations would no longer represent the specified parameter space. ### `cores` (optional) @@ -73,7 +236,7 @@ unc <- uncertainty( ipm = ipm, samples = samples, pars = pars, - kernels = kernels, + mega_mat = mega_mat, vr_table = vr_table, cores = parallelly::availableCores(omit = 1) ) @@ -95,7 +258,7 @@ A numeric vector containing the population growth rates (λ) calculated from eac ### `mod_uncert` -A single numeric value giving the estimated model uncertainty, calculated as the variance of the sampled λ values (`var(lambdas))`). +A single numeric value giving the estimated model uncertainty, calculated as the variance of the sampled λ values (`var(lambdas)`). ### `params_uncert` @@ -142,7 +305,7 @@ Allowing the model structure to change among bootstrap replicates introduces mod #### Consider the sampling design -Bootstrap resampling should reflect the biological and statistical structure of the original study whenever possible/ Depending on the study design, it may be appropriate to stratify sampling by factors such as year, site, or experimental block so that each bootstrap replicate preserves important features of the data. However, it is important to consider the tradeoff between stratification structure and the resulting sample size of each resampled stratified unit. +Bootstrap resampling should reflect the biological and statistical structure of the original study whenever possible. Depending on the study design, it may be appropriate to stratify sampling by factors such as year, site, or experimental block so that each bootstrap replicate preserves important features of the data. However, it is important to consider the tradeoff between stratification structure and the resulting sample size of each resampled stratified unit. The appropriate resampling strategy depends on the size and structure of the dataset, and no single approach is appropriate for all applications. Users should select a strategy that is consistent with the assumptions of their demographic models and the biological processes being studied. From 6da7eb4eb775abc5b37378c1283e5e4169da74ca Mon Sep 17 00:00:00 2001 From: anon <285962824+anon-lupineuncertainty@users.noreply.github.com> Date: Mon, 24 Aug 2026 15:57:26 +0200 Subject: [PATCH 36/36] Change values plotted in plot type "param" for uncertainty objects --- DESCRIPTION | 2 +- NEWS.md | 12 ++++++++++++ R/utils-uncertainty.R | 13 +++++++------ vignettes/uncertainty-analysis.Rmd | 4 ++-- 4 files changed, 22 insertions(+), 9 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 117fcd4..b7aeb74 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ipmr Title: Integral Projection Models -Version: 0.0.7 +Version: 0.0.8 Authors@R: c( person(given = "Sam", family = "Levin", diff --git a/NEWS.md b/NEWS.md index 562cd39..2773e75 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,15 @@ +# ipmr 0.0.8 + +This release adds methods for quantifying how parameter uncertainty propagates through IPMs. + +## Features + +- Adds 'uncertainty()' to quantify the contribution of parameter uncertainty to variance in population growth rate ('lambda'). + +- Adds 'sensitivity()' to calculate model sensitivity and elasticity to changes in parameter values. + +- Adds a new vignette demonstrating uncertainty analysis in 'ipmr'. + # ipmr 0.0.7 Fixes a test that was failing for `r-devel` on CRAN. There are no other changes. diff --git a/R/utils-uncertainty.R b/R/utils-uncertainty.R index bc0eaff..5bca904 100644 --- a/R/utils-uncertainty.R +++ b/R/utils-uncertainty.R @@ -711,25 +711,26 @@ validate_ipm_uncertainty <- function(ipm, pars, samples, mega_mat, vr_table, title_wrap <- function(x) stringr::str_wrap(x, width = 28) - p1 <- ggplot2::ggplot(df, ggplot2::aes(x = .data$cv, y = .data$parameter, + p1 <- ggplot2::ggplot(df, ggplot2::aes(x = .data$variance, y = .data$parameter, color = .data$vital_rate)) + ggplot2::geom_vline(xintercept = 0) + ggplot2::geom_point(size = 2.5) + - ggplot2::geom_segment(ggplot2::aes(x = 0, xend = .data$cv, + ggplot2::geom_segment(ggplot2::aes(x = 0, xend = .data$variance, yend = .data$parameter), linewidth = 1) + - ggplot2::labs(x = "Coefficient of variation", + ggplot2::scale_x_sqrt() + + ggplot2::labs(x = "Variance", y = "Parameter", title = title_wrap("Parameter uncertainty")) + ggplot2::scale_color_manual(values = vr_colors, guide = "none") + base_theme - p2 <- ggplot2::ggplot(df, ggplot2::aes(x = .data$elasticity, y = .data$parameter, + p2 <- ggplot2::ggplot(df, ggplot2::aes(x = .data$sensitivity, y = .data$parameter, color = .data$vital_rate)) + ggplot2::geom_vline(xintercept = 0) + ggplot2::geom_point(size = 2.5) + - ggplot2::geom_segment(ggplot2::aes(x = 0, xend = .data$elasticity, + ggplot2::geom_segment(ggplot2::aes(x = 0, xend = .data$sensitivity, yend = .data$parameter), linewidth = 1) + - ggplot2::labs(x = "Elasticity", + ggplot2::labs(x = "Sensitivity", y = NULL, title = title_wrap("Model sensitivity to parameter")) + ggplot2::scale_color_manual(values = vr_colors, guide = "none") + diff --git a/vignettes/uncertainty-analysis.Rmd b/vignettes/uncertainty-analysis.Rmd index 1e5f411..ffb2f36 100644 --- a/vignettes/uncertainty-analysis.Rmd +++ b/vignettes/uncertainty-analysis.Rmd @@ -331,8 +331,8 @@ plot(uncert_out, type = "param") ``` The parameter-level plot summarizes the uncertainty analysis for every parameter specified in the analysis. Three panels are displayed: -- **Parameter uncertainty**, represented as the coefficient of variation (CV) of the sampled parameter values -- **Model sensitivity**, shown as the mean elasticity of λ to each parameter across all sampled models +- **Parameter uncertainty**, represented as the variance of the sampled parameter values +- **Model sensitivity**, shown as the mean sensitivity of λ to each parameter across all sampled models - **Parameter contribution to uncertainty**, showing the estimated proportion of model uncertainty attributable to each parameter Parameters are grouped and colored according to the vital rate to which they correspond, when supplied, making it easy to distinguish whether uncertainty arises from variation within a vital rate or differences among vital rates.