## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>",
  fig.width = 7, fig.height = 4.5
)
library(PITS)

## ----load-data----------------------------------------------------------------
data("example_cfr_data")
head(example_cfr_data)

## ----plot-predata, fig.cap = "Monthly CFR over 24 pre-intervention months."----
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))

## ----estimate-params----------------------------------------------------------
params <- estimate_its_params(example_cfr_data)

## ----hypothesis---------------------------------------------------------------
level_change <- -3.0   # 3 pp reduction in CFR

## ----single-design------------------------------------------------------------
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)

## ----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
)

## ----plot-sweep, fig.cap = "Power curve. Dashed red = 80% target; green = minimum adequate n_post."----
plot_power_curve(
  sweep,
  main = "Power vs follow-up duration — CDSS/CFR example",
  xlab = "Post-implementation months (n_post)",
  ylab = "Estimated power (%)"
)

## ----its-example, fig.cap = "Simulated ITS under the alternative hypothesis with fitted segmented regression."----
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"
)

## ----grid---------------------------------------------------------------------
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)

## ----multi-site---------------------------------------------------------------
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)

## ----summary-table, echo = FALSE----------------------------------------------
knitr::kable(
  sweep[, c("n_post", "power_pct", "interpretation")],
  col.names = c("n_post (months)", "Power (%)", "Interpretation"),
  caption = sprintf(
    "Power by follow-up duration. baseline = %.1f%%, sigma = %.2f, rho = %.2f, level_change = %g pp.",
    params$baseline, params$sigma, params$rho, level_change
  )
)

## ----one-call-----------------------------------------------------------------
result_quick <- estimate_and_calculate(
  data = example_cfr_data,
  level_change = -3, n_post = 30,
  n_sim = 200, seed = 123
)
cat(sprintf("Power: %.1f%%  (%s)\n",
            result_quick$power_pct, result_quick$interpretation))