## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5
)
options(hcinfer.use_emoji = FALSE)

## ----run-hcbeta---------------------------------------------------------------
library(hcinfer)

fit <- lm(murder ~ hs_grad + poverty + single, data = Crime2009)
result <- hcinfer(fit, type = "hcbeta")

summary(result)

## ----hcbeta-method-params-----------------------------------------------------
result$method_params

## ----hcbeta-leverage-top5-----------------------------------------------------
diagnostics <- data.frame(
  state = Crime2009$state[as.integer(result$observation)],
  leverage = unname(result$leverage),
  weight = unname(result$weights),
  residual = unname(result$residuals)
)

head(diagnostics[order(-diagnostics$leverage), ], 5)

## ----hcbeta-weight-top5-------------------------------------------------------
head(diagnostics[order(-diagnostics$weight), ], 5)

## ----hcbeta-weight-plot, fig.alt = "Scatterplot of HCbeta adjustment factors against leverage values for the Crime2009 model."----
plot(vcov_hc(fit, type = "hcbeta"))

## ----hcbeta-sensitivity-------------------------------------------------------
settings <- list(
  default            = list(),
  stronger_exponent  = list(c1 = 10),
  faster_decay       = list(c2 = 1.0),
  tighter_truncation = list(lower = 0.05, upper = 0.90),
  capped_shapes      = list(a_max = 50, b_max = 50)
)

sensitivity <- lapply(names(settings), function(setting) {
  res <- do.call(hcinfer, c(list(fit, type = "hcbeta"), settings[[setting]]))
  row <- tests(res, parm = "single")
  ci <- confint(res, parm = "single")
  data.frame(
    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 <- do.call(rbind, sensitivity)
sensitivity