R/SLIC.R

In SLIC: LIC for Distributed Skewed Regression

#' SLIC function based on LIC with skewed error distributions
#'
#' The SLIC function extends the LIC method by assuming that the error term
#' follows a skewed distribution (Skew-Normal, Skew-t, or Skew-Laplace),
#' thereby improving the length and information optimisation criterion.
#'
#' @encoding UTF-8
#' @param X is a design matrix
#' @param Y is a random response vector of observed values
#' @param alpha is the significance level
#' @param K is the number of subsets
#' @param nk is the sample size of subsets
#' @param dist_type is the type of skewed error distribution:
#'   "skew_normal", "skew_t", or "skew_laplace"
#'
#' @return MUopt, Bopt, MAEMUopt, MSEMUopt, opt, Yopt
#' @export
#'
#' @examples
#' set.seed(123)
#' n <- 1000
#' p <- 5
#' X <- matrix(rnorm(n * p), ncol = p)
#' beta <- runif(p, 1, 2)
#' e <- sn::rsn(n = n, xi = 0, omega = 1, alpha = 5)
#' Y <- X %*% beta + e
#' SLIC(X, Y, alpha = 0.05, K = 10, dist_type = "skew_normal")
#'
#' @importFrom LaplacesDemon rslaplace
#' @importFrom LaplacesDemon rst
#' @importFrom sn rsn
#' @importFrom stats qt
#' @importFrom stats runif
SLIC <- function(X, Y, alpha = 0.05, K = 10, nk = NULL, dist_type = "skew_normal") {
  if (is.null(nk)) nk <- nrow(X) / K
  n <- nrow(X); p <- ncol(X)

  N <- L1 <- rep(NA, K)
  Rm <- matrix(0, nrow = nk, ncol = K)
  mr <- matrix(0, nrow = K, ncol = nk)

  for (i in 1:K) {
    mr[i, ] <- sample(1:n, nk, replace = FALSE)
    r <- matrix(c(1:nk, mr[i, ]), ncol = nk, byrow = TRUE)
    Rm[, i] <- r[2, ]
    R <- matrix(0, nrow = nk, ncol = n)
    R[t(r)] <- 1

    X1 <- R %*% X
    Y1 <- R %*% Y

    Hr <- X1 %*% solve(crossprod(X1)) %*% t(X1)
    I1 <- diag(rep(1, nk))
    SY <- sqrt(t(Y1) %*% (I1 - Hr) %*% Y1) / (nk - p)
    C1 <- sum(diag(X1 %*% solve(crossprod(X1)) %*% t(X1))) / nk
    L1[i] <- 2 * SY * C1 * qt(1 - alpha / 2, nk - p)
    N[i] <- det(t(X1) %*% X1)
  }

  opt1 <- Rm[, which.min(L1)]
  opt2 <- Rm[, which.max(N)]
  opt <- intersect(opt1, opt2)
  Yopt <- Y[opt]
  Xopt <- X[opt, , drop = FALSE]
  Bopt <- solve(crossprod(Xopt)) %*% t(Xopt) %*% Yopt
  MUopt <- Xopt %*% Bopt
  Nopt <- length(Yopt)
  E5 <- (t(Yopt - MUopt) %*% (Yopt - MUopt)) / Nopt
  A5 <- sum(abs(Yopt - MUopt)) / Nopt

  return(list(
    MUopt = MUopt, Bopt = Bopt, MAEMUopt = A5,
    MSEMUopt = E5, opt = opt, Yopt = Yopt
  ))
}