---
title: "Benchmark Validation for Classic and Batched Stepped-Wedge Designs"
author: "Amanda Li"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Benchmark Validation for Classic and Batched Stepped-Wedge Designs}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  fig.width = 7,
  fig.height = 4.8
)

benchmark_dir <- system.file("benchmarks", package = "stepwedgepower")
if (!nzchar(benchmark_dir)) {
  candidates <- c(
    file.path("..", "inst", "benchmarks"),
    file.path("inst", "benchmarks")
  )
  benchmark_dir <- candidates[dir.exists(candidates)][1L]
}
benchmark_dir <- normalizePath(benchmark_dir, mustWork = TRUE)
old_benchmark_options <- options(stepwedgepower.benchmark_dir = benchmark_dir)

for (file in c(
  "benchmark_helpers.R",
  "table2_software_comparison.R",
  "figure3_calendar_delay.R",
  "table3_separate_time_grid.R",
  "figure5_limited_clusters.R"
)) {
  source(file.path(benchmark_dir, file), local = knitr::knit_global())
}

run_external <- tolower(Sys.getenv(
  "STEPWEDGEPOWER_RUN_EXTERNAL",
  unset = "false"
)) %in% c("1", "true", "yes", "on")

benchmark_output <- file.path(
  tempdir(), "stepwedgepower-benchmark-vignette"
)
```

## Purpose

This vignette reproduces four numerical artifacts from Lin et al. (2026),
*Power and Sample Size Implications of Delayed Cluster Initiation in Stepped
Wedge Trials*:

- the software comparison in Table 2;
- the calendar-time delay curves in Figure 3;
- the separate-time power grid in Table 3;
- the limited-cluster variance curves in Figure 5.

The code is stored under `inst/benchmarks/`, so it does not add functions to
the package namespace. Reference values, tolerances, provenance notes, and
documented source discrepancies are stored in machine-readable CSV files.

The benchmark separates manuscript reproduction from validation of the main
`stepwedgepower` engine. The manuscript calculations include continuous LMM
approximations and marginal GEE power, whereas `stepwedgepower` primarily uses
simulation with a binomial GLMM. Exact agreement is required only when the same
method and parameterization are being reproduced.

## Provenance and prespecified tolerances

```{r manifest}
manifest <- .sp_bench_read_reference(
  "benchmark_manifest.csv", benchmark_dir
)
knitr::kable(manifest)
```

```{r reference-environment}
reference_environment <- .sp_bench_read_reference(
  "reference_environment.csv", benchmark_dir
)
knitr::kable(reference_environment)
```

The source materials contain several inconsistencies that are retained in the
benchmark record rather than silently resolved.

```{r discrepancies}
source_notes <- .sp_bench_read_reference(
  "source_discrepancies.csv", benchmark_dir
)
knitr::kable(source_notes)
```

Most importantly, the prose for Figure 3 states \(p_1=0.3\), while the public R
code and plotted values use \(p_1=0.2\). The numerical reproduction below uses
\(p_1=0.2\).

## Table 2: software comparison

PASS and the Shiny CRT Calculator are reference-only because they are not
invoked by this R package. Calls to `swCRTdesign` and `swdpwr` are rerun only
when the optional packages are installed and the environment variable
`STEPWEDGEPOWER_RUN_EXTERNAL=true` is set.

```{r table2}
table2 <- run_table2_benchmark(
  output_dir = benchmark_output,
  benchmark_dir = benchmark_dir,
  run_external = run_external,
  strict = FALSE
)

table2_display <- table2$validation[, c(
  "scenario", "method", "reference_power", "reproduced_power",
  "absolute_difference", "status", "package_version", "pass"
)]
knitr::kable(table2_display, digits = 3)
```

When external execution is disabled, the reference values remain visible and
the reproducible rows are marked as skipped.

## Figure 3: power under delayed batch initiation

For a single batch with \(s\) steps and \(s+1\) periods, two copies of the
within-batch design are placed on a common calendar-time axis. The second batch
is shifted by \(d\) periods. The benchmark forms the calendar-time design
matrix, the exchangeable covariance matrix, and obtains the treatment variance
from

\[
  \left(Z^\mathsf{T}V^{-1}Z\right)^{-1}.
\]

The pooled binary variance and the normal-approximation power formula match the
manuscript computation. The maximal non-overlap delay is included; its power
equals the no-delay value under the stated equivalence.

```{r figure3-results}
figure3 <- run_figure3_benchmark(
  output_dir = benchmark_output,
  benchmark_dir = benchmark_dir,
  strict = TRUE,
  make_plot = FALSE
)

knitr::kable(
  figure3$results[, c("design", "delay", "variance", "power")],
  digits = 4
)
```

```{r figure3-plot}
.sp_bench_plot_figure3(figure3$results)
```

```{r figure3-endpoints}
knitr::kable(figure3$endpoint_check, digits = 12)
```

## Table 3: separate-time GEE power grid

Table 3 fixes the total sample size at 1,200 and evaluates four combinations of
steps, periods, and batches. The public analysis uses `swdpwr` with a marginal
GEE model, logit link, and a very small control-profile increment to include
time effects. Although the manuscript caption calls these simulated results,
the public code computes analytical GEE power.

```{r table3}
table3 <- run_table3_benchmark(
  output_dir = benchmark_output,
  benchmark_dir = benchmark_dir,
  run_external = run_external,
  strict = FALSE
)

table3_display <- table3$validation[, c(
  "design", "p0", "p1", "icc", "reference_power",
  "reproduced_power", "absolute_difference", "status",
  "package_version", "pass"
)]
knitr::kable(table3_display, digits = 3)
```

## Figure 5: limited-cluster variance curves

After division by the common multiplicative factor \(2a\), Appendix C gives

\[
\begin{aligned}
V_{\mathrm{classic}}(\phi)
  &= \frac{(1-\phi)(2\phi+1)}{\phi+1},\\
V_{\mathrm{delay\,1}}(\phi)
  &= \frac{1+\phi}{2},\\
V_{\mathrm{delay\,2}}(\phi)
  &= 1.
\end{aligned}
\]

```{r figure5-results}
figure5 <- run_figure5_benchmark(
  output_dir = benchmark_output,
  benchmark_dir = benchmark_dir,
  strict = TRUE,
  make_plot = FALSE
)

knitr::kable(figure5$crossing_check, digits = 12)
```

The closed forms are also checked against an independent WLS matrix-inversion
calculation for every grid value with \(0 \leq \phi < 1\). The boundary
\(\phi=1\) is retained for the analytical curves but excluded from matrix
validation because the covariance matrix is singular there.

```{r figure5-matrix-check}
figure5_matrix_differences <- figure5$validation[
  figure5$validation$matrix_check_applicable,
  c(
    "classic_matrix_absolute_difference",
    "one_unit_delay_matrix_absolute_difference",
    "two_unit_delay_matrix_absolute_difference"
  ),
  drop = FALSE
]
knitr::kable(
  data.frame(
    check = names(figure5_matrix_differences),
    maximum_absolute_difference = vapply(
      figure5_matrix_differences,
      max,
      numeric(1),
      na.rm = TRUE
    )
  ),
  digits = 12
)
```

The classic and one-unit-delay curves cross at
\(\phi=1/\sqrt{5}\), approximately 0.4472.

```{r figure5-plot}
.sp_bench_plot_figure5(figure5$results)
```

## Running the complete suite

The installed benchmark directory contains a single runner:

```r
benchmark_dir <- system.file(
  "benchmarks",
  package = "stepwedgepower"
)

Sys.setenv(
  STEPWEDGEPOWER_RUN_EXTERNAL = "true",
  STEPWEDGEPOWER_BENCHMARK_STRICT = "true",
  STEPWEDGEPOWER_BENCHMARK_OUTPUT =
    file.path(getwd(), "stepwedgepower-benchmark-results")
)

source(file.path(benchmark_dir, "run_all.R"))
```

The runner writes CSV files, PDF figures, an RDS object, package-version
metadata, and `sessionInfo()`. External package versions are recorded because
software updates may alter rounded benchmark values.

## Relationship to simulation-based validation

These manuscript benchmarks answer whether the numerical settings and
analytical calculations have been reproduced. They do not imply that
conditional GLMM power, marginal GEE power, and linear-approximation power
should be identical. The package function `validate_engine()` separately
cross-checks the simulation engine against an independently written simulator
using the same data-generating and fitted models.

## References

Lin T, Bailey BE, Staras S, Bylund C, Lee J-H (2026). *Power and Sample Size
Implications of Delayed Cluster Initiation in Stepped Wedge Trials*.
Manuscript dated April 20, 2026.

Companion code: <https://github.com/tuolin123/BatchedSW>.

Kasza J, Bowden R, Hooper R, Forbes AB (2022). The batched stepped wedge
design: a design robust to delays in cluster recruitment. *Statistics in
Medicine*, 41, 3627--3641.

```{r cleanup, include=FALSE}
options(old_benchmark_options)
```