diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9117220dd..b2e9cb1db 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,7 @@ All notable changes to TEPP are documented here. The format follows Keep a Chang
 
 ## [Unreleased]
 
+- `network_analysis` posterior edge estimation now carries exact two-sided Fisher z-transform p-values against `rho = 0` (erfc evaluated by an all-positive confluent series plus the Laplace continued fraction, locked to libm reference values at 1e-13 relative tolerance), percentile-bootstrap credible intervals and selection fractions over posterior draws, Benjamini–Hochberg step-up admission on those exact p-values instead of the complement of a thresholded fraction, an explicit fail-closed `edge_drop_probability` parameter for consensus co-assignment resampling (replacing a hardcoded 0.1), and honest per-replicate stability admission; the greedy partition helper is renamed to state that it makes no modularity-optimization claim. APA 7 entries added: Benjamini & Hochberg (1995), Efron (1979), Fisher (1921), Hennig (2007), Monti (2003).
 - Restored protected-main gate integrity after the consolidation merges: hourly-scheduler prompt-contract tests now assert the gap-baseline-derived task contract (Gap ID naming, no invented weights) instead of stale increment-specific tokens, the operator-gap register inventory matches the live 33-PR queue, `evidence_core::image_unit` non-image/empty-subtype refusals and `load_union_branch_totals` valid-record accumulation have exact coverage, and the README crate fence plus duplicate registry entries stay deduped.
 - Branch-coverage diagnostics on the post-consolidation head exposed two uncovered outcomes in `evidence_core::image_unit` (`is_image_media_type_token` non-image prefix and empty-subtype refusals), one uncovered authored line (the strip-prefix refusal), and lost valid-record coverage for `load_union_branch_totals`; exact red-to-green cases now cover the non-image/empty-subtype data URIs and per-coordinate True/False accumulation.
 - Repaired post-consolidation merge fallout that left protected `main` red: restored the lost `return True` in the `check_coverage.py` match-guard branch, removed the shadowed duplicate `load_union_branch_totals` and `_is_multiline_match_guard` definitions plus duplicate workspace-crate entries (`episode_membership`, `analysis_engine`) from the contract tuple and Cargo member arrays, split two union-fused four-tuples back into `(variant, message)` pairs in the `event_core` error table, repaired the fused `identity_recovery_rate` body in `episode_membership::window`, deduplicated the checked-arithmetic eligible-count block in `analysis_engine`, fixed four-argument `unit()` test call sites, rebalanced the README crate-list fence around all 54 unique crates, and deduplicated the `location_membership`/`validation_core`/`tepp_api` architecture-table rows. Also documents private `PLAUSIBLE_IMAGE_MEDIA_TYPES` so `cargo doc -D warnings` passes.
diff --git a/crates/network_analysis/src/consensus.rs b/crates/network_analysis/src/consensus.rs
new file mode 100644
index 000000000..2678ee9be
--- /dev/null
+++ b/crates/network_analysis/src/consensus.rs
@@ -0,0 +1,301 @@
+//! Co-assignment consensus clustering from repeated label-invariant partitions.
+//!
+//! Runs multiple rounds of deterministic greedy partitioning over
+//! randomly perturbed admitted positive edges, builds a co-assignment
+//! matrix across replicates, and derives consensus clusters by
+//! thresholding that matrix. The resampling-based consensus view and the
+//! stability rationale follow Monti (2003) and Hennig (2007); the edge
+//! perturbation probability is an explicit parameter with provenance,
+//! never an implicit constant.
+
+#![forbid(unsafe_code)]
+#![allow(
+    clippy::doc_markdown,
+    clippy::items_after_statements,
+    clippy::must_use_candidate,
+    clippy::needless_for_each,
+    clippy::cast_precision_loss,
+    clippy::missing_errors_doc
+)]
+#![deny(missing_docs)]
+
+use crate::edges::NetworkEdge;
+use crate::error::NetworkEstimatorError;
+use std::collections::HashMap;
+
+/// Consensus clustering output.
+///
+/// * `assignments` – one entry per topic: `Some(cluster_id)` or `None` for unclustered.
+/// * `co_assignment` – K × K symmetric matrix of co-assignment frequencies.
+#[derive(Debug, Clone, PartialEq)]
+pub struct ConsensusClusterOutput {
+    /// Per-topic cluster assignment (`None` = unclustered).
+    pub assignments: Vec<Option<usize>>,
+    /// K × K co-assignment frequency matrix.
+    pub co_assignment: Vec<Vec<f64>>,
+}
+
+/// Derive consensus clusters from repeatedly perturbed partitions.
+///
+/// Each replicate independently drops every admitted positive edge with
+/// probability `edge_drop_probability`, repartitions the surviving edges,
+/// and accumulates label-invariant co-assignment counts. The final
+/// assignment thresholds the co-assignment frequency at
+/// `consensus_threshold`.
+///
+/// # Arguments
+///
+/// * `edges` – admitted positive edges (source, target, effect).
+/// * `k_topics` – total number of topics.
+/// * `n_replicates` – number of perturbed partitions to generate; at least 1.
+/// * `consensus_threshold` – minimum co-assignment fraction for
+///   same-cluster; must lie inside [0, 1].
+/// * `edge_drop_probability` – per-edge independent drop probability in
+///   each replicate; must be finite and inside [0, 1). The value is an
+///   explicit design parameter of the resampling scheme (Monti, 2003;
+///   Hennig, 2007), not an internal constant.
+/// * `seed` – deterministic seed.
+///
+/// # Errors
+///
+/// Fails closed on zero topics or replicates, a threshold outside
+/// [0, 1], or a non-finite drop probability at or above 1.
+pub fn consensus_clusters(
+    edges: &[NetworkEdge],
+    k_topics: usize,
+    n_replicates: usize,
+    consensus_threshold: f64,
+    edge_drop_probability: f64,
+    seed: u64,
+) -> Result<ConsensusClusterOutput, NetworkEstimatorError> {
+    if k_topics == 0 {
+        return Err(NetworkEstimatorError::DimensionMismatch);
+    }
+    if n_replicates == 0 {
+        return Err(NetworkEstimatorError::ZeroReplicates);
+    }
+    if !consensus_threshold.is_finite() || !(0.0..=1.0).contains(&consensus_threshold) {
+        return Err(NetworkEstimatorError::InvalidProbability);
+    }
+    if !edge_drop_probability.is_finite() || !(0.0..1.0).contains(&edge_drop_probability) {
+        return Err(NetworkEstimatorError::InvalidProbability);
+    }
+
+    // Build adjacency from positive edges only.
+    let mut adjacency: HashMap<usize, Vec<usize>> = HashMap::new();
+    for edge in edges {
+        if edge.effect > 0.0 {
+            adjacency.entry(edge.source).or_default().push(edge.target);
+            adjacency.entry(edge.target).or_default().push(edge.source);
+        }
+    }
+
+    // Deterministic LCG shared across the crate.
+    let mut state = seed ^ 0x9E37_79B9_7F4A_7C15;
+
+    // Accumulate co-assignment counts.
+    let mut co_count = vec![vec![0_u64; k_topics]; k_topics];
+
+    for _ in 0..n_replicates {
+        // Perturb: drop each edge independently with the caller-supplied
+        // probability so cluster recovery is stress-tested by resampling.
+        let perturbed: Vec<&NetworkEdge> = edges
+            .iter()
+            .filter(|_| {
+                state = state
+                    .wrapping_mul(6_364_136_223_846_793_005)
+                    .wrapping_add(1_442_695_040_888_963_407);
+                let draw = ((state >> 33) as f64 / (u64::MAX >> 33) as f64).min(1.0);
+                draw >= edge_drop_probability
+            })
+            .collect();
+
+        let partition = greedy_union_partition(&perturbed, k_topics);
+
+        for i in 0..k_topics {
+            for j in 0..k_topics {
+                if partition[i] == partition[j] {
+                    co_count[i][j] += 1;
+                }
+            }
+        }
+    }
+
+    // Build co-assignment frequency matrix and derive final assignment.
+    let mut co_freq = vec![vec![0.0_f64; k_topics]; k_topics];
+    let mut assignments: Vec<Option<usize>> = vec![None; k_topics];
+    let mut next_cluster = 0_usize;
+
+    for i in 0..k_topics {
+        if assignments[i].is_some() {
+            continue;
+        }
+        // Find all j ≥ i still unassigned whose co-assignment with i
+        // reaches the threshold.
+        let members: Vec<usize> = (i..k_topics)
+            .filter(|&j| {
+                assignments[j].is_none()
+                    && co_count[i][j] as f64 / n_replicates as f64 >= consensus_threshold
+            })
+            .collect();
+
+        if members.len() < 2 {
+            // Singletons stay unclustered unless they pair strongly.
+            if !members.is_empty()
+                && co_count[i][i] as f64 / n_replicates as f64 >= consensus_threshold
+                && adjacency.contains_key(&i)
+            {
+                assignments[i] = Some(next_cluster);
+                next_cluster += 1;
+            }
+            continue;
+        }
+
+        for &member in &members {
+            assignments[member] = Some(next_cluster);
+        }
+        next_cluster += 1;
+    }
+
+    for i in 0..k_topics {
+        for j in 0..k_topics {
+            co_freq[i][j] = co_count[i][j] as f64 / n_replicates as f64;
+        }
+    }
+
+    Ok(ConsensusClusterOutput {
+        assignments,
+        co_assignment: co_freq,
+    })
+}
+
+/// Greedy union-find partition over edges sorted by descending effect.
+///
+/// This is a deterministic single-pass stand-in for Leiden until a vetted
+/// implementation is adopted; it makes no modularity-optimization claim.
+fn greedy_union_partition(edges: &[&NetworkEdge], k: usize) -> Vec<usize> {
+    let mut parent: Vec<usize> = (0..k).collect();
+
+    fn find(parent: &mut [usize], mut x: usize) -> usize {
+        while parent[x] != x {
+            parent[x] = parent[parent[x]];
+            x = parent[x];
+        }
+        x
+    }
+
+    fn union(parent: &mut [usize], a: usize, b: usize) {
+        let root_a = find(parent, a);
+        let root_b = find(parent, b);
+        if root_a != root_b {
+            parent[root_b] = root_a;
+        }
+    }
+
+    // Sort edges by descending effect.
+    let mut sorted: Vec<&NetworkEdge> = edges.to_vec();
+    sorted.sort_by(|a, b| b.effect.total_cmp(&a.effect));
+
+    for edge in &sorted {
+        union(&mut parent, edge.source, edge.target);
+    }
+
+    // Normalise labels to 0..n_clusters.
+    let mut label_map: HashMap<usize, usize> = HashMap::new();
+    parent
+        .iter()
+        .map(|&root| {
+            let len = label_map.len();
+            *label_map.entry(root).or_insert(len)
+        })
+        .collect()
+}
+
+#[cfg(test)]
+#[allow(clippy::float_cmp)]
+mod tests {
+    use super::*;
+
+    fn edge(source: usize, target: usize, effect: f64) -> NetworkEdge {
+        NetworkEdge {
+            source,
+            target,
+            effect,
+            lower: effect,
+            upper: effect,
+            p_value: 0.0,
+            selection_probability: 1.0,
+        }
+    }
+
+    #[test]
+    fn invalid_parameters_fail_closed() {
+        let chain = vec![edge(0, 1, 0.9), edge(1, 2, 0.8)];
+        assert!(matches!(
+            consensus_clusters(&chain, 0, 5, 0.5, 0.1, 1),
+            Err(NetworkEstimatorError::DimensionMismatch)
+        ));
+        assert!(matches!(
+            consensus_clusters(&chain, 3, 0, 0.5, 0.1, 1),
+            Err(NetworkEstimatorError::ZeroReplicates)
+        ));
+        assert!(matches!(
+            consensus_clusters(&chain, 3, 5, 1.5, 0.1, 1),
+            Err(NetworkEstimatorError::InvalidProbability)
+        ));
+        assert!(matches!(
+            consensus_clusters(&chain, 3, 5, -0.1, 0.1, 1),
+            Err(NetworkEstimatorError::InvalidProbability)
+        ));
+        assert!(matches!(
+            consensus_clusters(&chain, 3, 5, 0.5, 1.0, 1),
+            Err(NetworkEstimatorError::InvalidProbability)
+        ));
+        assert!(matches!(
+            consensus_clusters(&chain, 3, 5, 0.5, f64::NAN, 1),
+            Err(NetworkEstimatorError::InvalidProbability)
+        ));
+    }
+
+    #[test]
+    fn zero_drop_probability_is_fully_deterministic() {
+        // With no perturbation every replicate sees the identical chain,
+        // so all three topics co-assign always and land in one cluster.
+        let chain = vec![edge(0, 1, 0.95), edge(1, 2, 0.85)];
+        let output = consensus_clusters(&chain, 3, 25, 0.99, 0.0, 11).unwrap();
+        assert_eq!(output.assignments[0], output.assignments[1]);
+        assert_eq!(output.assignments[1], output.assignments[2]);
+        assert!((output.co_assignment[0][2] - 1.0).abs() < 1e-12);
+    }
+
+    #[test]
+    fn weak_edges_do_not_merge_under_perturbation() {
+        // Two strongly-bound pairs of topics plus one isolated topic:
+        // within-pair co-assignment stays high across perturbed
+        // replicates, while the isolated topic never joins any cluster.
+        let pairs = vec![edge(0, 1, 0.97), edge(2, 3, 0.96)];
+        let output = consensus_clusters(&pairs, 5, 200, 0.6, 0.2, 21).unwrap();
+        assert_eq!(output.assignments[0], output.assignments[1]);
+        assert!(output.assignments[0].is_some());
+        assert_eq!(output.assignments[2], output.assignments[3]);
+        assert!(output.assignments[2].is_some());
+        assert_ne!(output.assignments[0], output.assignments[2]);
+        assert_eq!(output.assignments[4], None);
+        assert!((output.co_assignment[4][4] - 1.0).abs() < 1e-12);
+        // The within-pair edge survives roughly 80% of perturbations, so
+        // its co-assignment stays well above the consensus threshold
+        // while never reaching a deterministic 1.0.
+        let within_pair = output.co_assignment[0][1];
+        assert!(within_pair > 0.6 && within_pair < 1.0, "co = {within_pair}");
+        // Cross-pair topics never co-assign: they share no edge.
+        assert!((output.co_assignment[0][2]).abs() < 1e-12);
+    }
+
+    #[test]
+    fn identical_inputs_and_seed_reproduce_output_exactly() {
+        let chain = vec![edge(0, 1, 0.9), edge(1, 2, 0.8)];
+        let first = consensus_clusters(&chain, 3, 15, 0.6, 0.1, 7).unwrap();
+        let second = consensus_clusters(&chain, 3, 15, 0.6, 0.1, 7).unwrap();
+        assert_eq!(first, second);
+    }
+}
diff --git a/crates/network_analysis/src/edges.rs b/crates/network_analysis/src/edges.rs
new file mode 100644
index 000000000..c2dedf676
--- /dev/null
+++ b/crates/network_analysis/src/edges.rs
@@ -0,0 +1,590 @@
+//! Posterior log-ratio edge estimation for topic–topic networks.
+//!
+//! Topic proportions are compositional. This module transforms posterior
+//! draws into isometric log-ratio (ILR) coordinates via the sequential
+//! Egozcue basis supplied by `topic_measurement`, then computes draw-level
+//! Pearson correlations in that orthonormal space.
+//!
+//! Every reported quantity traces to an authoritative primary source:
+//!
+//! - The two-sided p-value against `rho = 0` uses the Fisher z-transform
+//!   `z = atanh(r) * sqrt(n - 3)` with the exact normal-tail identity
+//!   `p = erfc(|z| / sqrt(2))` (Fisher, 1921).
+//! - The complementary error function is evaluated with an all-positive
+//!   confluent series for small arguments and the Laplace continued
+//!   fraction for the tail; both branches are locked by known-truth
+//!   reference values with strict tolerances.
+//! - Credible intervals and selection probabilities are percentile
+//!   bootstrap quantities over posterior draws (Efron, 1979).
+//! - Edge admission controls the false discovery rate with the
+//!   Benjamini–Hochberg step-up procedure applied to the exact p-values
+//!   (Benjamini & Hochberg, 1995). No heuristic constants are used.
+
+#![forbid(unsafe_code)]
+#![allow(
+    clippy::doc_markdown,
+    clippy::items_after_statements,
+    clippy::must_use_candidate,
+    clippy::needless_for_each,
+    clippy::cast_precision_loss,
+    clippy::missing_errors_doc
+)]
+#![deny(missing_docs)]
+
+use crate::error::NetworkEstimatorError;
+
+/// One estimated topic–topic association edge.
+///
+/// `effect` is the full-sample Pearson correlation across draws in ILR
+/// space. `lower` / `upper` bound the percentile-bootstrap credible
+/// interval at the requested level. `p_value` is the exact two-sided
+/// Fisher z-test p-value against `rho = 0`. `selection_probability` is
+/// the fraction of bootstrap replicates whose resampled correlation
+/// reaches at least `admission_threshold` in absolute value.
+#[derive(Debug, Clone, PartialEq)]
+pub struct NetworkEdge {
+    /// Zero-based index of the first topic.
+    pub source: usize,
+    /// Zero-based index of the second topic (`source` < `target`).
+    pub target: usize,
+    /// Full-sample correlation in ILR space.
+    pub effect: f64,
+    /// Lower bound of the percentile-bootstrap credible interval.
+    pub lower: f64,
+    /// Upper bound of the percentile-bootstrap credible interval.
+    pub upper: f64,
+    /// Exact two-sided Fisher z-test p-value against `rho = 0`.
+    pub p_value: f64,
+    /// Bootstrap fraction of replicates reaching the admission threshold.
+    pub selection_probability: f64,
+}
+
+/// Complementary error function for non-negative arguments.
+///
+/// Small and moderate arguments use the confluent all-positive series
+/// `erf(x) = (2x/sqrt(pi)) e^{-x^2} sum_k (2x^2)^k / (2k+1)!!`; the
+/// subtraction from one is bounded away from cancellation because the
+/// series branch is only used while `erfc(x) >= 1e-11` relative error
+/// remains negligible. Tail arguments use the Laplace continued fraction
+/// for `int_x_inf e^{-t^2} dt`, evaluated with the modified Lentz
+/// algorithm so every partial denominator stays positive.
+pub(crate) fn erfc_nonnegative(x: f64) -> f64 {
+    const PI_SQRT: f64 = 1.772_453_850_905_516;
+    debug_assert!(x >= 0.0, "erfc_nonnegative requires x >= 0");
+    if x == 0.0 {
+        return 1.0;
+    }
+    if x < 3.0 {
+        // Confluent all-positive series for erf(x); the final subtraction
+        // loses at most ~1e-11 relative accuracy at the branch point,
+        // which known-truth tests lock below.
+        let xx = x * x;
+        let mut term = 1.0_f64;
+        let mut sum = 1.0_f64;
+        let mut k = 0.0_f64;
+        loop {
+            k += 1.0;
+            term *= 2.0 * xx / (2.0 * k + 1.0);
+            sum += term;
+            if term <= sum * 1e-18 || k >= 500.0 {
+                break;
+            }
+        }
+        let erf_x = 2.0 * x * sum * (-xx).exp() / PI_SQRT;
+        return (1.0 - erf_x).max(0.0);
+    }
+    // Laplace continued fraction:
+    // erfc(x) = (e^{-x^2} / sqrt(pi)) / (x + 1/(2x + 2/(x + 3/(2x + ...))))
+    // evaluated bottom-up with a fixed number of levels chosen from x.
+    let xx = x * x;
+    let mut cf = 0.0_f64;
+    // Deeper levels converge geometrically faster as x grows; 200 levels
+    // is far beyond the point where the remainder is below f64 epsilon
+    // even at the smallest tail argument entering this branch.
+    for level in (1..=200usize).rev() {
+        let numerator = level as f64;
+        let denominator = if level % 2 == 1 { 2.0 * x } else { x };
+        cf = numerator / (denominator + cf);
+    }
+    (-xx).exp() / (PI_SQRT * (x + cf))
+}
+
+/// Exact two-sided p-value of the Fisher z-transform test against rho=0.
+///
+/// With `n_obs` paired observations and sample correlation `r`, the
+/// statistic `z = atanh(r) sqrt(n_obs - 3)` is treated as standard
+/// normal under the null (Fisher, 1921), so the two-sided p-value is
+/// `erfc(|z| / sqrt(2))`.
+pub(crate) fn fisher_two_sided_p_value(r: f64, n_obs: usize) -> f64 {
+    let magnitude = r.abs().min(1.0);
+    if magnitude >= 1.0 {
+        return 0.0;
+    }
+    if n_obs < 4 || !magnitude.is_finite() {
+        return 1.0;
+    }
+    let z = magnitude.atanh() * ((n_obs - 3) as f64).sqrt();
+    erfc_nonnegative(z / std::f64::consts::SQRT_2)
+}
+
+/// Type-7 linear-interpolation quantile of an ascending-sorted sample.
+fn sorted_quantile(sorted: &[f64], probability: f64) -> f64 {
+    let last = sorted.len() - 1;
+    let position = probability * last as f64;
+    #[allow(
+        clippy::cast_possible_truncation,
+        clippy::cast_sign_loss,
+        clippy::cast_precision_loss
+    )]
+    let lower_index = position.floor() as usize;
+    let upper_index = (lower_index + 1).min(last);
+    let weight = position - lower_index as f64;
+    sorted[lower_index] * (1.0 - weight) + sorted[upper_index] * weight
+}
+
+/// Compute per-pair posterior edges with bootstrap uncertainty.
+///
+/// For every ILR coordinate pair the function reports the full-sample
+/// correlation, a two-sided exact p-value against `rho = 0`, a
+/// percentile-bootstrap credible interval at `ci_level`, and the
+/// bootstrap selection fraction at `admission_threshold`
+/// (Efron, 1979; Benjamini & Hochberg, 1995 govern later admission).
+///
+/// # Arguments
+///
+/// * `draws` – one row per posterior draw; each row holds the ILR
+///   coordinates for all K topics (length K − 1 after ILR).
+/// * `admission_threshold` – minimum absolute resampled correlation for
+///   a replicate to count toward the selection fraction; must be finite
+///   and non-negative.
+/// * `ci_level` – central credible-interval level in (0, 1), for
+///   example 0.95.
+/// * `n_resamples` – number of bootstrap resamples; at least 1.
+/// * `seed` – deterministic seed for the resampling generator.
+///
+/// # Errors
+///
+/// Fails closed on empty draws, inconsistent dimensions, fewer than 3
+/// observations per coordinate, a non-finite or negative admission
+/// threshold, a `ci_level` outside the open unit interval, or zero
+/// resamples.
+pub fn posterior_correlation_matrix(
+    draws: &[Vec<f64>],
+    admission_threshold: f64,
+    ci_level: f64,
+    n_resamples: usize,
+    seed: u64,
+) -> Result<Vec<NetworkEdge>, NetworkEstimatorError> {
+    if draws.is_empty() {
+        return Err(NetworkEstimatorError::EmptyDraws);
+    }
+    if !admission_threshold.is_finite() || admission_threshold < 0.0 {
+        return Err(NetworkEstimatorError::InvalidThreshold);
+    }
+    if !(ci_level.is_finite() && ci_level > 0.0 && ci_level < 1.0) {
+        return Err(NetworkEstimatorError::InvalidConfidenceLevel);
+    }
+    if n_resamples == 0 {
+        return Err(NetworkEstimatorError::ZeroReplicates);
+    }
+    let dim = draws[0].len();
+    if dim < 2 || draws.iter().any(|row| row.len() != dim) {
+        return Err(NetworkEstimatorError::DimensionMismatch);
+    }
+
+    let n_draws = draws.len();
+    let k = dim;
+    let mut cols = vec![Vec::with_capacity(n_draws); k];
+    for row in draws {
+        for (coordinate, &value) in row.iter().enumerate() {
+            cols[coordinate].push(value);
+        }
+    }
+
+    let full_rs = cross_draw_correlations(&cols)?;
+    let n_pairs = full_rs.len();
+
+    // Deterministic resampling of draw indices with replacement. The
+    // LCG constants are the widely published PCG-style multiplier and
+    // increment used elsewhere in this crate for reproducibility.
+    let mut state = seed ^ 0x9E37_79B9_7F4A_7C15;
+    let mut resample = |out: &mut Vec<usize>| {
+        for slot in out.iter_mut() {
+            state = state
+                .wrapping_mul(6_364_136_223_846_793_005)
+                .wrapping_add(1_442_695_040_888_963_407);
+            *slot = ((state >> 33) as usize) % n_draws;
+        }
+    };
+
+    let mut boot_matrix = vec![Vec::<f64>::with_capacity(n_resamples); n_pairs];
+    let mut indices = vec![0_usize; n_draws];
+    for _ in 0..n_resamples {
+        resample(&mut indices);
+        let boot_cols: Vec<Vec<f64>> = cols
+            .iter()
+            .map(|column| indices.iter().map(|&index| column[index]).collect())
+            .collect();
+        let rs = cross_draw_correlations(&boot_cols)?;
+        for (pair, value) in rs.into_iter().enumerate() {
+            boot_matrix[pair].push(value);
+        }
+    }
+
+    let alpha_tail = (1.0 - ci_level) / 2.0;
+    let mut edges = Vec::with_capacity(n_pairs);
+    let mut pair = 0_usize;
+    for i in 0..k {
+        for j in (i + 1)..k {
+            let mut sorted = boot_matrix[pair].clone();
+            sorted.sort_by(f64::total_cmp);
+            let selected = sorted
+                .iter()
+                .filter(|&&value| value.abs() >= admission_threshold)
+                .count();
+            edges.push(NetworkEdge {
+                source: i,
+                target: j,
+                effect: full_rs[pair],
+                lower: sorted_quantile(&sorted, alpha_tail),
+                upper: sorted_quantile(&sorted, 1.0 - alpha_tail),
+                p_value: fisher_two_sided_p_value(full_rs[pair], n_draws),
+                selection_probability: selected as f64 / n_resamples as f64,
+            });
+            pair += 1;
+        }
+    }
+    Ok(edges)
+}
+
+/// Apply the primary multiplicity-corrected admission policy.
+///
+/// An edge is admitted when its Benjamini–Hochberg step-up criterion
+/// passes at `fdr_alpha` over the exact per-edge p-values, its
+/// credible interval excludes zero, and its selection fraction reaches
+/// `min_selection_probability` (Benjamini & Hochberg, 1995).
+pub fn admit_edges(
+    mut edges: Vec<NetworkEdge>,
+    min_selection_probability: f64,
+    fdr_alpha: f64,
+) -> Vec<NetworkEdge> {
+    edges.sort_by(|a, b| a.p_value.total_cmp(&b.p_value));
+    let total = edges.len();
+    if total > 0 {
+        let mut largest_passing_rank = 0_usize;
+        for (rank_zero, edge) in edges.iter().enumerate() {
+            let rank = rank_zero + 1;
+            let critical = fdr_alpha * rank as f64 / total as f64;
+            if edge.p_value <= critical {
+                largest_passing_rank = rank;
+            }
+        }
+        edges.truncate(largest_passing_rank);
+    }
+    edges.retain(|edge| {
+        edge.lower * edge.upper > 0.0 && edge.selection_probability >= min_selection_probability
+    });
+    edges.sort_by_key(|edge| (edge.source, edge.target));
+    edges
+}
+
+/// Apply the per-replicate admission rule used inside stability scoring.
+///
+/// A replicate admits an edge when its absolute correlation reaches
+/// `admission_threshold` and its exact Fisher p-value survives the
+/// Benjamini–Hochberg step-up at `fdr_alpha`. Credible intervals are not
+/// defined inside a single replicate, so they play no role here.
+pub fn admit_edges_within_replicate(
+    mut edges: Vec<NetworkEdge>,
+    admission_threshold: f64,
+    fdr_alpha: f64,
+) -> Vec<NetworkEdge> {
+    edges.retain(|edge| edge.effect.abs() >= admission_threshold);
+    edges.sort_by(|a, b| a.p_value.total_cmp(&b.p_value));
+    let total = edges.len();
+    if total > 0 {
+        let mut largest_passing_rank = 0_usize;
+        for (rank_zero, edge) in edges.iter().enumerate() {
+            let rank = rank_zero + 1;
+            let critical = fdr_alpha * rank as f64 / total as f64;
+            if edge.p_value <= critical {
+                largest_passing_rank = rank;
+            }
+        }
+        edges.truncate(largest_passing_rank);
+    }
+    edges
+}
+
+/// Compute cross-draw Pearson correlations for every coordinate pair.
+///
+/// Called once per replicate; produces one correlation per pair.
+pub(crate) fn cross_draw_correlations(
+    ilr_columns: &[Vec<f64>],
+) -> Result<Vec<f64>, NetworkEstimatorError> {
+    let k = ilr_columns.len();
+    if k < 2 {
+        return Err(NetworkEstimatorError::DimensionMismatch);
+    }
+    let n_obs = ilr_columns[0].len();
+    if n_obs < 3 {
+        return Err(NetworkEstimatorError::InsufficientObservations);
+    }
+    for column in ilr_columns {
+        if column.len() != n_obs {
+            return Err(NetworkEstimatorError::DimensionMismatch);
+        }
+    }
+    let mut out = Vec::with_capacity(k * (k - 1) / 2);
+    for i in 0..k {
+        for j in (i + 1)..k {
+            out.push(pearson(&ilr_columns[i], &ilr_columns[j]));
+        }
+    }
+    Ok(out)
+}
+
+/// Plain Pearson correlation between two equal-length slices.
+fn pearson(a: &[f64], b: &[f64]) -> f64 {
+    let n = a.len() as f64;
+    let mean_a = a.iter().sum::<f64>() / n;
+    let mean_b = b.iter().sum::<f64>() / n;
+    let mut covariance = 0.0;
+    let mut variance_a = 0.0;
+    let mut variance_b = 0.0;
+    for index in 0..a.len() {
+        let da = a[index] - mean_a;
+        let db = b[index] - mean_b;
+        covariance += da * db;
+        variance_a += da * da;
+        variance_b += db * db;
+    }
+    let denominator = (variance_a * variance_b).sqrt();
+    if denominator < f64::EPSILON {
+        0.0
+    } else {
+        covariance / denominator
+    }
+}
+
+#[cfg(test)]
+#[allow(clippy::float_cmp)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn empty_draws_fail_closed() {
+        let result = posterior_correlation_matrix(&[], 0.3, 0.95, 10, 1);
+        assert!(matches!(result, Err(NetworkEstimatorError::EmptyDraws)));
+    }
+
+    #[test]
+    fn dimension_mismatch_detected() {
+        let draws = vec![vec![1.0, 2.0], vec![1.0]];
+        let result = posterior_correlation_matrix(&draws, 0.3, 0.95, 10, 1);
+        assert!(matches!(
+            result,
+            Err(NetworkEstimatorError::DimensionMismatch)
+        ));
+    }
+
+    #[test]
+    fn invalid_inputs_fail_closed() {
+        let draws = vec![vec![1.0, 2.0]; 12];
+        assert!(matches!(
+            posterior_correlation_matrix(&draws, -0.1, 0.95, 10, 1),
+            Err(NetworkEstimatorError::InvalidThreshold)
+        ));
+        assert!(matches!(
+            posterior_correlation_matrix(&draws, f64::NAN, 0.95, 10, 1),
+            Err(NetworkEstimatorError::InvalidThreshold)
+        ));
+        assert!(matches!(
+            posterior_correlation_matrix(&draws, 0.3, 1.0, 10, 1),
+            Err(NetworkEstimatorError::InvalidConfidenceLevel)
+        ));
+        assert!(matches!(
+            posterior_correlation_matrix(&draws, 0.3, 0.0, 10, 1),
+            Err(NetworkEstimatorError::InvalidConfidenceLevel)
+        ));
+        assert!(matches!(
+            posterior_correlation_matrix(&draws, 0.3, 0.95, 0, 1),
+            Err(NetworkEstimatorError::ZeroReplicates)
+        ));
+        let degenerate = vec![vec![1.0]];
+        assert!(matches!(
+            posterior_correlation_matrix(&degenerate, 0.3, 0.95, 10, 1),
+            Err(NetworkEstimatorError::DimensionMismatch)
+        ));
+    }
+
+    #[test]
+    fn erfc_matches_known_truth_reference_values() {
+        // Reference values from the platform libm double-precision erfc
+        // (the authoritative standard implementation); each is locked
+        // with a strict relative tolerance so any wrong constant or
+        // truncated series fails loudly instead of silently shifting
+        // p-values. Values below 3e-16 relative to one are not usable
+        // anchors in f64 and start from 6.0 downward only as smoke.
+        let references = [
+            (0.25_f64, 0.723_673_609_831_763_1),
+            (0.5, 0.479_500_122_186_953_5),
+            (1.0, 0.157_299_207_050_285_16),
+            (1.5, 0.033_894_853_524_689_274),
+            (2.0, 0.004_677_734_981_047_264_5),
+            (3.0, 2.209_049_699_858_543_8e-5),
+            (4.0, 1.541_725_790_028_002e-8),
+            (5.0, 1.537_459_794_428_035e-12),
+            (6.0, 2.151_973_671_249_891_3e-17),
+        ];
+        for (argument, expected) in references {
+            let computed = erfc_nonnegative(argument);
+            let relative_error = if expected == 0.0 {
+                computed.abs()
+            } else {
+                ((computed - expected) / expected).abs()
+            };
+            assert!(
+                relative_error < 1e-13,
+                "erfc({argument}) = {computed}, expected {expected}"
+            );
+        }
+        assert_eq!(
+            erfc_nonnegative(0.0).total_cmp(&1.0),
+            std::cmp::Ordering::Equal
+        );
+    }
+
+    #[test]
+    fn fisher_p_value_hits_exact_normal_anchor() {
+        // Choose n so sqrt(n-3) = 10 and r such that atanh(r)*10 equals
+        // the exact 97.5% normal quantile 1.959963985..., whose two-sided
+        // p-value is 0.05 by construction of that quantile.
+        let quantile_975 = 1.959_963_984_540_054_f64;
+        let r = (quantile_975 / 10.0).tanh();
+        let p = fisher_two_sided_p_value(r, 103);
+        assert!((p - 0.05).abs() < 1e-9, "p = {p}");
+        assert_eq!(
+            fisher_two_sided_p_value(0.0, 30).total_cmp(&1.0),
+            std::cmp::Ordering::Equal
+        );
+        assert_eq!(
+            fisher_two_sided_p_value(1.0, 30),
+            fisher_two_sided_p_value(-1.0, 30)
+        );
+        assert_eq!(
+            fisher_two_sided_p_value(1.0, 30).total_cmp(&f64::EPSILON),
+            std::cmp::Ordering::Less
+        );
+        let positive = fisher_two_sided_p_value(0.4, 50);
+        let negative = fisher_two_sided_p_value(-0.4, 50);
+        assert_eq!(positive.total_cmp(&negative), std::cmp::Ordering::Equal);
+        assert!(positive < fisher_two_sided_p_value(0.2, 50));
+    }
+
+    #[test]
+    fn known_positive_and_negative_correlations_are_recovered() {
+        let mut draws = Vec::with_capacity(60);
+        for i in 0..60_i32 {
+            let x = f64::from(i) / 60.0;
+            draws.push(vec![x, x + 0.01 * (f64::from(i).sin()), -x]);
+        }
+        let edges = posterior_correlation_matrix(&draws, 0.5, 0.95, 40, 7).unwrap();
+        assert!(edges[0].effect > 0.99, "r(0,1) = {}", edges[0].effect);
+        assert!(edges[1].effect < -0.99, "r(0,2) = {}", edges[1].effect);
+        // Strong signal must carry decisive evidence on every axis.
+        assert!(edges[0].p_value < 1e-20, "p(0,1) = {}", edges[0].p_value);
+        assert!((edges[0].selection_probability - 1.0).abs() < f64::EPSILON);
+        assert!(edges[0].lower > 0.9 && edges[0].upper > edges[0].lower);
+        // Independent coordinates keep large p-values.
+        let independent = posterior_correlation_matrix(
+            &(0..80)
+                .map(|i| {
+                    let x = f64::from(i);
+                    vec![x.sin(), (2.0 * x).cos(), 3.0 * x + 0.001 * (7.0 * x).sin()]
+                })
+                .collect::<Vec<_>>(),
+            0.5,
+            0.95,
+            20,
+            11,
+        )
+        .unwrap();
+        assert!(
+            independent[0].p_value > 1e-3,
+            "p = {}",
+            independent[0].p_value
+        );
+    }
+
+    #[test]
+    fn bootstrap_is_deterministic_per_seed() {
+        let draws: Vec<Vec<f64>> = (0..40)
+            .map(|i| {
+                let x = f64::from(i) / 40.0;
+                vec![x, x + 1e-3, -x]
+            })
+            .collect();
+        let first = posterior_correlation_matrix(&draws, 0.5, 0.95, 25, 42).unwrap();
+        let second = posterior_correlation_matrix(&draws, 0.5, 0.95, 25, 42).unwrap();
+        assert_eq!(first, second);
+        let other_seed = posterior_correlation_matrix(&draws, 0.5, 0.95, 25, 43).unwrap();
+        // Different seeds may perturb bootstrap bounds but never the
+        // deterministic point estimate or exact p-value.
+        assert_eq!(first[0].effect, other_seed[0].effect);
+        assert_eq!(first[0].p_value, other_seed[0].p_value);
+    }
+
+    #[test]
+    fn benjamin_hochberg_admits_only_surviving_prefix() {
+        let make_edge = |source: usize, p: f64| NetworkEdge {
+            source,
+            target: source + 10,
+            effect: 0.9,
+            lower: 0.5,
+            upper: 0.99,
+            p_value: p,
+            selection_probability: 1.0,
+        };
+        let edges = vec![
+            make_edge(0, 0.001),
+            make_edge(1, 0.008),
+            make_edge(2, 0.039),
+            make_edge(3, 0.041),
+            make_edge(4, 0.2),
+        ];
+        // Sorted p-values [0.001, 0.008, 0.039, 0.041, 0.2] against
+        // critical values alpha*i/5 at alpha = 0.05: rank 1 passes
+        // (0.001 <= 0.01), rank 2 passes (0.008 <= 0.02), ranks 3-5 all
+        // fail (0.039 > 0.03), so the step-up keeps exactly the first
+        // two — the largest prefix with every member passing.
+        let admitted = admit_edges(edges, 0.5, 0.05);
+        assert_eq!(admitted.len(), 2);
+        assert_eq!(admitted[0].source, 0);
+        assert_eq!(admitted[1].source, 1);
+        // A stricter selection floor removes edges independently.
+        let strict = admit_edges(vec![make_edge(0, 0.001), make_edge(1, 0.008)], 1.01, 0.05);
+        assert!(strict.is_empty());
+    }
+
+    #[test]
+    fn within_replicate_rule_combines_threshold_and_fdr() {
+        let make_edge = |source: usize, effect: f64, p: f64| NetworkEdge {
+            source,
+            target: source + 10,
+            effect,
+            lower: effect,
+            upper: effect,
+            p_value: p,
+            selection_probability: 1.0,
+        };
+        let edges = vec![
+            make_edge(0, 0.9, 0.0001),
+            make_edge(1, 0.2, 0.0001),
+            make_edge(2, 0.8, 0.9),
+        ];
+        let admitted = admit_edges_within_replicate(edges, 0.5, 0.05);
+        assert_eq!(admitted.len(), 1);
+        assert_eq!(admitted[0].source, 0);
+    }
+}
diff --git a/crates/network_analysis/src/error.rs b/crates/network_analysis/src/error.rs
index ab4256482..0336ba27c 100644
--- a/crates/network_analysis/src/error.rs
+++ b/crates/network_analysis/src/error.rs
@@ -27,9 +27,49 @@ impl fmt::Display for NetworkError {
 
 impl std::error::Error for NetworkError {}
 
+/// Fail-closed errors from the posterior network estimator.
+#[derive(Clone, Copy, Debug, Eq, PartialEq)]
+#[non_exhaustive]
+pub enum NetworkEstimatorError {
+    /// Posterior draw set is empty.
+    EmptyDraws,
+    /// Coordinate dimensions are inconsistent across draws or too small.
+    DimensionMismatch,
+    /// Fewer than 3 ILR coordinates available for pairwise correlation.
+    InsufficientCoordinates,
+    /// Fewer than 3 observations per coordinate for correlation.
+    InsufficientObservations,
+    /// Bootstrap replicates must be at least 1.
+    ZeroReplicates,
+    /// An admission threshold must be finite and non-negative.
+    InvalidThreshold,
+    /// A credible-interval level must lie strictly inside (0, 1).
+    InvalidConfidenceLevel,
+    /// A probability parameter must be finite and inside [0, 1).
+    InvalidProbability,
+}
+
+impl fmt::Display for NetworkEstimatorError {
+    fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
+        let message = match self {
+            Self::EmptyDraws => "posterior draw set is empty",
+            Self::DimensionMismatch => "coordinate dimensions are inconsistent",
+            Self::InsufficientCoordinates => "fewer than 3 ILR coordinates for correlation",
+            Self::InsufficientObservations => "fewer than 3 observations per coordinate",
+            Self::ZeroReplicates => "bootstrap replicates must be at least 1",
+            Self::InvalidThreshold => "admission threshold must be finite and non-negative",
+            Self::InvalidConfidenceLevel => "credible-interval level must be inside (0, 1)",
+            Self::InvalidProbability => "probability parameter must be inside [0, 1)",
+        };
+        formatter.write_str(message)
+    }
+}
+
+impl std::error::Error for NetworkEstimatorError {}
+
 #[cfg(test)]
 mod tests {
-    use super::NetworkError;
+    use super::{NetworkError, NetworkEstimatorError};
 
     #[test]
     fn error_messages_are_stable() {
@@ -50,4 +90,44 @@ mod tests {
             assert_eq!(error.to_string(), message);
         }
     }
+
+    #[test]
+    fn estimator_error_messages_are_stable() {
+        for (error, message) in [
+            (
+                NetworkEstimatorError::EmptyDraws,
+                "posterior draw set is empty",
+            ),
+            (
+                NetworkEstimatorError::DimensionMismatch,
+                "coordinate dimensions are inconsistent",
+            ),
+            (
+                NetworkEstimatorError::InsufficientCoordinates,
+                "fewer than 3 ILR coordinates for correlation",
+            ),
+            (
+                NetworkEstimatorError::InsufficientObservations,
+                "fewer than 3 observations per coordinate",
+            ),
+            (
+                NetworkEstimatorError::ZeroReplicates,
+                "bootstrap replicates must be at least 1",
+            ),
+            (
+                NetworkEstimatorError::InvalidThreshold,
+                "admission threshold must be finite and non-negative",
+            ),
+            (
+                NetworkEstimatorError::InvalidConfidenceLevel,
+                "credible-interval level must be inside (0, 1)",
+            ),
+            (
+                NetworkEstimatorError::InvalidProbability,
+                "probability parameter must be inside [0, 1)",
+            ),
+        ] {
+            assert_eq!(error.to_string(), message);
+        }
+    }
 }
diff --git a/crates/network_analysis/src/lib.rs b/crates/network_analysis/src/lib.rs
index 37a8d7fed..dfa99347b 100644
--- a/crates/network_analysis/src/lib.rs
+++ b/crates/network_analysis/src/lib.rs
@@ -1,14 +1,21 @@
 #![forbid(unsafe_code)]
 #![deny(missing_docs)]
-//! Compositional network and clustering gates for TRSL-TM outputs.
+//! Compositional network, clustering, and posterior-association analysis
+//! for TRSL-TM outputs.
 //!
 //! Raw topic proportions are not ordinary Euclidean coordinates. Cluster
 //! recovery is scored with label-invariant pair precision and recall against
-//! known truth (ADR 0005/0012).
+//! known truth (ADR 0005/0012). The posterior estimator transforms draws into
+//! isometric log-ratio coordinates, estimates topic–topic associations with
+//! uncertainty and bootstrap stability, applies multiplicity-corrected edge
+//! admission, and derives co-assignment consensus clusters.
 
 mod cluster;
+mod consensus;
+mod edges;
 mod error;
 mod simplex;
+mod stability;
 
 /// Opaque cluster identity.
 pub use cluster::ClusterLabel;
@@ -16,7 +23,25 @@ pub use cluster::ClusterLabel;
 pub use cluster::cluster_pair_precision;
 /// Pair recall of recovered clusters.
 pub use cluster::cluster_pair_recall;
+/// Consensus clustering output.
+pub use consensus::ConsensusClusterOutput;
+/// Derive consensus clusters from repeated partitions.
+pub use consensus::consensus_clusters;
+/// One estimated topic–topic association edge.
+pub use edges::NetworkEdge;
+/// Apply multiplicity-corrected edge-admission policy.
+pub use edges::admit_edges;
+/// Apply the per-replicate admission rule used inside stability scoring.
+pub use edges::admit_edges_within_replicate;
+/// Compute posterior correlation matrix from ILR draws.
+pub use edges::posterior_correlation_matrix;
 /// Fail-closed network-analysis errors.
 pub use error::NetworkError;
+/// Fail-closed posterior-network-estimator errors.
+pub use error::NetworkEstimatorError;
 /// Refuse raw simplex proportions as Euclidean coordinates.
 pub use simplex::refuse_raw_simplex_as_euclidean;
+/// Per-edge bootstrap stability score.
+pub use stability::BootstrapEdgeStability;
+/// Run bootstrap replicates and return stability scores.
+pub use stability::bootstrap_edge_stability;
diff --git a/crates/network_analysis/src/stability.rs b/crates/network_analysis/src/stability.rs
new file mode 100644
index 000000000..0a49672a5
--- /dev/null
+++ b/crates/network_analysis/src/stability.rs
@@ -0,0 +1,269 @@
+//! Bootstrap edge-stability scoring.
+//!
+//! Resamples posterior draws with replacement, recomputes cross-draw
+//! correlations, admits edges inside each replicate with the exact
+//! Fisher p-value under Benjamini–Hochberg control plus the absolute
+//! threshold rule, and reports the fraction of replicates in which each
+//! edge is admitted with a sign consistent with the full-sample
+//! estimate. Resampling-based stability assessment follows the
+//! cluster-stability framework of Hennig (2007) and the consensus
+//! resampling view of Monti (2003).
+
+#![forbid(unsafe_code)]
+#![allow(
+    clippy::doc_markdown,
+    clippy::items_after_statements,
+    clippy::must_use_candidate,
+    clippy::needless_for_each,
+    clippy::cast_precision_loss,
+    clippy::missing_errors_doc
+)]
+#![deny(missing_docs)]
+
+use crate::edges::{
+    NetworkEdge, admit_edges_within_replicate, cross_draw_correlations, fisher_two_sided_p_value,
+};
+use crate::error::NetworkEstimatorError;
+
+/// Bootstrap stability result for one candidate edge.
+#[derive(Debug, Clone, PartialEq)]
+pub struct BootstrapEdgeStability {
+    /// Zero-based index of the first topic.
+    pub source: usize,
+    /// Zero-based index of the second topic.
+    pub target: usize,
+    /// Fraction of bootstrap replicates where the edge is admitted and
+    /// its correlation sign matches the full-sample estimate.
+    pub stability: f64,
+}
+
+/// Run `n_replicates` nonparametric bootstrap replicates over posterior
+/// ILR draws and return per-edge stability scores.
+///
+/// # Arguments
+///
+/// * `ilr_draws` – one row per posterior draw; columns are ILR coordinates.
+/// * `admission_threshold` – minimum absolute resampled correlation for a
+///   replicate to admit an edge; must be finite and non-negative.
+/// * `fdr_alpha` – Benjamini–Hochberg level applied to the exact per-pair
+///   p-values inside every replicate.
+/// * `n_replicates` – number of bootstrap resamples; at least 1.
+/// * `seed` – deterministic seed for reproducibility.
+///
+/// # Errors
+///
+/// Fails closed on empty draws, inconsistent dimensions, fewer than 3
+/// observations, an invalid threshold, or zero replicates.
+pub fn bootstrap_edge_stability(
+    ilr_draws: &[Vec<f64>],
+    admission_threshold: f64,
+    fdr_alpha: f64,
+    n_replicates: usize,
+    seed: u64,
+) -> Result<Vec<BootstrapEdgeStability>, NetworkEstimatorError> {
+    if ilr_draws.is_empty() || ilr_draws[0].is_empty() {
+        return Err(NetworkEstimatorError::EmptyDraws);
+    }
+    if !admission_threshold.is_finite() || admission_threshold < 0.0 {
+        return Err(NetworkEstimatorError::InvalidThreshold);
+    }
+    if n_replicates == 0 {
+        return Err(NetworkEstimatorError::ZeroReplicates);
+    }
+    let k = ilr_draws[0].len();
+    let d = ilr_draws.len();
+
+    // Transpose to column-major: columns[c] = values of coordinate c across draws.
+    let mut columns = vec![Vec::with_capacity(d); k];
+    for row in ilr_draws {
+        for (coordinate, &value) in row.iter().enumerate() {
+            columns[coordinate].push(value);
+        }
+    }
+
+    let full_rs = cross_draw_correlations(&columns)?;
+
+    // Simple deterministic LCG PRNG shared across the crate.
+    let mut state = seed ^ 0x9E37_79B9_7F4A_7C15;
+
+    let mut admitted_count = vec![0_u64; k * (k - 1) / 2];
+    let mut sign_match_count = vec![0_u64; k * (k - 1) / 2];
+
+    for _replicate in 0..n_replicates {
+        let indices: Vec<usize> = (0..d)
+            .map(|_| {
+                state = state
+                    .wrapping_mul(6_364_136_223_846_793_005)
+                    .wrapping_add(1_442_695_040_888_963_407);
+                ((state >> 33) as usize) % d
+            })
+            .collect();
+
+        let boot_columns: Vec<Vec<f64>> = columns
+            .iter()
+            .map(|column| indices.iter().map(|&index| column[index]).collect())
+            .collect();
+
+        let rs = cross_draw_correlations(&boot_columns)?;
+        let replicate_edges = build_replicate_edges(&rs, d);
+        let admitted =
+            admit_edges_within_replicate(replicate_edges, admission_threshold, fdr_alpha);
+        let admitted_set: std::collections::HashSet<(usize, usize)> = admitted
+            .iter()
+            .map(|edge| (edge.source, edge.target))
+            .collect();
+
+        let mut pair = 0_usize;
+        for i in 0..k {
+            for j in (i + 1)..k {
+                if admitted_set.contains(&(i, j)) {
+                    admitted_count[pair] += 1;
+                    let full_sign = full_rs[pair].signum();
+                    let replicate_sign = rs[pair].signum();
+                    #[allow(clippy::float_cmp)]
+                    if full_sign != 0.0 && replicate_sign == full_sign {
+                        sign_match_count[pair] += 1;
+                    }
+                }
+                pair += 1;
+            }
+        }
+    }
+
+    let mut out = Vec::with_capacity(k * (k - 1) / 2);
+    let mut pair = 0_usize;
+    for i in 0..k {
+        for j in (i + 1)..k {
+            out.push(BootstrapEdgeStability {
+                source: i,
+                target: j,
+                stability: if admitted_count[pair] > 0 {
+                    sign_match_count[pair] as f64 / admitted_count[pair] as f64
+                } else {
+                    0.0
+                },
+            });
+            pair += 1;
+        }
+    }
+    Ok(out)
+}
+
+/// Construct candidate edges from one bootstrap replicate's correlations.
+///
+/// Each edge carries its exact Fisher p-value against `rho = 0`
+/// computed at the draw count of the original sample; interval bounds
+/// are undefined inside a replicate, so they mirror the effect and are
+/// never consulted by [`admit_edges_within_replicate`].
+fn build_replicate_edges(rs: &[f64], n_obs: usize) -> Vec<NetworkEdge> {
+    // Smallest k with k(k-1)/2 >= pair count; exact for well-formed input
+    // via the quadratic formula k = (1 + sqrt(1 + 8 * pairs)) / 2.
+    #[allow(clippy::manual_midpoint)]
+    let approximated = ((8.0 * rs.len() as f64 + 1.0).sqrt() + 1.0) / 2.0;
+    #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
+    let truncated = approximated as usize;
+    let base = truncated.max(2);
+    let k = if base * (base - 1) / 2 < rs.len() {
+        base + 1
+    } else {
+        base
+    };
+    let mut edges = Vec::with_capacity(rs.len());
+    let mut pair = 0_usize;
+    for i in 0..k {
+        for j in (i + 1)..k {
+            edges.push(NetworkEdge {
+                source: i,
+                target: j,
+                effect: rs[pair],
+                lower: rs[pair],
+                upper: rs[pair],
+                p_value: fisher_two_sided_p_value(rs[pair], n_obs),
+                selection_probability: 1.0,
+            });
+            pair += 1;
+        }
+    }
+    edges
+}
+
+#[cfg(test)]
+#[allow(clippy::float_cmp)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn invalid_inputs_fail_closed() {
+        assert!(matches!(
+            bootstrap_edge_stability(&[], 0.5, 0.05, 10, 1),
+            Err(NetworkEstimatorError::EmptyDraws)
+        ));
+        assert!(matches!(
+            bootstrap_edge_stability(&[vec![1.0]], 0.5, 0.05, 10, 1),
+            Err(NetworkEstimatorError::DimensionMismatch)
+        ));
+        let draws = vec![vec![1.0, 2.0]; 12];
+        assert!(matches!(
+            bootstrap_edge_stability(&draws, -0.5, 0.05, 10, 1),
+            Err(NetworkEstimatorError::InvalidThreshold)
+        ));
+        assert!(matches!(
+            bootstrap_edge_stability(&draws, 0.5, 0.05, 0, 1),
+            Err(NetworkEstimatorError::ZeroReplicates)
+        ));
+    }
+
+    #[test]
+    fn strong_edges_are_perfectly_stable_and_weak_edges_are_not_admitted() {
+        let mut draws = Vec::with_capacity(60);
+        for i in 0..60_i32 {
+            let x = f64::from(i) / 60.0;
+            draws.push(vec![x, x + 1e-9, -x]);
+        }
+        let scores = bootstrap_edge_stability(&draws, 0.9, 0.05, 30, 5).unwrap();
+        assert_eq!(scores.len(), 3);
+        // Pair (0,1): r ≈ +1 survives every replicate with matching sign.
+        assert!(
+            scores[0].stability > 0.99,
+            "stability(0,1) = {}",
+            scores[0].stability
+        );
+        // Pair (0,2): r ≈ −1 also admitted; sign consistency counts both.
+        assert!(scores[1].stability > 0.99);
+        // Pair (1,2): r ≈ −1 as well under this construction? x vs −(x+eps).
+        assert!(scores[2].stability > 0.99);
+    }
+
+    #[test]
+    fn independent_noise_never_reaches_stable_admission() {
+        let draws: Vec<Vec<f64>> = (0..50)
+            .map(|i| {
+                let x = f64::from(i);
+                vec![x.sin(), (3.7 * x).cos(), (9.1 * x).sin()]
+            })
+            .collect();
+        let scores = bootstrap_edge_stability(&draws, 0.95, 0.01, 25, 9).unwrap();
+        for score in &scores {
+            assert!(
+                score.stability < 0.5,
+                "unexpected stable noise edge {} -> {}: {}",
+                score.source,
+                score.target,
+                score.stability
+            );
+        }
+    }
+
+    #[test]
+    fn bootstrap_is_deterministic_per_seed() {
+        let draws: Vec<Vec<f64>> = (0..40)
+            .map(|i| {
+                let x = f64::from(i) / 40.0;
+                vec![x, x + 1e-6, -x]
+            })
+            .collect();
+        let first = bootstrap_edge_stability(&draws, 0.9, 0.05, 20, 42).unwrap();
+        let second = bootstrap_edge_stability(&draws, 0.9, 0.05, 20, 42).unwrap();
+        assert_eq!(first, second);
+    }
+}
diff --git a/crates/role_contradiction/src/lib.rs b/crates/role_contradiction/src/lib.rs
index 808a9c108..67cdf1dc9 100644
--- a/crates/role_contradiction/src/lib.rs
+++ b/crates/role_contradiction/src/lib.rs
@@ -12,6 +12,8 @@ mod role;
 
 /// Fail-closed role-contradiction errors.
 pub use error::RoleContradictionError;
+/// Closed vocabulary of commercial roles that can change over time.
+pub use role::ContextualRole;
 /// Fraction of recovered contextual roles that match known truth.
 pub use role::identity_recovery_rate;
 /// Refuse a contradictory customer/competitor pair in one group.
@@ -20,5 +22,3 @@ pub use role::refuse_contradictory_roles;
 pub use role::refuse_role_as_entity_class;
 /// Return whether two roles contradict in the same group.
 pub use role::roles_contradict;
-/// Closed vocabulary of commercial roles that can change over time.
-pub use role::ContextualRole;
diff --git a/crates/role_contradiction/src/role.rs b/crates/role_contradiction/src/role.rs
index 8d37287b4..a893b32d5 100644
--- a/crates/role_contradiction/src/role.rs
+++ b/crates/role_contradiction/src/role.rs
@@ -104,8 +104,8 @@ pub fn identity_recovery_rate(
 #[cfg(test)]
 mod tests {
     use super::{
-        identity_recovery_rate, refuse_contradictory_roles, refuse_role_as_entity_class,
-        roles_contradict, ContextualRole,
+        ContextualRole, identity_recovery_rate, refuse_contradictory_roles,
+        refuse_role_as_entity_class, roles_contradict,
     };
     use crate::RoleContradictionError;
 
diff --git a/crates/role_contradiction/tests/role_contradiction_contract.rs b/crates/role_contradiction/tests/role_contradiction_contract.rs
index 0cf83f46c..7bf4628e4 100644
--- a/crates/role_contradiction/tests/role_contradiction_contract.rs
+++ b/crates/role_contradiction/tests/role_contradiction_contract.rs
@@ -1,8 +1,8 @@
 //! Customer and competitor cannot occupy the same group at once.
 
 use role_contradiction::{
-    identity_recovery_rate, refuse_contradictory_roles, refuse_role_as_entity_class,
-    roles_contradict, ContextualRole, RoleContradictionError,
+    ContextualRole, RoleContradictionError, identity_recovery_rate, refuse_contradictory_roles,
+    refuse_role_as_entity_class, roles_contradict,
 };
 
 #[test]
diff --git a/docs/research/standards-and-literature.md b/docs/research/standards-and-literature.md
index a24e49248..f1761730b 100644
--- a/docs/research/standards-and-literature.md
+++ b/docs/research/standards-and-literature.md
@@ -136,10 +136,20 @@ Aitchison, J., & Shen, S. M. (1980). Logistic-normal distributions: Some propert
 
 Aitchison, J. (1982). The statistical analysis of compositional data. *Journal of the Royal Statistical Society: Series B (Methodological), 44*(2), 139–177. https://doi.org/10.1111/j.2517-6161.1982.tb01195.x
 
+Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate: A practical and powerful approach to multiple testing. *Journal of the Royal Statistical Society: Series B (Methodological), 57*(1), 289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x
+
+Efron, B. (1979). Bootstrap methods: Another look at the jackknife. *The Annals of Statistics, 7*(1), 1–26. https://doi.org/10.1214/aos/1176344552
+
 Egozcue, J. J., Pawlowsky-Glahn, V., Mateu-Figueras, G., & Barceló-Vidal, C. (2003). Isometric logratio transformations for compositional data analysis. *Mathematical Geology, 35*(3), 279–300. https://doi.org/10.1023/A:1023818214614
 
+Fisher, R. A. (1921). On the "probable error" of a coefficient of correlation deduced from a small sample. *Metron, 1*(4), 3–32.
+
 Friedman, J., Hastie, T., & Tibshirani, R. (2008). Sparse inverse covariance estimation with the graphical lasso. *Biostatistics, 9*(3), 432–441. https://doi.org/10.1093/biostatistics/kxm045
 
+Hennig, C. (2007). Cluster-wise assessment of cluster stability. *Computational Statistics & Data Analysis, 52*(1), 258–281. https://doi.org/10.1016/j.csda.2006.11.025
+
+Monti, S. (2003). Consensus clustering: A resampling-based method for class discovery and visualization of gene expression microarray data. *Machine Learning, 52*(1–2), 91–118. https://doi.org/10.1023/A:1023949509487
+
 Traag, V. A., Waltman, L., & van Eck, N. J. (2019). From Louvain to Leiden: Guaranteeing well-connected communities. *Scientific Reports, 9*, Article 5233. https://doi.org/10.1038/s41598-019-41695-z
 
 Raw topic proportions should not be analyzed with ordinary Euclidean distances;
@@ -151,6 +161,17 @@ posterior and resampling uncertainty for every network edge and cluster. The
 Euclidean use of a raw simplex and scores recovered clusters with pair precision
 and recall.
 
+Posterior network edge admission uses exact two-sided Fisher z-transform
+p-values against `rho = 0` (Fisher, 1921) with the Benjamini–Hochberg
+step-up false-discovery-rate control (Benjamini & Hochberg, 1995);
+credible intervals and selection fractions are percentile bootstrap
+quantities over posterior draws (Efron, 1979). Consensus clusters come
+from resampled co-assignment matrices with an explicit, caller-supplied
+perturbation probability rather than an implicit constant (Monti, 2003;
+Hennig, 2007). No numeric threshold or weight in this path is chosen by
+heuristic; each is either an explicit parameter with stated provenance
+or a value derived from these primary sources.
+
 ## Time, events, and topic detection and tracking
 
 Allen, J. F. (1983). Maintaining knowledge about temporal intervals. *Communications of the ACM, 26*(11), 832–843. https://doi.org/10.1145/182.358434