--- title: "HRRI: Holobiont Redox Resilience Index — End-to-End Workflow" author: - name: Mitra Ghotbi email: mitra.ghotbi@gmail.com date: "`r Sys.Date()`" package: HRRI output: rmarkdown::html_vignette: toc: true toc_depth: 3 number_sections: true css: css/hrri.css fig_width: 7 fig_height: 4.2 dev: png vignette: > %\VignetteIndexEntry{HRRI: Holobiont Redox Resilience Index — End-to-End Workflow} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.2, fig.align = "center", dpi = 150, out.width = "100%", message = FALSE, warning = FALSE ) ## Consistent figure styling across the vignette ------------------------- if (requireNamespace("ggplot2", quietly = TRUE)) { old_theme <- ggplot2::theme_set( ggplot2::theme_minimal(base_size = 11) + ggplot2::theme( panel.grid.minor = ggplot2::element_blank(), panel.grid.major = ggplot2::element_line(linewidth = 0.3, colour = "#dde3e1"), axis.title = ggplot2::element_text(colour = "#4a5451"), axis.text = ggplot2::element_text(colour = "#4a5451"), strip.text = ggplot2::element_text(face = "bold", colour = "#1c2321"), plot.title = ggplot2::element_text(face = "bold", colour = "#1c2321"), legend.position = "bottom" ) ) } ## Vignette-local adapter: retain every required design identifier even when ## a pipeline release returns only scores or a subset of identifier columns. attach_hrri_ids <- function(scores, id) { scores <- as.data.frame(scores) id <- as.data.frame(id) keys <- c("plot", "depth", "plant_id", "time") if (anyDuplicated(names(scores)) || anyDuplicated(names(id))) { stop("Score and identifier tables must have unique column names.") } missing_keys <- setdiff(keys, names(id)) if (length(missing_keys)) { stop("sim$id is missing: ", paste(missing_keys, collapse = ", ")) } if (!nrow(id) || nrow(scores) != nrow(id)) { stop("Expected one pipeline score row per sim$id row.") } if (!"RRI" %in% names(scores) || !is.numeric(scores[["RRI"]])) { stop("Pipeline scores must contain a numeric RRI column.") } make_key <- function(x) { if (anyNA(x[keys])) stop("Observation identifiers cannot be missing.") values <- lapply(x[keys], as.character) if (any(vapply(values, function(v) any(grepl("\034", v, fixed = TRUE)), logical(1)))) { stop("Observation identifiers contain the reserved key separator.") } do.call(paste, c(values, list(sep = "\034"))) } id_key <- make_key(id) if (anyDuplicated(id_key)) { stop("sim$id must uniquely identify plot x depth x plant x time rows.") } ## Put scores in the original simulator order whenever usable keys exist. if (all(keys %in% names(scores))) { score_key <- make_key(scores) if (anyDuplicated(score_key)) stop("Duplicate observation keys in scores.") idx <- match(id_key, score_key) alignment <- "observation keys" } else if ("row_id" %in% names(id) && "row_id" %in% names(scores)) { id_row <- as.character(id[["row_id"]]) score_row <- as.character(scores[["row_id"]]) if (anyNA(id_row) || anyNA(score_row) || anyDuplicated(id_row) || anyDuplicated(score_row)) { stop("row_id must be unique and nonmissing for key-based alignment.") } idx <- match(id_row, score_row) alignment <- "row_id" } else { ## Scores-only releases must preserve their input row order. Row count ## equality alone cannot prove this; do not sort either table beforehand. idx <- seq_len(nrow(id)) alignment <- "input row order (pipeline contract)" } if (anyNA(idx)) stop("Some simulator identifiers have no matching scores.") scores <- scores[idx, , drop = FALSE] shared <- intersect(c(keys, "row_id"), intersect(names(scores), names(id))) for (nm in shared) { if (!identical(as.character(scores[[nm]]), as.character(id[[nm]]))) { stop("Conflicting or misaligned identifier column: ", nm) } } ## Rebuild identifiers once; never create depth.x/depth.y pairs. out <- cbind(id[keys], scores[setdiff(names(scores), keys)]) rownames(out) <- NULL attr(out, "id_alignment") <- alignment out } finite_mean <- function(x) { x <- x[is.finite(x)] if (length(x)) mean(x) else NA_real_ } ``` # Introduction The **HRRI** package implements exploratory, multi-domain diagnostics for describing how soil–plant–microbiome systems buffer and recover from hydroclimatic redox disturbances (Ghotbi *et al.*, 2026). The theoretical framing distinguishes four properties. They are not all identifiable from a single observation curve: | Property | Symbol | Interpretation | |---|---|---| | Capacity | $Q$ | Electron-accepting and electron-donating inventory available within the system (mmol e⁻ kg⁻¹) | | Connectivity | $\alpha$ | Fraction of $Q$ electrochemically accessible to porewater | | Kinetics | $k$ | Characteristic rate of electron exchange under physicochemical and biological constraints (h⁻¹) | | Memory | $M$ | Legacy of prior disturbances retained through persistent biogeochemical, microbial and physiological states that influence future system responses | These combine through the accessible-capacity formula: $$C_{\rm acc} = \sum_i Q_i \cdot \alpha_i \cdot \left(1 - e^{-k_i \tau}\right)$$ and aggregate into the Holobiont Redox Resilience Index: $$\mathrm{RRI}_{it} = w_P \cdot P_{it} + w_S \cdot S_{it} + w_M \cdot M_{it}$$ where $P$ (Physiology), $S$ (Soil), and $M$ (Microbial) are domain scores. This vignette walks through the full workflow from in-silico data generation to RRI computation, accessible-capacity estimation, and recovery-signature metrics — all driven by `simulate_redox_holobiont()`. # In-silico Data Generation `simulate_redox_holobiont()` is the package's master simulator. It generates synthetic longitudinal observations across all three holobiont domains.
Only Fe and Mn inventories carry closed-balance checks. Carbon, nitrogen, sulfur and oxygen budgets are **not** closed, and all rate parameters are illustrative rather than field-calibrated. Simulated output supports software demonstration and falsifiable model checks, not empirical ecological inference.
```{r simulate} library(HRRI) packageVersion("HRRI") ## Compatibility shim ----------------------------------------------------- ## rri_pipeline() is the convenience wrapper around rri_pipeline_st(). ## If the *installed* HRRI predates the wrapper, define an equivalent local ## version so this vignette knits against either release. Reinstall the ## package (see README) to use the exported function directly. if (!exists("rri_pipeline", mode = "function")) { message("Installed HRRI has no rri_pipeline(); using a vignette-local wrapper.") rri_pipeline <- function(dat = NULL, soil = NULL, plant = NULL, micro = NULL, id = NULL, domain_weights = c(Physio = 0.4, Soil = 0.35, Micro = 0.25), ...) { stopifnot(setequal(names(domain_weights), c("Physio", "Soil", "Micro"))) w <- domain_weights[c("Physio", "Soil", "Micro")] w <- w / sum(w) res <- rri_pipeline_st( ROS_flux = plant, Eh_stability = soil, micro_data = micro, id = id, w1 = unname(w[1]), w2 = unname(w[2]), w3 = unname(w[3]), ... ) res$scores <- res$row_scores res } } ## Reproducible 1-cycle flood-drain experiment ## n_plot=2, n_depth=2, n_plant=3, n_time=30 -> 360 rows sim <- simulate_redox_holobiont( n_plot = 2, n_depth = 2, n_plant = 3, n_time = 30, p_micro = 20, seed = 42, scenario = "flood_drain", n_cycles = 1, disturbance_strength = 0.70, history_strength = 0.55, decoupling = 0.20 ) ## Top-level structure names(sim) nrow(sim$id) # one row per plot × depth × plant × time ``` ## Design identifiers ```{r design} head(sim$id[, c("plot","depth","plant_id","time","cycle","phase","WFPS")]) ``` ## Soil geochemical outputs ```{r soil} head(sim$soil_data[, c("EAC","EDC","Cacc_EAC","Cacc_total","Cacc_fraction", "FeIII_poor_crystalline_mmol_kg", "FeII_mmol_kg","Eh","pH")]) ``` ## Fe mass-balance verification ```{r conservation} ## Maximum absolute error should be < 0.01 mmol kg-1 sim$conservation_checks ``` ## Plant physiology ```{r plant} head(sim$plant_data[, c("SPAD","FvFm","ROL","ROS_load","aerenchyma")]) ``` ## Microbial functional genes ```{r genes} ## 18 genes spanning Fe-cycling, denitrification, nitrification, ## methanogenesis, and sulfur cycling colnames(sim$micro_gene_abundance) summary(sim$micro_gene_abundance[, "mcrA"]) # methanogenesis gene ``` # Accessible-Capacity Estimation `rri_accessible_capacity()` computes $C_{\rm acc}$ for arbitrary mineralogical reservoirs with explicitly supplied accessibility and exchange-rate parameters. The values below are illustrative model inputs, not estimates or validated literature defaults. ```{r cacc} ## Subset one plot-depth unit for illustration idx <- sim$id$plot == "P1" & sim$id$depth == "D1" & sim$id$plant_id == "Plant1" sdf <- sim$soil_data[idx, ] ## Define reservoir specifications ## (Q_col names must match columns in sdf) res_spec <- list( reactive_FeIII = list( Q_col = "FeIII_poor_crystalline_mmol_kg", alpha = "alpha_accept", # column name: per-row connectivity k = "k_accept_h", # column name: per-row kinetics type = "EAC" ), crystalline_FeIII = list( Q_col = "FeIII_crystalline_mmol_kg", alpha = 0.20, # attenuated connectivity for crystalline phases k = 0.008, # h-1: slow exchange (goethite/hematite) type = "EAC" ), FeII_pool = list( Q_col = "FeII_mmol_kg", alpha = "alpha_donate", k = "k_donate_h", type = "EDC" ) ) ## tau = 24 h (diurnal event timescale) cap <- rri_accessible_capacity(sdf, res_spec, tau = 24, normalise = FALSE, return_components = TRUE) ## Per-component summary (returned because return_components = TRUE) cap$components ## Mean accessible vs. total inventory cat("Mean Cacc_raw:", mean(cap$cacc_raw, na.rm=TRUE), "mmol e- kg-1\n") cat("Mean fraction :", mean(cap$cacc_fraction, na.rm=TRUE), "\n") if ("ck_limited" %in% names(cap) && length(cap$ck_limited)) { cat("CK-limited rows:", sum(cap$ck_limited, na.rm=TRUE), "/", sum(!is.na(cap$ck_limited)), "classified rows\n") } else { cat("CK-limited classification is not returned by this HRRI version.\n") } ``` ## Effect of event timescale τ ```{r tau_sweep} tau_vals <- c(1, 6, 24, 72, 168, 720) # 1 h to 30 d cacc_tau <- sapply(tau_vals, function(tt) { r <- rri_accessible_capacity(sdf, res_spec, tau = tt, normalise = FALSE) mean(r$cacc_raw, na.rm = TRUE) }) data.frame(tau_h = tau_vals, Cacc_mean = round(cacc_tau, 2)) ``` # RRI Pipeline `rri_pipeline()` integrates available observed domains into an exploratory composite. The simulator's hidden architecture columns are excluded. We use explicit measured anchors to orient otherwise arbitrary PCA axes. ```{r rri_pipeline} rri_out <- rri_pipeline( plant = sim$ROS_flux, soil = sim$Eh_stability, micro = log1p(sim$micro_gene_abundance), id = sim$id, mode = "snapshot", scaling = "pnorm", direction_anchor_phys = "FvFm", direction_anchor_soil = "Eh", direction_anchor_micro = "mtrA", domain_weights = c(Physio=0.35, Soil=0.40, Micro=0.25) ) ## Align once and reuse this identifier-complete table downstream. rri_scored <- attach_hrri_ids(rri_out$row_scores, sim$id) attr(rri_scored, "id_alignment") summary(rri_scored$RRI) head(rri_scored[, c("plot", "depth", "plant_id", "time", "RRI", "Physio", "Soil", "Micro")]) ```
Alignment uses complete observation keys, or a shared unique `row_id` when available. If neither is returned the pipeline must preserve input row order — matching row counts alone do **not** establish alignment. Shared identifier columns are checked for conflicts rather than silently overwritten.
## Agreement with the prescribed target A single pooled correlation is the wrong summary here, for two reasons. First, these 360 rows are **12 trajectories observed at 30 time points**, not 360 independent observations. Rows within a trajectory are strongly dependent, so an interval computed from the row count is far too narrow. Second, Pearson's $r$ measures *association*, not *agreement*. A score equal to twice the target plus a constant correlates with it perfectly while matching it nowhere. Lin's concordance correlation coefficient penalises departure from the 1:1 line and is the quantity that belongs beside it. `rri_accuracy()` reports both, with intervals obtained by resampling whole trajectories rather than rows. ```{r validation} ## rri_scored is aligned to sim$id, and hence to its latent_truth vector. truth <- sim$latent_truth if (!is.numeric(truth) || length(truth) != nrow(rri_scored)) { stop("latent_truth must be a numeric vector with one value per sim$id row.") } ## One independent experimental unit = one plot x depth x plant trajectory. traj <- interaction(rri_scored$plot, rri_scored$depth, rri_scored$plant_id, drop = TRUE) acc <- rri_accuracy( score = rri_scored$RRI, target = truth, cluster = traj, n_boot = 500, n_perm = 500, seed = 42 ) acc ```
`effective_n` in the dependence table, not the row count, is what governs precision. Where the design effect is well above one, the naive interval should not be quoted: the function prints both widths so the difference is visible rather than asserted.
Splitting the error says which kind of disagreement is present, and the three components sum to the mean squared error exactly. ```{r validation_decomp} acc$decomposition[, c("component", "percent")] ## Exactness check: the residual is numerical noise, not a rounding allowance. c(mse = attr(acc$decomposition, "mse"), residual = attr(acc$decomposition, "residual")) ``` Large squared bias is a systematic offset, removable by recentring. Large variance mismatch means the score is flatter or more volatile than the target. Large lack of correlation means the score does not track the target's pattern, and no rescaling will repair it. `plot_rri_accuracy()` draws the same four questions as one figure. ```{r validation_figure, fig.width=9.5, fig.height=7.5, out.width="100%"} plot_rri_accuracy(acc, score_label = "RRI", target_label = "Prescribed target") ```
**Reading it.** **A** puts the fitted line against the dashed 1:1 line; a flatter fit means the score compresses the target's range. Open points are trajectory means, the level at which the units are independent. **B** is a Bland-Altman plot: a scatter that slopes or fans out shows disagreement that depends on level, which no correlation coefficient can reveal. **C** is the headline: the violet distribution resamples rows and is too narrow, the teal one resamples trajectories and is honest; the bars beneath give both widths. **D** partitions the mean squared error exactly. **What it does not show.** None of the four panels speaks to out-of-sample performance. The target is prescribed by the same simulator that produced the inputs, so a tight panel A means the estimator is self-consistent, not that it would recover an unobserved field quantity.
`latent_truth` and `RRI` are produced by the same generator. Everything above therefore quantifies **internal consistency** — whether the estimator recovers the target its own simulator prescribed. It is not independent predictive validation, and it is not evidence that the individual latent parameters have been identified. An empirical claim requires a target measured independently of the score, replicated across independent experimental units.
# Recovery Signatures `rri_recovery_metrics()` summarizes the RRI trajectory around a specified disturbance window. The output fields depend on the installed package version; the complete returned table is displayed below without assuming legacy names. A score minimum or persistent departure alone does not establish biochemical pathway truncation or alternative routing. Those interpretations require independent process evidence. The function requires **one row per time point per group**. Because the pipeline returns one row per plant, we must first average RRI across plants within each plot × depth × time cell before calling recovery metrics. ```{r recovery} ## Step 1 — use the aligned score table created in the pipeline chunk. ## Step 2 — aggregate to one row per plot × depth × time (mean over plants). ## The data-frame method avoids formula-level complete-case filtering. ## An all-missing group remains NA; it is not replaced by zero. rri_agg <- stats::aggregate( x = rri_scored["RRI"], by = rri_scored[c("plot", "depth", "time")], FUN = finite_mean ) rri_agg <- rri_agg[order(rri_agg$plot, rri_agg$depth, rri_agg$time), , drop = FALSE] rownames(rri_agg) <- NULL stopifnot(!anyDuplicated(rri_agg[c("plot", "depth", "time")])) ## Step 3 — extract recovery signatures. ## Pass rri_agg directly (group columns are already inside it; no id= needed). metrics <- rri_recovery_metrics( res = rri_agg, time_col = "time", group_cols = c("plot","depth"), perturb_start = 8, perturb_end = 18, rri_col = "RRI", forcing_col = NULL # no measured forcing supplied ) if (!is.data.frame(metrics) || nrow(metrics) == 0L) { stop("rri_recovery_metrics() returned no nonempty recovery data frame.") } names(metrics) metrics ``` With `forcing_col = NULL`, a returned hysteresis-related field must be interpreted according to the installed function's documented definition; it does not demonstrate a measured forcing–response loop. The limits 8 and 18 are example analysis boundaries, in the units of `time`; verify them against the simulator's event schedule before interpreting recovery rates. # Disturbance-History Sensitivity Repeated-cycle simulations can test consequences of this simulator's stated rules. They do not establish a universal mineralogical ratchet. In this version, crystallisation is continuous rather than restricted to reoxidation events, and end-state EAC is not constrained to decline monotonically with cycle count. ```{r history_sensitivity} history <- do.call(rbind, lapply(1:4, function(nc) { z <- simulate_redox_holobiont(n_plot=1, n_depth=1, n_plant=2, n_time=30, p_micro=5, seed=99, n_cycles=nc, disturbance_strength=0.70) keep <- z$id$plant_id=="Plant1" data.frame(n_cycles=nc, EAC_end=tail(z$soil_data$EAC[keep],1), memory_end=tail(z$latent_state$memory[keep],1)) })) history ``` Interpret the direction and magnitude as a model sensitivity result. An empirical claim about hydrological memory requires independent observations. ```{r restore-theme, include=FALSE} ## theme_set() changes state that persists for the rest of the session. ## Vignettes build in their own process so nothing outside is affected, but ## restoring is the same courtesy CRAN asks for with par() and options(). if (exists("old_theme")) ggplot2::theme_set(old_theme) ``` # Session Information ```{r sessionInfo} sessionInfo() ``` # References ## Published methods Keiluweit, M., Wanzek, T., Kleber, M., Nico, P., & Fendorf, S. (2017). Anaerobic microsites have an unaccounted role in soil carbon stabilization. *Nature Communications*, **8**, 1771. Klüpfel, L., Piepenbrock, A., Kappler, A., & Sander, M. (2014). Humic substances as fully regenerable electron acceptors in recurrently anoxic environments. *Nature Geoscience*, **7**, 195–200. Kobayashi, K., & Salam, M. U. (2000). Comparing simulated and measured values using mean squared deviation and its components. *Agronomy Journal*, **92**, 345–352. Lin, L. I. (1989). A concordance correlation coefficient to evaluate reproducibility. *Biometrics*, **45**, 255–268. Sander, M., Hofstetter, T. B., & Gorski, C. A. (2015). Electrochemical analyses of redox-active iron minerals: a review of nonmediated and mediated approaches. *Environmental Science & Technology*, **49**, 5862–5878. Thompson, A., Chadwick, O. A., Rancourt, D. G., & Chorover, J. (2006). Iron-oxide crystallinity increases during soil redox oscillations. *Geochimica et Cosmochimica Acta*, **70**, 1710–1727. ## Companion manuscripts These describe the framework this package implements. None is published and two are under review; the entries are provisional and should be replaced with the published versions. Ghotbi, M., Ghotbi, M., Komluski, J., & Holtgrewe-Stukenbrock, E. H. HRRI: direction-aware diagnostics for soil–plant–microbiome redox recovery across hydroclimatic disturbances. *In preparation.* Ghotbi, M., Kolody, B. C., Ghotbi, M., & Holtgrewe-Stukenbrock, E. A Theory of Hydroclimatic Redox Resilience. *Submitted to Communications Earth & Environment.* — the source of the capacity, connectivity, kinetics and memory decomposition used throughout this vignette. Ghotbi, M., Ghotbi, M., Mühling, K. H., & Stukenbrock, E. H. Rhizosphere redox recovery after hydrological disturbances: mechanisms across the soil–plant–microbiome continuum. *Submitted to Soil Biology & Biochemistry.*