library(here)
library(tidyverse)
library(gt)
source(here("04-topics/rep-chv2011/Rcode/chv2011-data-prep.R"))
source(here("04-topics/rep-chv2011/Rcode/chv2011-mte-core.R"))
source(here("04-topics/rep-chv2011/Rcode/chv2011-gt-quarto.R"))

B <- as.integer(Sys.getenv("CHV2011_BOOT", "50"))
ana <- load_chv2011_analysis(trim = TRUE)
df <- ana$data
mte <- estimate_mte_polynomial(df)
pairs <- compute_late_pairs(mte)

boot <- chv2011_mte_bootstrap_indices(nrow(df), B = B)
raw_full <- load_chv2011_raw() |> add_chv2011_derived()

boot_diff <- matrix(NA_real_, nrow = nrow(pairs), ncol = B)
for (b in seq_len(B)) {
  tryCatch({
    dfb <- if (boot$mode == "bootdata") {
      load_chv2011_boot_sample(b, for_analysis = TRUE) |>
        chv2011_attach_propensity(trim = TRUE)
    } else {
      bootstrap_analysis_draw(raw_full, boot$idx_mat[, b])
    }
    boot_diff[, b] <- compute_late_pairs(estimate_mte_polynomial(dfb))$diff
  }, error = function(e) NULL)
}

table_data <- pairs |>
  mutate(
    se = apply(boot_diff, 1, sd, na.rm = TRUE),
    p = 2 * pnorm(-abs(diff / se)),
    cell = chv2011_coef_cell(diff, se, chv2011_stars(p))
  ) |>
  select(
    `Interval 1` = interval1,
    `Interval 2` = interval2,
    `LATE diff` = cell
  )

joint_p <- tryCatch(stats::t.test(pairs$diff)$p.value, error = function(e) NA_real_)

gt_tbl <- table_data |>
  chv2011_quarto_blank_df() |>
  gt() |>
  tab_header(
    title = "Table 4(b)",
    subtitle = sprintf(
      "Test of equality of adjacent LATEs (joint p = %.3f; bootstrap %s, B = %d)",
      joint_p, boot$mode, B
    )
  )

save(table_data, gt_tbl, joint_p, boot, file = here("04-topics/rep-chv2011/Rcode/Table_4b.RData"))