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-gt-quarto.R"))

B_BOOT <- as.integer(Sys.getenv("CHV2011_BOOT", "250"))

raw <- load_chv2011_raw() |> add_chv2011_derived()
fit <- glm(choice_formula("baseline"), data = raw, family = binomial(link = "logit"))
avder <- average_marginal_derivatives(fit, raw)
se <- bootstrap_avg_derivatives(raw, B = B_BOOT)

show_vars <- c(
  "cafqt", "mhgc", "numsibs", "urban14", "lavlocwage17", "avurate",
  "pub4", "lwage5_17", "lurate_17", "tuition"
)
labels <- c(
  "Corrected AFQT", "Mother's Years of Schooling", "Number of Siblings",
  "Urban Residence at 14", "\"Permanent\" Local Log Earnings at 17",
  "\"Permanent\" State Unemployment Rate at 17",
  "Presence of a College at 14", "Local Log Earnings at 17",
  "Local Unemployment Rate at 17 (in %)",
  "Tuition in 4 Year Public Colleges at 17 (in $100)"
)

coef_tidy <- broom::tidy(fit)

table_data <- tibble(variable = show_vars, label = labels) |>
  left_join(avder, by = "variable") |>
  left_join(coef_tidy, by = c("variable" = "term")) |>
  mutate(
    se = se[match(variable, avder$variable)],
    coef_cell = if_else(
      !is.na(estimate),
      chv2011_coef_cell(estimate, std.error, chv2011_stars(p.value)),
      " "
    ),
    deriv_cell = chv2011_coef_cell(
      avg_deriv, se, chv2011_stars(2 * pnorm(-abs(avg_deriv / se)))
    )
  ) |>
  select(label, coef_cell, deriv_cell)

gt_tbl <- table_data |>
  chv2011_quarto_blank_df() |>
  gt() |>
  tab_header(
    title = "Table A-4",
    subtitle = "College decision model — coefficients and average derivatives"
  ) |>
  cols_label(
    label = " ",
    coef_cell = "Coefficient",
    deriv_cell = "Average Derivative"
  )

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