R/AIC_Calc.R

In spind: Spatial Methods and Indices

#'Akaike Information Criterion with correction for sample size
#'
#' @description Calculates AIC and AICc
#'
#' @param formula A model formula
#' @param family Family used to fit the model. \code{gaussian}, \code{binomial}, or \code{poisson} are supported
#' @param data A data frame
#' @param mu Fitted values from a model
#' @param n.eff Effective number of observations. Default is NULL
#'
#' @return A list with the following components
#' \describe{
#'   \item{\code{loglik}}{Log likelihood of the model}
#'   \item{\code{df}}{Degrees of freedom}
#'   \item{\code{AIC}}{AIC score for the specified model}
#'   \item{\code{AICc}}{AIC score corrected for small sample sizes}
#' }
#'
#' @author Gudrun Carl, Sam Levin
#'
#' @examples
#' data(musdata)
#' coords <- musdata[ ,4:5]
#' mglm <- glm(musculus ~ pollution + exposure, "poisson", musdata)
#'
#' aic <- aic.calc(musculus ~ pollution + exposure, "poisson", musdata,
#'                 mglm$fitted)
#' aic$AIC
#'
#' @importFrom stats model.frame model.matrix
#' @export

aic.calc <- function(formula, family, data, mu, n.eff = NULL){

  X <- stats::model.matrix(formula, data)
  if(is.vector(stats::model.frame(formula, data)[[1]])){
    y <- stats::model.frame(formula, data)[[1]]
    ntr <- 1
  }
  if(family == "binomial" & is.matrix(stats::model.frame(formula, data)[[1]])){
    y <- stats::model.frame(formula, data)[[1]][ ,1]
    ntr <- stats::model.frame(formula, data)[[1]][ ,1] +
           stats::model.frame(formula, data)[[1]][ ,2]
  }

  n <- dim(X)[1]
  nvar <- dim(X)[2]
  if(family == "gaussian"){
    sos <- sum((y - mu)^2)                    # sum of squares
    sigma2 <- sos / n                        # variance
    #loglik<- -n/2*log(2*pi) - n*log(sigma2^(1/2)) - 1/(2*sigma2)*sos
    #loglik<- -n/2*log(2*pi) - n*log(sigma2^(1/2)) - n/2
    #loglik<- -n/2*log(2*pi) - n/2*log(sigma2) - n/2
    #loglik<- -n/2*(log(2*pi) + log(sigma2) +1)
    loglik <- -n / 2 * (log(2 * pi * sigma2) + 1)   # log likelihood
    if(!is.null(n.eff)) loglik <- -n.eff / 2 * (log(2 * pi * sigma2) + 1)
    K <- nvar + 1            # nvar (number of pred.+interc.) +1 (for variance)
  }
  if(family == "binomial"){                         # choose= Binomialkoeff.
    loglik <- sum(y * log(mu / (1 - mu)) +
                  ntr * log(1 - mu) + log(choose(ntr, y)))

    K <- nvar              # without variance (since var. is function of mean)
  }
  if(family == "poisson"){
    # loglik<- sum(y*log(mu)-mu)  # useful for delta in multimodel inference
    loglik <- sum(y * log(mu) - mu) - sum(log(factorial(y)))
    K <- nvar              # without variance (since var. is function of mean)

  }

  AIC  <- -2 * loglik + 2 * K              # AIC
  AICc <- -2 * loglik + 2 * K + 2 * K * (K + 1) / (n - K - 1)    # AICc


  return(list(loglik = loglik, df = K, AIC = AIC, AICc = AICc))
}