The vignettes so far audit evidence that is already in hand. This one covers the other half of the package: ranking tissue-to-blood candidates from stage-resolved effect tables, and then committing to a decision rule before the destination cohort is opened.
The motivation is a negative result. Two PDAC serum audits were
prespecified, hashed, externally timestamped, and opened exactly once.
Both terminated without producing a performance statistic: one stopped
at a label-blind technical QC gate because too few frozen candidates
were quantifiable in the public serum deposit, and the other lock of
that dataset and the second dataset stopped because the deposited
outcome labels did not match what had been frozen. The record is in
audit/PROVENANCE.md in the source repository. That is why
the runner below has a terminal INCONCLUSIVE state that is
deliberately not a failure of the marker.
Everything in this vignette runs on simulated data, so no download is needed and the numbers are reproducible from the seeds shown.
sb_simulate() generates a stage-resolved effect table, a
blood detection table, and the feature-level truth used only for
evaluation. The translation scenario plants coherent
tissue-only decoys and inflammation-associated blood decoys, so a method
that merely rewards strong tissue signal can be caught out.
sim <- sb_simulate(n_features = 400L, scenario = "translation", seed = 1L)
vapply(sim, nrow, integer(1))
#> effects detection truth
#> 8141 1600 400sb_validate_effects() is the input contract. It
normalises the five required columns, drops non-finite estimates, and
refuses duplicated feature/cohort/modality/contrast rows rather than
silently averaging them.
effects <- sb_validate_effects(sim$effects)
str(effects[, c("feature", "cohort", "modality", "contrast", "estimate")])
#> 'data.frame': 8141 obs. of 5 variables:
#> $ feature : chr "G00002" "G00003" "G00004" "G00006" ...
#> $ cohort : chr "R1" "R1" "R1" "R1" ...
#> $ modality: chr "rna" "rna" "rna" "rna" ...
#> $ contrast: chr "early_vs_reference" "early_vs_reference" "early_vs_reference" "early_vs_reference" ...
#> $ estimate: num 0.109 0.3686 0.2381 0.0237 -0.3891 ...
unique(effects$modality)
#> [1] "rna" "protein" "plasma"A detection in a shallow cohort is more informative than the same detection in a cohort where nearly everything is seen, so raw detection frequency is not used directly.
stagebridge() combines eight components into a weighted
geometric score. The geometric form matters: a candidate with one
near-zero component cannot be rescued by the other seven, which is the
intended behaviour for a translation claim.
target_relation must be declared before
looking at the destination data. "concordant" asks for
tissue and blood to move the same way; "inverse" asks for
the opposite; "agnostic" scores magnitude only.
ranked <- stagebridge(
sim$effects,
detection = sim$detection,
target_relation = "concordant"
)
head(ranked[, c(
"feature", "stagebridge_score", "score_rank",
"pareto_front", "pareto_tier"
)], 5)
#> feature stagebridge_score score_rank pareto_front pareto_tier
#> 1 G00299 0.7940782 1 1 A
#> 2 G00290 0.7656246 2 1 A
#> 3 G00330 0.7535901 3 1 A
#> 4 G00289 0.7477191 4 1 A
#> 5 G00329 0.6742572 6 1 ATwo tiers are reported and they are deliberately separate.
pareto_front is exact non-dominance over the eight
components; score_percentile_tier is a score-only cut. A
candidate can lead on score while being dominated, and that disagreement
is information.
loco <- sb_loco_stability(sim$effects, sim$detection, top_k = 20L)
head(loco, 3)
#> feature median_loco_rank worst_loco_rank loco_rank_iqr top_k_frequency n_loco
#> 1 G00299 1 5 0.00 1 8
#> 2 G00290 2 5 0.25 1 8
#> 3 G00330 3 4 0.25 1 8worst_loco_rank is the number to read. A candidate whose
median rank is 2 but whose worst leave-one-cohort-out rank is 150 is
carried by a single cohort.
The package ships the comparators rather than only its own score, so that claims of improvement can be checked instead of asserted.
scores <- sb_method_scores(sim$effects, sim$detection)
truth <- sim$truth$truth[match(scores$feature, sim$truth$feature)]
methods <- setdiff(names(scores), "feature")
evaluation <- data.frame(
method = methods,
auroc = vapply(methods, function(m) {
sb_binary_auc(truth, scores[[m]])
}, numeric(1)),
auprc = vapply(methods, function(m) {
sb_average_precision(truth, scores[[m]])
}, numeric(1)),
recall_top20 = vapply(methods, function(m) {
sb_topk_recall(truth, scores[[m]], k = 20L)
}, numeric(1)),
row.names = NULL
)
evaluation[order(-evaluation$auroc), ]
#> method auroc auprc recall_top20
#> 3 signed_stouffer 0.9742908 0.6715781 0.5000000
#> 5 rank_product 0.9678635 0.6242474 0.5000000
#> 4 intersection 0.9676418 0.6115087 0.4583333
#> 2 vote_count 0.9675310 0.6057679 0.4166667
#> 1 mean_abs_z 0.9629876 0.5864892 0.3750000
#> 6 svd_loading 0.9621011 0.5809682 0.4166667
#> 7 observed_rank_product 0.9580009 0.8005503 0.7083333
#> 8 StageBridge 0.9132314 0.7795034 0.6666667On this scenario StageBridge has the lowest AUROC of the
eight methods while ranking second on AUPRC and on top-20
recall. That is the expected shape of the trade-off, not a bug: the
score is tuned to put translatable candidates at the very top of a
shortlist, and it pays for that with global ranking quality.
observed_rank_product beats it on all three summaries here.
Read this as a demonstration of the evaluation harness, not as a claim
that one method dominates – a single simulated scenario cannot settle
that, and the scenario argument exists precisely so that assumptions can
be violated on purpose.
sb_multiblock_plsda_scores() is the matched-sample
comparator. Unlike everything above it needs the same observations, in
the same order, in every block, so it does not apply to unpaired
cohort-level effect tables. It is shown here for completeness.
outcome <- factor(rep(c("Normal", "Tumor"), each = 12))
dn <- list(paste0("S", seq_len(24)), paste0("F", seq_len(8)))
rna <- matrix(stats::rnorm(24 * 8), nrow = 24, dimnames = dn)
protein <- rna + stats::rnorm(24 * 8, sd = 0.5)
head(sb_multiblock_plsda_scores(list(RNA = rna, Protein = protein), outcome), 3)
#> feature MBPLSDA MBPLSDA_RNA MBPLSDA_Protein
#> F1 F1 1.00000 1.000 1.000
#> F5 F5 0.73951 0.875 0.625
#> F6 F6 0.73951 0.625 0.875Everything above is exploratory and re-runnable.
sb_one_shot() is not. It takes a frozen panel, a frozen
direction per feature, and a rule fixed in advance, then returns exactly
one of PASS, FAIL or
INCONCLUSIVE.
str(sb_default_one_shot_rule())
#> List of 13
#> $ minimum_broad_candidate_coverage : int 200
#> $ minimum_top10_coverage : int 5
#> $ minimum_top20_coverage : int 10
#> $ minimum_per_group_detections_for_panel_gene : int 14
#> $ panel_size : int 10
#> $ permutation_replicates : int 2000
#> $ permutation_seed : int 20260718
#> $ bootstrap_replicates : int 1000
#> $ bootstrap_seed : int 20260719
#> $ minimum_pass_auc : num 0.7
#> $ maximum_pass_permutation_p : num 0.05
#> $ minimum_direction_consistent_nominal_feature_replications: int 2
#> $ feature_replication_nominal_p : num 0.05The gate stack runs first: group ontology and counts, broad candidate coverage, top-10 and top-20 panel coverage, and complete-panel quantifiability. The statistical endpoint runs only if every gate passes.
make_cohort <- function(delta, n_covered = 320L, seed = 20260719L) {
frozen <- data.frame(
feature = sprintf("G%03d", 1:320),
frozen_direction = rep(c(1L, -1L), length.out = 320L)
)
panel <- data.frame(
feature = frozen$feature[1:20],
frozen_direction = frozen$frozen_direction[1:20],
rank_weight = 1 / log2(2:21)
)
samples <- c(
sprintf("CP_%02d", 1:20), sprintf("CN_%02d", 1:20),
sprintf("HC_%02d", 1:20)
)
grid <- expand.grid(
gene = frozen$feature[seq_len(n_covered)], sample = samples,
KEEP.OUT.ATTRS = FALSE, stringsAsFactors = FALSE
)
grid$group <- sub("_.*$", "", grid$sample)
grid$log2_abundance <- withr::with_seed(
seed, stats::rnorm(nrow(grid), 8, 1)
)
if (delta != 0) {
hit <- grid$group == "CP" & grid$gene %in% panel$feature[1:10]
direction <- panel$frozen_direction[match(grid$gene, panel$feature)]
grid$log2_abundance[hit] <-
grid$log2_abundance[hit] + delta * direction[hit]
}
list(
abundance = grid, frozen = frozen, panel = panel,
groups = list(
case = "CP", primary_reference = "CN", secondary_reference = "HC"
),
counts = c(CP = 20L, CN = 20L, HC = 20L)
)
}
# Fewer resamples than the locked rule, to keep the vignette quick.
quick_rule <- function() {
rule <- sb_default_one_shot_rule()
rule$permutation_replicates <- 200L
rule$bootstrap_replicates <- 200L
rule
}
audit <- function(...) {
cohort <- make_cohort(...)
sb_one_shot(
cohort$abundance, cohort$frozen, cohort$panel,
cohort$groups, cohort$counts, quick_rule()
)
}A cohort with a real planted signal reaches PASS:
passed <- audit(delta = 1.2)
c(status = passed$status, auc = round(passed$primary$auc, 3))
#> status auc
#> "PASS" "0.99"The same cohort with the signal removed reaches FAIL.
The endpoint ran and did not meet the criteria – this is a statement
about the biology:
failed <- audit(delta = 0)
c(status = failed$status, runner_entered = failed$statistical_runner_entered)
#> status runner_entered
#> "FAIL" "TRUE"Now a coverage failure, the case behind the first real PDAC lock. The
signal is still present, but only 50 of the 320 frozen candidates are
quantifiable, so the coverage gate halts execution before the endpoint.
The result is INCONCLUSIVE: the test never ran, and nothing
has been learned about the marker:
halted <- audit(delta = 1.2, n_covered = 50L)
c(status = halted$status, runner_entered = halted$statistical_runner_entered)
#> status runner_entered
#> "INCONCLUSIVE" "FALSE"
halted$failed_gates
#> [1] "broad_coverage"Collapsing that third outcome into FAIL would have
reported a biological conclusion that the data cannot support. Keeping
it separate is the point of the design.
A gate stack strict enough never to pass anything would produce the
same INCONCLUSIVE results on real data while proving
nothing. The package’s test suite therefore drives the runner to all
three terminal states and checks that AUROC rises monotonically with the
size of the planted signal
(tests/testthat/test-sb-one-shot.R). PASS is
reachable, so the real-data terminations are attributable to assay
coverage and to outcome-label mismatches rather than to an
over-conservative rule.
The two recorded PDAC audits, their accessions and their terminal
states are documented in audit/PROVENANCE.md. The lock
manifests hash paths relative to the root of the archived working tree
and are covered by external RFC 3161 timestamps, so they verify only
inside that archive; the copies under audit/locks/ are a
disclosure record. Regenerating a manifest would destroy the property it
exists to establish.
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] data.table_1.18.6.1 pdactrace_0.99.25 rmarkdown_2.32
#>
#> loaded via a namespace (and not attached):
#> [1] SummarizedExperiment_1.43.0 gtable_0.3.6
#> [3] ellipse_0.5.0 xfun_0.61
#> [5] bslib_0.12.0 ggplot2_4.0.3
#> [7] ggrepel_0.9.8 Biobase_2.73.2
#> [9] lattice_0.23-1 vctrs_0.7.3
#> [11] tools_4.6.1 generics_0.1.4
#> [13] stats4_4.6.1 parallel_4.6.1
#> [15] tibble_3.3.1 rARPACK_0.11-0
#> [17] pkgconfig_2.0.3 Matrix_1.7-6
#> [19] RColorBrewer_1.1-3 S7_0.2.2
#> [21] S4Vectors_0.51.10 mixOmics_6.37.0
#> [23] lifecycle_1.0.5 stringr_1.6.0
#> [25] compiler_4.6.1 farver_2.1.2
#> [27] statmod_1.5.2 Seqinfo_1.3.2
#> [29] codetools_0.2-20 htmltools_0.5.9
#> [31] sys_3.4.3 buildtools_1.0.0
#> [33] sass_0.4.10 yaml_2.3.12
#> [35] tidyr_1.3.2 pillar_1.11.1
#> [37] jquerylib_0.1.4 MASS_7.3-66
#> [39] BiocParallel_1.47.0 limma_3.99.0
#> [41] DelayedArray_0.39.7 cachem_1.1.0
#> [43] abind_1.4-8 RSpectra_0.16-2
#> [45] tidyselect_1.2.1 digest_0.6.39
#> [47] stringi_1.8.9 purrr_1.2.2
#> [49] reshape2_1.4.5 dplyr_1.2.1
#> [51] maketools_1.3.2 labeling_0.4.3
#> [53] fastmap_1.2.0 grid_4.6.1
#> [55] cli_3.6.6 SparseArray_1.13.3
#> [57] magrittr_2.0.5 S4Arrays_1.13.1
#> [59] corpcor_1.6.10 withr_3.0.3
#> [61] scales_1.4.0 XVector_0.53.0
#> [63] matrixStats_1.5.0 igraph_2.3.3
#> [65] otel_0.2.0 gridExtra_2.3.1
#> [67] evaluate_1.0.5 knitr_1.52
#> [69] GenomicRanges_1.65.4 IRanges_2.47.5
#> [71] rlang_1.3.0 Rcpp_1.1.2
#> [73] glue_1.8.1 BiocGenerics_0.59.12
#> [75] jsonlite_2.0.0 plyr_1.8.9
#> [77] R6_2.6.1 MatrixGenerics_1.25.0