R/PLSdeviance.R

In abnormally-distributed/cvreg: Cross Validation and Robust Estimation Utilities

#' Calculate Akaike Information Criterion for gpls models
#'
#' @param model a fitted model
#' @param correction whether or not to calculate the small sample-bias corrected AIC (AICc).
#'
#' @return numeric value
#' @export
#' @keywords internal
#' @method AIC gpls
#' @references Sugiura, N. (1978). Further analysis of the data by Akaike's information criterion and the finite corrections. Comm. Statist. A 7, 13-26.
AIC.gpls = function(model, correction = F){
  logLik = logLik.gpls(model)
  k = model$ncomp + 1
  n = length(model$y_obs)
  aic <- (-2*logLik) + (2 * k)
  if (correction){
    aic <- aic + (((2*k)*(k + 1)) / (n - k - 1))
  }
  return(aic)
}

#' Calculate Bayesian Information Criterion for gpls models
#'
#' @param model a fitted model
#' @param correction whether or not to calculate the small sample-bias corrected BIC (BICc).
#' @return numeric value
#' @export
#' @keywords internal
#' @method BIC gpls
#'
#' @references McQuarrie, A.D. (1999). A small-sample correction for the Schwarz BIC model selection criterion, Statistics & Probability Letters 44 pp.79-86. doi: 10.1016/S0167-7152(98)00294-6
#'
BIC.gpls = function(model, correction = F){
  logLik = logLik.gpls(model)
  k = model$ncomp + 1
  n = nrow(model$model.matrix)

  if (correction){
    p =  (k * log(n) * n) / (n - k - 1)
  } else {
    p = k * log(n)
  }
  (-2*logLik) + p
}

#' Calculate model log likelihood for gpls models
#'
#' @param model a fitted model
#'
#' @return numeric value
#' @export
#' @keywords internal
#' @method logLik gpls
#'
logLik.gpls <- function(model) {

  y <- model$y_obs
  mu <- model$fitted
  family = model$family
  wts <- rep(1, length(y))

  if (family == "binomial" | family == "bernoulli") {
    if (!is.numeric(y)){
      y <- as.numeric(y) - 1
      y <- y == max(y)
      y <- as.numeric(y)
    }
    logLik <- dbinom(y, size = 1, prob = mu, log = TRUE)

  } else if (family == "poisson") {

    logLik <- dpois(y, mu, log = TRUE)

  } else if (family == "gaussian") {

    logLik <- dnorm(y, mu, sd(y - mu), log = TRUE)

  } else if (family == "Gamma"){
    sigma2 <- var(y - mu)
    logLik <- dgamma(y, mu^2/sigma2, mu/sigma2, log = TRUE)

  } else if (family == "inverse.gaussian"){

    dinvgauss <- function (x, mu = 1, sigma = 1, log = FALSE){
      if (any(mu < 0))
        stop(paste("mu must be positive", "\n", ""))
      if (any(sigma < 0))
        stop(paste("sigma must be positive", "\n", ""))
      if (any(x < 0))
        stop(paste("x must be positive", "\n", ""))
      log.lik <- (-0.5 * log(2 * pi) - log(sigma) - (3/2) * log(x) -
                    ((x - mu)^2)/(2 * sigma^2 * (mu^2) * x))
      if (log == FALSE)
        fy <- exp(log.lik)
      else fy <- log.lik
      fy
    }

    logLik <- dinvgauss(y, mu, 1, log = TRUE)

  } else if (family == "negative.binomial"){
    theta = model$theta

    dnegbinom <- function (x, mu = 1, theta = 1, log = FALSE) {
      if (any(mu <= 0))
        stop(paste("mu must be greater than 0 ", "\n", ""))
      if (any(theta <= 0))
        stop(paste("theta must be greater than 0 ", "\n", ""))
      if (length(theta) > 1)
        fy <- ifelse(theta > 1e-04,
                     dnbinom(x, size = mu/theta, mu = mu, log = log),
                     dpois(x = x, lambda = mu, log = log))
      else fy <- if (theta < 1e-04)
        dpois(x = x, lambda = mu, log = log)
      else dnbinom(x, size = mu/theta, mu = mu, log = log)
      fy
    }

    logLik <- dnegbinom(round(x), round(mu), theta, log = TRUE)
  }

  return(sum(logLik))
}


#' Calculate observation-wise log likelihood values
#'
#' @title Calculate observation-wise log likelihood values
#' @param x argument
#' @param ... other arguments
#' @examples
#' logLikFun()
#'
#' @rdname logLikFun
#' @export logLikFun
logLikFun <- function(x, ...){
  UseMethod("logLikFun")
}


#' Calculate observation-wise log likelihood values for gpls models
#'
#' @param model a fitted model
#'
#' @return numeric value
#' @export
#' @method logLikFun gpls
#' @keywords internal
logLikFun.gpls <- function(model) {

  y <- model$y_obs
  mu <- model$fitted
  family = model$family
  wts <- rep(1, length(y))

  if (family == "binomial") {
    if (!is.numeric(y)){
      y <- as.numeric(y) - 1
      y <- y == max(y)
      y <- as.numeric(y)
    }
    logLik <- dbinom(y, size = 1, prob = mu, log = TRUE)

  } else if (family == "poisson") {

    logLik <- dpois(y, mu, log = TRUE)

  } else if (family == "gaussian") {

    logLik <- dnorm(y, mu, sd(y - mu), log = TRUE)

  } else if (family == "Gamma"){
    sigma2 <- var(y - mu)
    logLik <- dgamma(y, mu^2/sigma2, mu/sigma2, log = TRUE)

  } else if (family == "inverse.gaussian"){

    dinvgauss <- function (x, mu = 1, sigma = 1, log = FALSE){
      if (any(mu < 0))
        stop(paste("mu must be positive", "\n", ""))
      if (any(sigma < 0))
        stop(paste("sigma must be positive", "\n", ""))
      if (any(x < 0))
        stop(paste("x must be positive", "\n", ""))
      log.lik <- (-0.5 * log(2 * pi) - log(sigma) - (3/2) * log(x) - ((x - mu)^2)/(2 * sigma^2 * (mu^2) * x))
      if (log == FALSE)
        fy <- exp(log.lik)
      else fy <- log.lik
      fy
    }

    logLik <- dinvgauss(y, mu, mu^3, log = TRUE)

  } else if (family == "negative.binomial"){

    theta = model$theta
    dnegbinom <- function (x, mu = 1, theta = 1, log = FALSE) {
      if (any(mu <= 0))
        stop(paste("mu must be greater than 0 ", "\n", ""))
      if (any(theta <= 0))
        stop(paste("theta must be greater than 0 ", "\n", ""))
      if (length(theta) > 1)
        fy <- ifelse(theta > 1e-04,
                     dnbinom(x, size = mu/theta, mu = mu, log = log),
                     dpois(x = x, lambda = mu, log = log)
                     )
      else fy <- if (theta < 1e-04)
        dpois(x = x, lambda = mu, log = log)
      else dnbinom(x, size = mu/theta, mu = mu, log = log)
      fy
    }

    logLik <- dnegbinom(round(x), round(mu), theta, log = TRUE)
  }

  return(logLik)
}