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

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

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

fit <- lm(expenditure ~ income_scaled + income_scaled_sq, data = schools)
result <- hcinfer(fit, type = "hcbeta")

summary(result)

## -----------------------------------------------------------------------------
tests(result)

## -----------------------------------------------------------------------------
dplyr::select(tests(result, parm = "income_scaled_sq"), term, std_error, p_value)

## -----------------------------------------------------------------------------
result$method_params

## -----------------------------------------------------------------------------
diagnostics <- tibble::tibble(
  observation = result$observation,
  leverage = result$leverage,
  weight = result$weights,
  residual = result$residuals
) |>
  dplyr::mutate(
    state = schools$state[as.integer(observation)],
    .before = leverage
  )

diagnostics |>
  dplyr::arrange(dplyr::desc(leverage)) |>
  dplyr::slice_head(n = 5)

## -----------------------------------------------------------------------------
diagnostics |>
  dplyr::arrange(dplyr::desc(weight)) |>
  dplyr::slice_head(n = 5)

## ----fig.width = 7, fig.height = 4, fig.alt = "Scatterplot of HCbeta adjustment factors against leverage values for the public-schools model."----
plot(vcov_hc(fit, type = "hcbeta"))

## -----------------------------------------------------------------------------
hcbeta_cov <- vcov_hc(fit, type = "hcbeta")
hcbeta_cov
vcov(hcbeta_cov)

## -----------------------------------------------------------------------------
hcbeta_cov$method_params
summary(hcbeta_cov)

## -----------------------------------------------------------------------------
sensitivity_results <- list(
  default = hcinfer(fit, type = "hcbeta"),
  c1_5 = hcinfer(fit, type = "hcbeta", c1 = 5),
  tighter_truncation = hcinfer(fit, type = "hcbeta", lower = 0.02, upper = 0.98)
)

sensitivity <- purrr::imap(sensitivity_results, \(res, setting) {
  row <- tests(res, parm = "income_scaled_sq")
  ci <- confint(res, parm = "income_scaled_sq")

  tibble::tibble(
    setting = setting,
    std_error = row$std_error,
    p_value = row$p_value,
    conf_low = ci$conf_low,
    conf_high = ci$conf_high,
    max_weight = max(res$weights)
  )
})
sensitivity <- dplyr::bind_rows(sensitivity)
sensitivity

## -----------------------------------------------------------------------------
hc3 <- hcinfer(fit, type = "hc3")

extract_income <- function(res, method) {
  row <- tests(res, parm = "income_scaled_sq")
  ci <- confint(res, parm = "income_scaled_sq")
  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
  )
}

compare_hc3 <- dplyr::bind_rows(
  extract_income(result, "hcbeta"),
  extract_income(hc3, "hc3")
)
compare_hc3