Motivating example
A hospital is planning to deploy a clinical decision support
system (CDSS) intended to reduce its case fatality rate (CFR).
The research team wants to answer, prospectively:
“How many months of post-implementation follow-up do we need to
reliably detect a meaningful reduction in CFR?”
Step 1 — Load pre-intervention data
data("example_cfr_data")
head(example_cfr_data)
#> time outcome
#> 1 1 12.79
#> 2 2 14.58
#> 3 3 16.44
#> 4 4 18.41
#> 5 5 14.47
#> 6 6 13.21
plot(example_cfr_data$time, example_cfr_data$outcome,
type = "o", pch = 16, col = "steelblue", lwd = 1.2,
xlab = "Month", ylab = "CFR (%)",
main = "Pre-intervention case fatality rate",
las = 1, bty = "l")
abline(lm(outcome ~ time, data = example_cfr_data),
col = "firebrick", lwd = 2, lty = 2)
legend("topright", legend = c("Observed", "Linear trend"),
col = c("steelblue", "firebrick"), lty = 1,
pch = c(16, NA), bty = "n", cex = 0.85)
grid(NULL, NULL, lwd = 0.6, col = rgb(0, 0, 0, 0.1))
Monthly CFR over 24 pre-intervention
months.
Step 2 — Estimate nuisance parameters
params <- estimate_its_params(example_cfr_data)
#>
#> ============================================================
#> PITS - ITS Parameter Estimation
#> ============================================================
#>
#> Estimation method: GLS-ML
#> Observations (n_pre): 24
#> Baseline: 14.5853
#> Sigma (residual SD): 1.4876 (10.2% of baseline)
#> Rho (AR1): 0.3832
#> Pre-trend (per unit): 0.0039
#>
#> ------------------------------------------------------------
#> Copy these into calculate_power():
#>
#> n_pre = 24
#> baseline = 14.5853
#> sigma = 1.4876
#> rho = 0.3832
#> level_change = ??? # YOUR clinical hypothesis
#> slope_change = 0 # set if testing slope
#>
#> ============================================================
{r diag-plots, fig.width = 8, fig.height = 6, fig.cap = "Diagnostic plots for pre-intervention residuals."} diagnose_params(example_cfr_data)
The residual ACF confirms moderate positive autocorrelation (rho ≈
0.37), typical for monthly hospital data. Ignoring this would overstate
power.
Step 4 — Calculate power for a single design
Note on simulation size. To keep this vignette fast
to build, the examples below use a small number of Monte Carlo
replications (n_sim = 100–200). For a reportable analysis,
use n_sim = 1000 or more for stable power estimates.
result_24 <- calculate_power(
n_pre = params$n_pre, n_post = 24,
baseline = params$baseline, level_change = level_change,
sigma = params$sigma, rho = params$rho,
test = "level", n_sim = 200, seed = 123
)
print(result_24)
#>
#> ============================================================
#> PITS - ITS Power Estimate
#> ============================================================
#>
#> Design: Single site
#> Test: level
#> n_pre: 24
#> n_post: 24
#> Baseline: 14.585
#> Level change: -3.000
#> Sigma: 1.488
#> Rho: 0.383
#> Alpha: 0.050
#> N simulations: 200
#> N converged: 200 (100.0%)
#>
#> ------------------------------------------------------------
#> POWER = 77.5% - Borderline (60-79%)
#> ------------------------------------------------------------
24 months is borderline. The sweep below identifies the minimum
adequate duration.
Step 5 — Optimise with a design sweep
sweep <- power_sweep(
sweep_post = c(12, 18, 24, 30, 36, 42, 48),
n_pre = params$n_pre,
baseline = params$baseline,
level_change = level_change,
sigma = params$sigma,
rho = params$rho,
test = "level",
n_sim = 200, seed = 123
)
#>
#> =======================================================
#> PITS - Design optimisation sweep
#> =======================================================
#> n_post Power Interpretation
#> -------------------------------------------------------
#> 12 70.5% Borderline (60-79%)
#> 18 73.5% Borderline (60-79%)
#> 24 77.5% Borderline (60-79%)
#> 30 80.5% Adequate (>= 80%) <-- adequate
#> 36 87.5% Adequate (>= 80%) <-- adequate
#> 42 87.0% Adequate (>= 80%) <-- adequate
#> 48 91.5% Adequate (>= 80%) <-- adequate
#> =======================================================
plot_power_curve(
sweep,
main = "Power vs follow-up duration — CDSS/CFR example",
xlab = "Post-implementation months (n_post)",
ylab = "Estimated power (%)"
)
Power curve. Dashed red = 80% target; green =
minimum adequate n_post.
Step 6 — Plot a simulated ITS realisation
plot_its_example(
n_pre = params$n_pre, n_post = 30,
baseline = params$baseline, level_change = level_change,
sigma = params$sigma, rho = params$rho,
seed = 7,
xlab = "Month", ylab = "CFR (%)",
main = "Simulated ITS — CDSS effect on case fatality rate"
)
Simulated ITS under the alternative hypothesis
with fitted segmented regression.
Step 7 — Sensitivity analysis
How sensitive is the required follow-up to uncertainty in
sigma and rho?
grid <- build_param_grid(
n_post = c(18, 24, 30, 36),
level_change = level_change,
sigma = c(1.0, round(params$sigma, 1), 2.0),
rho = c(0.2, 0.4, 0.6),
n_pre = params$n_pre,
baseline = params$baseline
)
grid_results <- run_power_grid(grid, n_sim = 100, seed = 42, verbose = FALSE)
{r heatmap-sigma, fig.height = 5, fig.cap = "Power by n_post and sigma (rho = 0.4). Dashed line = first adequate column."} plot_power_heatmap( grid_results[grid_results$rho == 0.4, ], x = "n_post", y = "sigma", main = "Power heatmap: n_post × sigma (rho = 0.4)", xlab = "Post-implementation months", ylab = "Residual SD (sigma)" )
{r heatmap-rho, fig.height = 5, fig.cap = "Power by n_post and rho (sigma fixed at estimated value)."} sigma_val <- round(params$sigma, 1) plot_power_heatmap( grid_results[round(grid_results$sigma, 1) == sigma_val, ], x = "n_post", y = "rho", main = sprintf("Power heatmap: n_post × rho (sigma = %.1f)", sigma_val), xlab = "Post-implementation months", ylab = "Autocorrelation (rho)" )
Step 8 — Multi-site analysis
If the CDSS is deployed simultaneously across three hospitals,
pooling evidence increases power substantially.
sites <- list(
list(name = "Hospital A", n_pre = 24, n_post = 24,
baseline = params$baseline, level_change = level_change,
slope_change = 0, sigma = params$sigma, rho = params$rho),
list(name = "Hospital B", n_pre = 24, n_post = 24,
baseline = 17.5, level_change = level_change,
slope_change = 0, sigma = 2.0, rho = 0.35),
list(name = "Hospital C", n_pre = 24, n_post = 24,
baseline = 13.0, level_change = level_change,
slope_change = 0, sigma = 1.8, rho = 0.42)
)
result_multi <- calculate_power_multi(
sites, test = "level", n_sim = 200, seed = 123
)
print(result_multi)
#>
#> ============================================================
#> PITS - ITS Power Estimate
#> ============================================================
#>
#> Design: Multi-site (3 sites)
#> Test: level
#> Alpha: 0.050
#> N simulations: 200
#> N converged: 200 (100.0%)
#>
#> ------------------------------------------------------------
#> POWER = 96.0% - Adequate (>= 80%)
#> ------------------------------------------------------------
One-call shortcut
result_quick <- estimate_and_calculate(
data = example_cfr_data,
level_change = -3, n_post = 30,
n_sim = 200, seed = 123
)
#>
#> ============================================================
#> PITS - ITS Parameter Estimation
#> ============================================================
#>
#> Estimation method: GLS-ML
#> Observations (n_pre): 24
#> Baseline: 14.5853
#> Sigma (residual SD): 1.4876 (10.2% of baseline)
#> Rho (AR1): 0.3832
#> Pre-trend (per unit): 0.0039
#>
#> ------------------------------------------------------------
#> Copy these into calculate_power():
#>
#> n_pre = 24
#> baseline = 14.5853
#> sigma = 1.4876
#> rho = 0.3832
#> level_change = ??? # YOUR clinical hypothesis
#> slope_change = 0 # set if testing slope
#>
#> ============================================================
cat(sprintf("Power: %.1f%% (%s)\n",
result_quick$power_pct, result_quick$interpretation))
#> Power: 86.0% (Adequate (>= 80%))