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"))
source(here("04-topics/rep-chv2011/Rcode/chv2011-mte-core.R"))

ana <- load_chv2011_analysis(trim = TRUE)
df <- ana$data

# Normal switching regression: separate wage equations by S (proxy for A-6).
y1 <- lm(as.formula(paste("wage ~", paste(chv2011_x, collapse = " + "))), data = df |> filter(state == 1))
y0 <- lm(as.formula(paste("wage ~", paste(chv2011_x, collapse = " + "))), data = df |> filter(state == 0))

show <- c("exp", "expsq", "cafqt", "mhgc", "numsibs", "urban14", "lavlocwage17", "avurate", "lwage5", "lurate")
labels <- c(
  "Experience", "Experience Squared", "Corrected AFQT", "Mother's Years of Schooling",
  "Number of Siblings", "Urban Residence at 14", "Permanent Local Log Earnings at 17",
  "Permanent State Unemployment at 17", "Local Log Earnings 1991", "Local Unemployment 1991"
)

t1 <- broom::tidy(y1) |> filter(term %in% show)
t0 <- broom::tidy(y0) |> filter(term %in% show)

table_data <- tibble(term = show, label = labels) |>
  left_join(t1 |> rename(mu1 = estimate, se1 = std.error), by = c("term" = "term")) |>
  left_join(t0 |> rename(mu0 = estimate, se0 = std.error), by = c("term" = "term")) |>
  mutate(
    diff = mu1 - mu0,
    cell_mu1 = chv2011_coef_cell(mu1, se1),
    cell_mu0 = chv2011_coef_cell(mu0, se0),
    cell_diff = chv2011_coef_cell(diff, sqrt(se1^2 + se0^2))
  ) |>
  select(label, cell_mu1, cell_mu0, cell_diff)

gt_tbl <- table_data |>
  chv2011_quarto_blank_df() |>
  gt() |>
  tab_header(
    title = "Table A-6",
    subtitle = "Maximum likelihood estimates, normal switching regression (OLS-by-S proxy)"
  ) |>
  cols_label(
    label = "Variable",
    cell_mu1 = "mu_1(X)",
    cell_mu0 = "mu_0(X)",
    cell_diff = "mu_1(X) - mu_0(X)"
  )

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