---
title: "Using HCbeta"
vignette: >
  %\VignetteIndexEntry{Using HCbeta}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

HCbeta is the default estimator in `hcinfer()`. This vignette shows how to use
it, inspect the stored quantities, and run small sensitivity checks with its
method-specific arguments.

## Run HCbeta

Fit an OLS model and request HCbeta explicitly. This is equivalent to the
default `hcinfer(fit)` call.

```{r}
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)
```

## Extract coefficient test results

Use `tests()` to extract the coefficient-level Wald results as a tibble.

```{r}
tests(result)
```

For example, extract the robust standard error and p-value for the quadratic
income term.

```{r}
dplyr::select(tests(result, parm = "income_scaled_sq"), term, std_error, p_value)
```

## Inspect HCbeta parameters

HCbeta stores both user-facing constants and estimated quantities in
`method_params`.

```{r}
result$method_params
```

The most useful entries for routine inspection are:

| Entry | Meaning |
|-------|---------|
| `c1`, `c2` | Constants controlling the HCbeta exponent. |
| `lower`, `upper` | Truncation limits for leverage complements. |
| `a_tilde`, `b_tilde` | Adjusted Beta shape parameters. |
| `zeta` | Shrinkage weight used in the Beta parameter adjustment. |

These values are diagnostics. In routine use, the defaults should usually be
left unchanged.

## Inspect leverage and weights

HCbeta, like the other estimators, stores leverage values and robust weights.
This table shows the observations with the largest leverages.

```{r}
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)
```

You can also sort by robust weight.

```{r}
diagnostics |>
  dplyr::arrange(dplyr::desc(weight)) |>
  dplyr::slice_head(n = 5)
```

The covariance object can be plotted directly to display adjustment factors
against leverages.

```{r, 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"))
```

## Use the covariance-only interface

Use `vcov_hc()` when you only need the HCbeta covariance matrix and diagnostics.

```{r}
hcbeta_cov <- vcov_hc(fit, type = "hcbeta")
hcbeta_cov
vcov(hcbeta_cov)
```

The covariance object stores the same method parameters and diagnostics.

```{r}
hcbeta_cov$method_params
summary(hcbeta_cov)
```

## Run a sensitivity check

The HCbeta constants can be passed through `...`. A sensitivity check compares
the default result with a small set of alternative settings.

```{r}
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
```

This table helps identify whether the reported inference is sensitive to small
changes in HCbeta tuning constants. The defaults remain the recommended starting
point.

## Compare HCbeta with one classical estimator

For a quick comparison, put HCbeta next to a classical estimator such as HC3.

```{r}
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
```

## Practical workflow

For routine analyses:

1. Fit an OLS model with `lm()`.
2. Run `hcinfer(fit)` or `hcinfer(fit, type = "hcbeta")`.
3. Use `summary(result)` for readable output.
4. Use `tests(result)`, `confint(result)`, `coef(result)`, and `vcov(result)` for extracted values.
5. Inspect `result$weights`, `result$leverage`, and `result$method_params` when diagnostics are needed.

Use `vignette("hcinfer-comparison")` when you want to compare HCbeta with other
estimators implemented in the package.