---
title: "Pairs bootstrap inference with hcinfer"
vignette: >
  %\VignetteIndexEntry{Pairs bootstrap inference with 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)
```

The pairs (case) bootstrap is a resampling method for the coefficients of an
ordinary least squares model. It resamples whole observations, refits the model
many times, and summarizes the resulting spread of the estimates. Because it
makes no assumption about the form of the error variance, it provides an
empirical cross-check on the analytic heteroskedasticity-consistent (HC)
standard errors that `hcinfer` computes, and this vignette shows how to run it,
read its output, and compare it with the HC estimators.

## The method

For each replicate the bootstrap draws `n` rows with replacement from the
original data and refits OLS on that resample. Writing the estimate on replicate
`r` as the vector of refitted coefficients, the bootstrap standard error of a
coefficient is the standard deviation of its replicate values. `boot_pairs()`
offers three interval types: `"percentile"` uses the empirical quantiles of the
replicates, `"basic"` (reverse percentile) reflects those quantiles about the
original estimate, and `"normal"` uses the estimate plus or minus a normal
quantile times the bootstrap standard error.

## Data and model

The model below uses the `PublicSchools2` data, which has complete
observations for all 51 states, and includes a `south` indicator that the
original `PublicSchools` data lack.

```{r}
library(hcinfer)

schools <- PublicSchools2
schools$income_scaled <- schools$income / 10000

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

## Running the bootstrap

A single call fits the resamples and stores the estimates, standard errors,
bias, and intervals. Supplying `seed` makes the run reproducible.

```{r}
boot <- boot_pairs(fit, B = 2000, seed = 123)
boot
```

## Extractors

The `coef()`, `vcov()`, and `confint()` methods pull out the pieces you need.
`coef()` returns the original OLS estimates, `vcov()` the bootstrap covariance
matrix, and `confint()` the interval table. Note the argument-name asymmetry:
`boot_pairs()` sets the default interval type via `ci_type`, while
`confint()` overrides it via `type`.

```{r}
coef(boot)
vcov(boot)
confint(boot)
```

`confint()` can recompute intervals at a different level or type directly from
the stored replicates, without rerunning the bootstrap.

```{r}
confint(boot, level = 0.99, type = "basic")
confint(boot, parm = "south", level = 0.90)
```

## Visualizing the intervals

`plot()` draws each coefficient as its estimate with its bootstrap interval,
colored by whether the interval excludes zero.

```{r bootstrap-ci-plot, fig.alt = "Pairs bootstrap confidence intervals for the public-schools regression coefficients."}
plot(boot)
```

## An empirical reference for HC standard errors

Because the pairs bootstrap assumes nothing about the error variance, its
standard errors are a useful cross-check on the analytic HC estimators. The
following table places the OLS, bootstrap, HCbeta, and HC3 standard errors side
by side.

```{r}
data.frame(
  term = boot$table$term,
  ols = sqrt(diag(vcov(fit))),
  bootstrap = boot$table$std_error,
  hcbeta = sqrt(diag(vcov(hcinfer(fit, type = "hcbeta")))),
  hc3 = sqrt(diag(vcov(hcinfer(fit, type = "hc3"))))
)
```

The bootstrap and HC standard errors should broadly agree; large disagreements
are worth investigating, often at high-leverage points.

## Reproducibility

With a fixed `seed`, two runs are identical, and the call restores the caller's
random-number stream so it does not disturb a surrounding analysis.

```{r}
a <- boot_pairs(fit, B = 1000, seed = 7)
b <- boot_pairs(fit, B = 1000, seed = 7)
identical(a$replicates, b$replicates)
```

## Running in parallel

For large `B` or large `n`, set `cores` to `2` or more to distribute the
replicate fits across worker processes (requires the mirai and carrier
packages). The numeric result is identical to a sequential run with the same
seed; parallelism only changes the speed. The default `cores = 1` runs
sequentially.

```{r, eval = FALSE}
boot_pairs(fit, B = 10000, cores = 4, seed = 1)
```

## Practical guidance

Use a few thousand replicates for stable standard errors and more for stable
tail quantiles of the percentile and basic intervals. Choose `ci_type` to match
your needs: percentile and basic intervals adapt to skewness in the replicate
distribution, while normal intervals are symmetric. If a resample is rank
deficient it is dropped with a warning, and the summaries use the remaining
replicates.