## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(hcinfer.use_emoji = TRUE)

## -----------------------------------------------------------------------------
library(hcinfer)

hc_methods()

## -----------------------------------------------------------------------------
schools <- PublicSchools |>
  dplyr::mutate(
    income_scaled = income / 10000,
    income_scaled_sq = income_scaled^2
  )

fit <- lm(expenditure ~ income_scaled + income_scaled_sq, data = schools)

## -----------------------------------------------------------------------------
methods <- c("hc0", "hc3", "hc4", "hc4m", "hcbeta")

results <- purrr::map(methods, \(method) {
  hcinfer(fit, type = method)
})
names(results) <- methods

## -----------------------------------------------------------------------------
extract_term <- function(result, method, term = "income_scaled_sq") {
  row <- tests(result, parm = term)
  ci <- confint(result, parm = term)

  tibble::tibble(
    method = method,
    estimate = row$estimate,
    std_error = row$std_error,
    p_value = row$p_value,
    conf_low = ci$conf_low,
    conf_high = ci$conf_high,
    reject = row$reject
  )
}

comparison <- purrr::imap(results, extract_term)
comparison <- dplyr::bind_rows(comparison)
comparison

## -----------------------------------------------------------------------------
comparison <- comparison |>
  dplyr::mutate(interval_width = conf_high - conf_low)

comparison

## ----fig.width = 7, fig.height = 4, fig.alt = "Bar chart comparing robust standard errors for the quadratic income term across HC estimators."----
ggplot2::ggplot(comparison, ggplot2::aes(x = method, y = std_error)) +
  ggplot2::geom_col(fill = "#305c8a") +
  ggplot2::labs(
    x = "Estimator",
    y = "Robust standard error for the quadratic income term"
  ) +
  ggplot2::theme_minimal(base_size = 12)

## -----------------------------------------------------------------------------
comparison |>
  dplyr::select(method, conf_low, conf_high, interval_width)

## -----------------------------------------------------------------------------
cov_hc3 <- vcov_hc(fit, type = "hc3")
cov_hc3
vcov(cov_hc3)

## ----fig.width = 7, fig.height = 4, fig.alt = "Scatterplot of HC3 adjustment factors against leverage values for the public-schools model."----
plot(cov_hc3)

## -----------------------------------------------------------------------------
summary(cov_hc3)

## -----------------------------------------------------------------------------
covariances <- purrr::map(methods, \(method) {
  vcov_hc(fit, type = method)
})
names(covariances) <- methods

diagnostic_comparison <- purrr::imap(covariances, \(cov, method) {
  tibble::tibble(
    method = method,
    max_leverage = max(cov$leverage),
    max_weight = max(cov$weights),
    median_weight = stats::median(cov$weights)
  )
})
diagnostic_comparison <- dplyr::bind_rows(diagnostic_comparison)
diagnostic_comparison

## ----fig.width = 8, fig.height = 6, fig.alt = "Faceted scatterplot of HC adjustment factors against leverage values for HC3, HC4, HC4m, and HCbeta."----
figure_methods <- c("hc3", "hc4", "hc4m", "hcbeta")
figure_covariances <- purrr::map(figure_methods, \(method) {
  vcov_hc(fit, type = method)
})
names(figure_covariances) <- figure_methods

weight_comparison <- purrr::imap(figure_covariances, \(cov, method) {
  tibble::tibble(
    method = cov$label,
    leverage = cov$leverage,
    weight = cov$weights,
    high_leverage = cov$leverage > 3 * cov$p / cov$n
  )
}) |>
  dplyr::bind_rows() |>
  dplyr::mutate(
    method = factor(method, levels = c("HC3", "HC4", "HC4m", "HCbeta"))
  )

ggplot2::ggplot(weight_comparison, ggplot2::aes(x = leverage, y = weight)) +
  ggplot2::geom_point(
    ggplot2::aes(color = high_leverage),
    size = 1.8,
    alpha = 0.85
  ) +
  ggplot2::facet_wrap(~method, scales = "free_y", ncol = 2) +
  ggplot2::scale_color_manual(
    values = c(`FALSE` = "#2c5f8a", `TRUE` = "#c0392b"),
    guide = "none"
  ) +
  ggplot2::labs(
    x = expression(h[t]),
    y = expression(g[t])
  ) +
  ggplot2::theme_minimal(base_size = 12) +
  ggplot2::theme(
    strip.text = ggplot2::element_text(face = "bold"),
    panel.grid.minor = ggplot2::element_blank()
  )

## -----------------------------------------------------------------------------
comparison |>
  dplyr::select(method, estimate, std_error, p_value, conf_low, conf_high)