---
title: "Introduction to hcinfer"
vignette: >
  %\VignetteIndexEntry{Introduction to hcinfer}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

This vignette gives a compact overview of `hcinfer`; no prior reading is
needed.

## Data and model

The examples use the `PublicSchools` data included with the package. The model
is the quadratic income model used in the HCbeta application; the `lm()` call
omits observations with missing model variables (Wisconsin expenditure is
missing, so the model fits 50 of the 51 observations).

```{r}
library(hcinfer)

schools <- PublicSchools
schools$income_scaled <- schools$income / 10000
schools$income_scaled_sq <- schools$income_scaled^2

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

## Available estimators

Use `hc_methods()` to list the estimators and their default arguments.

```{r}
hc_methods()
```

## Robust inference

Use `hcinfer()` to perform heteroskedasticity-robust inference for an `lm()`
object. The default estimator is HCbeta.

```{r}
result <- hcinfer(fit)
result
```

Use `summary()` for a readable report with tests on regression coefficients,
confidence intervals, leverage diagnostics, and heteroskedasticity-robust
diagnostics.

```{r}
summary(result)
```

You can choose another covariance matrix estimator with `type`.

```{r}
result_hc3 <- hcinfer(fit, type = "hc3")
result_hc3
```

## Test extraction

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

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

Select a single coefficient by name or position.

```{r}
tests(result, parm = "income_scaled_sq")
tests(result, parm = 3)
```

Change the significance level (`alpha`) to update the rejection decision. The
test statistic and p-value are not recomputed.

```{r}
tests(result, alpha = 0.10)
```

## Confidence interval extraction

Use `confint()` to extract heteroskedasticity-robust confidence intervals.

```{r}
confint(result)
```

You can select a coefficient and change the confidence level.

```{r}
confint(result, parm = "income_scaled_sq")
confint(result, parm = "income_scaled_sq", level = 0.90)
```

## Coefficients and covariance matrices

The `coef()` method returns the OLS estimates stored in `result`.

```{r}
coef(result)
```

The `vcov()` method extracts the robust covariance matrix stored in the
`hcinfer` object.

```{r}
robust_vcov <- vcov(result)
robust_vcov
```

Heteroskedasticity-robust standard errors are the square roots of the diagonal
entries.

```{r}
sqrt(diag(robust_vcov))
```

## Covariance-only workflow

Use `vcov_hc()` when you only need the heteroskedasticity-robust covariance
matrix and diagnostics.

```{r}
cov_hcbeta <- vcov_hc(fit)
cov_hcbeta
```

The same `summary()` and `vcov()` generics work for covariance objects.

```{r}
summary(cov_hcbeta)
vcov(cov_hcbeta)
```

You can also choose another estimator and pass its method constants.

```{r}
cov_hc5 <- vcov_hc(fit, type = "hc5", k = 0.7)
cov_hc5
```

## Plots

Use `plot()` on an `hcinfer` object to display robust confidence intervals.

```{r intro-ci-plot, fig.alt = "Robust confidence intervals for the public-schools regression coefficients."}
plot(result)
```

Select one coefficient with `parm`.

```{r intro-ci-single-plot, fig.alt = "Robust confidence interval for the quadratic income coefficient."}
plot(result, parm = "income_scaled_sq")
```

Use `plot()` on a `vcov_hc()` object to show leverage values and HC adjustment
factors.

```{r intro-vcov-plot, fig.alt = "HCbeta adjustment factors plotted against leverage values for the public-schools regression."}
plot(cov_hcbeta)
```

Set `label_top` to control how many observations with the largest adjustment
factors are labeled.

```{r intro-hc3-plot, fig.alt = "HC3 adjustment factors plotted against leverage values with the two largest weights labeled."}
plot(vcov_hc(fit, type = "hc3"), label_top = 2)
```

## Typical workflow

The routine analysis fits a model, passes it to `hcinfer()`, and inspects the
results with `summary()`, `tests()`, `confint()`, and `plot()`.

```{r, eval = FALSE}
fit <- lm(y ~ x1 + x2, data = data)
result <- hcinfer(fit, type = "hcbeta")

summary(result)
tests(result)
confint(result)
plot(result)
```