Looverage

looverage.R

# Projection matrix updates ---
set.seed(753)
prec <- .Machine$double.eps

N <- 54 # Simulate data
x <- as.matrix(c(rnorm(N), rnorm(3, 6, 0.25), rnorm(3, 8, 0.25)))
y <- c(
  x[seq(N)] * -0.5 + rnorm(N, 0, 1),
  x[seq(N + 1, N + 3)] * 0.1 + rnorm(3, 0, 0.1),
  x[seq(N + 4, N + 6)] * 0.4 + rnorm(3, 0, 0.1)
)
m <- lm(y ~ x - 1)

get_H <- \(x) x %*% solve(crossprod(x), t(x))

loo_verage <- function(i, r) # Looverage of after removing r
  H[i, i] + exp(log(H[i, r]^2) - log(1 - H[r, r]))

loo_H <- function(i, j = i, r) # H[i, j] after removing r
  H[i, j] + H[i, r] * H[j, r] / (1 - H[r, r])

lko_H <- function(i, j = i, R) { # H[i, j] after removing set R
  if(length(R) == 0) return(H[i, j])
  Recall(i, j, R[-1]) +
    Recall(i, R[1], R[-1]) * Recall(j, R[1], R[-1]) /
    (1 - Recall(R[1], R[1], R[-1]))
}

H <- get_H(x) # Projection matrix

# Leave-one-out leverage in row i for removals of column j
LOO_verage <- diag(H) + t(exp(log(H^2) - log(1 - diag(H))))

# Updating formula for the inverse to get the coefficients
update_inv <- function(XX_inv, X_rm) {
  XX_inv + (XX_inv %*% crossprod(X_rm) %*% XX_inv) /
    as.numeric(1 - X_rm %*% tcrossprod(XX_inv, X_rm))
}

# Leave-one-out DFBETA
Sx_inv <- chol2inv(chol(crossprod(x)))
dfe <- resid(m) + t(H * resid(m) / (1 - diag(H))) # Updated error
LOO_dfbeta <- matrix(NA_real_, NROW(x), NROW(x))
for(j in seq_len(NROW(x))) {
  LOO_dfbeta[, j] <- (t(update_inv(Sx_inv, x[j, , drop = FALSE]) %*%
    t(x * dfe[, j])) / (1 - LOO_verage[, j]))[, 1L]
}

# Check equivalence
abs(loo_verage(1, r = 2) -
  get_H(x[-2, ])[1, 1]) < prec
abs(loo_H(1, 2, r = N) -
  get_H(x[-N, ])[1, 2]) < prec
abs(lko_H(1, 2, R = seq(3, 7)) -
  get_H(x[-seq(3, 7), ])[1, 2]) < prec
all(abs(LOO_verage[-N, N] -
  diag(get_H(x[-N, ]))) < prec)
all(abs(LOO_dfbeta[-1, 1] -
  dfbeta(lm(y[-1] ~ x[-1, ] - 1))[, 1]) < prec)
abs(coef(m)[1] - dfbeta(m)[1, 1] - LOO_dfbeta[2, 1] -
  coef(lm(y[-1:-2] ~ x[-1:-2, ] - 1))[1]) < prec