# Table 4(a): Test of linearity of E(Y|X,P) - port of bootstrap_MT.R
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-table4a-core.R"))
source(here("04-topics/rep-chv2011/Rcode/chv2011-gt-quarto.R"))

data_path <- chv2011_table4a_data_path()
if (!file.exists(data_path)) {
  export_chv2011_table4a_data(data_path)
}

env <- new.env()
load(data_path, envir = env)
data <- env$data

B <- as.integer(Sys.getenv("CHV2011_BOOT_4A", "1000"))
res <- chv2011_table4a_run(data, B = B)

table_data <- tibble(
  polynomial_order = res$degrees,
  p_value_boot = res$p_boot,
  p_value_asymp_mt = res$p_asymp_mt,
  reject = if_else(res$p_boot < res$critical_rw, "Reject", "Fail to reject")
) |>
  mutate(
    p_asymp = chv2011_fmt_num(p_value_asymp_mt),
    cell = chv2011_fmt_num(p_value_boot)
  ) |>
  select(polynomial_order, p_asymp, cell, reject)

gt_tbl <- table_data |>
  chv2011_quarto_blank_df() |>
  gt() |>
  tab_header(
    title = "Table 4(a)",
    subtitle = sprintf(
      "Linearity of E(Y|X,P); Romano-Wolf critical value (10%%) = %.3f; joint: %s",
      res$critical_rw,
      if (res$reject_rw) "Reject" else "Fail to reject"
    )
  ) |>
  cols_label(
    polynomial_order = "Highest order of polynomial in P",
    p_asymp = "Asymptotic p-value (Murphy-Topel)",
    cell = "Bootstrap p-value",
    reject = "Decision (RW)"
  )

critical_rw <- res$critical_rw

save(
  table_data,
  gt_tbl,
  res,
  critical_rw,
  file = here("04-topics/rep-chv2011/Rcode/Table_4a.RData")
)