R/train_LightRuleFit.R

In rtemis: Schema-Driven Machine Learning and Visualization

# train_LightRuleFit.R
# ::rtemis::
# 2025- EDG rtemis.org

# %% .rule_coefs ----
#' Extract rule coefficients from the GLMNET step of LightRuleFit
#'
#' `coef()` on a `glmnet`/`cv.glmnet` model returns a sparse matrix, except
#' for multinomial (multiclass) models, where it returns a *list* of sparse
#' matrices, one per outcome class. This normalizes both to a dense numeric
#' matrix of rule coefficients (rules in rows, one column per coefficient
#' set), with the intercept row dropped.
#'
#' @param model `glmnet` or `cv.glmnet` model.
#'
#' @return Numeric matrix: rules x coefficient sets. Columns are named by
#' outcome class for multiclass models, "Coefficient" otherwise.
#'
#' @author EDG
#' @keywords internal
#' @noRd
.rule_coefs <- function(model) {
  coefs <- stats::coef(model)
  if (is.list(coefs)) {
    # Multinomial: one sparse matrix per class, each with the intercept first.
    out <- do.call(
      cbind,
      lapply(coefs, function(m) as.matrix(m)[-1L, 1L])
    )
    colnames(out) <- names(coefs)
  } else {
    out <- as.matrix(coefs)[-1L, 1L, drop = FALSE]
    colnames(out) <- "Coefficient"
  }
  out
} # /rtemis::.rule_coefs


# %% .rule_importance ----
#' Reduce a rule's coefficient set to a single importance value
#'
#' Single-coefficient models (binary classification, regression) have one
#' signed coefficient per rule, which is returned as-is: the sign is
#' meaningful (direction of the rule's effect) and drives the descending
#' ranking of rules.
#'
#' Multinomial (multiclass) models have one coefficient per outcome class per
#' rule, so importance is the sum of their absolute values (the L1 norm of the
#' coefficient row). `varimp_super()` also reports each class's coefficients
#' as their own measures.
#'
#' @param coef_matrix Numeric matrix: rules x coefficient sets, from
#' `.rule_coefs`.
#'
#' @return Numeric vector, one value per rule.
#'
#' @author EDG
#' @keywords internal
#' @noRd
.rule_importance <- function(coef_matrix) {
  if (NCOL(coef_matrix) == 1L) {
    coef_matrix[, 1L]
  } else {
    rowSums(abs(coef_matrix))
  }
} # /rtemis::.rule_importance


# %% train_.LightRuleFitHyperparameters ----
#' Train a LightRuleFit model
#'
#' Train a LightRuleFit model using LightGBM and GLMNET.
#'
#' @param hyperparameters `LightRuleFitHyperparameters` object: make using [setup_LightRuleFit].
#' @param x tabular data: Training set.
#' @param weights Numeric vector: Case weights.
#' @param dat_validation tabular data: Validation set.
#' @param verbosity Integer: If > 0, print messages.
#'
#' @author EDG
#' @keywords internal
#' @noRd
method(train_, LightRuleFitHyperparameters) <- function(
  hyperparameters,
  x,
  weights = NULL,
  dat_validation = NULL,
  execution_config = setup_FutureExecution(),
  verbosity = 1L
) {
  # Dependencies ----
  check_dependencies("lightgbm", "glmnet", "matrixStats", "gsubfn")

  # Hyperparameters ----
  # Hyperparameters must be either untunable or frozen by `train`.
  if (needs_tuning(hyperparameters)) {
    rtemis.core::abort(
      "Hyperparameters must be fixed - use train() instead.",
      class = c("rtemis_value_error", "rtemis_input_error")
    )
  }

  # Data ----
  check_supervised(
    x = x,
    dat_validation = dat_validation,
    allow_missing = TRUE,
    verbosity = verbosity
  )
  type <- supervised_type(x)
  nclasses <- if (type == "Classification") nlevels(x[[ncol(x)]]) else 1L

  # IFW for LightGBM ----
  # See setup_LightRuleFit: You can choose to use IFW for both steps with `ifw = TRUE` OR control each steps individually using `ifw_lightgbm` and `ifw_glmnet`.
  lightgbm_weights <- if (hyperparameters[["ifw_lightgbm"]]) {
    ifw(x[[ncol(x)]], verbosity = verbosity)
  } else {
    weights
  }

  # Train Gradient Boosting using LightGBM ----
  # Every hyperparameter this class shares with LightGBM, derived by
  # `LightRuleFit_lightgbm_params()`.
  lgbm_parameters <- update(
    setup_LightGBM(),
    get_hyperparams(hyperparameters, LightRuleFit_lightgbm_params())
  )
  lgbm_parameters@ifw <- hyperparameters[["ifw_lightgbm"]]
  mod_lgbm <- train(
    x = x,
    dat_validation = dat_validation,
    weights = lightgbm_weights,
    hyperparameters = lgbm_parameters,
    # tuner_config = tuner_config, # ? add tuner_config to LightRuleFitHyperparameters
    outer_resampling_config = NULL,
    execution_config = execution_config,
    verbosity = verbosity
  )

  # Extract Rules from Boosted Trees ----
  lgbm_rules <- extract_rules(
    mod_lgbm@model,
    n_iter = NULL,
    xnames = names(x),
    factor_levels = get_factor_levels(x),
    verbosity = verbosity
  )
  # A tree with no split yields an empty rule, which selects no subset.
  lgbm_rules <- lgbm_rules[nzchar(lgbm_rules)]
  if (length(lgbm_rules) == 0L) {
    rtemis.core::abort(
      "The boosting stage made no split on these ",
      NROW(x),
      " cases, so there are no rules to select. Lower `min_data_in_leaf` ",
      "or `min_sum_hessian_in_leaf`, or provide more cases.",
      class = c("rtemis_value_error", "rtemis_input_error")
    )
  }

  # Match cases x rules ----
  cases_by_rules <- match_cases_by_rules(x, lgbm_rules, verbosity = verbosity)

  # IFW for LASSO ----
  glmnet_weights <- if (hyperparameters[["ifw_glmnet"]]) {
    ifw(x[[ncol(x)]], verbosity = verbosity)
  } else {
    weights
  }

  # LASSO: Select Rules ----
  lasso_hyperparameters <- setup_GLMNET(
    alpha = hyperparameters[["alpha_glmnet"]],
    lambda = hyperparameters[["lambda_glmnet"]]
  )
  dat_rules <- data.frame(cases_by_rules, y = x[[ncol(x)]])
  colnames(dat_rules)[ncol(dat_rules)] <- colnames(x)[ncol(x)]
  mod_glmnet <- train(
    dat_rules,
    hyperparameters = lasso_hyperparameters,
    weights = glmnet_weights,
    execution_config = execution_config,
    verbosity = verbosity
  )

  # Rule coefficients ----
  # Multiclass models carry one coefficient per class per rule: rules are
  # selected if nonzero for any class, and reduced to a single importance
  # value for ranking (see `.rule_importance`). The per-class coefficients
  # are kept verbatim as extra columns of `rules_coefs`.
  coef_matrix <- .rule_coefs(mod_glmnet@model)
  rules_coefs <- data.frame(
    Rule = lgbm_rules,
    Coefficient = .rule_importance(coef_matrix)
  )
  if (NCOL(coef_matrix) > 1L) {
    rules_coefs <- cbind(rules_coefs, as.data.frame(coef_matrix))
  }
  nonzero_index <- which(rowSums(abs(coef_matrix)) > 0)
  rules_selected <- lgbm_rules[nonzero_index]
  cases_by_rules_selected <- cases_by_rules[, nonzero_index]
  Ncases_by_rules <- matrixStats::colSums2(cases_by_rules_selected)

  # Empirical risk ----
  if (type == "Classification" && nclasses == 2) {
    x <- as.data.table(x)
    empirical_risk <- vector("numeric", length(rules_selected))
    for (i in seq_along(rules_selected)) {
      match <- x[eval(parse(text = rules_selected[i])), ]
      freq <- table(match[[ncol(match)]])
      empirical_risk[i] <- freq[mod_glmnet@binclasspos] / sum(freq)
    }
  } else {
    empirical_risk <- NULL
  }

  # Format Rules ----
  # => Check format_LightRuleFit_rules' use of gsubfn::gsubfn
  rules_selected_formatted <- format_LightRuleFit_rules(
    rules_selected,
    decimal_places = 2
  )
  # appease R CMD check
  Coefficient <- NULL
  rules_selected_formatted_coefs <- data.table(
    Rule_ID = seq(rules_selected_formatted),
    Rule = rules_selected_formatted,
    N_Cases = Ncases_by_rules,
    Coefficient = rules_coefs[["Coefficient"]][nonzero_index]
  )
  # For multiclass, `Coefficient` is the aggregate importance used for
  # ranking; break out the signed per-class coefficients so the listing is
  # complete. (setorder below reorders these columns in lockstep.)
  if (NCOL(coef_matrix) > 1L) {
    rules_selected_formatted_coefs <- cbind(
      rules_selected_formatted_coefs,
      as.data.table(coef_matrix[nonzero_index, , drop = FALSE])
    )
  }
  if (type == "Classification" && nclasses == 2) {
    # appease R CMD check
    Empirical_Risk <- NULL
    rules_selected_formatted_coefs[, Empirical_Risk := empirical_risk]
  }
  setorder(rules_selected_formatted_coefs, -Coefficient)

  # LightRuleFit ----
  model <- LightRuleFit(
    model_lightgbm = mod_lgbm,
    model_glmnet = mod_glmnet,
    rules = lgbm_rules,
    rules_coefs = rules_coefs,
    rules_index = nonzero_index,
    rules_selected = rules_selected,
    rules_selected_formatted = rules_selected_formatted,
    rules_selected_formatted_coefs = rules_selected_formatted_coefs,
    y_levels = levels(x[[ncol(x)]]),
    xnames = names(x)[-ncol(x)],
    complexity_metrics = data.frame(
      n_rules_total = length(lgbm_rules),
      n_nonzero_rules = length(nonzero_index)
    )
  )
  # The boosting and lasso stages are fitted with their own hyperparameters,
  # which record the values they resolved; this records them on the
  # LightRuleFit settings they came from.
  lgbm_resolved <- mod_lgbm@hyperparameters
  hyperparameters <- record_backend_values(
    hyperparameters,
    c(
      lgbm_resolved@hyperparameters[LightRuleFit_lightgbm_params()],
      list(lambda_glmnet = mod_glmnet@hyperparameters[["lambda"]])
    )
  )
  list(model = model, preprocessor = NULL, hyperparameters = hyperparameters)
} # /rtemis::train_.LightRuleFitHyperparameters


# %% predict_super.LightRuleFitHyperparameters ----
#' Predict from LightRuleFit LightGBM model
#'
#' @param model LightRuleFit object trained using `train_LightRuleFit`.
#' @param newdata data.frame or similar: Data to predict on.
#'
#' @keywords internal
#' @noRd
method(predict_super, LightRuleFit) <- function(
  model,
  newdata,
  type = NULL,
  execution_config = NULL,
  verbosity = 0L
) {
  check_inherits(newdata, "data.frame")

  rules <- model@rules
  cases_by_rules <- match_cases_by_rules(newdata, rules, verbosity = verbosity)
  datm <- data.matrix(cases_by_rules)
  if (model@model_lightgbm@type == "Classification") {
    predicted <- predict(
      model@model_glmnet@model,
      newx = datm,
      type = "response"
    )
    if (length(model@y_levels) == 2) {
      predicted[, 1]
    } else {
      predicted
    }
  } else {
    as.numeric(predict(model@model_glmnet@model, newx = datm))
  }
} # /rtemis::predict_super.LightRuleFit


# %% varimp_super.LightRuleFit ----
#' Get variable importance from LightRuleFit model
#'
#' @param model LightRuleFit object trained using `train_LightRuleFit`.
#'
#' @keywords internal
#' @noRd
method(varimp_super, LightRuleFit) <- function(model) {
  coef_matrix <- .rule_coefs(model@model_glmnet@model)
  rules <- rownames(coef_matrix)
  multiclass <- NCOL(coef_matrix) > 1L
  alpha <- model@model_glmnet@hyperparameters[["alpha"]]
  penalty <- if (alpha == 1) {
    "lasso"
  } else if (alpha == 0) {
    "ridge"
  } else {
    "elastic net"
  }
  indicator <- paste0(
    "a rule is an indicator, so its coefficient does not depend on predictor ",
    "units."
  )
  measures <- list(
    Coefficient = importance_measure(
      rules,
      .rule_importance(coef_matrix),
      kind = if (multiclass) "coefficient_magnitude" else "coefficient",
      signed = !multiclass,
      direction = if (multiclass) "larger" else "absolute",
      description = if (multiclass) {
        paste0(
          "Sum over classes of the absolute ",
          penalty,
          " coefficients of each rule extracted from the boosted trees; ",
          indicator
        )
      } else {
        paste0(
          "The ",
          penalty,
          " coefficient of each rule extracted from the boosted trees, on the ",
          "scale of the linear predictor; ",
          indicator
        )
      }
    )
  )
  if (multiclass) {
    # A class label is kept as the measure name unless it repeats a name
    # already used, such as the aggregate's.
    for (level in colnames(coef_matrix)) {
      key <- level
      if (key %in% names(measures)) {
        key <- paste0(level, " (class)")
      }
      measures[[key]] <- importance_measure(
        rules,
        coef_matrix[, level],
        kind = "coefficient",
        signed = TRUE,
        direction = "absolute",
        description = paste0(
          "The ",
          penalty,
          " coefficient of each rule for class ",
          level,
          ", on the scale of the linear predictor."
        )
      )
    }
  }
  VariableImportance(measures = measures)
} # /rtemis::varimp_super.LightRuleFit