--- title: "Introduction to gkrreg: Gaussian kernel robust regression" author: "Eufrásio de Andrade Lima Neto, Marcelo Rodrigo Portela Ferreira" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Introduction to gkrreg: Gaussian kernel robust regression} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} bibliography: references.bib --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4, fig.align = "center" ) set.seed(1) ``` ## Overview The **gkrreg** package implements the *Gaussian Kernel Robust Regression* (GKRReg, or GKRR) method proposed by @decarvalho2017. 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: - `gkrr()` — model fitting with four width hyper-parameter estimators (`"s1"`, `"s2"`, `"s3"`, `"s4"`) plus automatic data-driven selection (`"auto"`); - `gkrr_boot()` — bootstrap inference (standard errors, confidence intervals and p-values); - `summary.gkrr()` — a coefficient table modelled after `summary.lm()`, with sandwich inference by default and bootstrap inference when available; - `plot.gkrr()` — six diagnostic panels emphasising kernel weights; - Six real datasets from the robust regression literature. --- ## 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 @decarvalho2017). ### 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 ```{r basic} library(gkrreg) data(belgium_calls) fit <- gkrr(calls ~ year, data = belgium_calls, sigma_method = "s3") print(fit) ``` `summary()` uses the **sandwich variance estimator** by default, producing standard errors, confidence intervals and Wald z-test p-values with no additional computation: ```{r summary-sandwich} summary(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: ```{r vcov} round(vcov(fit), 6) round(vcov(fit, type = "HC3"), 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 [@white1982]. 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 [@mackinnon1985; @cribarineto2004], which rescale the residuals using the leverage $h_{ii}$. `gkrr_vcov()` provides the analogues for GKRReg, selected through `type`: ```{r hc-family} 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) ``` 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: ```{r neff} c(n = nrow(stars_cyg), neff_sum = sum(fit_stars$weights), neff_kish = sum(fit_stars$weights)^2 / sum(fit_stars$weights^2)) ``` 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 [@kish1965], `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: ```{r leverage-diagnostics} V <- gkrr_vcov(fit_stars, type = "HC3") c(n = nrow(stars_cyg), n_wneg = attr(V, "n_wneg"), n_fallback = attr(V, "n_fallback")) ``` 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: ```{r jackknife} jk <- gkrr_jackknife(fit_stars) ok <- is.finite(jk$rel_error) & !jk$fallback summary(abs(jk$rel_error[ok])) ``` #### `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: ```{r kernel-update} 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) ``` #### 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: ```{r vcov-all} all_v <- gkrr_vcov_all(fit_stars) head(all_v$se, 8) unlist(all_v$diagnostics[c("neff_sum", "neff_kish", "n_wneg", "n_fallback")]) ``` 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, @efron1987; 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()`** ```{r boot-true, eval = FALSE} 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** ```{r boot-separate, eval = FALSE} 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. ```{r plots, fig.width = 7, fig.height = 3} 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 ```{r comparison, message = FALSE, fig.width = 6, fig.height = 4} 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)) 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 ```{r belgium, fig.width = 6, fig.height = 4} 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") ``` ```{r belgium-diag, fig.width = 7, fig.height = 3} 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 ```{r delivery} 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) ``` ```{r delivery-diag, fig.width = 7, fig.height = 3} 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 ```{r mammals, fig.width = 6, fig.height = 4} 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 ```{r session} sessionInfo() ``` ## References