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

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

This vignette shows how to run HCbeta on the `Crime2009` dataset, the 2009 U.S.
crime application from the HCbeta paper, which regresses the murder rate on
high-school graduation, poverty, and single-parent household rates. It inspects
HCbeta's method parameters and diagnostic quantities, and runs a small
sensitivity check on its tuning controls. It assumes the introduction
(`vignette("introduction", package = "hcinfer")`).

## Run HCbeta

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

```{r run-hcbeta}
library(hcinfer)

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

summary(result)
```

## Inspect HCbeta parameters

HCbeta stores its six user-facing controls and eight estimated quantities in
`method_params`. The adjustable controls are `c1`, `c2`, `lower`, `upper`,
`a_max`, and `b_max`. The remaining entries are computed from the fitted design
by method of moments and shrinkage.

```{r hcbeta-method-params}
result$method_params
```

The table below maps every printed name to its role and, where applicable, to
the corresponding mathematical symbol.

| Entry | Symbol | Meaning |
|-------|--------|---------|
| `c1` | $c_1$ | Exponent constant (default 7) |
| `c2` | $c_2$ | Exponent decay rate (default 0.75) |
| `lower` | | Lower truncation limit for $w_t$ (default 0.01) |
| `upper` | | Upper truncation limit for $w_t$ (default 0.99) |
| `a_max` | $A_{\max}$ | Upper cap for $\tilde a$ (default 10000, valid range $[50,\;25000]$) |
| `b_max` | $B_{\max}$ | Upper cap for $\tilde b$ (default 10000, valid range $[50,\;25000]$) |
| `mu_hat` | $\hat\mu$ | Mean of the truncated leverage complements $w_t$ |
| `s2_w` | $s_w^2$ | Variance of $w_t$ |
| `phi_hat` | $\hat\phi$ | Estimated dispersion |
| `a_hat` | $\hat a$ | Raw moment shape for the Beta family |
| `b_hat` | $\hat b$ | Raw moment shape for the Beta family |
| `zeta` | $\zeta$ | Shrinkage weight toward $a = b = 1$ |
| `a_tilde` | $\tilde a$ | Adjusted shape after shrinkage and floor |
| `b_tilde` | $\tilde b$ | Adjusted shape after shrinkage and floor |

A fixed shape floor $\varepsilon = 0.01$ is applied via $\max(\cdot,\;\varepsilon)$
after shrinkage and before the $\min(\cdot,\;A_{\max})$ caps. It is part of the
HCbeta definition, not a user argument, and is not accepted through `...`. It is
distinct from `lower`, which truncates the leverage complements $w_t$ before the
Beta CDF is evaluated.

## 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 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)
```

You can also sort by robust weight to see which observations contribute most to
the variance estimate.

```{r hcbeta-weight-top5}
head(diagnostics[order(-diagnostics$weight), ], 5)
```

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

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

## Run a sensitivity check

All six adjustable HCbeta controls can be passed through `...`. The sensitivity
check below compares the default result with a small set of alternative
settings, each varying only declared controls so that the interpretation
remains tied to HCbeta. For every setting it reports the robust standard error,
p-value, and confidence interval for the focus coefficient `single`, together
with the largest adjustment factor.

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

In this model the exponent constants `c1` and `c2` and the truncation window
drive the HCbeta correction: raising `c1` to 10 pushes the robust SE upward,
while increasing `c2` to 1.0 pulls it downward. Tightening the truncation
bounds to $[0.05,\;0.90]$ also changes the result by restricting the range of
leverage complements fed to the Beta CDF. In contrast, the shape caps `a_max`
and `b_max` do not change the output when set to their minimum admissible value
of 50, because the adjusted shapes $\tilde a \approx 3.4$ and $\tilde b \approx
0.7$ sit far below that floor. This illustrates a practical guardrail: when the
design lacks extreme leverage complements, the caps remain inactive and the
inference is driven primarily by the exponent and truncation settings.