## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse   = TRUE,
  comment    = "#>",
  fig.width  = 6,
  fig.height = 4,
  fig.align  = "center"
)
set.seed(1)

## ----basic--------------------------------------------------------------------
library(gkrreg)

data(belgium_calls)

fit <- gkrr(calls ~ year, data = belgium_calls, sigma_method = "s3")
print(fit)

## ----summary-sandwich---------------------------------------------------------
summary(fit)

## ----vcov---------------------------------------------------------------------
round(vcov(fit), 6)
round(vcov(fit, type = "HC3"), 6)

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

## ----neff---------------------------------------------------------------------
c(n         = nrow(stars_cyg),
  neff_sum  = sum(fit_stars$weights),
  neff_kish = sum(fit_stars$weights)^2 / sum(fit_stars$weights^2))

## ----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"))

## ----jackknife----------------------------------------------------------------
jk <- gkrr_jackknife(fit_stars)
ok <- is.finite(jk$rel_error) & !jk$fallback
summary(abs(jk$rel_error[ok]))

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

## ----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")])

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

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

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

## ----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")
}

## ----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")

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

## ----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, 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------------------------------------------------------------------
sessionInfo()