Introduction to gkrreg: Gaussian kernel robust regression

Eufrásio de Andrade Lima Neto, Marcelo Rodrigo Portela Ferreira

2026-10-14

Overview

The gkrreg package implements the Gaussian Kernel Robust Regression (GKRReg, or GKRR) method proposed by De Carvalho et al. (2017). The method fits a linear regression model by iteratively re-weighting observations using the Gaussian kernel, so that outliers and leverage points are automatically down-weighted. Convergence is theoretically guaranteed (Propositions 4.1 and 4.2 of the paper).

The package provides:

The GKRR method

Model and objective function

GKRReg minimises

\[S(\boldsymbol{\beta}) = 2\sum_{i=1}^n \bigl[1 - G(y_i,\, \hat{\mu}_i)\bigr], \qquad G(a,b) = \exp\!\bigl(-{(a-b)^2}/{\gamma^2}\bigr),\]

where \(\hat{\mu}_i = \mathbf{x}_i^\top\boldsymbol{\beta}\) and \(\gamma^2 > 0\) is the kernel width hyper-parameter. An observation with a large residual contributes close to 2 to \(S\); a perfectly fitted observation contributes 0.

Estimation algorithm

Differentiating \(S\) yields an IRWLS problem with kernel weights \(k_{ii}^{(t)} = G(y_i, \hat{\mu}_i^{(t-1)}) \in (0,1]\). The algorithm starts from OLS and alternates between WLS and weight updates until convergence (Propositions 4.1–4.2 of De Carvalho et al. (2017)).

Width hyper-parameter \(\gamma^2\)

Five options are available via the sigma_method argument:

Method Description Recommended scenario
"s1" Mean of the 0.1 and 0.9 quantiles of OLS squared residuals on a sub-sample (Caputo estimator) Clean data
"s2" Pairwise median of \((y_i - \hat{\mu}_j^{\text{OLS}})^2\), \(i \neq j\) Y-space outliers \(\leq 10\%\); X-space outliers \(\leq 15\%\)
"s3" Residual variance \(\sum(y_i - \hat{\mu}_i)^2 / (n - p - 1)\) Y-space outliers \(\geq 15\%\); leverage points
"s4" AICc-selected bandwidth \(h^2\) via sm::h.select() Large samples (\(n \geq 200\))
"auto" Selects among S1–S4 by out-of-bag bootstrap MSE When the outlier scenario is unknown

When "auto" is used, a message is emitted showing which method was selected and the comparative OOB MSEs. The selection uses \(B = 99\) replicates by default, configurable via auto_args = list(B = ...).

Basic usage

library(gkrreg)

data(belgium_calls)

fit <- gkrr(calls ~ year, data = belgium_calls, sigma_method = "s3")
print(fit)
#> 
#> GKRReg -- Gaussian Kernel Robust Regression
#> --------------------------------------------
#> gamma^2 : 31.61  (method: s3)
#> Iterations: 9  [converged]
#> 
#> Coefficients:
#> (Intercept)        year 
#>     -6.5764      0.1351 
#> 
#> (Sandwich inference available -- use summary() for the full table)

summary() uses the sandwich variance estimator by default, producing standard errors, confidence intervals and Wald z-test p-values with no additional computation:

summary(fit)
#> 
#> Call:
#> gkrr(formula = calls ~ year, data = belgium_calls, sigma_method = "s3")
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -0.6165 -0.1917 -0.0271  3.5213 18.4537 
#> 
#> Coefficients:
#>             Estimate Std. Error CI 95% lower CI 95% upper   p-value    
#> (Intercept)  -6.5764     1.1145      -8.7609      -4.3920 3.621e-09 ***
#> year          0.1351     0.0203       0.0954       0.1748 2.529e-11 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> (Sandwich SE (asymptotic Wald z-test))
#> 
#> gamma^2: 31.61  |  method: s3  |  iterations: 9  |  converged: TRUE
#> R-squared: -0.1219  |  Weighted R-squared: 0.5499
#> Kernel weights -- min: 0.0000  mean: 0.7501  max: 1.0000
#> Note: sandwich inference may be less reliable here (n = 24 (small sample); 12% of observations received near-zero kernel weight (< 0.01)).
#>   Consider bootstrap inference via boot = TRUE in gkrr() or
#>   summary(fit, boot = gkrr_boot(fit)).

The sandwich covariance matrix is accessible via vcov(), which also takes a type argument selecting a finite-sample correction (see Finite-sample corrections below); the default type = "HC0" returns the matrix computed during fitting:

round(vcov(fit), 6)
#>             (Intercept)      year
#> (Intercept)    1.242183 -0.022546
#> year          -0.022546  0.000410
round(vcov(fit, type = "HC3"), 6)
#>             (Intercept)      year
#> (Intercept)    1.805784 -0.032744
#> year          -0.032744  0.000595
#> attr(,"type")
#> [1] "HC3"
#> attr(,"leverage")
#> [1] "bread"
#> attr(,"neff")
#> [1] "sum"
#> attr(,"neff_value")
#> [1] 18.00199
#> attr(,"n_wneg")
#> [1] 6
#> attr(,"n_fallback")
#> [1] 6

Statistical inference

Sandwich variance estimator (default)

The original paper does not provide a sampling distribution for \(\hat{\boldsymbol{\beta}}\). The WLS covariance matrix \((X^\top \hat{K} X)^{-1}\) from the final IRWLS step treats \(\hat{K}\) as fixed and therefore underestimates the true variance.

gkrreg implements an analytic sandwich variance estimator based on the theory of generalised \(M\)-estimators. At convergence the GKRReg estimating equations are:

\[\sum_{i=1}^n k_{ii}(\boldsymbol{\beta})\,\mathbf{x}_i (y_i - \mathbf{x}_i^\top\boldsymbol{\beta}) = \mathbf{0},\]

and the asymptotic sandwich covariance is \(n^{-1}\hat{A}^{-1}\hat{B}\hat{A}^{-1}\), where

\[\hat{A} = \frac{1}{n}\,\mathbf{X}^\top \operatorname{diag}\!\left[k_{ii}\!\left(1 - \frac{2\hat{e}_i^2}{\hat{\gamma}^2}\right)\right] \mathbf{X}, \qquad \hat{B} = \frac{1}{n}\,\mathbf{X}^\top \operatorname{diag}\!\left(k_{ii}^2\,\hat{e}_i^2\right)\mathbf{X}.\]

This corresponds to the HC0 class of heteroskedasticity-robust covariance estimators (White 1982). The correction term \((1 - 2\hat{e}_i^2/\hat{\gamma}^2)\) in \(\hat{A}\) arises directly from differentiating the Gaussian kernel with respect to \(\boldsymbol{\beta}\).

When is the sandwich reliable? The sandwich is consistent and efficient asymptotically, but it has two distinct limitations, and they call for different remedies. Its finite-sample behaviour can be improved while keeping the computation analytic, which is what the HC corrections of the next subsection do. Its treatment of \(\hat{\gamma}^2\) as fixed cannot: propagating that variability requires the bootstrap. summary() emits a note when the sample is small or a large share of the observations has been down-weighted to near zero.

Finite-sample corrections: the GKRReg-HC family

The estimator above is of the HC0 class: it uses the squared residuals in \(\hat{B}\) with no finite-sample correction. In the classical linear model this is improved by the HC1–HC4 corrections (MacKinnon and White 1985; Cribari-Neto 2004), which rescale the residuals using the leverage \(h_{ii}\). gkrr_vcov() provides the analogues for GKRReg, selected through type:

fit_stars <- gkrr(log_light ~ log_temp, data = stars_cyg,
                  sigma_method = "s4")
round(sapply(c("HC0", "HC1", "HC2", "HC3", "HC4"),
             function(ty) sqrt(diag(gkrr_vcov(fit_stars, type = ty)))), 5)
#>                 HC0     HC1     HC2     HC3     HC4
#> (Intercept) 2.24600 2.41405 2.39076 2.29924 2.11415
#> log_temp    0.55623 0.59785 0.59106 0.56643 0.52452

All five reduce to their classical counterparts as \(\gamma^2 \to \infty\), the limit in which the Gaussian kernel becomes flat and GKRReg reduces to OLS. Two features of the method shape them.

The relevant sample size is the effective one

An observation down-weighted to near zero contributes neither information nor a degree of freedom, so the nominal \(n\) overstates the sample:

c(n         = nrow(stars_cyg),
  neff_sum  = sum(fit_stars$weights),
  neff_kish = sum(fit_stars$weights)^2 / sum(fit_stars$weights^2))
#>         n  neff_sum neff_kish 
#>  47.00000  14.88378  19.88266

HC1 and HC4 are calibrated by this quantity rather than by \(n\), through neff = "sum" (\(\sum_i k_{ii}\)) or Kish’s effective sample size (Kish 1965), neff = "kish", which penalises the dispersion of the weights and not only their sum.

The weighted leverage is not confined to \([0,1]\)

The leave-one-out expansion that motivates HC2–HC4 delivers the diagonal of \(\mathbf{X}(\mathbf{X}^\top\hat{W}\mathbf{X})^{-1}\mathbf{X}^\top\hat{W}\) with \(w_i = k_{ii}(1 - 2\hat{e}_i^2/\hat{\gamma}^2)\) — the weights of the bread, not the kernel weights. Because \(w_i < 0\) whenever \(|\hat{e}_i| > \hat{\gamma}/\sqrt{2}\), which is where the influence function peaks, \(\hat{W}\) need not be positive definite. This has no analogue in the classical family, where \(h_{ii} \in [0,1]\) by construction. Observations whose leverage leaves the admissible range fall back to the always-admissible kernel-weighted leverage, and both counts are reported as attributes:

V <- gkrr_vcov(fit_stars, type = "HC3")
c(n = nrow(stars_cyg),
  n_wneg     = attr(V, "n_wneg"),
  n_fallback = attr(V, "n_fallback"))
#>          n     n_wneg n_fallback 
#>         47         34         35

Negative bread weights are the normal regime of a redescending estimator rather than an edge case, so these two numbers are worth reading alongside any HC2–HC4 standard error: they bound the extent to which the correction rests on the one-step approximation. gkrr_jackknife() checks that approximation directly, by refitting the model with each observation deleted in turn:

jk <- gkrr_jackknife(fit_stars)
ok <- is.finite(jk$rel_error) & !jk$fallback
summary(abs(jk$rel_error[ok]))
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> 0.002115 0.020396 0.064049 0.224700 0.230524 1.037622

kernel_update is the most consequential choice

Because the kernel is redescending, inflating a residual reduces the weight attached to it. By default gkrr_vcov() recomputes the weight at the inflated residual, \(\tilde{k}_{ii} = \exp(-\tilde{e}_i^2/\hat{\gamma}^2)\), which keeps the meat internally consistent but also makes the correction partly self-limiting. Holding the weight at \(k_{ii}\) instead gives a markedly larger correction, and the gap widens the smaller \(\hat{\gamma}\) is relative to the residual scale — that is, exactly where the finite-sample correction is most needed:

cmp <- function(formula, data, sigma_method) {
  f <- gkrr(formula, data = data, sigma_method = sigma_method, seed = 1)
  j <- length(coef(f))                      # slope
  se <- function(...) unname(sqrt(diag(gkrr_vcov(f, ...)))[j])
  c(gamma      = sqrt(f$gamma2),
    HC0        = se("HC0"),
    HC3_update = se("HC3", kernel_update = TRUE),
    HC3_fixed  = se("HC3", kernel_update = FALSE))
}
tab <- rbind(
  belgium_calls = cmp(calls     ~ year,     belgium_calls, "s3"),
  mammals       = cmp(log_brain ~ log_body, mammals,       "s3"),
  stars_cyg     = cmp(log_light ~ log_temp, stars_cyg,     "s4"))
round(cbind(tab, ratio = tab[, "HC3_fixed"] / tab[, "HC3_update"]), 5)
#>                 gamma     HC0 HC3_update HC3_fixed   ratio
#> belgium_calls 5.62234 0.02025    0.02439   0.02485 1.01887
#> mammals       0.69429 0.01857    0.01893   0.01998 1.05579
#> stars_cyg     0.23018 0.55623    0.56643   1.01091 1.78471

Comparing every variant at once

gkrr_vcov_all() evaluates the whole family — each correction crossed with both leverages and both effective sample sizes — and returns the standard errors together with the diagnostics. It is meant for simulation studies and for inspecting a single fit:

all_v <- gkrr_vcov_all(fit_stars)
head(all_v$se, 8)
#>   type leverage neff kernel_update (Intercept)  log_temp
#> 1  HC0    bread  sum          TRUE    2.246004 0.5562296
#> 2  HC1    bread  sum          TRUE    2.414045 0.5978457
#> 3  HC2    bread  sum          TRUE    2.390764 0.5910582
#> 4  HC3    bread  sum          TRUE    2.299240 0.5664274
#> 5  HC4    bread  sum          TRUE    2.114151 0.5245157
#> 6  HC2   kernel  sum          TRUE    2.261746 0.5593232
#> 7  HC3   kernel  sum          TRUE    2.332883 0.5742915
#> 8  HC4   kernel  sum          TRUE    2.244117 0.5559112
unlist(all_v$diagnostics[c("neff_sum", "neff_kish", "n_wneg", "n_fallback")])
#>   neff_sum  neff_kish     n_wneg n_fallback 
#>   14.88378   19.88266   34.00000   35.00000

A default correction is deliberately not recommended here. HC3 is the usual choice in the classical setting, but for GKRReg the ordering of the family under contamination has yet to be established by simulation, and HC2 — which uses \((1-m_i)^{-1/2}\) in place of \((1-m_i)^{-1}\) and is correspondingly less sensitive as \(m_i \to 1\) — may prove the better default. Until that evidence is in, type = "HC0" remains the package default, so that results from earlier versions are reproduced exactly.

Bootstrap inference

gkrr_boot() runs a pairs bootstrap, re-fitting the complete GKRReg algorithm — including re-estimating \(\gamma^2\) — on each replicate. This captures all sources of variability and is more reliable than the sandwich for small samples or heavy contamination.

Three CI types are available: "percentile", "normal", and "bca" (bias-corrected and accelerated, Efron (1987); the default and recommended). Bootstrap p-values use the centred-t approach:

\[p_j = \frac{2}{B}\,\#\!\left\{|\hat\beta_j^* - \hat\beta_j| \geq |\hat\beta_j|\right\}.\]

Option A — boot = TRUE inside gkrr()

fit_b <- gkrr(calls ~ year, data = belgium_calls,
              sigma_method = "s3",
              boot      = TRUE,
              boot_args = list(B = 999, type = "bca", seed = 1))
summary(fit_b)

Option B — separate gkrr_boot() call

boot <- gkrr_boot(fit, B = 999, type = "bca", seed = 1)
summary(fit, boot = boot)
plot(boot, which = 1)   # histogram + shaded CI per coefficient
plot(boot, which = 2)   # bootstrap scatter-plot matrix

Choosing among the procedures

The analytic corrections and the bootstrap fix different deficiencies of HC0, so they are complementary rather than competing.

Analytic (HC0–HC4) Pairs bootstrap
Cost one fit \(B\) fits
Finite-sample correction explicit implicit, via resampling
Variability of \(\hat{\gamma}^2\) not propagated propagated
Selection of \(\hat{\gamma}^2\) not propagated propagated
Skewness of \(\hat{\boldsymbol{\beta}}\) symmetric Wald intervals BCa adjusts
Needs \(\hat{W}\) positive definite HC2–HC4 only no
Deterministic yes no

Every member of the analytic family is computed at a fixed \(\hat{\gamma}^2\). If estimating the kernel width contributes substantially to the variance of \(\hat{\boldsymbol{\beta}}\), all of HC0–HC4 understate it by a common factor, and no choice of correction recovers it. That is the situation the bootstrap exists for.

Two quantities available from a single fit anticipate it at no cost: the effective sample size, and the share of observations with near-zero kernel weight. summary() emits a note when \(n < 50\) or more than 10% of the observations fall below \(k_{ii} = 0.01\) — a weight that corresponds to a residual beyond \(2.15\hat{\gamma}\), a large deviation on any scale. When both procedures have been computed, it also warns if their standard errors differ by more than 25% for any coefficient; when the analytic one is the smaller of the two, the neglected variability of \(\hat{\gamma}^2\) is the likely reason.

Diagnostic plots

plot.gkrr() provides six panels. Point size is inversely proportional to \(k_{ii}\): outliers appear large and red; well-fitted observations appear small and blue.

oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
par(mfrow = c(1, 3))
plot(fit, which = 1, ask = FALSE)   # residuals vs. fitted
plot(fit, which = 3, ask = FALSE)   # weight vs. residual + kernel curve
plot(fit, which = 4, ask = FALSE)   # weight vs. index

par(oldpar)

Panel 3 overlays the theoretical curve \(G(e) = \exp(-e^2/\hat{\gamma}^2)\), making the down-weighting mechanism directly visible. Observations 15–20 in belgium_calls (a recording error) receive near-zero weights and are effectively excluded from the fit.

which Description
1 Residuals vs. fitted values
2 Observed vs. fitted (\(y\) vs. \(\hat{\mu}\))
3 Kernel weight vs. residual with theoretical curve
4 Kernel weight vs. observation index
5 Normal QQ-plot of residuals, coloured by weight
6 Convergence of \(S(\boldsymbol{\beta})\)

Comparison with other methods

data(kootenay)

fit_ols  <- lm(newgate ~ libby, data = kootenay)
fit_gkrr <- gkrr(newgate ~ libby, data = kootenay, sigma_method = "s1")

if (requireNamespace("MASS", quietly = TRUE)) {
  fit_m  <- MASS::rlm(newgate ~ libby, data = kootenay, method = "M")
  fit_mm <- MASS::rlm(newgate ~ libby, data = kootenay, method = "MM")
} else { fit_m <- fit_mm <- NULL }

tab <- rbind(OLS  = coef(fit_ols),
             GKRR = coef(fit_gkrr))
if (!is.null(fit_m))
  tab <- rbind(tab, M = coef(fit_m), MM = coef(fit_mm))
print(round(tab, 4))
#>      (Intercept)   libby
#> OLS      23.1649 -0.0138
#> GKRR      5.4561  0.6216
#> M        23.1942 -0.0186
#> MM        5.4667  0.6209

plot(kootenay$libby, kootenay$newgate,
     xlab = "Libby flow", ylab = "Newgate flow",
     main = "Kootenay River -- X-space outlier (1934)",
     pch = 19, col = "grey60")
points(kootenay["1934","libby"], kootenay["1934","newgate"],
       col = "red", pch = 17, cex = 1.6)
abline(fit_ols,  col = "black",     lwd = 2, lty = 2)
abline(fit_gkrr, col = "firebrick", lwd = 2)
if (!is.null(fit_m)) {
  abline(fit_m,  col = "darkorange", lwd = 2, lty = 3)
  abline(fit_mm, col = "purple",     lwd = 2, lty = 4)
  legend("topleft", c("OLS","GKRR","M","MM","1934"),
         col = c("black","firebrick","darkorange","purple","red"),
         lty = c(2,1,3,4,NA), pch = c(NA,NA,NA,NA,17), lwd = 2, bty = "n")
} else {
  legend("topleft", c("OLS","GKRR","1934"),
         col = c("black","firebrick","red"),
         lty = c(2,1,NA), pch = c(NA,NA,17), lwd = 2, bty = "n")
}

Real-data applications

Belgium international calls — Y-space outliers

fit_ols  <- lm(calls ~ year, data = belgium_calls)
fit_gkrr <- gkrr(calls ~ year, data = belgium_calls, sigma_method = "s3")

plot(belgium_calls$year + 1900, belgium_calls$calls,
     xlab = "Year", ylab = "Calls (tens of millions)",
     main = "Belgium International Calls",
     pch = 19, col = "grey60")
points(belgium_calls$year[15:20] + 1900, belgium_calls$calls[15:20],
       col = "red", pch = 17, cex = 1.4)
abline(fit_ols,  col = "black",     lwd = 2, lty = 2)
abline(fit_gkrr, col = "firebrick", lwd = 2)
legend("topleft", c("OLS","GKRR (S3)","Outliers (1964-69)"),
       col = c("black","firebrick","red"),
       lty = c(2,1,NA), pch = c(NA,NA,17), lwd = 2, bty = "n")

oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
par(mfrow = c(1, 3))
plot(fit_gkrr, which = 1, ask = FALSE)
plot(fit_gkrr, which = 3, ask = FALSE)
plot(fit_gkrr, which = 4, ask = FALSE)

par(oldpar)

Delivery time — leverage points in multiple regression

data(delivery)
fit_ols  <- lm(delivery_time ~ n_products + distance, data = delivery)
fit_gkrr <- gkrr(delivery_time ~ n_products + distance, data = delivery,
                 sigma_method = "s3")
round(rbind(OLS  = coef(fit_ols),
            GKRR = coef(fit_gkrr)), 4)
#>      (Intercept) n_products distance
#> OLS       2.3412     1.6159   0.0144
#> GKRR      4.0607     1.3891   0.0145
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
par(mfrow = c(1, 3))
plot(fit_gkrr, which = 2, ask = FALSE)
plot(fit_gkrr, which = 4, ask = FALSE)
plot(fit_gkrr, which = 5, ask = FALSE)

par(oldpar)

Mammals — leverage on the log scale

data(mammals)
fit_ols  <- lm(log_brain ~ log_body, data = mammals)
fit_gkrr <- gkrr(log_brain ~ log_body, data = mammals, sigma_method = "s3")

plot(mammals$log_body, mammals$log_brain,
     xlab = "log(body mass, kg)", ylab = "log(brain mass, g)",
     main = "Brain vs. Body Mass (62 Mammal Species)",
     pch = 19, col = "grey60", cex = 0.8)
elephants <- mammals$species %in% c("African elephant","Asian elephant")
points(mammals$log_body[elephants], mammals$log_brain[elephants],
       col = "red", pch = 17, cex = 1.5)
abline(fit_ols,  col = "black",     lwd = 2, lty = 2)
abline(fit_gkrr, col = "firebrick", lwd = 2)
legend("topleft", c("OLS","GKRR (S3)","Elephants"),
       col = c("black","firebrick","red"),
       lty = c(2,1,NA), pch = c(NA,NA,17), lwd = 2, bty = "n")

Session information

sessionInfo()
#> R version 4.6.0 (2026-04-24)
#> Platform: aarch64-apple-darwin23
#> Running under: macOS Tahoe 26.6.2
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: America/Fortaleza
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] gkrreg_0.5.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] digest_0.6.39     R6_2.6.1          fastmap_1.2.0     xfun_0.57        
#>  [5] cachem_1.1.0      knitr_1.51        htmltools_0.5.9   rmarkdown_2.31   
#>  [9] lifecycle_1.0.5   cli_3.6.6         sass_0.4.10       jquerylib_0.1.4  
#> [13] sm_2.2-6.0        compiler_4.6.0    rstudioapi_0.18.0 tools_4.6.0      
#> [17] evaluate_1.0.5    bslib_0.11.0      yaml_2.3.12       otel_0.2.0       
#> [21] jsonlite_2.0.0    rlang_1.2.0       MASS_7.3-65

References

Cribari-Neto, Francisco. 2004. “Asymptotic Inference Under Heteroskedasticity of Unknown Form.” Computational Statistics & Data Analysis 45 (2): 215–33. https://doi.org/10.1016/S0167-9473(02)00366-3.
De Carvalho, Francisco de A. T., Eufrásio de A. Lima Neto, and Marcelo R. P. Ferreira. 2017. “A Robust Regression Method Based on Exponential-Type Kernel Functions.” Neurocomputing 234: 58–74. https://doi.org/10.1016/j.neucom.2016.12.035.
Efron, Bradley. 1987. “Better Bootstrap Confidence Intervals.” Journal of the American Statistical Association 82 (397): 171–85. https://doi.org/10.1080/01621459.1987.10478410.
Kish, Leslie. 1965. Survey Sampling. John Wiley & Sons.
MacKinnon, James G., and Halbert White. 1985. “Some Heteroskedasticity-Consistent Covariance Matrix Estimators with Improved Finite Sample Properties.” Journal of Econometrics 29 (3): 305–25. https://doi.org/10.1016/0304-4076(85)90158-7.
White, Halbert. 1982. “Maximum Likelihood Estimation of Misspecified Models.” Econometrica 50 (1): 1–25. https://doi.org/10.2307/1912526.