Study design, inference and reproducibility
A reproducible experiment needs a declared endpoint, a comparison design, a replication schedule and an explicit multiplicity family. This tutorial shows those decisions using synthetic logistics items and dynamic network challenges. It does not provide calibrated application-specific assumptions.
Define the experiment before looking at outcomes
The unit of analysis is an independently simulated batch. Items within a batch share queues, vehicles, observations and disruption risk history, so they are not independent statistical replicates. Use batch-level outcomes for the functions below.
master_seed <- 20260512L
scenario <- example_scenario("constrained", N = 6L)
policies <- c("P0", "P1", "P3")
arm_seeds <- setNames(master_seed + seq_along(policies) - 1L, policies)
n_rep <- 6L
bootstrap_resamples <- 10000L
endpoint <- "viable_stage2"
comparison_family <- "teaching-policy-comparisons-at-Stage2"
data.frame(policy = policies, seed = unname(arm_seeds), n_rep = n_rep)
#> policy seed n_rep
#> 1 P0 20260512 6
#> 2 P1 20260513 6
#> 3 P3 20260514 6The master seed is 20260512. Different arms receive documented, distinct seeds derived from it. Six replications keep this vignette small; choose the study replication count using a prospective Monte Carlo precision assessment. The bootstrap uses the full reporting default of 10,000 resamples. A smaller bootstrap count may be useful for an interactive preview but must be identified as such and replaced for a final analysis.
The primary logistics endpoint is viable delivery, operationalized here by arrival viable at Stage 2 by the administrative horizon. This is a proportion of the original batch, not a Kaplan-Meier viability probability at an arbitrary time. A item still viable elsewhere at the horizon is an endpoint failure. Arrival is recorded before unloading and ends item-condition hazard; transport can continue to Stage 3. Downstream disruptions can still affect the Stage 3 endpoint.
Run independent policy arms
arms <- setNames(lapply(policies, function(policy) {
simulate_shield(scenario, policy = policy, disruption_model = "static",
n_rep = n_rep, seed = arm_seeds[[policy]])
}), policies)
summary_P0 <- summarize_shield(arms$P0)
summary_P0[summary_P0$endpoint %in%
c("viable_stage2", "t_stage2_med", "util_L4", "n_pending"), ]
#> endpoint n missing median q1 q3
#> 2 viable_stage2 6 0 0.3333333 0.3333333 0.3333333
#> 4 t_stage2_med 6 0 11.0000000 11.0000000 11.0000000
#> 9 util_L4 6 0 0.2500000 0.2500000 0.2500000
#> 14 n_pending 6 0 6.0000000 6.0000000 6.0000000Each row in an arm is one replication; each row above describes an endpoint across those replications. The medians and quartiles do not describe the parameter distribution or uncertainty about the model’s scientific validity.
plot(stats::ecdf(arms$P0[[endpoint]]), verticals = TRUE, do.points = FALSE,
xlim = c(0, 1), ylim = c(0, 1),
xlab = "Proportion arriving viable at Stage 2",
ylab = "Cumulative proportion of batch replications", main = "")
lines(stats::ecdf(arms$P1[[endpoint]]), verticals = TRUE, do.points = FALSE,
lty = 2)
lines(stats::ecdf(arms$P3[[endpoint]]), verticals = TRUE, do.points = FALSE,
lty = 3)
legend("bottomright", legend = policies, lty = 1:3, bty = "n")
Neutral teaching run: batch-level Stage 2 arrival-viable proportions by policy.
The neutral scenario has zero exposure and item-condition hazards. Differences in arrival outcomes can reflect congestion or the administrative horizon. The plot is a teaching output and establishes no policy advantage in practice.
Report nonparametric comparisons in full
compare_policies() applies a two-sided, tie-corrected
Mann-Whitney U approximation to independent arms. The test compares
distributions; it is not specifically a test of medians. It also returns
Cliff’s delta and the difference in medians with percentile bootstrap
95% intervals. This is a primary nonparametric framework, with no
preliminary normality test.
comparisons <- do.call(rbind, lapply(c("P1", "P3"), function(policy) {
compare_policies(
arms[[policy]][[endpoint]], arms$P0[[endpoint]],
label_x = policy, label_y = "P0", family = comparison_family,
B = bootstrap_resamples, seed = master_seed
)
}))
comparisons[c("arm_x", "arm_y", "n_x", "n_y", "U", "p_raw_log10",
"cliff_delta", "cliff_lwr", "cliff_upr",
"med_diff", "med_diff_lwr", "med_diff_upr",
"med_diff_ci_degenerate")]
#> arm_x arm_y n_x n_y U p_raw_log10 cliff_delta cliff_lwr cliff_upr
#> 1 P1 P0 6 6 18 0.000000 0.0000000 0 0.0000000
#> 2 P3 P0 6 6 6 -1.608912 -0.6666667 -1 -0.3333333
#> med_diff med_diff_lwr med_diff_upr med_diff_ci_degenerate
#> 1 0.0000000 0.0000000 0 TRUE
#> 2 -0.1666667 -0.1666667 0 FALSEEffects are oriented as x minus y. Positive
values favor x for viability; for delay or resource use,
whether larger values are desirable depends on the question. Do not
infer direction from the policy labels alone. Cliff’s delta contrasts
the probabilities of a randomly drawn x exceeding or
falling below a randomly drawn y, accounting for ties.
Retain sample sizes, median/IQR descriptives, the U statistic, raw and adjusted p-values, effect estimates and intervals in the analysis output. The bootstrap intervals quantify Monte Carlo sampling uncertainty conditional on the model. They do not include parameter uncertainty, structural uncertainty, measurement validity or uncertainty from choosing a model after seeing the result.
A bootstrap interval can collapse on this discrete endpoint. Inspect
the med_diff_ci_degenerate flag and the interval endpoints;
a point interval is not evidence of empirical certainty.
Declare the complete FDR family
The family below contains both planned policy contrasts at the same endpoint. Adjust its complete set of p-values together. If a study defines additional endpoints or scenarios in the same family, include them before adjustment. Selecting only significant or interesting contrasts changes the family.
comparisons$p_bh_log10 <- adjust_fdr(comparisons$p_raw_log10)
comparisons[c("family", "arm_x", "arm_y", "p_raw_log10", "p_bh_log10")]
#> family arm_x arm_y p_raw_log10 p_bh_log10
#> 1 teaching-policy-comparisons-at-Stage2 P1 P0 0.000000 0.000000
#> 2 teaching-policy-comparisons-at-Stage2 P3 P0 -1.608912 -1.307882The helpers retain base-10 log p-values so extreme tails remain
distinguishable. p_raw_log10 is authoritative when a linear
p-value underflows; use log probabilities or a properly formatted bound
instead of reporting a computed probability as zero.
adjust_fdr() implements Benjamini-Hochberg adjustment for
one supplied family and leaves missing entries missing. It does not
infer a hypothesis grid or silently pool separate families.
More than two arms
For a prespecified omnibus comparison, use Kruskal-Wallis and retain its H statistic and degrees of freedom. Any decision to gate pairwise contrasts on that omnibus result belongs to the declared study analysis plan.
values <- unlist(lapply(arms, `[[`, endpoint), use.names = FALSE)
groups <- rep(names(arms), vapply(arms, nrow, integer(1)))
compare_multiple(values, groups, family = "teaching-omnibus-at-Stage2")
#> family test arms H df p_raw
#> 1 teaching-omnibus-at-Stage2 Kruskal-Wallis P0|P1|P3 9.714286 2 0.00777266
#> p_raw_log10 p_underflow
#> 1 -2.10943 FALSEThe omnibus helper reports the test statistic and p-value. Pairwise effect sizes and their intervals come from explicitly chosen contrasts; the function does not select post-hoc comparisons or manufacture a multi-arm effect size. For association questions, the study’s primary nonparametric choice is Spearman; specify the variables, independent unit and any bootstrap interval separately.
Paired disruptions-disabled intervention
A paired intervention is a different design. It requires identical parameters, policy, disruption-response family, replication count and seed schedule, with disruption enablement as the only changed setting. Item draws are indexed by item and minute and use a stream separate from policy and disruption draws.
with_disruptions <- simulate_shield(
scenario, policy = "P3", disruption_model = "adaptive", n_rep = n_rep,
seed = master_seed, disruptions_enabled = TRUE
)
without_disruptions <- simulate_shield(
scenario, policy = "P3", disruption_model = "adaptive", n_rep = n_rep,
seed = master_seed, disruptions_enabled = FALSE
)
compare_intervention(with_disruptions, without_disruptions,
B = bootstrap_resamples, seed = master_seed)
#> endpoint n_pairs contrast mean_difference lower
#> 1 viable_stage2 6 no_disruptions minus disruptions 0 0
#> upper point_interval
#> 1 0 TRUEThe contrast is no disruptions minus disruptions, estimated as a paired mean difference of batch proportions with a paired percentile bootstrap interval. Pair indices are resampled jointly. The helper checks matching identities and aligns rows by replication identifier.
The neutral exposure is zero, so the intervention has no effect here.
In a user-specified disruption risk experiment, this is a total modelled
intervention effect, including downstream decisions and
administrative-horizon effects. It is not identical to
n_disruption_loss or a count of direct disruption losses.
Matched arms must not be passed to the independent-arm Mann-Whitney U
helper.
Preserve missing outcomes and replay the experiment
Arrival-time medians are conditional on arrival. If an arm contains
an undefined median, compare_policies() rejects that
nonfinite input. Decide prospectively whether a comparison of defined
batch medians is meaningful, report the missing counts by arm, and keep
the viability/arrival endpoint alongside it. Replacing a missing time by
zero or the horizon changes the estimand.
Simulation and bootstrap calls restore the caller’s RNG state and kind, while recording their own reproducible configuration. Repeating a run with the same parameters, software and seeds reproduces its batch results:
state_before <- .Random.seed
replayed <- simulate_shield(scenario, policy = "P0", disruption_model = "static",
n_rep = n_rep, seed = arm_seeds[["P0"]])
stopifnot(identical(arms$P0, replayed),
identical(state_before, .Random.seed))Replay one batch or distribute a fixed schedule
Use replication_seeds() to construct the complete cell
schedule once, then send individual seeds to
simulate_replication() from your analysis scheduler. The
latter takes a replication seed directly; it does not sample a new cell
schedule or attach study identifiers. Checkpoint naming and workload
allocation remain the analysis project’s responsibility.
schedule <- replication_seeds(n_rep = n_rep, seed = arm_seeds[["P0"]])
stopifnot(identical(schedule, arms$P0$replicate_seed))
first_batch <- simulate_replication(
scenario, policy = "P0", disruption_model = "static",
replication_seed = schedule[1]
)
stopifnot(all(vapply(names(first_batch), function(column) {
identical(first_batch[[column]], arms$P0[[column]][1])
}, logical(1))))Generate the declared full schedule before splitting into chunks. Restarting a master seed independently for each chunk repeats the beginning of the schedule and creates duplicate replications. For an disruptions-disabled pair, use the same schedule in both configurations. Preserve the total replication count and the mapping from each seed to its replication identifier in durable checkpoints.
For analysis pipelines that already manage their own data frames,
mann_whitney_u() exposes the tested U calculation and log
tail probability, bootstrap_statistic() exposes named
independent or paired interval estimators, and
effect_magnitude() provides descriptive Cliff’s delta
labels. These primitives do not establish independence or pairing for
the caller; prefer the higher-level comparison functions when using
complete run objects.
Save the complete request
Save the request, full result object, analysis settings and session information. An RDS retains attributes required for pairing, whereas CSV retains only table columns.
request <- list(parameters = scenario, policy = "P0", disruption_model = "static",
n_rep = n_rep, seed = arm_seeds[["P0"]])
bundle <- list(request = request, run = arms$P0,
analysis = list(endpoint = endpoint,
family = comparison_family,
bootstrap_resamples = bootstrap_resamples,
master_seed = master_seed),
session_info = sessionInfo())
bundle_file <- tempfile("shieldR-study-", fileext = ".rds")
saveRDS(bundle, bundle_file)
restored <- readRDS(bundle_file)
stopifnot(identical(do.call(simulate_shield, restored$request), restored$run))For a real analysis, use a named output directory and retain an immutable copy of the exact source archive or package commit, the scenario input hash, seed schedule, replication counts and an environment lockfile. A version number alone cannot distinguish two uncommitted source states with the same version. Study scheduling, frozen input/output manifests and publication exports remain in the analysis project. Method changes require explicit validation before adoption by a study.
sessionInfo()
#> R version 4.6.2 (2026-10-29)
#> 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] shieldR_0.6.0 rmarkdown_2.32
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 jsonlite_2.0.0 compiler_4.6.2 promises_1.5.0
#> [5] Rcpp_1.1.2 parallel_4.6.2 later_1.4.8 jquerylib_0.1.4
#> [9] scales_1.4.0 yaml_2.3.12 fastmap_1.2.0 mime_0.13
#> [13] ggplot2_4.0.3 R6_2.6.1 knitr_1.52 htmlwidgets_1.6.4
#> [17] maketools_1.3.2 shiny_1.14.0 RColorBrewer_1.1-3 bslib_0.12.0
#> [21] rlang_1.3.0 DT_0.34.0 cachem_1.1.0 httpuv_1.6.17
#> [25] xfun_0.62 S7_0.2.2 sass_0.4.10 sys_3.4.3
#> [29] otel_0.2.0 memoise_2.0.1 cli_3.6.6 magrittr_2.0.5
#> [33] crosstalk_1.2.2 digest_0.6.39 grid_4.6.2 fontawesome_0.5.3
#> [37] xtable_1.8-8 lifecycle_1.0.5 vctrs_0.7.3 evaluate_1.0.5
#> [41] glue_1.8.1 farver_2.1.2 buildtools_1.0.0 tools_4.6.2
#> [45] htmltools_0.5.9