## ----setup, message = FALSE, warning = FALSE---------------------------------- knitr::opts_chunk$set(echo = TRUE, tidy = FALSE, collapse = TRUE, comment = "#>") library(raretrans) ## ----define-matrices---------------------------------------------------------- # ── Transition matrix T ──────────────────────────────────────────────────────── # Each column must sum to ≤ 1 (the remainder is mortality). # Rows = destination stage; Columns = source stage. T_mat <- matrix( c(0.2368, 0.0000, 0.0000, # seed → seed, nonrep, rep 0.0010, 0.0100, 0.0200, # nonrep → seed, nonrep, rep 0.0000, 0.1429, 0.1364), # rep → seed, nonrep, rep nrow = 3, ncol = 3, byrow = TRUE, dimnames = list(c("seed", "nonrep", "rep"), c("seed", "nonrep", "rep")) ) # ── Fecundity matrix F ───────────────────────────────────────────────────────── # F[i,j] = average number of stage-i individuals produced by one stage-j parent. F_mat <- matrix( c(0, 0, 1.9638, # reproductive adults produce seeds 0, 0, 8.3850, # reproductive adults produce nonreproductive offspring 0, 0, 0.0000), # no fecundity contributions to the reproductive stage nrow = 3, ncol = 3, byrow = TRUE, dimnames = list(c("seed", "nonrep", "rep"), c("seed", "nonrep", "rep")) ) # ── Wrap in a named list — the format expected by raretrans ──────────────────── TF <- list(T = T_mat, F = F_mat) TF ## ----define-N----------------------------------------------------------------- N <- c(seed = 80, nonrep = 25, rep = 12) N ## ----observed-lambda---------------------------------------------------------- # Dominant eigenvalue = asymptotic growth rate λ A_obs <- TF$T + TF$F lambda_obs <- Re(eigen(A_obs)$values[1]) cat("Observed lambda:", round(lambda_obs, 3), "\n") ## ----col-sums----------------------------------------------------------------- colSums(TF$T) ## ----uniform-prior------------------------------------------------------------ # Uniform prior: equal weight to all fates (including death) P_uniform <- matrix(0.25, nrow = 4, ncol = 3) fill_transitions(TF, N, P = P_uniform) ## ----compare-uniform---------------------------------------------------------- TF$T ## ----informative-prior-------------------------------------------------------- P_info <- matrix( c(0.25, 0.05, 0.00, # → seed 0.05, 0.55, 0.15, # → nonrep 0.00, 0.15, 0.70, # → rep 0.70, 0.25, 0.15), # → dead nrow = 4, ncol = 3, byrow = TRUE ) # Verify columns sum to 1 colSums(P_info) ## ----informative-fill--------------------------------------------------------- T_posterior <- fill_transitions(TF, N, P = P_info) T_posterior ## ----priorweight-------------------------------------------------------------- fill_transitions(TF, N, P = P_info, priorweight = 0.5) ## ----fertility-prior---------------------------------------------------------- # alpha = prior "offspring counts"; beta = prior "adult counts" # Only the (seed, rep) and (nonrep, rep) entries are reproductive alpha <- matrix( c(NA, NA, 1e-5, NA, NA, 1e-5, NA, NA, NA), nrow = 3, ncol = 3, byrow = TRUE ) beta <- matrix( c(NA, NA, 1e-5, NA, NA, 1e-5, NA, NA, NA), nrow = 3, ncol = 3, byrow = TRUE ) F_posterior <- fill_fertility(TF, N, alpha = alpha, beta = beta) F_posterior ## ----combined-posterior------------------------------------------------------- posterior <- list( T = fill_transitions(TF, N, P = P_info), F = fill_fertility(TF, N, alpha = alpha, beta = beta) ) posterior # Asymptotic growth rate of the posterior matrix lambda_post <- Re(eigen(posterior$T + posterior$F)$values[1]) cat("Posterior lambda:", round(lambda_post, 3), "\n") cat("Observed lambda:", round(lambda_obs, 3), "\n") ## ----returnType-TN------------------------------------------------------------ TN <- fill_transitions(TF, N, P = P_info, returnType = "TN") TN ## ----returnType-A------------------------------------------------------------- fill_transitions(TF, N, P = P_info, returnType = "A") ## ----returnType-ab------------------------------------------------------------ fill_fertility(TF, N, alpha = alpha, beta = beta, returnType = "ab") ## ----sim-lambda, fig.width=6, fig.height=4, fig.cap="Posterior distribution of the asymptotic growth rate λ for *Chamaedorea elegans*. The dashed vertical line marks λ = 1 (stable population). Values to the right indicate growth; values to the left indicate decline."---- set.seed(20240301) # for reproducibility # Simulate 1000 projection matrices from the posterior sims <- sim_transitions(TF, N, P = P_info, alpha = alpha, beta = beta, priorweight = 0.5, samples = 1000) # Extract λ from each simulated matrix lambdas <- sapply(sims, function(mat) Re(eigen(mat)$values[1])) # Summarise cat("Posterior median λ:", round(median(lambdas), 3), "\n") cat("95% credible interval: [", round(quantile(lambdas, 0.025), 3), ",", round(quantile(lambdas, 0.975), 3), "]\n") cat("P(λ > 1):", round(mean(lambdas > 1), 3), "\n") # Plot hist(lambdas, breaks = 30, main = expression("Posterior distribution of " * lambda), xlab = expression(lambda), col = "steelblue", border = "white") abline(v = 1, lty = 2, lwd = 2, col = "firebrick") ## ----cri---------------------------------------------------------------------- cri <- transition_CrI(TF, N, P = P_info, priorweight = 0.5, ci = 0.95, stage_names = c("seed", "nonrep", "rep")) cri ## ----plot-cri, fig.width=7, fig.height=4, fig.cap="Posterior mean transition probabilities (points) and 95% credible intervals (lines) for *Chamaedorea elegans*. Each panel shows the fate distribution from one source stage."---- plot_transition_CrI(cri) ## ----plot-cri-no-dead, fig.width=7, fig.height=4, fig.cap="As above but excluding the dead fate."---- plot_transition_CrI(cri, include_dead = FALSE) ## ----plot-density, fig.width=7, fig.height=7, fig.cap="Posterior beta density for each transition in the *C. elegans* matrix. Columns = source stage (from); rows = destination stage (to). Shaded region = 95% credible interval. Zero-probability transitions show a degenerate spike at 0."---- plot_transition_density(TF, N, P = P_info, priorweight = 0.5, stage_names = c("seed", "nonrep", "rep")) ## ----plot-density-no-dead, fig.width=7, fig.height=5, fig.cap="As above but excluding the dead fate row."---- plot_transition_density(TF, N, P = P_info, priorweight = 0.5, stage_names = c("seed", "nonrep", "rep"), include_dead = FALSE)