R/03_obtaining_gof.R

In Stan125/ghp: Generalisation of hierarchical partitioning

#' Obtain Goodness of Fit Figures
#'
#' Step Three: Obtain Goodness of Fit Figures for a set of computed models
#'
#' This function takes a list generated by \code{\link{mfit}} to obtain all
#' possible models and then exports a 'goodfit' object, used for further
#' partitioning.
#'
#' @param mfits A list with five components: \enumerate{ \item \strong{models}:
#'   A list with one or two elements. Each element has all possible model fits
#'   in it. Each element represents one distributional parameter. \item
#'   \strong{modelids}: A character vector with model id's corresponding to the
#'   models. in the models element \item \strong{expl_names}: names of
#'   explanatory variables (ungrouped case) or group names (grouped case). \item
#'   \strong{npar}: Number of parameters which are modeled. \item
#'   \strong{method}: One of "lm", "gamlss". } Best created by
#'   \code{\link{mfit}}.
#' @param gof Goodness of fit which should be obtained for all models. Which is
#'   available depends on the method with which the models were computed.
#'   Currently: \itemize{\item \strong{lm method}: "AIC", "r.squared", "loglik",
#'   "deviance" \item \strong{gamlss method}: "AIC", "deviance", "loglik", "R2e}
#' @return A gof object (type list), which has the following elements:
#'   \enumerate{\item \strong{gofs}: A list with \code{npar} elements, each
#'   being a vector with the goodness of fits of the models. \item
#'   \strong{model_ids}: A character vector with id's of the models: e.g. "x0",
#'   "x1", ... \item \strong{expl_names}: Names of explanatory variables
#'   (grouped) or names of groups (grouped). \item \strong{npar} Number of
#'   modeled paramateters. Can be 1 or 2. \item \strong{method}: The method used
#'   to compute models. Can be one of "lm" or "gamlss". \item \strong{gof}: A
#'   single character depicting the goodness of fit that was extracted.}
#' @export

gof <- function(mfits, gof) {
  # Obtain gofs
  if (mfits$method == "lm")
    gofs <- gof_lm(mfits, gof)
  if (mfits$method == "gamlss")
    gofs <- gof_gamlss(mfits, gof)

  # Assemble gof list
  gofs_list <- list(gofs = gofs, model_ids = mfits$model_ids,
                    expl_names = mfits$expl_names, npar = mfits$npar,
                    method = mfits$method, gof = gof)

  # Make GOF class
  class(gofs_list) <- "goodfit"

  # Return it
  return(gofs_list)
}

#' Internal: Get GOF's for models computed with "lm" method
#'
#' @keywords internal

gof_lm <- function(mfits, gof) {
  # Stop if goodness of fit not implemented
  available_gofs <- c("AIC", "r.squared", "loglik", "deviance")
  if (!gof %in% available_gofs)
    stop("Goodness of Fit not implemented")

  # Get gofs for mu parameter
  gofs_mu <- gof_getter(mfits$models$mu, gof)

  # Return gofs
  return(list(mu = gofs_mu))
}

#' Internal: Get GOF's for models computed with "gamlss" method
#'
#' @keywords internal

gof_gamlss <- function(mfits, gof) {
  # Stop if goodness of fit not implemented
  available_gofs <- c("AIC", "deviance", "loglik", "R2e")
  if (!gof %in% available_gofs)
    stop("Goodness of Fit not implemented")

  # Get gofs for mu parameter
  gofs_mu <- gof_getter(mfits$models$mu, gof)

  # Get GOF's for sigma parameter
  if (mfits$npar == 2)
    gofs_sigma <- gof_getter(mfits$models$sigma, gof)

  # Return list depending on npar
  if (mfits$npar == 1)
    return(list(mu = gofs_mu))
  if (mfits$npar == 2)
    return(list(mu = gofs_mu, sigma = gofs_sigma))
}

#' Function to obtain specific GOF for all given models
#'
#' @keywords internal

gof_getter <- function(models, gof) {
  if (gof == "r.squared")
    gofs <- sapply(models, FUN = function(x) return(summary(x)$r.squared))
  if (gof == "AIC")
    gofs <- sapply(models, FUN = function(x) return(AIC(x)))
  if (gof == "loglik")
    gofs <- sapply(models, FUN = function(x) return(logLik(x)))
  if (gof == "deviance")
    gofs <- sapply(models, FUN = function(x) return(deviance(x)))
  if (gof == "R2e") {
    L0 <- logLik(models[[1]])
    n <- models[[1]]$N
    gofs <- sapply(models, FUN = function(x) {
      Lm <- logLik(x)
      return(1 - (Lm / L0)^(-(2 / n) * L0))
    })
  }
  return(gofs)
}