# Pairs bootstrap inference with hcinfer

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.

    library(hcinfer)

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

    fit <- lm(expenditure ~ income_scaled + south, data = schools)
    fit
    #>
    #> Call:
    #> lm(formula = expenditure ~ income_scaled + south, data = schools)
    #>
    #> Coefficients:
    #>   (Intercept)  income_scaled          south
    #>         -3111           4763          -1112

## Running the bootstrap

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

    boot <- boot_pairs(fit, B = 2000, seed = 123)
    boot
    #>
    #> ── Pairs bootstrap inference ───────────────────────────────────────────────────
    #> Model: `expenditure ~ income_scaled + south`
    #> Observations: 51 | Parameters: 3
    #> Replicates: 2000 of 2000 valid | Interval: percentile at 95.0%
    #> Execution: sequential | Seed: 123
    #> # A tibble: 3 × 5
    #>   term          estimate bias   boot_se ci
    #>   <chr>         <chr>    <chr>  <chr>   <chr>
    #> 1 (Intercept)   -3111    -461.9 2801    [-1.041e+04, 979.8]
    #> 2 income_scaled 4763     103    627.5   [3865, 6415]
    #> 3 south         -1112    52.32  902.2   [-2851, 764.1]

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

    coef(boot)
    #>   (Intercept) income_scaled         south
    #>     -3110.908      4762.569     -1112.153
    vcov(boot)
    #>               (Intercept) income_scaled      south
    #> (Intercept)       7845298    -1711565.8 -1140786.7
    #> income_scaled    -1711566      393699.5   168161.0
    #> south            -1140787      168161.0   813921.3
    confint(boot)
    #> # A tibble: 3 × 4
    #>   term          conf_low conf_high level
    #>   <chr>            <dbl>     <dbl> <dbl>
    #> 1 (Intercept)    -10413.      980.  0.95
    #> 2 income_scaled    3865.     6415.  0.95
    #> 3 south           -2851.      764.  0.95

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

    confint(boot, level = 0.99, type = "basic")
    #> # A tibble: 3 × 4
    #>   term          conf_low conf_high level
    #>   <chr>            <dbl>     <dbl> <dbl>
    #> 1 (Intercept)     -9108.     8298.  0.99
    #> 2 income_scaled    2206.     6102.  0.99
    #> 3 south           -3858.      984.  0.99
    confint(boot, parm = "south", level = 0.90)
    #> # A tibble: 1 × 4
    #>   term  conf_low conf_high level
    #>   <chr>    <dbl>     <dbl> <dbl>
    #> 1 south   -2525.      446.   0.9

## Visualizing the intervals

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

    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.

    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"))))
    )
    #>                        term       ols bootstrap    hcbeta       hc3
    #> (Intercept)     (Intercept) 2932.8749 2800.9458 2431.8340 2230.6434
    #> income_scaled income_scaled  636.0126  627.4548  543.1465  498.2300
    #> south                 south 1030.4299  902.1759  943.9061  887.2695

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.

    a <- boot_pairs(fit, B = 1000, seed = 7)
    b <- boot_pairs(fit, B = 1000, seed = 7)
    identical(a$replicates, b$replicates)
    #> [1] TRUE

## 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.

    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.