R/gam_vibes.R

In Stan125/vibe: Variable Importance Beyond Linear Models

#' @title Variable Importance calculation for a `gam` object
#'
#' @description `vibe.gam` takes a fitted [mgcv::gam()] object and calculates
#' variable importance metrics by fitting the submodels required, extracting the
#' desired goodness-of-fit metric and applying variable importance metrics to
#' it.
#'
#' @param object A `gam` object, typically result of [mgcv::gam()].
#' @param varimp One of `c("hp", "relweights")`, which stand for hierarchical
#'   partitioning and relative weights.
#' @param gof Goodness-of-fit metric, the changes of which shall be
#'   analysed
#' @param ncores Number of cores used for the model fitting process, happening
#'   in [`part_core`].
#' @param progress Boolean. Do you want to see a progress bar?
#' @param ... Other arguments
#'
#' @importFrom stats model.frame
#' @importFrom stats family
#' @export
#' @examples
#' library("mgcv")
#' gam_ocat <- gam(
#'   satisfaction ~ admin + hygiene + time_appointment +
#'     quality_dr + diagnosis_exactness + equipment_modern +
#'     friendly_workers + parking_playingrooms_cafes,
#'   data = vibe::sat, family = ocat(R = 3)
#' )
#' hp_gam <- vibe(gam_ocat, varimp = "hp", gof = "R2e", progress = FALSE)
#' rw_gam <- vibe(gam_ocat, varimp = "relweights", gof = "R2e")
#' print(hp_gam)
#' print(rw_gam)
vibe.gam <- function(object,
                     varimp = "hp",
                     gof = "R2e",
                     ncores = 1,
                     progress = TRUE,
                     ...) {
  # Defensive Programming - is everything supplied the way it should be?
  args_supported(
    object = object,
    varimp = varimp,
    gof = gof,
    progress = progress
  )

  # Obtain data
  base_df <- model.frame(object)
  depvar <- base_df[, 1]
  depvar_name <- names(base_df)[1]
  expl_df <- base_df[, -c(1)]

  # Obtain family
  fam <- family(object)

  # Model Class - MC with added EE since it sounds cool
  mcee <- class_finder(object)

  if (varimp == "hp") {
    ## Obtain model ids
    model_ids <- mids(ncol(expl_df))

    ## Fit models and get goodness-es of fit
    gofs <- fit_and_gof(
      depvar = depvar,
      expl_df = expl_df,
      fam = fam,
      ncores = ncores,
      progress = progress,
      gof = gof,
      class = mcee,
      depvar_name = depvar_name,
      base_df = base_df
    )

    # Name vector with ID's
    names(gofs) <- model_ids
    gof_list <- list(
      gofs = list(mu = gofs),
      model_ids = model_ids,
      expl_names = list(mu = colnames(expl_df)),
      npar = 1,
      gof = gof,
      varimp = varimp
    )

    # Do hierarchical partitioning
    gof_res <- part(gof_list)

    # Summarize into nice format
    result <- make_vibe(
      results = gof_res,
      depvar_name = depvar_name,
      varimp = varimp,
      class = mcee
    )

    # Return
    return(result)
  } else if (varimp == "relweights") {
    # Relative Weights
    relweight_res <- rel_weights(
      expl_df = expl_df,
      fam = fam,
      depvar = depvar,
      gof = gof,
      class = mcee
    )

    # Summarize into nice format
    result <- make_vibe(
      results = relweight_res,
      depvar_name = depvar_name,
      varimp = varimp,
      class = mcee
    )

    # Return
    return(result)
  }
}