R/eel.lm.R

In markean/bayesEL: Bayesian Empicial Likelihood

#' Extended empirical likelihood for linear regression
#'
#' Extended empirical likelihood for linear regression
#'
#' @param b a vector of parameters to be tested.
#' @param x a matrix or vector of data. Each row is an observation vector.
#' @param y a matrix or vector of responses(numeric).
#' @param abstol an optional value for the absolute convergence tolerance. Defaults to 1e-8.
#' @param maxit an optional value for the maximum number of iterations. Defaults to 100.
#'
#' @return An object with S3 class "eel"
#'
#' @examples
#' b <- c(1, 1)
#' x <- matrix(rnorm(100), nrow = 50)
#' y <- rnorm(50)
#' eel.lm(b, x, y)
#'
#' @import Matrix
#' @importFrom stats coef lm
#' @export

eel.lm <- function(b, x, y, abstol = 1e-8, maxit = 100) {
  ## Data
  x <- as.matrix(x)
  if (length(b) != ncol(x))
    stop("length(b) != ncol(x)")
  if (rankMatrix(x) != ncol(x))
    stop("model matrix must have full column rank")
  colnames(x) <- NULL
  y <- drop(y)
  n <- nrow(x)
  bhat <- unname(coef(lm(y ~ ., data = data.frame(x[, 2], y))))
  # distance between xbar and theta
  distance <- sapply(dist(rbind(b, bhat), method = "euclidean"),
                     function(x) {
                       attributes(x) <- NULL
                       x
                     })
  # nothing to extend when theta == xbar
  if (distance == 0) {
    el_result <- el.lm(b, x, y, abstol, maxit)
    result <- list()
    class(result) <- "eel"
    result$logLR <- el_result$logLR
    result$logL <- el_result$logL
    result$w <- el_result$w
    result$lambda <- el_result$lambda
    result$iterations <- el_result$iterations
    return(result)
  }
  direction_unit_vec <- (b - bhat) / distance


  ## EEL
  f <- function(alpha) {
    alpha *
      (1 - el2.lm(bhat + alpha * direction_unit_vec, x, y)$logLR / n) - distance
  }
  # wrapper for f to remove Inf value
  g <- function(alpha) {
    ifelse(f(alpha) == Inf, .Machine$double.xmax, f(alpha))
  }
  # compute the expansion distance
  alpha <- uniroot(
    g,
    lower = 0,
    upper = distance,
    extendInt = "upX",
    tol = .Machine$double.eps^.25,
    check.conv = F
  )$root
  # preimage of the expansion
  b_prime <- bhat + alpha * direction_unit_vec


  ## Result
  el.result <- el.lm(b_prime, x, y, abstol, maxit)
  result <- list()
  class(result) <- "eel"
  result$logLR <- el.result$logLR
  result$logL <- el.result$logL
  result$w <- el.result$w
  result$lambda <- el.result$lambda
  result$iterations <- el.result$iterations
  result$message <- ifelse(all.equal(sum(result$w), 1) == T,
                           "preimage is in the convex hull",
                           "something wrong")
  result
}