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/DESCRIPTION b/DESCRIPTION
index 0a0d6d1..8966889 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",
@@ -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
diff --git a/NAMESPACE b/NAMESPACE
index 3d6b3a9..f774dec 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -56,6 +56,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)
@@ -73,6 +74,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)
@@ -138,80 +140,97 @@ 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)
-importFrom(Rcpp,
-  evalCpp,
-  sourceCpp
-)
-importFrom(grDevices,
-  grey,
-  rainbow
-)
-importFrom(graphics,
-  abline,
-  axis,
-  contour,
-  image,
-  layout,
-  par
-)
+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)
+importFrom(graphics,axis)
+importFrom(graphics,contour)
+importFrom(graphics,image)
+importFrom(graphics,layout)
+importFrom(graphics,par)
 importFrom(magrittr,"%>%")
 importFrom(methods,hasArg)
-importFrom(purrr,
-  flatten,
-  map,
-  map2,
-  map_chr,
-  map_dbl
-)
-importFrom(rlang,
-  ":=",
-  .data,
-  as_environment,
-  call_args,
-  caller_env,
-  child_env,
-  empty_env,
-  enexpr,
-  enexprs,
-  enquo,
-  enquos,
-  ensym,
-  ensyms,
-  env_names,
-  env_parents,
-  eval_bare,
-  expr,
-  expr_text,
-  exprs,
-  is_empty,
-  is_quosure,
-  list2,
-  new_function,
-  parse_expr,
-  parse_exprs,
-  quo,
-  quo_is_null,
-  quo_name,
-  quo_text,
-  quos,
-  sym,
-  syms
-)
-importFrom(stats,
-  bw.nrd0,
-  dnorm,
-  na.omit,
-  runif,
-  setNames
-)
-importFrom(utils,
-  globalVariables,
-  packageDescription
-)
+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)
+importFrom(rlang,call_args)
+importFrom(rlang,caller_env)
+importFrom(rlang,child_env)
+importFrom(rlang,empty_env)
+importFrom(rlang,enexpr)
+importFrom(rlang,enexprs)
+importFrom(rlang,enquo)
+importFrom(rlang,enquos)
+importFrom(rlang,ensym)
+importFrom(rlang,ensyms)
+importFrom(rlang,env_names)
+importFrom(rlang,env_parents)
+importFrom(rlang,eval_bare)
+importFrom(rlang,expr)
+importFrom(rlang,expr_text)
+importFrom(rlang,exprs)
+importFrom(rlang,is_empty)
+importFrom(rlang,is_quosure)
+importFrom(rlang,list2)
+importFrom(rlang,new_function)
+importFrom(rlang,parse_expr)
+importFrom(rlang,parse_exprs)
+importFrom(rlang,quo)
+importFrom(rlang,quo_is_null)
+importFrom(rlang,quo_name)
+importFrom(rlang,quo_text)
+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)
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/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
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"
+  )
+)
diff --git a/R/sensitivity.R b/R/sensitivity.R
new file mode 100644
index 0000000..69c1060
--- /dev/null
+++ b/R/sensitivity.R
@@ -0,0 +1,205 @@
+#' @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 using
+#' finite differences.
+#'
+#' @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 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.
+#'
+#'@details
+#' 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.
+#'
+#' 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
+#' \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}: 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, mega_mat = NULL, delta = 1e-4,
+                        bounds = NULL) {
+
+  ## ---- Input checks
+
+  val <- validate_ipm_sensitivity(ipm, pars, mega_mat, delta, bounds)
+
+  ipm      <- val$ipm
+  pars     <- val$pars
+  pars_all <- val$pars_all
+  mega_mat <- val$mega_mat
+  delta    <- val$delta
+  bounds   <- val$bounds
+
+  ## ---- Setup
+
+  npar        <- length(pars)
+  spar        <- numeric(npar)
+  diff_method <- character(npar)
+
+  kernel <- make_iter_kernel(ipm, mega_mat = mega_mat)
+
+  if (!is.list(kernel) || is.null(kernel[[1]])) {
+    stop("Kernel construction failed.")
+  }
+
+  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 = mega_mat)
+    dominant_eigen(ker[[1]])
+  }
+
+  ## ---- Loop over each parameter
+
+  for (i in seq_len(npar)) {
+
+    par_now <- pars[i]
+    theta <- pars_all[[par_now]]
+
+    ## ---- Define bounds (default none)
+
+    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
+
+    ## ---- Forward difference
+
+    if (use_forward) {
+      lambda_up <- perturb_lambda(par_now, theta + delta)
+
+      spar[i] <- (lambda_up - lambda_all)/delta
+      diff_method[i] <- "forward"
+
+      ## ---- 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"
+
+      ## ---- Central difference
+    } else {
+      lambda_down <- perturb_lambda(par_now, theta - delta)
+      lambda_up <- perturb_lambda(par_now, theta + delta)
+
+      spar[i] <- (lambda_up - lambda_down)/(2*delta)
+      diff_method[i] <- "central"
+    }
+  }
+
+  ## ---- Output
+
+  sens <- spar
+  elas <- spar * abs(as.numeric(pars_all[pars]) / lambda_all)
+
+  names(sens) <- pars
+  names(elas) <- pars
+  names(diff_method) <- pars
+
+  return(list(
+    sensitivity = sens,
+    elasticity = elas,
+    difference_method = diff_method
+  ))
+}
diff --git a/R/uncertainty.R b/R/uncertainty.R
new file mode 100644
index 0000000..43f7ede
--- /dev/null
+++ b/R/uncertainty.R
@@ -0,0 +1,311 @@
+#' @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 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
+#'  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.
+#'
+#'@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"}.
+#'
+#' 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
+#' 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
+#' 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
+#' # 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
+#' # 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, mega_mat = NULL, vr_table,
+                        delta = 1e-4, bounds = NULL, 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, mega_mat,
+                                  vr_table, delta, bounds, cores)
+
+  ipm      <- val$ipm
+  pars     <- val$pars
+  samples  <- val$samples
+  mega_mat <- val$mega_mat
+  vr_table <- val$vr_table
+  delta    <- val$delta
+  bounds   <- val$bounds
+  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", "mega_mat", "delta", "bounds"),
+    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 = mega_mat)
+    mat_base <- ker_base[[1]]
+
+    lambda_base <- Re(eigen(mat_base, only.values = TRUE)$values[1])
+
+    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 = mega_mat)
+      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]
+      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)),
+      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)
+
+  ## ---- 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"))
+  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),
+    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)) {
+    dat_temp <- samples[, par_idx[i]]
+    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 <- stats::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 %>%
+    dplyr::left_join(
+      sens_lookup %>% dplyr::rename(sens1 = sensitivity),
+      by = c("par1" = "parameter")
+    ) %>%
+    dplyr::left_join(
+      sens_lookup %>% dplyr::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(.data$par1) |>
+    dplyr::summarize(variance_sum = sum(.data$variance_part), .groups = "drop")
+
+  var_cont$variance_prop <- var_cont$variance_sum / var_lambda
+
+  stats_temp <- dplyr::left_join(stats_temp, var_cont,
+                                 by = c("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(.data$vital_rate) %>%
+    dplyr::summarize(
+      variance_sum = sum(.data$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
+  )
+}
diff --git a/R/utils-uncertainty.R b/R/utils-uncertainty.R
new file mode 100644
index 0000000..5bca904
--- /dev/null
+++ b/R/utils-uncertainty.R
@@ -0,0 +1,818 @@
+#' 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
+  )
+}
+
+#' Warn if finite-difference perturbation may cross parameter boundaries
+#'
+#' 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.
+#' @param bounds Named list whose elements are numeric vectors of length two
+#'  defining parameter boundary constraints.
+#'
+#' @noRd
+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)
+    }
+
+  }
+
+  warn_pars <- unique(warn_pars)
+
+  if (length(warn_pars) > 0) {
+
+    warning(
+      "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)
+
+}
+
+
+
+#' Validate inputs to sensitivity function
+#'
+#' Returns usable objects for the sensitivity function.
+#'
+#' @noRd
+validate_ipm_sensitivity <- function(ipm, pars, mega_mat, delta, bounds) {
+
+  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.")
+  }
+
+  ## ---- mega_mat
+  is_simple <- grepl("simple", class(ipm)[1])
+
+  if (is_simple) {
+
+    # 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.")
+    }
+
+  } else {
+
+    # General IPMs require an explicitly specified mega_mat
+    if (is.null(mega_mat)) {
+      stop("`mega_mat` must be provided for general IPMs.")
+    }
+
+    if (!is.character(mega_mat) || length(mega_mat) == 0) {
+      stop("`mega_mat` must be a non-empty character vector or string.")
+    }
+  }
+
+  ## ---- 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.")
+  }
+
+  ## ---- 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)
+  )
+
+  if (any(outside)) {
+    stop(
+      "Current parameter value(s) fall outside the supplied bounds: ",
+      paste(names(bounds)[outside], collapse = ", ")
+    )
+  }
+
+  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) {
+      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 = pars_all[pars],
+    delta = delta,
+    bounds = bounds
+  )
+
+  ## return cleaned objects for main function
+  list(
+    ipm       = ipm,
+    pars      = pars,
+    pars_all  = pars_all,
+    mega_mat  = mega_mat,
+    delta     = delta,
+    bounds    = bounds
+  )
+}
+
+
+#' 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, mega_mat, vr_table,
+                                     delta, bounds, 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 = ", "))
+  }
+
+  ## ---- mega_mat
+  is_simple <- grepl("simple", class(ipm)[1])
+
+  if (is_simple) {
+
+    # 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.")
+    }
+
+  } else {
+
+    # General IPMs require an explicitly specified mega_mat
+    if (is.null(mega_mat)) {
+      stop("`mega_mat` must be provided for general IPMs.")
+    }
+
+    if (!is.character(mega_mat) || length(mega_mat) == 0) {
+      stop("`mega_mat` must be a non-empty character vector or string.")
+    }
+  }
+
+  ## ---- 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.")
+  }
+
+  ## ---- 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)
+  )
+
+  if (any(outside)) {
+    stop(
+      "Current parameter value(s) fall outside the supplied bounds: ",
+      paste(names(bounds)[outside], collapse = ", ")
+    )
+  }
+
+  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) {
+      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."
+    )
+  }
+
+  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,
+    delta = delta,
+    bounds = bounds
+  )
+
+  ## ---- 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,
+    mega_mat = mega_mat,
+    vr_table = vr_table,
+    delta    = delta,
+    bounds   = bounds,
+    cores    = as.integer(cores)
+  )
+}
+
+#' 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
+
+  # Order parameters by contribution
+  df$parameter <- factor(
+    df$parameter,
+    levels = rev(df$parameter[order(df$vital_rate)])
+  )
+
+  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) stringr::str_wrap(x, width = 28)
+
+  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$variance,
+                                       yend = .data$parameter), linewidth = 1) +
+    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$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$sensitivity,
+                                       yend = .data$parameter), linewidth = 1) +
+    ggplot2::labs(x = "Sensitivity",
+                  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 = .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 = .data$variance_prop,
+                                       yend = .data$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 + patchwork::plot_layout(ncol = 3, axes = 'collect_y')
+}
+
+#' 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
+
+  df$vital_rate <- factor(
+    df$vital_rate,
+    levels = df$vital_rate[order(df$variance_sum)]
+  )
+
+  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" ) +
+    ggplot2::scale_fill_manual(values = vr_colors) +
+    ggplot2::theme_minimal() +
+    ggplot2::guides(fill = "none")
+}
+
+#' 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
+
+  total_var <- x$mod_uncert
+
+
+  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(
+      ggplot2::aes( xmin = 0.7,
+                    xmax = 1.3,
+                    ymin = 0,
+                    ymax = total_var,
+                    color = "Total variance in lambda"),
+      inherit.aes = FALSE,
+      fill = NA,
+      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 lambda" = "black") ) +
+    ggplot2::guides( fill = ggplot2::guide_legend(order = 1),
+                     color = ggplot2::guide_legend(order = 2) ) +
+    ggplot2::theme_minimal() +
+    ggplot2::theme( axis.text.x = ggplot2::element_blank(),
+                    axis.ticks.x = ggplot2::element_blank() )
+
+  p
+}
diff --git a/man/plot.ipmr_uncertainty.Rd b/man/plot.ipmr_uncertainty.Rd
new file mode 100644
index 0000000..f2f59d0
--- /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 \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{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..0b0ef02
--- /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 \code{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..b15f362
--- /dev/null
+++ b/man/sensitivity.Rd
@@ -0,0 +1,117 @@
+% 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, mega_mat = NULL, delta = 1e-04, bounds = NULL)
+}
+\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{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.}
+
+\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"}.
+}
+}
+\description{
+Computes numerical sensitivities and elasticities of the
+dominant eigenvalue (\eqn{\lambda}) with respect to model parameters using
+finite differences.
+}
+\details{
+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.
+
+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
+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
new file mode 100644
index 0000000..845c02a
--- /dev/null
+++ b/man/uncertainty.Rd
@@ -0,0 +1,123 @@
+% 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,
+  mega_mat = NULL,
+  vr_table,
+  delta = 1e-04,
+  bounds = NULL,
+  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 dataframe containing sampled parameter values.}
+
+\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}{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.}
+
+\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.
+}
+\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"}.
+
+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
+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
+# 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
new file mode 100644
index 0000000..6a28d55
--- /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, mega_mat = 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, mega_mat = 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, mega_mat = 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, mega_mat = 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, 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
new file mode 100644
index 0000000..a8f4c4a
--- /dev/null
+++ b/tests/testthat/test-sensitivity.R
@@ -0,0 +1,600 @@
+# 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("sensitivity errors on missing ipm", {
+  expect_error(
+    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, mega_mat = 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, mega_mat = 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"),
+    mega_mat = 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, mega_mat = c("P", "F")),
+    "Only density-independent IPMs are currently supported."
+  )
+})
+
+test_that("sensitivity rejects stochastic IPMs", {
+
+  ipm <- make_test_ipm("simple_di_stoch")
+
+  expect_error(
+    sensitivity(ipm, mega_mat = c("P", "F")),
+    "Only deterministic IPMs are currently supported."
+  )
+})
+
+test_that("sensitivity errors on unknown parameter", {
+
+  ipm <- make_test_ipm("simple_di_det")
+
+  expect_error(
+    sensitivity(ipm, pars = "not_a_param", mega_mat = c("P", "F")),
+    "Unknown parameter"
+  )
+})
+
+test_that("sensitivity rejects non-list bounds", {
+
+  ipm <- make_test_ipm("simple_di_det")
+
+  expect_error(
+    sensitivity(
+      ipm,
+      mega_mat = 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,
+      mega_mat = 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,
+      mega_mat = 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,
+      mega_mat = 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,
+      mega_mat = 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,
+      mega_mat = 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,
+      mega_mat = 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,
+      mega_mat = 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,
+      mega_mat = 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,
+      mega_mat = c("P", "F"),
+      bounds = list(
+        s_int = c(-Inf, Inf),
+        f_s_int = c(0, Inf)
+      )
+    )
+
+  )
+
+})
+
+test_that("sensitivity errors when delta exceeds feasible interval", {
+
+  ipm <- make_test_ipm("simple_di_det")
+
+  pars <- parameters(ipm)
+
+  expect_error(
+
+    sensitivity(
+      ipm,
+      mega_mat = c("P", "F"),
+      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("sensitivity returns difference method", {
+
+  ipm <- make_test_ipm("simple_di_det")
+
+  res <- sensitivity(
+    ipm,
+    mega_mat = 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,
+    mega_mat = c("P", "F"),
+    bounds = list(s_int = c(pars$s_int, Inf))
+  )
+
+  expect_warning(
+    sensitivity(
+      ipm,
+      mega_mat = 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,
+    mega_mat = c("P", "F"),
+    bounds = list(s_int = c(-Inf, pars$s_int))
+  )
+
+  expect_warning(
+    sensitivity(
+      ipm,
+      mega_mat = 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")
+
+  res <- sensitivity(ipm, mega_mat = 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("bounded sensitivity matches unbounded sensitivity when parameter is away from bounds", {
+
+  ipm <- make_test_ipm("simple_di_det")
+
+  pars <- parameters(ipm)
+
+  res1 <- sensitivity(
+    ipm,
+    mega_mat = c("P", "F")
+  )
+
+  res2 <- sensitivity(
+    ipm,
+    mega_mat = 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
new file mode 100644
index 0000000..5139d27
--- /dev/null
+++ b/tests/testthat/test-uncertainty.R
@@ -0,0 +1,892 @@
+# 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")
+
+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", {
+
+  expect_error(
+    uncertainty(pars = NULL, samples = NULL, mega_mat = NULL, vr_table = NULL),
+    "`ipm` must be provided"
+  )
+})
+
+test_that("uncertainty rejects non-list IPM", {
+
+  expect_error(
+    uncertainty(ipm = 5, samples = data.frame(), mega_mat = 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(), mega_mat = 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(), mega_mat = 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,
+    mega_mat = 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,
+    mega_mat = c("P","F"),
+    vr_table = vr
+  )
+
+  expect_true(all(res$params_uncert$parameter %in%
+                    c("s_int", "g_int")))
+})
+
+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,
+      mega_mat = c("P", "F"),
+      vr_table = vr,
+      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")
+  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,
+      mega_mat = c("P", "F"),
+      bounds = list(c(0, 1)),
+      vr_table = vr
+    ),
+    "must be a named list"
+  )
+
+})
+
+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,
+      mega_mat = c("P", "F"),
+      bounds = list(
+        fake_parameter = c(0, 1)
+      ),
+      vr_table = vr
+    ),
+    "Unknown parameter"
+  )
+
+})
+
+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,
+      mega_mat = c("P", "F"),
+      bounds = list(
+        s_int = c(0)
+      ),
+      vr_table = vr
+    ),
+    "length 2"
+  )
+
+})
+
+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,
+      mega_mat = c("P", "F"),
+      bounds = list(
+        s_int = c(0, 1, 2)
+      ),
+      vr_table = vr
+    ),
+    "length 2"
+  )
+
+})
+
+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,
+      mega_mat = c("P", "F"),
+      bounds = list(
+        s_int = c("a", "b")
+      ),
+      vr_table = vr
+    ),
+    "must be numeric"
+  )
+
+})
+
+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,
+      mega_mat = c("P", "F"),
+      bounds = list(
+        s_int = c(NA, 1)
+      ),
+      vr_table = vr
+    ),
+    "may not contain NA"
+  )
+
+})
+
+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,
+      mega_mat = c("P", "F"),
+      bounds = list(
+        s_int = c(1, 0)
+      ),
+      vr_table = vr
+    ),
+    "lower bound"
+  )
+
+})
+
+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,
+      mega_mat = c("P", "F"),
+      bounds = list(
+        s_int = c(10, 20)
+      ),
+      vr_table = vr
+    ),
+    "outside the 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"
+  samples <- as.data.frame(lapply(parameters(ipm), function(x) rep(x, 50)))
+
+  expect_no_error(
+
+    uncertainty(
+      ipm,
+      samples = samples,
+      mega_mat = c("P", "F"),
+      bounds = list(
+        s_int = c(-Inf, Inf),
+        f_s_int = c(0, Inf)
+      ),
+      vr_table = vr
+    )
+
+  )
+
+})
+
+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)
+
+  expect_error(
+
+    uncertainty(
+      ipm,
+      samples = samples,
+      mega_mat = c("P", "F"),
+      delta = 5,
+      bounds = list(
+        s_int = c(pars$s_int - 1,
+                  pars$s_int + 1)
+      ),
+      vr_table = vr
+    ),
+
+    "delta"
+
+  )
+
+})
+
+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,
+    mega_mat = 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,
+    mega_mat = 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,
+    mega_mat = 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,
+    mega_mat = 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,
+      mega_mat = 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,
+      mega_mat = 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,
+      mega_mat = 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 ----------------------------------------------
+
+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,
+    mega_mat = 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,
+    mega_mat = c("P","F"),
+    vr_table = vr
+  )
+
+  sens <- sensitivity(ipm, mega_mat = c("P","F"))$sensitivity
+
+  expect_equal(
+    unname(unc$params_uncert$sensitivity),
+    unname(sens[unc$params_uncert$parameter]),
+    tolerance = 1e-6
+  )
+})
+
+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,
+    mega_mat = c("P", "F"),
+    vr_table = vr,
+    bounds = bounds
+  )
+
+  sens <- sensitivity(
+    ipm,
+    mega_mat = 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")
+
+  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,
+    mega_mat = c("P","F"),
+    vr_table = vr
+  )
+
+  expect_true("vital_rate" %in% names(res$params_uncert))
+})
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)
+})
diff --git a/vignettes/uncertainty-analysis.Rmd b/vignettes/uncertainty-analysis.Rmd
new file mode 100644
index 0000000..ffb2f36
--- /dev/null
+++ b/vignettes/uncertainty-analysis.Rmd
@@ -0,0 +1,365 @@
+---
+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 `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 
+
+## 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 four 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.  
+
+### `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.  
+
+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. 
+
+### `delta` (optional)  
+
+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)  
+
+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, eval = F, echo = T}
+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, eval = F, echo = T}
+unc <- uncertainty(
+  ipm = ipm,
+  samples = samples,
+  pars = pars,
+  mega_mat = mega_mat,
+  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, 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
+- 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 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.  
+
+#### 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 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.  
+
+### 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.