STD all washouts

# ===== Packages (minimal, consistent) =====
library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(lubridate)

Attaching package: 'lubridate'

The following objects are masked from 'package:base':

    date, intersect, setdiff, union

library(purrr)
library(ggplot2)
library(lme4)

Loading required package: Matrix

library(lmerTest)

Attaching package: 'lmerTest'

The following object is masked from 'package:lme4':

    lmer

The following object is masked from 'package:stats':

    step

library(tibble)
library(scales)

Attaching package: 'scales'

The following object is masked from 'package:purrr':

    discard

library(gt)

# Optional safeguards against masking
select <- dplyr::select
filter <- dplyr::filter
lag <- dplyr::lag

# ===== Load + merge =====
cpa  <- read.csv("~/AHA paper/distribution/data/cpa1.csv")
hosp <- read.csv("~/AHA paper/distribution/data/hosp_characteristics.csv")
df0  <- merge(cpa, hosp, by = "UNIQUE_SITE_ID")

# ===== Cohort restrictions (CONFIRMED FINAL) =====
df_cohort <- df0 %>%
  filter(CPAAGE_YEARS_D >= 0.08, CPAAGE_YEARS_D <= 18) %>%
  filter(CPA_EVENT == "Y") %>%
  mutate(INDEXEVENT = ifelse(is.na(INDEXEVENT), 1L, as.integer(INDEXEVENT))) %>%
  filter(INDEXEVENT == 1L) %>%
  mutate(
    EVENTLOCATION = as.numeric(EVENTLOCATION),
    date = as.Date(DISCHARGEDATE, format = "%Y-%m-%d"),
    year = as.integer(format(date, "%Y"))
  ) %>%
  filter(EVENTLOCATION %in% c(12, 7, 19)) %>%
  filter(!is.na(date), year <= 2024)

nrow(df_cohort)

[1] 10368

# ===== Denominator checks =====
denoms <- tibble(
  cohort_n = nrow(df_cohort),
  surv_nonmissing = df_cohort %>% filter(!is.na(DISCHARGEDISPOSITION)) %>% nrow(),
  epi_nonmissing_120 = df_cohort %>%
    filter(!is.na(TIMEFROMPULSELESSTOEPI_D), TIMEFROMPULSELESSTOEPI_D <= 120) %>%
    nrow()
)
denoms

# A tibble: 1 × 3
  cohort_n surv_nonmissing epi_nonmissing_120
     <int>           <int>              <int>
1    10368           10367               4883

# ===== Fixed guideline cut dates (do not change) =====
cuts <- tibble::tribble(
  ~Update, ~cut,
  "2005",  as.Date("2006-05-01"),
  "2010",  as.Date("2011-05-01"),
  "2015",  as.Date("2016-05-01"),
  "2020",  as.Date("2021-05-01")
)

# ===== Make windows for arbitrary washout months =====
# For washout > 0:
#   wash_start = cut - wash_months
#   wash_end   = cut - 1 day
#   start_pre  = previous wash_start (first fixed at 2000-01-01)
#   end_post   = day before next wash_start (last fixed at 2023-12-31)
#
# For washout == 0:
#   No months are excluded; windows become contiguous cut-to-cut segments
make_windows <- function(wash_months = 6,
                         first_start_pre = as.Date("2000-01-01"),
                         final_end_post  = as.Date("2023-12-31")) {

  w <- cuts %>%
    arrange(cut) %>%
    mutate(
      wash_start = if (wash_months == 0) as.Date(NA) else (cut %m-% months(wash_months)),
      wash_end   = if (wash_months == 0) as.Date(NA) else (cut - days(1))
    )

  # start_pre
  w$start_pre <- first_start_pre
  if (nrow(w) > 1) {
    for (i in 2:nrow(w)) w$start_pre[i] <- w$wash_start[i - 1]
  }

  # end_post
  w$end_post <- final_end_post
  if (nrow(w) > 1) {
    for (i in 1:(nrow(w) - 1)) w$end_post[i] <- w$wash_start[i + 1] - days(1)
  }

  # Special-case: no washout (exclude nothing; define contiguous windows)
  if (wash_months == 0) {
    # Make the washout filter never exclude anything
    w$wash_start <- as.Date("9999-12-31")
    w$wash_end   <- as.Date("1900-01-01")

    # Define start/end windows by cut-to-cut
    w$start_pre <- c(first_start_pre, cuts$cut[-nrow(cuts)])
    w$end_post  <- c(cuts$cut[-1] - days(1), final_end_post)
  }

  w %>% mutate(across(start_pre:end_post, as.Date))
}

windows_0  <- make_windows(0)
windows_6  <- make_windows(6)
windows_12 <- make_windows(12)
windows_18 <- make_windows(18)

windows_0

# A tibble: 4 × 6
  Update cut        wash_start wash_end   start_pre  end_post  
  <chr>  <date>     <date>     <date>     <date>     <date>    
1 2005   2006-05-01 9999-12-31 1900-01-01 2000-01-01 2011-04-30
2 2010   2011-05-01 9999-12-31 1900-01-01 2006-05-01 2016-04-30
3 2015   2016-05-01 9999-12-31 1900-01-01 2011-05-01 2021-04-30
4 2020   2021-05-01 9999-12-31 1900-01-01 2016-05-01 2023-12-31

windows_6

# A tibble: 4 × 6
  Update cut        wash_start wash_end   start_pre  end_post  
  <chr>  <date>     <date>     <date>     <date>     <date>    
1 2005   2006-05-01 2005-11-01 2006-04-30 2000-01-01 2010-10-31
2 2010   2011-05-01 2010-11-01 2011-04-30 2005-11-01 2015-10-31
3 2015   2016-05-01 2015-11-01 2016-04-30 2010-11-01 2020-10-31
4 2020   2021-05-01 2020-11-01 2021-04-30 2015-11-01 2023-12-31

windows_12

# A tibble: 4 × 6
  Update cut        wash_start wash_end   start_pre  end_post  
  <chr>  <date>     <date>     <date>     <date>     <date>    
1 2005   2006-05-01 2005-05-01 2006-04-30 2000-01-01 2010-04-30
2 2010   2011-05-01 2010-05-01 2011-04-30 2005-05-01 2015-04-30
3 2015   2016-05-01 2015-05-01 2016-04-30 2010-05-01 2020-04-30
4 2020   2021-05-01 2020-05-01 2021-04-30 2015-05-01 2023-12-31

windows_18

# A tibble: 4 × 6
  Update cut        wash_start wash_end   start_pre  end_post  
  <chr>  <date>     <date>     <date>     <date>     <date>    
1 2005   2006-05-01 2004-11-01 2006-04-30 2000-01-01 2009-10-31
2 2010   2011-05-01 2009-11-01 2011-04-30 2004-11-01 2014-10-31
3 2015   2016-05-01 2014-11-01 2016-04-30 2009-11-01 2019-10-31
4 2020   2021-05-01 2019-11-01 2021-04-30 2014-11-01 2023-12-31

# ===== Survival recode =====
surv_df <- df_cohort %>%
  filter(!is.na(DISCHARGEDISPOSITION)) %>%
  mutate(survived = ifelse(as.character(DISCHARGEDISPOSITION) == "2", 1L, 0L))

table(surv_df$survived, useNA = "ifany")

   0    1 
5847 4520 

# ===== Survival: site × quarter table =====
quarterly_data <- surv_df %>%
  mutate(quarter = floor_date(date, "quarter")) %>%
  group_by(UNIQUE_SITE_ID, quarter) %>%
  summarise(
    n         = n(),
    successes = sum(survived == 1L, na.rm = TRUE),
    failures  = n - successes,
    surv_rate = successes / n,
    .groups   = "drop"
  )

# ===== Epi: site × month table (mean time) =====
epi_monthly <- df_cohort %>%
  filter(!is.na(TIMEFROMPULSELESSTOEPI_D), TIMEFROMPULSELESSTOEPI_D <= 120) %>%
  mutate(month_date = floor_date(date, "month")) %>%
  group_by(UNIQUE_SITE_ID, month_date) %>%
  summarise(
    n = n(),
    mean_epi = mean(TIMEFROMPULSELESSTOEPI_D, na.rm = TRUE),
    .groups = "drop"
  )

# ===== Epi: annual patient-level summary (for annual plot) =====
time_to_epi_year <- df_cohort %>%
  filter(!is.na(TIMEFROMPULSELESSTOEPI_D), TIMEFROMPULSELESSTOEPI_D <= 120) %>%
  mutate(year = year(date)) %>%
  group_by(year) %>%
  summarise(
    Mean = mean(TIMEFROMPULSELESSTOEPI_D, na.rm = TRUE),
    SEM  = sd(TIMEFROMPULSELESSTOEPI_D, na.rm = TRUE) / sqrt(n()),
    n    = n(),
    .groups = "drop"
  ) %>%
  filter(year <= 2024)

# ===== Survival ITS: quarter-level =====
build_its_quarter <- function(df, start_pre, wash_start, wash_end, cut, end_post) {
  df %>%
    filter(quarter >= start_pre, quarter <= end_post) %>%
    filter(!(quarter >= wash_start & quarter <= wash_end)) %>%
    mutate(
      post       = as.integer(quarter >= cut),
      time_after = if_else(quarter >= cut, as.numeric(quarter - cut) / 365.25, 0)
    )
}

analyze_its_quarter <- function(df) {

  # Risk difference (pp): weighted by n
  lin <- lmer(
    surv_rate ~ post + time_after + (1 | UNIQUE_SITE_ID),
    data = df, weights = n
  )
  bl <- fixef(lin); Vl <- as.matrix(vcov(lin))
  rd_post  <- bl["post"];       se_rd_post  <- sqrt(Vl["post","post"])
  rd_slope <- bl["time_after"]; se_rd_slope <- sqrt(Vl["time_after","time_after"])

  # Odds ratio: binomial GLMM using exact successes/failures (no explicit weights)
  glm <- glmer(
    cbind(successes, failures) ~ post + time_after + (1 | UNIQUE_SITE_ID),
    data = df, family = binomial(link = "logit"),
    control = glmerControl(optimizer = "bobyqa")
  )
  bp <- fixef(glm); Vp <- as.matrix(vcov(glm))
  logor_post  <- bp["post"];       se_logor_post  <- sqrt(Vp["post","post"])
  logor_slope <- bp["time_after"]; se_logor_slope <- sqrt(Vp["time_after","time_after"])

  z  <- function(e,s) ifelse(s > 0, e/s, NA_real_)
  pv <- function(e,s) 2 * pnorm(abs(z(e,s)), lower.tail = FALSE)

  tibble(
    Metric = c("Immediate change", "Slope difference"),
    `Risk Difference (pp)` = c(100*rd_post, 100*rd_slope),
    `RD 95% CI` = c(
      sprintf("(%.2f, %.2f)", 100*(rd_post - 1.96*se_rd_post),   100*(rd_post + 1.96*se_rd_post)),
      sprintf("(%.2f, %.2f)", 100*(rd_slope - 1.96*se_rd_slope), 100*(rd_slope + 1.96*se_rd_slope))
    ),
    `RD p` = c(pv(rd_post, se_rd_post), pv(rd_slope, se_rd_slope)),
    `Odds Ratio` = c(exp(logor_post), exp(logor_slope)),
    `OR 95% CI` = c(
      sprintf("(%.2f, %.2f)", exp(logor_post - 1.96*se_logor_post),  exp(logor_post + 1.96*se_logor_post)),
      sprintf("(%.2f, %.2f)", exp(logor_slope - 1.96*se_logor_slope),exp(logor_slope + 1.96*se_logor_slope))
    ),
    `OR p` = c(pv(logor_post, se_logor_post), pv(logor_slope, se_logor_slope))
  )
}

# ===== Epi ITS: monthly =====
build_its_month <- function(df, start_pre, wash_start, wash_end, cut, end_post) {
  df %>%
    filter(month_date >= start_pre, month_date <= end_post) %>%
    filter(!(month_date >= wash_start & month_date <= wash_end)) %>%
    mutate(
      post       = as.integer(month_date >= cut),
      time_after = if_else(month_date >= cut, as.numeric(month_date - cut) / 365.25, 0)
    )
}

analyze_its_epi <- function(df) {
  model <- lmer(mean_epi ~ post + time_after + (1 | UNIQUE_SITE_ID),
                data = df, weights = n)
  b <- fixef(model); V <- as.matrix(vcov(model))
  e_post  <- b["post"];       se_post  <- sqrt(V["post","post"])
  e_slope <- b["time_after"]; se_slope <- sqrt(V["time_after","time_after"])

  z  <- function(e,s) ifelse(s > 0, e/s, NA_real_)
  pv <- function(e,s) 2 * pnorm(abs(z(e,s)), lower.tail = FALSE)

  tibble(
    Metric = c("Immediate change", "Slope difference"),
    Estimate = c(e_post, e_slope),
    `95% CI` = c(
      sprintf("(%.2f, %.2f)", e_post  - 1.96*se_post,  e_post  + 1.96*se_post),
      sprintf("(%.2f, %.2f)", e_slope - 1.96*se_slope, e_slope + 1.96*se_slope)
    ),
    p_value = c(pv(e_post, se_post), pv(e_slope, se_slope))
  )
}

# ===== Run all 4 washout scenarios =====
run_scenario <- function(windows_dates, label) {

  # Survival quarter ITS
  surv_res <- purrr::pmap_dfr(
    as.list(as.data.frame(windows_dates, stringsAsFactors = FALSE)),
    function(Update, start_pre, wash_start, wash_end, cut, end_post) {
      d <- build_its_quarter(quarterly_data, start_pre, wash_start, wash_end, cut, end_post)
      message(sprintf("[SURV-Quarter | %s | %s] site-quarters=%d (pre=%d, post=%d)",
                      label, Update, nrow(d), sum(d$post == 0), sum(d$post == 1)))
      analyze_its_quarter(d) %>% mutate(Update = Update, .before = 1)
    }
  ) %>%
    mutate(
      Scenario = label,
      `Risk Difference (pp)` = round(`Risk Difference (pp)`, 2),
      `RD p` = signif(`RD p`, 3),
      `Odds Ratio` = round(`Odds Ratio`, 3),
      `OR p` = signif(`OR p`, 3)
    ) %>%
    dplyr::select(
      Scenario, Update, Metric,
      `Risk Difference (pp)`, `RD 95% CI`, `RD p`,
      `Odds Ratio`, `OR 95% CI`, `OR p`
    ) %>%
    arrange(Update, match(Metric, c("Immediate change", "Slope difference")))

  # Epi monthly ITS
  epi_res <- purrr::pmap_dfr(
    as.list(as.data.frame(windows_dates, stringsAsFactors = FALSE)),
    function(Update, start_pre, wash_start, wash_end, cut, end_post) {
      d <- build_its_month(epi_monthly, start_pre, wash_start, wash_end, cut, end_post)
      message(sprintf("[EPI-Month | %s | %s] site-months=%d (pre=%d, post=%d)",
                      label, Update, nrow(d), sum(d$post == 0), sum(d$post == 1)))
      analyze_its_epi(d) %>% mutate(Update = Update, .before = 1)
    }
  ) %>%
    mutate(
      Scenario = label,
      Estimate = round(Estimate, 2),
      p_value  = signif(p_value, 3)
    ) %>%
    dplyr::select(Scenario, Update, Metric, Estimate, `95% CI`, p_value) %>%
    arrange(Update, match(Metric, c("Immediate change", "Slope difference")))

  list(survival_quarter = surv_res, epi_month = epi_res)
}

scenario_list <- list(
  "No washout"    = windows_0,
  "Washout 6 mo"  = windows_6,
  "Washout 12 mo" = windows_12,
  "Washout 18 mo" = windows_18
)

all_results <- purrr::imap(scenario_list, run_scenario)

[SURV-Quarter | No washout | 2005] site-quarters=884 (pre=241, post=643)

[SURV-Quarter | No washout | 2010] site-quarters=1250 (pre=643, post=607)

[SURV-Quarter | No washout | 2015] site-quarters=1422 (pre=607, post=815)

[SURV-Quarter | No washout | 2020] site-quarters=1308 (pre=815, post=493)

[EPI-Month | No washout | 2005] site-months=862 (pre=184, post=678)

[EPI-Month | No washout | 2010] site-months=1324 (pre=678, post=646)

[EPI-Month | No washout | 2015] site-months=1610 (pre=646, post=964)

[EPI-Month | No washout | 2020] site-months=1589 (pre=964, post=625)

[SURV-Quarter | Washout 6 mo | 2005] site-quarters=780 (pre=186, post=594)

[SURV-Quarter | Washout 6 mo | 2010] site-quarters=1190 (pre=649, post=541)

[SURV-Quarter | Washout 6 mo | 2015] site-quarters=1328 (pre=590, post=738)

[SURV-Quarter | Washout 6 mo | 2020] site-quarters=1297 (pre=804, post=493)

[EPI-Month | Washout 6 mo | 2005] site-months=761 (pre=136, post=625)

[EPI-Month | Washout 6 mo | 2010] site-months=1244 (pre=673, post=571)

[EPI-Month | Washout 6 mo | 2015] site-months=1506 (pre=624, post=882)

[EPI-Month | Washout 6 mo | 2020] site-months=1582 (pre=957, post=625)

[SURV-Quarter | Washout 12 mo | 2005] site-quarters=698 (pre=159, post=539)

[SURV-Quarter | Washout 12 mo | 2010] site-quarters=1109 (pre=621, post=488)

[SURV-Quarter | Washout 12 mo | 2015] site-quarters=1253 (pre=592, post=661)

[SURV-Quarter | Washout 12 mo | 2020] site-quarters=1273 (pre=780, post=493)

[EPI-Month | Washout 12 mo | 2005] site-months=671 (pre=117, post=554)

[EPI-Month | Washout 12 mo | 2010] site-months=1125 (pre=621, post=504)

[EPI-Month | Washout 12 mo | 2015] site-months=1406 (pre=628, post=778)

[EPI-Month | Washout 12 mo | 2020] site-months=1545 (pre=920, post=625)

[SURV-Quarter | Washout 18 mo | 2005] site-quarters=595 (pre=132, post=463)

[SURV-Quarter | Washout 18 mo | 2010] site-quarters=993 (pre=572, post=421)

[SURV-Quarter | Washout 18 mo | 2015] site-quarters=1174 (pre=601, post=573)

[SURV-Quarter | Washout 18 mo | 2020] site-quarters=1252 (pre=759, post=493)

[EPI-Month | Washout 18 mo | 2005] site-months=572 (pre=98, post=474)

[EPI-Month | Washout 18 mo | 2010] site-months=993 (pre=560, post=433)

[EPI-Month | Washout 18 mo | 2015] site-months=1320 (pre=637, post=683)

[EPI-Month | Washout 18 mo | 2020] site-months=1521 (pre=896, post=625)

# ===== Combined survival table across all scenarios =====
survival_all <- bind_rows(lapply(all_results, `[[`, "survival_quarter")) %>%
  arrange(
    factor(Scenario, levels = c("No washout", "Washout 6 mo", "Washout 12 mo", "Washout 18 mo")),
    Update,
    match(Metric, c("Immediate change", "Slope difference"))
  )

print(survival_all, n = Inf, width = Inf)

# A tibble: 32 × 9
   Scenario      Update Metric           `Risk Difference (pp)` `RD 95% CI`   
   <chr>         <chr>  <chr>                             <dbl> <chr>         
 1 No washout    2005   Immediate change                   9.27 (2.47, 16.06) 
 2 No washout    2005   Slope difference                   0.96 (-0.65, 2.58) 
 3 No washout    2010   Immediate change                   0.72 (-4.50, 5.94) 
 4 No washout    2010   Slope difference                   1.22 (-0.33, 2.76) 
 5 No washout    2015   Immediate change                   5.57 (1.12, 10.01) 
 6 No washout    2015   Slope difference                   0.13 (-1.25, 1.51) 
 7 No washout    2020   Immediate change                  -1.5  (-7.00, 3.99) 
 8 No washout    2020   Slope difference                   1    (-2.45, 4.45) 
 9 Washout 6 mo  2005   Immediate change                  10.5  (1.89, 19.04) 
10 Washout 6 mo  2005   Slope difference                   0.49 (-1.36, 2.34) 
11 Washout 6 mo  2010   Immediate change                   1.01 (-4.44, 6.45) 
12 Washout 6 mo  2010   Slope difference                   1.64 (-0.19, 3.48) 
13 Washout 6 mo  2015   Immediate change                   5.4  (0.73, 10.08) 
14 Washout 6 mo  2015   Slope difference                   0.39 (-1.20, 1.98) 
15 Washout 6 mo  2020   Immediate change                  -1.41 (-6.92, 4.10) 
16 Washout 6 mo  2020   Slope difference                   1.16 (-2.29, 4.60) 
17 Washout 12 mo 2005   Immediate change                   9.01 (-0.25, 18.28)
18 Washout 12 mo 2005   Slope difference                   0.8  (-1.29, 2.90) 
19 Washout 12 mo 2010   Immediate change                   2.19 (-3.54, 7.92) 
20 Washout 12 mo 2010   Slope difference                   1.27 (-0.87, 3.42) 
21 Washout 12 mo 2015   Immediate change                   5.38 (0.49, 10.27) 
22 Washout 12 mo 2015   Slope difference                   1.08 (-0.76, 2.93) 
23 Washout 12 mo 2020   Immediate change                  -1.71 (-7.22, 3.80) 
24 Washout 12 mo 2020   Slope difference                   1.26 (-2.19, 4.70) 
25 Washout 18 mo 2005   Immediate change                  10.8  (0.85, 20.77) 
26 Washout 18 mo 2005   Slope difference                  -0.71 (-3.29, 1.87) 
27 Washout 18 mo 2010   Immediate change                   2.16 (-3.89, 8.21) 
28 Washout 18 mo 2010   Slope difference                   1.72 (-0.86, 4.30) 
29 Washout 18 mo 2015   Immediate change                   4.75 (-0.48, 9.98) 
30 Washout 18 mo 2015   Slope difference                   1.51 (-0.74, 3.75) 
31 Washout 18 mo 2020   Immediate change                  -0.88 (-6.46, 4.70) 
32 Washout 18 mo 2020   Slope difference                   1.12 (-2.36, 4.59) 
    `RD p` `Odds Ratio` `OR 95% CI`   `OR p`
     <dbl>        <dbl> <chr>          <dbl>
 1 0.00752        1.55  (1.17, 2.06) 0.00256
 2 0.243          1.04  (0.98, 1.11) 0.181  
 3 0.786          1.02  (0.84, 1.24) 0.852  
 4 0.123          1.06  (1.00, 1.12) 0.056  
 5 0.0142         1.27  (1.07, 1.50) 0.00678
 6 0.856          1.01  (0.95, 1.06) 0.834  
 7 0.592          0.938 (0.76, 1.16) 0.552  
 8 0.569          1.04  (0.91, 1.19) 0.537  
 9 0.0167         1.64  (1.13, 2.36) 0.00894
10 0.603          1.03  (0.95, 1.10) 0.491  
11 0.717          1.03  (0.84, 1.27) 0.787  
12 0.0788         1.08  (1.01, 1.16) 0.0343 
13 0.0235         1.26  (1.05, 1.51) 0.0118 
14 0.63           1.02  (0.96, 1.08) 0.586  
15 0.616          0.941 (0.76, 1.16) 0.575  
16 0.511          1.05  (0.92, 1.20) 0.474  
17 0.0565         1.51  (1.01, 2.24) 0.043  
18 0.453          1.04  (0.96, 1.13) 0.313  
19 0.453          1.08  (0.87, 1.34) 0.491  
20 0.244          1.06  (0.98, 1.15) 0.145  
21 0.0312         1.26  (1.04, 1.52) 0.017  
22 0.249          1.05  (0.98, 1.12) 0.204  
23 0.543          0.929 (0.75, 1.15) 0.502  
24 0.474          1.05  (0.92, 1.21) 0.436  
25 0.0333         1.65  (1.07, 2.55) 0.0241 
26 0.59           0.977 (0.88, 1.08) 0.662  
27 0.484          1.07  (0.84, 1.36) 0.568  
28 0.192          1.09  (0.98, 1.20) 0.11   
29 0.0748         1.23  (1.00, 1.51) 0.0449 
30 0.189          1.07  (0.98, 1.16) 0.143  
31 0.758          0.962 (0.78, 1.19) 0.722  
32 0.53           1.05  (0.92, 1.20) 0.489  

# Optional prettier HTML table
gt::gt(survival_all)

Scenario	Update	Metric	Risk Difference (pp)	RD 95% CI	RD p	Odds Ratio	OR 95% CI	OR p
No washout	2005	Immediate change	9.27	(2.47, 16.06)	0.00752	1.551	(1.17, 2.06)	0.00256
No washout	2005	Slope difference	0.96	(-0.65, 2.58)	0.24300	1.045	(0.98, 1.11)	0.18100
No washout	2010	Immediate change	0.72	(-4.50, 5.94)	0.78600	1.019	(0.84, 1.24)	0.85200
No washout	2010	Slope difference	1.22	(-0.33, 2.76)	0.12300	1.058	(1.00, 1.12)	0.05600
No washout	2015	Immediate change	5.57	(1.12, 10.01)	0.01420	1.266	(1.07, 1.50)	0.00678
No washout	2015	Slope difference	0.13	(-1.25, 1.51)	0.85600	1.006	(0.95, 1.06)	0.83400
No washout	2020	Immediate change	-1.50	(-7.00, 3.99)	0.59200	0.938	(0.76, 1.16)	0.55200
No washout	2020	Slope difference	1.00	(-2.45, 4.45)	0.56900	1.043	(0.91, 1.19)	0.53700
Washout 6 mo	2005	Immediate change	10.47	(1.89, 19.04)	0.01670	1.635	(1.13, 2.36)	0.00894
Washout 6 mo	2005	Slope difference	0.49	(-1.36, 2.34)	0.60300	1.026	(0.95, 1.10)	0.49100
Washout 6 mo	2010	Immediate change	1.01	(-4.44, 6.45)	0.71700	1.029	(0.84, 1.27)	0.78700
Washout 6 mo	2010	Slope difference	1.64	(-0.19, 3.48)	0.07880	1.078	(1.01, 1.16)	0.03430
Washout 6 mo	2015	Immediate change	5.40	(0.73, 10.08)	0.02350	1.260	(1.05, 1.51)	0.01180
Washout 6 mo	2015	Slope difference	0.39	(-1.20, 1.98)	0.63000	1.017	(0.96, 1.08)	0.58600
Washout 6 mo	2020	Immediate change	-1.41	(-6.92, 4.10)	0.61600	0.941	(0.76, 1.16)	0.57500
Washout 6 mo	2020	Slope difference	1.16	(-2.29, 4.60)	0.51100	1.050	(0.92, 1.20)	0.47400
Washout 12 mo	2005	Immediate change	9.01	(-0.25, 18.28)	0.05650	1.508	(1.01, 2.24)	0.04300
Washout 12 mo	2005	Slope difference	0.80	(-1.29, 2.90)	0.45300	1.044	(0.96, 1.13)	0.31300
Washout 12 mo	2010	Immediate change	2.19	(-3.54, 7.92)	0.45300	1.080	(0.87, 1.34)	0.49100
Washout 12 mo	2010	Slope difference	1.27	(-0.87, 3.42)	0.24400	1.062	(0.98, 1.15)	0.14500
Washout 12 mo	2015	Immediate change	5.38	(0.49, 10.27)	0.03120	1.260	(1.04, 1.52)	0.01700
Washout 12 mo	2015	Slope difference	1.08	(-0.76, 2.93)	0.24900	1.047	(0.98, 1.12)	0.20400
Washout 12 mo	2020	Immediate change	-1.71	(-7.22, 3.80)	0.54300	0.929	(0.75, 1.15)	0.50200
Washout 12 mo	2020	Slope difference	1.26	(-2.19, 4.70)	0.47400	1.054	(0.92, 1.21)	0.43600
Washout 18 mo	2005	Immediate change	10.81	(0.85, 20.77)	0.03330	1.652	(1.07, 2.55)	0.02410
Washout 18 mo	2005	Slope difference	-0.71	(-3.29, 1.87)	0.59000	0.977	(0.88, 1.08)	0.66200
Washout 18 mo	2010	Immediate change	2.16	(-3.89, 8.21)	0.48400	1.072	(0.84, 1.36)	0.56800
Washout 18 mo	2010	Slope difference	1.72	(-0.86, 4.30)	0.19200	1.086	(0.98, 1.20)	0.11000
Washout 18 mo	2015	Immediate change	4.75	(-0.48, 9.98)	0.07480	1.230	(1.00, 1.51)	0.04490
Washout 18 mo	2015	Slope difference	1.51	(-0.74, 3.75)	0.18900	1.067	(0.98, 1.16)	0.14300
Washout 18 mo	2020	Immediate change	-0.88	(-6.46, 4.70)	0.75800	0.962	(0.78, 1.19)	0.72200
Washout 18 mo	2020	Slope difference	1.12	(-2.36, 4.59)	0.53000	1.048	(0.92, 1.20)	0.48900

# ===== Combined epi table across all scenarios =====
epi_all <- bind_rows(lapply(all_results, `[[`, "epi_month")) %>%
  arrange(
    factor(Scenario, levels = c("No washout", "Washout 6 mo", "Washout 12 mo", "Washout 18 mo")),
    Update,
    match(Metric, c("Immediate change", "Slope difference"))
  )

print(epi_all, n = Inf, width = Inf)

# A tibble: 32 × 6
   Scenario      Update Metric           Estimate `95% CI`       p_value
   <chr>         <chr>  <chr>               <dbl> <chr>            <dbl>
 1 No washout    2005   Immediate change    -0.67 (-1.50, 0.17)  0.116  
 2 No washout    2005   Slope difference    -0.11 (-0.31, 0.10)  0.306  
 3 No washout    2010   Immediate change    -0.13 (-0.69, 0.43)  0.653  
 4 No washout    2010   Slope difference    -0.14 (-0.32, 0.03)  0.102  
 5 No washout    2015   Immediate change    -0.11 (-0.47, 0.26)  0.571  
 6 No washout    2015   Slope difference     0.07 (-0.04, 0.18)  0.239  
 7 No washout    2020   Immediate change     0.16 (-0.26, 0.59)  0.443  
 8 No washout    2020   Slope difference    -0.06 (-0.31, 0.20)  0.658  
 9 Washout 6 mo  2005   Immediate change    -1.24 (-2.25, -0.24) 0.015  
10 Washout 6 mo  2005   Slope difference    -0.04 (-0.28, 0.20)  0.751  
11 Washout 6 mo  2010   Immediate change    -0.21 (-0.82, 0.39)  0.485  
12 Washout 6 mo  2010   Slope difference    -0.14 (-0.35, 0.07)  0.192  
13 Washout 6 mo  2015   Immediate change    -0.17 (-0.54, 0.19)  0.356  
14 Washout 6 mo  2015   Slope difference     0.09 (-0.04, 0.21)  0.165  
15 Washout 6 mo  2020   Immediate change     0.2  (-0.18, 0.59)  0.303  
16 Washout 6 mo  2020   Slope difference    -0.03 (-0.26, 0.20)  0.802  
17 Washout 12 mo 2005   Immediate change    -1.45 (-2.54, -0.35) 0.00972
18 Washout 12 mo 2005   Slope difference    -0.07 (-0.36, 0.22)  0.62   
19 Washout 12 mo 2010   Immediate change    -0.18 (-0.82, 0.47)  0.595  
20 Washout 12 mo 2010   Slope difference    -0.16 (-0.41, 0.09)  0.218  
21 Washout 12 mo 2015   Immediate change    -0.29 (-0.69, 0.11)  0.154  
22 Washout 12 mo 2015   Slope difference     0.1  (-0.05, 0.25)  0.189  
23 Washout 12 mo 2020   Immediate change     0.22 (-0.16, 0.61)  0.255  
24 Washout 12 mo 2020   Slope difference    -0.02 (-0.25, 0.21)  0.864  
25 Washout 18 mo 2005   Immediate change    -1.57 (-2.81, -0.33) 0.013  
26 Washout 18 mo 2005   Slope difference     0.08 (-0.30, 0.45)  0.69   
27 Washout 18 mo 2010   Immediate change    -0.36 (-1.09, 0.36)  0.33   
28 Washout 18 mo 2010   Slope difference    -0.11 (-0.44, 0.22)  0.511  
29 Washout 18 mo 2015   Immediate change    -0.17 (-0.58, 0.24)  0.414  
30 Washout 18 mo 2015   Slope difference    -0.04 (-0.22, 0.14)  0.646  
31 Washout 18 mo 2020   Immediate change     0.3  (-0.06, 0.67)  0.0999 
32 Washout 18 mo 2020   Slope difference    -0.02 (-0.23, 0.20)  0.866  

# Optional prettier HTML table
gt::gt(epi_all)

Scenario	Update	Metric	Estimate	95% CI	p_value
No washout	2005	Immediate change	-0.67	(-1.50, 0.17)	0.11600
No washout	2005	Slope difference	-0.11	(-0.31, 0.10)	0.30600
No washout	2010	Immediate change	-0.13	(-0.69, 0.43)	0.65300
No washout	2010	Slope difference	-0.14	(-0.32, 0.03)	0.10200
No washout	2015	Immediate change	-0.11	(-0.47, 0.26)	0.57100
No washout	2015	Slope difference	0.07	(-0.04, 0.18)	0.23900
No washout	2020	Immediate change	0.16	(-0.26, 0.59)	0.44300
No washout	2020	Slope difference	-0.06	(-0.31, 0.20)	0.65800
Washout 6 mo	2005	Immediate change	-1.24	(-2.25, -0.24)	0.01500
Washout 6 mo	2005	Slope difference	-0.04	(-0.28, 0.20)	0.75100
Washout 6 mo	2010	Immediate change	-0.21	(-0.82, 0.39)	0.48500
Washout 6 mo	2010	Slope difference	-0.14	(-0.35, 0.07)	0.19200
Washout 6 mo	2015	Immediate change	-0.17	(-0.54, 0.19)	0.35600
Washout 6 mo	2015	Slope difference	0.09	(-0.04, 0.21)	0.16500
Washout 6 mo	2020	Immediate change	0.20	(-0.18, 0.59)	0.30300
Washout 6 mo	2020	Slope difference	-0.03	(-0.26, 0.20)	0.80200
Washout 12 mo	2005	Immediate change	-1.45	(-2.54, -0.35)	0.00972
Washout 12 mo	2005	Slope difference	-0.07	(-0.36, 0.22)	0.62000
Washout 12 mo	2010	Immediate change	-0.18	(-0.82, 0.47)	0.59500
Washout 12 mo	2010	Slope difference	-0.16	(-0.41, 0.09)	0.21800
Washout 12 mo	2015	Immediate change	-0.29	(-0.69, 0.11)	0.15400
Washout 12 mo	2015	Slope difference	0.10	(-0.05, 0.25)	0.18900
Washout 12 mo	2020	Immediate change	0.22	(-0.16, 0.61)	0.25500
Washout 12 mo	2020	Slope difference	-0.02	(-0.25, 0.21)	0.86400
Washout 18 mo	2005	Immediate change	-1.57	(-2.81, -0.33)	0.01300
Washout 18 mo	2005	Slope difference	0.08	(-0.30, 0.45)	0.69000
Washout 18 mo	2010	Immediate change	-0.36	(-1.09, 0.36)	0.33000
Washout 18 mo	2010	Slope difference	-0.11	(-0.44, 0.22)	0.51100
Washout 18 mo	2015	Immediate change	-0.17	(-0.58, 0.24)	0.41400
Washout 18 mo	2015	Slope difference	-0.04	(-0.22, 0.14)	0.64600
Washout 18 mo	2020	Immediate change	0.30	(-0.06, 0.67)	0.09990
Washout 18 mo	2020	Slope difference	-0.02	(-0.23, 0.20)	0.86600

# Same colors as your preferred figure
pals_cols <- c(
  "2005" = "#F15A5A",
  "2010" = "#7CB518",
  "2015" = "#1AA6B7",
  "2020" = "#9B5DE5"
)

plot_epi_annual_by_scenario <- function(windows_dates, label) {

  annual_base <- time_to_epi_year %>%
    mutate(year_date = as.Date(paste0(year, "-01-01")))

  annual_with_periods <- purrr::pmap_dfr(
    as.list(as.data.frame(windows_dates, stringsAsFactors = FALSE)),
    function(Update, start_pre, wash_start, wash_end, cut, end_post) {

      annual_base %>%
        filter(year_date >= start_pre, year_date <= end_post) %>%
        filter(!(year_date >= wash_start & year_date <= wash_end)) %>%
        mutate(
          Update  = as.character(Update),
          prepost = if_else(year_date < cut, "Pre", "Post"),
          Scenario = label
        )
    }
  ) %>%
    arrange(year_date, if_else(Update == "2015", 2L, 1L)) %>%
    mutate(
      Update  = factor(Update, levels = c("2005", "2010", "2015", "2020")),
      prepost = factor(prepost, levels = c("Pre", "Post"))
    )

  p <- ggplot(
    annual_with_periods,
    aes(x = year, y = Mean,
        color = Update, linetype = prepost,
        group = interaction(Update, prepost))
  )

  # Only add washout shading if this is NOT the no-washout scenario
  if (!grepl("No washout", label)) {
    p <- p +
      geom_rect(
        data = windows_dates,
        aes(xmin = year(wash_start), xmax = year(wash_end), ymin = -Inf, ymax = Inf),
        inherit.aes = FALSE,
        fill = "grey85", alpha = 0.18
      )
  }

  p +
    geom_vline(
      data = windows_dates,
      aes(xintercept = year(cut)),
      linetype = "dashed", color = "grey40", linewidth = 0.6
    ) +
    geom_line(linewidth = 0.9) +
    geom_point(size = 1.6, alpha = 0.7) +
    geom_errorbar(
      aes(ymin = Mean - SEM, ymax = Mean + SEM),
      width = 0.18, linewidth = 0.55, alpha = 0.9
    ) +
    scale_color_manual(values = pals_cols, drop = FALSE) +
    scale_x_continuous(
      breaks = seq(2000, 2024, 1),
      limits = c(2000, 2024)
    ) +
    guides(
      color = guide_legend(order = 1, override.aes = list(linetype = 1, linewidth = 1.2)),
      linetype = guide_legend(order = 2)
    ) +
    labs(
      title = paste0("Time to First Epinephrine (≤120 min) — ", label),
      subtitle = ifelse(
        grepl("No washout", label),
        "Annual mean ± SEM, colored by PALS update (no washout; contiguous pre/post windows)",
        "Annual mean ± SEM, colored by PALS update (pre/post windows; washout excluded per scenario)"
      ),
      x = "Year of Event",
      y = "Time to First Epinephrine (minutes)",
      color = "Guideline update",
      linetype = "Pre vs Post"
    ) +
    theme_minimal(base_size = 12) +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1),
      legend.position = "bottom",
      legend.box = "vertical",
      plot.title = element_text(face = "bold")
    )
}
plot_epi_annual_by_scenario(windows_0,  "No washout")

plot_epi_annual_by_scenario(windows_6,  "Washout 6 mo")

plot_epi_annual_by_scenario(windows_12, "Washout 12 mo")

plot_epi_annual_by_scenario(windows_18, "Washout 18 mo")

monthly_base <- surv_df %>%
  mutate(month_date = floor_date(date, "month")) %>%
  group_by(month_date) %>%
  summarise(
    total      = n(),
    survived   = sum(survived == 1L, na.rm = TRUE),
    proportion = survived / total,
    .groups    = "drop"
  )

pals_cols <- c(
  "2005" = "#F15A5A",
  "2010" = "#7CB518",
  "2015" = "#1AA6B7",
  "2020" = "#9B5DE5"
)

plot_surv_codeA_style <- function(windows_dates, label) {

  monthly_with_periods <- purrr::pmap_dfr(
    as.list(as.data.frame(windows_dates, stringsAsFactors = FALSE)),
    function(Update, start_pre, wash_start, wash_end, cut, end_post) {

      monthly_base %>%
        filter(month_date >= start_pre, month_date <= end_post) %>%
        filter(!(month_date >= wash_start & month_date <= wash_end)) %>%
        mutate(
          Update  = as.character(Update),
          prepost = if_else(month_date < cut, "Pre", "Post"),
          Scenario = label
        )
    }
  ) %>%
    filter(month_date >= as.Date("2000-01-01")) %>%
    arrange(month_date, if_else(Update == "2015", 2L, 1L)) %>%
    mutate(
      Update  = factor(Update, levels = c("2005", "2010", "2015", "2020")),
      prepost = factor(prepost, levels = c("Pre", "Post"))
    )

  p <- ggplot(
    monthly_with_periods,
    aes(x = month_date, y = proportion,
        color = Update, linetype = prepost)
  )

  # Only add washout shading when this is NOT the no-washout scenario
  if (!grepl("No washout", label)) {
    p <- p +
      geom_rect(
        data = windows_dates,
        aes(xmin = wash_start, xmax = wash_end, ymin = -Inf, ymax = Inf),
        inherit.aes = FALSE,
        fill = "grey85", alpha = 0.30
      )
  }

  p +
    geom_point(size = 0.9, stroke = 0.4, shape = 16, alpha = 0.85) +
    geom_smooth(method = "lm", se = TRUE, aes(weight = total), linewidth = 0.8) +
    geom_vline(
      data = windows_dates,
      aes(xintercept = cut),
      color = "grey30", linetype = "dashed", linewidth = 0.5
    ) +
    scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
    scale_x_date(
      date_labels = "%Y",
      breaks = seq(as.Date("2000-01-01"), as.Date("2024-01-01"), by = "1 year"),
      limits = c(as.Date("2000-01-01"), as.Date("2024-01-01"))
    ) +
    scale_color_manual(values = pals_cols, drop = FALSE) +
    labs(
      title = paste0("Survival to Discharge by PALS Guideline Periods — ", label),
      subtitle = ifelse(
        grepl("No washout", label),
        "Monthly proportions using contiguous pre/post windows with no washout",
        "Monthly proportions using same pre/post definitions as ITS (each update’s washout excluded)"
      ),
      x = "Date",
      y = "Survival to discharge (%)",
      color = "Guideline update",
      linetype = "Pre vs Post"
    ) +
    theme_minimal(base_size = 12) +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1),
      legend.position = "bottom",
      plot.title = element_text(face = "bold")
    )
}

plot_surv_codeA_style(windows_0,  "No washout")

`geom_smooth()` using formula = 'y ~ x'

plot_surv_codeA_style(windows_6,  "Washout 6 mo")

`geom_smooth()` using formula = 'y ~ x'

plot_surv_codeA_style(windows_12, "Washout 12 mo")

`geom_smooth()` using formula = 'y ~ x'

plot_surv_codeA_style(windows_18, "Washout 18 mo")

`geom_smooth()` using formula = 'y ~ x'

library(stringr)

survival_effects <- survival_all %>%
  mutate(
    washout_months = case_when(
      Scenario == "No washout"    ~ 0,
      Scenario == "Washout 6 mo"  ~ 6,
      Scenario == "Washout 12 mo" ~ 12,
      Scenario == "Washout 18 mo" ~ 18
    )
  )

plot_data <- survival_effects %>%
  select(Update, Metric, washout_months, `Risk Difference (pp)`)

ggplot(plot_data,
       aes(x = washout_months,
           y = `Risk Difference (pp)`,
           color = Metric,
           group = Metric)) +

  geom_line(linewidth = 1) +
  geom_point(size = 3) +

  facet_wrap(~ Update, nrow = 1) +

  scale_x_continuous(
    breaks = c(0, 6, 12, 18),
    labels = c("0\n(None)", "6", "12", "18")
  ) +

  labs(
    title = "Sensitivity of Survival Effect Size to Washout Duration",
    subtitle = "Interrupted time-series estimates across washout assumptions",
    x = "Washout duration (months)",
    y = "Risk Difference (percentage points)",
    color = "ITS Effect"
  ) +

  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "bottom"
  )

plot_data <- survival_effects %>%
  tidyr::separate(`RD 95% CI`,
                  into = c("low","high"),
                  sep = ",",
                  remove = FALSE) %>%
  mutate(
    low  = as.numeric(gsub("[()]", "", low)),
    high = as.numeric(gsub("[()]", "", high))
  )

ggplot(plot_data,
       aes(x = washout_months,
           y = `Risk Difference (pp)`,
           color = Metric,
           group = Metric)) +

  geom_line(linewidth = 1) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = low, ymax = high),
              width = 0.8, alpha = 0.6) +

  facet_wrap(~ Update, nrow = 1) +

  scale_x_continuous(
    breaks = c(0, 6, 12, 18),
    labels = c("0\n(None)", "6", "12", "18")
  ) +

  labs(
    title = "Sensitivity of Survival Effect Size to Washout Duration",
    subtitle = "Interrupted time-series estimates across washout assumptions",
    x = "Washout duration (months)",
    y = "Risk Difference (percentage points)",
    color = "ITS Effect"
  ) +

  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "bottom"
  )

library(dplyr)
library(tidyr)

Attaching package: 'tidyr'

The following objects are masked from 'package:Matrix':

    expand, pack, unpack

library(stringr)
library(ggplot2)

# ---- PALS palette (same as your epi plot) ----
pals_cols <- c(
  "2005" = "#F15A5A",
  "2010" = "#7CB518",
  "2015" = "#1AA6B7",
  "2020" = "#9B5DE5"
)

# ---- Build plot dataset from survival_all ----
sens_rd <- survival_all %>%
  mutate(
    washout_months = case_when(
      Scenario == "No washout"    ~ 0,
      Scenario == "Washout 6 mo"  ~ 6,
      Scenario == "Washout 12 mo" ~ 12,
      Scenario == "Washout 18 mo" ~ 18,
      TRUE ~ NA_real_
    ),
    sig = if_else(`RD p` < 0.05, "*", "")
  ) %>%
  # Parse "(low, high)" from RD 95% CI into numeric bounds
  separate(`RD 95% CI`, into = c("low", "high"), sep = ",\\s*", remove = FALSE) %>%
  mutate(
    low  = as.numeric(str_remove_all(low,  "[\\(\\)]")),
    high = as.numeric(str_remove_all(high, "[\\(\\)]")),
    Update = factor(Update, levels = c("2005", "2010", "2015", "2020")),
    Metric = factor(Metric, levels = c("Immediate change", "Slope difference"))
  ) %>%
  filter(!is.na(washout_months))

# ---- Plot: single comparison figure (facet by Metric; color by Update) ----
p_sens <- ggplot(sens_rd,
                 aes(x = washout_months,
                     y = `Risk Difference (pp)`,
                     color = Update,
                     group = Update)) +

  # 95% CI error bars
  geom_errorbar(aes(ymin = low, ymax = high), width = 0.8, linewidth = 0.6, alpha = 0.75) +

  # line + points
  geom_line(linewidth = 1.0, alpha = 0.9) +
  geom_point(size = 3.0) +

  # significance marker (RD p < 0.05)
  geom_text(aes(label = sig), vjust = -1.1, show.legend = FALSE, size = 5) +

  facet_wrap(~ Metric, ncol = 1, scales = "free_y") +

  scale_color_manual(values = pals_cols, drop = FALSE) +
  scale_x_continuous(
    breaks = c(0, 6, 12, 18),
    labels = c("0 (None)", "6", "12", "18"),
    limits = c(-0.5, 18.5)
  ) +

  labs(
    title = "Sensitivity Analysis: Survival Effect Size vs Washout Duration",
    subtitle = "Quarter-level ITS estimates (Risk Difference, percentage points). Error bars = 95% CI. * = p < 0.05",
    x = "Washout duration (months)",
    y = "Risk Difference (percentage points)",
    color = "PALS update"
  ) +

  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold"),
    strip.text = element_text(face = "bold"),
    legend.position = "bottom",
    panel.grid.minor = element_blank()
  )

p_sens

# PNG (journal-friendly raster)
ggsave("sens_survival_RD_vs_washout_300dpi.png",
       plot = p_sens, width = 10, height = 7, dpi = 300)

# Optional: PDF (vector, great for journals)
ggsave("sens_survival_RD_vs_washout.pdf",
       plot = p_sens, width = 10, height = 7)

library(ggplot2)
library(dplyr)

# ---- Conceptual curves (no real data; just shapes) ----
# Panel 1: System reset (immediate jump)
sys <- tibble(
  mech = "System reset\n(immediate level shift)",
  t = c(0, 4.9, 5.0, 10),
  y = c(0.20, 0.20, 0.30, 0.30)
)

# Panel 2: Behavioral adoption (gradual improvement / slope change)
beh <- tibble(
  mech = "Behavioral adoption\n(progressive slope change)",
  t = seq(0, 10, by = 0.2)
) %>%
  mutate(
    y = if_else(t < 5, 0.20, 0.20 + 0.02*(t - 5))  # gradual rise after implementation
  )

# Panel 3: Incremental refinement (small level shift then plateau)
ref <- tibble(
  mech = "Incremental refinement\n(small shift then stabilization)",
  t = c(0, 4.9, 5.0, 10),
  y = c(0.20, 0.20, 0.24, 0.24)
)

concept_df <- bind_rows(sys, beh, ref)

# Helper for the dashed implementation line
impl_line <- tibble(t = 5)

p_concept <- ggplot(concept_df, aes(x = t, y = y)) +
  geom_line(linewidth = 1.1) +
  geom_vline(data = impl_line, aes(xintercept = t),
             linetype = "dashed", linewidth = 0.7) +
  annotate("text",
           x = 5,
           y = 0.335,
           label = "Guideline update",
           vjust = 1,
           size = 3.5) +
  facet_wrap(~ mech, ncol = 1) +
  scale_x_continuous(
    breaks = c(0, 5, 10),
    labels = c("Pre", "Implementation", "Post")
  ) +
  scale_y_continuous(limits = c(0.16, 0.34)) +
  labs(
    title = "Conceptual Model of Guideline Implementation Effects in ITS",
    subtitle = "Stylized patterns illustrating immediate level shifts vs progressive adoption vs incremental refinement",
    x = NULL,
    y = "Outcome level (schematic)"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold"),
    strip.text = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  ) 

p_concept