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", "50"))

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

# Robinson (1988) partial linear: Y = X'b0 + g(P) + e; report average derivatives w.r.t. X.
rob <- function(data) {
  p <- data$phat
  y <- data$wage
  x <- model.matrix(as.formula(paste("~", paste(chv2011_x, collapse = " + "))), data = data)
  y_res <- y - stats::loess(y ~ p, span = 0.75)$fitted
  fit <- lm(y_res ~ x - 1)
  coef(fit)
}

est <- rob(df)
boot <- replicate(B_BOOT, {
  idx <- sample.int(nrow(df), nrow(df), replace = TRUE)
  rob(df[idx, , drop = FALSE])
})
se <- apply(boot, 1, sd)

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

table_data <- tibble(
  term = names(est),
  label = labels[names(est)],
  estimate = as.numeric(est),
  se = se
) |>
  filter(!is.na(label)) |>
  mutate(cell = chv2011_coef_cell(estimate, se)) |>
  select(label, cell)

gt_tbl <- table_data |>
  chv2011_quarto_blank_df() |>
  gt() |>
  tab_header(
    title = "Table A-7",
    subtitle = "Average derivatives, partially linear wage model (Robinson 1988)"
  ) |>
  cols_label(label = "Variable", cell = "Average derivative (SE)")

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