diff --git a/DESCRIPTION b/DESCRIPTION
index 9c9ceb4..b9dabd5 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -22,6 +22,7 @@ Imports:
 Suggests:
     aorsf,
     C50,
+    Cubist,
     dbarts,
     grf,
     knitr,
diff --git a/NAMESPACE b/NAMESPACE
index c42289f..638d0f7 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,6 +4,7 @@ S3method(active_predictors,C5.0)
 S3method(active_predictors,ObliqueForest)
 S3method(active_predictors,bart)
 S3method(active_predictors,cforest)
+S3method(active_predictors,cubist)
 S3method(active_predictors,grf)
 S3method(active_predictors,lgb.Booster)
 S3method(active_predictors,party)
@@ -22,6 +23,7 @@ S3method(extract_rules,C5.0)
 S3method(extract_rules,ObliqueForest)
 S3method(extract_rules,bart)
 S3method(extract_rules,cforest)
+S3method(extract_rules,cubist)
 S3method(extract_rules,grf)
 S3method(extract_rules,lgb.Booster)
 S3method(extract_rules,party)
diff --git a/R/C5.0.R b/R/C5.0.R
index f3a9ab6..3ee6cd0 100644
--- a/R/C5.0.R
+++ b/R/C5.0.R
@@ -817,7 +817,7 @@ c5_extract_one <- function(tree_num, tree_lines, num_trials) {
 # @param data Training data for factor level information
 #
 # @return Tibble with id, rules, and tree columns
-c5_extract_rules_from_rules_text <- function(x, tree, data) {
+c5_extract_rules_from_rules_text <- function(x, tree) {
   # Parse rules text
   rules_lines <- strsplit(x$rules, "\n")[[1]]
 
@@ -826,8 +826,7 @@ c5_extract_rules_from_rules_text <- function(x, tree, data) {
     tree,
     c5_extract_rules_one_trial,
     rules_lines = rules_lines,
-    num_trials = x$trials["Actual"],
-    data = data
+    num_trials = x$trials["Actual"]
   )
 
   # Combine and format
@@ -840,8 +839,7 @@ c5_extract_rules_from_rules_text <- function(x, tree, data) {
 c5_extract_rules_one_trial <- function(
   tree_num,
   rules_lines,
-  num_trials,
-  data
+  num_trials
 ) {
   # Find rule set boundaries for this trial
   if (num_trials == 1) {
@@ -903,14 +901,14 @@ c5_parse_rules_block <- function(rules_lines, data) {
     condition_lines <- rule_block[grepl('type="[123]"', rule_block)]
 
     # Parse conditions into expression
-    rules_list[[i]] <- c5_parse_rule_conditions(condition_lines, data)
+    rules_list[[i]] <- c5_parse_rule_conditions(condition_lines)
   }
 
   rules_list
 }
 
 # Parse condition lines into R expression
-c5_parse_rule_conditions <- function(condition_lines, data) {
+c5_parse_rule_conditions <- function(condition_lines) {
   if (length(condition_lines) == 0) {
     return(rlang::expr(TRUE))
   }
@@ -1029,20 +1027,10 @@ c5_extract_rules_one <- function(tree_num, x, data) {
 #'   extract rules from. Default is `1L` for the first tree or trial. Values must be
 #'   between 1 and the number of actual trials (`x$trials["Actual"]`).
 #'   For rule-based models, this parameter refers to the trial number.
-#' @param data Data.frame containing the training data. Required for C5.0
-#'   models to properly parse tree or rule structure with correct factor levels.
 #' @export
-extract_rules.C5.0 <- function(x, tree = 1L, data = NULL, ...) {
+extract_rules.C5.0 <- function(x, tree = 1L, ...) {
   rlang::check_installed("C50")
 
-  # Require data parameter
-  if (is.null(data)) {
-    cli::cli_abort(
-      "{.arg data} is required for {.fn extract_rules.C5.0}.",
-      "i" = "Provide the training data to extract rules correctly."
-    )
-  }
-
   # Detect model type
   is_rule_model <- !is.null(x$rules) && nchar(x$rules) > 0
 
@@ -1069,15 +1057,35 @@ extract_rules.C5.0 <- function(x, tree = 1L, data = NULL, ...) {
 
   if (is_rule_model) {
     # Parse rule-based models directly
-    return(c5_extract_rules_from_rules_text(x, tree, data))
+    return(c5_extract_rules_from_rules_text(x, tree))
   }
 
-  # Tree-based model: existing logic
-  results <- lapply(tree, c5_extract_rules_one, x = x, data = data)
+  # Tree-based model: Parse tree text directly without party conversion
+  # Extract rules for each requested tree
+  results <- lapply(tree, function(tree_num) {
+    # For now, parse the entire tree - multi-trial support can be added later
+    rules <- c5_extract_rules_from_tree(x$tree, trial = tree_num)
+
+    # Rename 'trial' column to 'tree' for consistency with existing API
+    if ("trial" %in% names(rules)) {
+      names(rules)[names(rules) == "trial"] <- "tree"
+    }
+
+    # Ensure tree column has correct value
+    rules$tree <- tree_num
+
+    # Ensure consistent column order: id, rules, tree
+    rules[c("id", "rules", "tree")]
+  })
 
   # Combine and sort by tree then id
-  dplyr::bind_rows(results) |>
+  result_df <- dplyr::bind_rows(results) |>
     dplyr::arrange(tree, id)
+
+  # Add the appropriate S3 classes
+  class(result_df) <- c("rule_set_C5.0", "rule_set", class(result_df))
+
+  result_df
 }
 
 # ------------------------------------------------------------------------------
@@ -1136,3 +1144,364 @@ active_predictors.C5.0 <- function(x, tree = 1L, ...) {
   dplyr::bind_rows(results) |>
     dplyr::arrange(tree)
 }
+
+# ------------------------------------------------------------------------------
+# Direct C5.0 tree parsing functions (without party conversion)
+# ------------------------------------------------------------------------------
+
+#' Parse a single line of C5.0 tree text
+#'
+#' @param line Character string containing one line of tree text in key="value" format
+#' @return Named list of parsed properties
+#' @noRd
+c5_parse_tree_line <- function(line) {
+  if (length(line) == 0 || nchar(line) == 0) {
+    return(list())
+  }
+
+  # Use the existing c5_parse_attributes function which handles key="value" format
+  props <- c5_parse_attributes(line)
+
+  # Convert certain fields to appropriate types
+  if ("type" %in% names(props)) {
+    props$type <- as.integer(props$type)
+  }
+  if ("forks" %in% names(props)) {
+    props$forks <- as.integer(props$forks)
+  }
+  if ("cut" %in% names(props)) {
+    props$cut <- as.numeric(props$cut)
+  }
+
+  props
+}
+
+#' Build hierarchical tree structure from flat list of parsed lines
+#'
+#' C5.0 trees are stored in pre-order traversal (parent before children).
+#' This function reconstructs the hierarchy using the forks count.
+#'
+#' @param parsed_lines List of parsed node properties
+#' @return Root node of hierarchical tree structure
+#' @noRd
+c5_build_tree_structure <- function(parsed_lines) {
+  if (length(parsed_lines) == 0) {
+    return(NULL)
+  }
+
+  # Skip header lines (id, entries)
+  start_idx <- 1
+  while (
+    start_idx <= length(parsed_lines) &&
+      !("type" %in% names(parsed_lines[[start_idx]]))
+  ) {
+    start_idx <- start_idx + 1
+  }
+
+  if (start_idx > length(parsed_lines)) {
+    return(NULL)
+  }
+
+  # Recursive function to build tree
+  build_node <- function(idx) {
+    if (idx > length(parsed_lines)) {
+      return(list(node = NULL, next_idx = idx))
+    }
+
+    node_props <- parsed_lines[[idx]]
+    node <- list(properties = node_props, children = list())
+
+    # If this is a leaf node (type=0 or no forks), return it
+    if (
+      node_props$type == 0 || is.null(node_props$forks) || node_props$forks == 0
+    ) {
+      return(list(node = node, next_idx = idx + 1))
+    }
+
+    # Otherwise, recursively build children
+    next_idx <- idx + 1
+    for (i in seq_len(node_props$forks)) {
+      result <- build_node(next_idx)
+      node$children[[i]] <- result$node
+      next_idx <- result$next_idx
+    }
+
+    return(list(node = node, next_idx = next_idx))
+  }
+
+  result <- build_node(start_idx)
+  result$node
+}
+
+#' Extract all root-to-leaf paths from tree structure
+#'
+#' @param tree_node Tree node (from c5_build_tree_structure)
+#' @param current_path List of conditions accumulated so far
+#' @return List of paths, each containing conditions, class, and freq
+#' @noRd
+c5_extract_tree_paths <- function(tree_node, current_path = list()) {
+  if (is.null(tree_node)) {
+    return(list())
+  }
+
+  # Check if this is a leaf node
+  if (
+    tree_node$properties$type == 0 ||
+      is.null(tree_node$children) ||
+      length(tree_node$children) == 0
+  ) {
+    # Return completed path
+    return(list(list(
+      conditions = current_path,
+      class = tree_node$properties$class,
+      freq = tree_node$properties$freq
+    )))
+  }
+
+  # For non-leaf nodes, recurse on each child
+  paths <- list()
+  for (i in seq_along(tree_node$children)) {
+    condition <- c5_create_branch_condition(
+      tree_node$properties,
+      i,
+      length(tree_node$children)
+    )
+    child_paths <- c5_extract_tree_paths(
+      tree_node$children[[i]],
+      c(current_path, list(condition))
+    )
+    paths <- c(paths, child_paths)
+  }
+
+  paths
+}
+
+#' Create condition for a specific branch
+#'
+#' @param node_props Node properties (parsed from tree line)
+#' @param branch_index Which branch (1-based)
+#' @param num_branches Total number of branches from this node
+#' @return List with att, op, and value fields
+#' @noRd
+c5_create_branch_condition <- function(node_props, branch_index, num_branches) {
+  type <- node_props$type
+  att <- node_props$att
+
+  if (type == 2) {
+    # Threshold split
+    cut <- node_props$cut
+
+    if (num_branches == 2) {
+      # Binary split (no missing values)
+      if (branch_index == 1) {
+        return(list(att = att, op = "<=", value = cut))
+      } else {
+        return(list(att = att, op = ">", value = cut))
+      }
+    } else if (num_branches == 3) {
+      # Three-way split (with missing value branch)
+      if (branch_index == 1) {
+        return(list(att = att, op = "is.na", value = NULL))
+      } else if (branch_index == 2) {
+        return(list(att = att, op = "<=", value = cut))
+      } else {
+        return(list(att = att, op = ">", value = cut))
+      }
+    }
+  } else if (type == 1) {
+    # Discrete split
+    # Look for val field in properties
+    # For discrete splits, each branch has its own value
+    # The first branch gets the first value, etc.
+    # We need to extract from the elts field if present, or val field
+
+    # Check if there's an elts field (for multiple values per branch)
+    if (!is.null(node_props$elts)) {
+      # Parse the elts field - it contains comma-separated quoted values
+      elts <- c5_parse_elts(node_props$elts)
+      if (branch_index <= length(elts)) {
+        val <- elts[branch_index]
+        return(list(att = att, op = "==", value = val))
+      }
+    }
+
+    # Otherwise look for val field
+    if (!is.null(node_props$val)) {
+      return(list(att = att, op = "==", value = node_props$val))
+    }
+
+    # Fallback - this shouldn't happen with valid trees
+    return(list(att = att, op = "==", value = paste0("branch", branch_index)))
+  } else if (type == 3) {
+    # Subset split
+    # For subset splits, each branch can match multiple values
+    # The elts field contains the values for each branch
+
+    if (!is.null(node_props$elts)) {
+      # The elts field may contain multiple sets of values
+      # One for each branch, separated by some delimiter
+      # Need to parse this carefully
+
+      # For now, parse as comma-separated and assign to branches
+      all_elts <- c5_parse_elts(node_props$elts)
+
+      # If there's only one value, use ==
+      if (length(all_elts) == 1) {
+        return(list(att = att, op = "==", value = all_elts[1]))
+      } else {
+        # Multiple values, use %in%
+        return(list(att = att, op = "%in%", value = all_elts))
+      }
+    }
+
+    # Fallback
+    return(list(att = att, op = "%in%", value = character(0)))
+  }
+
+  # Unknown type - shouldn't happen
+  return(list(att = att, op = "==", value = "unknown"))
+}
+
+#' Parse elts field containing categorical values
+#'
+#' The elts field contains comma-separated quoted values like:
+#' "val1","val2","val3"
+#'
+#' @param elts_str Character string with elts values
+#' @return Character vector of unquoted values
+#' @noRd
+c5_parse_elts <- function(elts_str) {
+  if (is.null(elts_str) || length(elts_str) == 0) {
+    return(character(0))
+  }
+
+  # If elts_str is a vector, process first element
+  if (length(elts_str) > 1) {
+    elts_str <- elts_str[1]
+  }
+
+  if (nchar(elts_str) == 0) {
+    return(character(0))
+  }
+
+  # Remove outer quotes if present
+  elts_str <- gsub('^"|"$', '', elts_str)
+
+  # Split on ","
+  parts <- strsplit(elts_str, '","', fixed = TRUE)[[1]]
+
+  # Clean up any remaining quotes
+  parts <- gsub('"', '', parts)
+
+  parts
+}
+
+#' Convert tree paths to rule format
+#'
+#' @param paths List of paths from c5_extract_tree_paths
+#' @param trial Which trial these rules are from (for multi-trial models)
+#' @return Tibble with rules in standard format
+#' @noRd
+c5_convert_paths_to_rules <- function(paths, trial = 1L) {
+  if (length(paths) == 0) {
+    return(tibble::tibble(
+      id = integer(),
+      rules = list(),
+      trial = integer()
+    ))
+  }
+
+  # Store trial in a local variable to ensure it's accessible in lapply
+  trial_num <- trial
+
+  # Convert each path to a rule
+  rules <- lapply(seq_along(paths), function(i) {
+    path <- paths[[i]]
+
+    # Convert conditions to R expressions
+    if (length(path$conditions) == 0) {
+      # No conditions - always true
+      rule_expr <- rlang::expr(TRUE)
+    } else {
+      exprs <- lapply(path$conditions, function(cond) {
+        att_sym <- rlang::sym(cond$att)
+
+        if (cond$op == "is.na") {
+          rlang::call2("is.na", att_sym)
+        } else if (cond$op == "%in%") {
+          rlang::call2("%in%", att_sym, rlang::call2("c", !!!cond$value))
+        } else {
+          rlang::call2(cond$op, att_sym, cond$value)
+        }
+      })
+
+      # Combine with AND
+      rule_expr <- combine_rule_elements(exprs)
+    }
+
+    list(
+      trial = trial_num,
+      id = i,
+      rule = rule_expr,
+      class = path$class,
+      freq = path$freq
+    )
+  })
+
+  # Extract components before creating tibble to avoid scoping issues
+  id_vec <- vapply(rules, function(r) r$id, integer(1))
+  rules_list <- lapply(rules, function(r) r$rule)
+  trial_vec <- vapply(rules, function(r) r$trial, integer(1))
+
+  # Convert to tibble with consistent naming and ordering
+  tibble::tibble(
+    id = id_vec,
+    rules = rules_list,
+    trial = trial_vec
+  )
+}
+
+#' Extract rules directly from C5.0 tree text
+#'
+#' This is the main entry point for direct tree parsing.
+#' It combines all the helper functions to extract rules from tree text.
+#'
+#' @param tree_text Character string containing the tree in C5.0 format
+#' @param trial Which trial to extract (0 for all, 1+ for specific trial)
+#' @return Tibble with extracted rules
+#' @noRd
+c5_extract_rules_from_tree <- function(tree_text, trial = 1L) {
+  # Split into lines
+  lines <- strsplit(tree_text, "\n")[[1]]
+
+  # Remove empty lines
+  lines <- lines[nchar(trimws(lines)) > 0]
+
+  if (length(lines) == 0) {
+    return(tibble::tibble(
+      id = integer(),
+      rules = list(),
+      trial = integer()
+    ))
+  }
+
+  # Parse each line
+  parsed_lines <- lapply(lines, c5_parse_tree_line)
+
+  # Build tree structure
+  tree_structure <- c5_build_tree_structure(parsed_lines)
+
+  if (is.null(tree_structure)) {
+    return(tibble::tibble(
+      id = integer(),
+      rules = list(),
+      trial = integer()
+    ))
+  }
+
+  # Extract paths
+  paths <- c5_extract_tree_paths(tree_structure)
+
+  # Convert to rules
+  c5_convert_paths_to_rules(paths, trial)
+}
diff --git a/R/cubist.R b/R/cubist.R
new file mode 100644
index 0000000..064e62e
--- /dev/null
+++ b/R/cubist.R
@@ -0,0 +1,242 @@
+#' Extract rules from a Cubist model
+#'
+#' @description
+#' Extracts rule conditions from a Cubist regression model as R expressions.
+#' Each rule consists of conditions that define when a linear model applies.
+#'
+#' @param x A `cubist` object from the Cubist package.
+#' @param committee An integer vector specifying which committee(s) to extract
+#'   rules from. Defaults to `1L` (first committee). Values must be between 1
+#'   and the total number of committees in the model.
+#' @param ... Not used.
+#'
+#' @return A tibble with columns:
+#'   - `committee`: Integer committee number
+#'   - `id`: Integer rule number within the committee
+#'   - `rules`: List column containing R expressions for each rule's conditions
+#'
+#' @details
+#' Cubist models use committees (similar to boosting iterations) where each
+#' committee contains multiple rules. Each rule has:
+#' - Conditions that determine when the rule applies (splits on predictors)
+#' - A linear model that makes predictions when conditions are met
+#'
+#' This function extracts the conditions as R expressions that can be evaluated
+#' on data. Rules with no conditions (applying to all data) return `TRUE`.
+#'
+#' The expressions use standard R operators:
+#' - Continuous splits: `>`, `<=`, etc.
+#' - Categorical single value: `==`
+#' - Categorical multiple values: `%in%`
+#' - Missing values: `is.na()`
+#'
+#' @examples
+#' \dontrun{
+#' library(Cubist)
+#' library(lorax)
+#'
+#' # Create sample data
+#' set.seed(1)
+#' n <- 100
+#' p <- 5
+#' X <- matrix(rnorm(n * p), n, p)
+#' colnames(X) <- paste0("x", 1:p)
+#' y <- X[, 1] + X[, 2]^2 + rnorm(n)
+#'
+#' # Fit Cubist model with multiple committees
+#' mod <- cubist(X, y, committees = 3)
+#'
+#' # Extract rules from first committee
+#' rules <- extract_rules(mod)
+#' rules
+#'
+#' # Extract from multiple committees
+#' rules_all <- extract_rules(mod, committee = 1:3)
+#'
+#' # Convert to readable text
+#' rule_text(rules$rules[[1]])
+#' }
+#'
+#' @seealso [rules::tidy.cubist()] for extracting rules as text strings
+#' @export
+extract_rules.cubist <- function(x, committee = 1L, ...) {
+  rlang::check_installed("Cubist")
+
+  # Validate committee parameter
+  if (!is.numeric(committee) || !all(committee == as.integer(committee))) {
+    cli::cli_abort(
+      "{.arg committee} must be an integer vector, not {.obj_type_friendly {committee}}.",
+      call = rlang::caller_env()
+    )
+  }
+
+  committee <- as.integer(committee)
+  committee <- unique(committee)
+
+  # Check bounds
+  num_committees <- x$committees
+  if (any(committee < 1L) || any(committee > num_committees)) {
+    cli::cli_abort(
+      "{.arg committee} values must be between 1 and {num_committees}, not {committee}.",
+      call = rlang::caller_env()
+    )
+  }
+
+  # Extract rules for each committee
+  results <- lapply(
+    committee,
+    cubist_extract_rules_from_committee,
+    x = x
+  )
+
+  # Combine results
+  result_df <- dplyr::bind_rows(results) |>
+    dplyr::arrange(committee, id)
+
+  # Add classes following lorax pattern
+  class(result_df) <- c("rule_set_cubist", "rule_set", class(result_df))
+
+  result_df
+}
+
+# Internal helper to extract rules from a single committee
+cubist_extract_rules_from_committee <- function(x, committee_num) {
+  # Get splits for this committee (may be NULL if no splits)
+  if (!is.null(x$splits)) {
+    # Note: splits$committee is integer, so we can compare directly
+    committee_splits <- x$splits[
+      x$splits$committee == committee_num,
+      ,
+      drop = FALSE
+    ]
+  } else {
+    committee_splits <- data.frame()
+  }
+
+  # Count rules in this committee
+  # Use coefficients to determine number of rules
+  # Note: committee and rule columns are character in Cubist
+  committee_coefs <- x$coefficients[
+    x$coefficients$committee == as.character(committee_num),
+    ,
+    drop = FALSE
+  ]
+  rule_nums <- as.integer(unique(committee_coefs$rule))
+  num_rules <- length(rule_nums)
+
+  if (num_rules == 0) {
+    return(tibble::tibble(
+      committee = integer(),
+      id = integer(),
+      rules = list()
+    ))
+  }
+
+  # Extract conditions for each rule
+  rule_exprs <- lapply(rule_nums, function(rule_num) {
+    # Get splits for this rule
+    rule_splits <- committee_splits[
+      committee_splits$rule == rule_num,
+      ,
+      drop = FALSE
+    ]
+
+    if (nrow(rule_splits) == 0) {
+      # No conditions - rule applies to all data
+      return(rlang::expr(TRUE))
+    }
+
+    # Convert each split to an expression
+    split_exprs <- lapply(seq_len(nrow(rule_splits)), function(i) {
+      cubist_split_to_expr(rule_splits[i, , drop = FALSE])
+    })
+
+    # Combine with AND
+    combine_rule_elements(split_exprs)
+  })
+
+  tibble::tibble(
+    committee = rep(committee_num, num_rules),
+    id = rule_nums,
+    rules = rule_exprs
+  )
+}
+
+# Internal helper to convert a single split to an R expression
+cubist_split_to_expr <- function(split_row) {
+  var_name <- split_row$variable
+  # Remove quotes that may be present for factor variables
+  var_name <- gsub('"', '', var_name)
+  var_sym <- rlang::sym(var_name)
+
+  if (split_row$type == "type2") {
+    # Continuous split
+    direction <- split_row$dir
+    value <- split_row$value
+
+    # Check for missing value condition
+    if (is.na(value) && direction == "=") {
+      return(rlang::call2("is.na", var_sym))
+    }
+
+    # Standard continuous split
+    return(rlang::call2(direction, var_sym, value))
+  } else if (split_row$type == "type3") {
+    # Categorical split
+    categories <- split_row$category
+
+    # Parse the categories (remove quotes, split on comma)
+    categories <- gsub('"', '', categories)
+    categories <- trimws(strsplit(categories, ",")[[1]])
+
+    if (length(categories) == 1) {
+      # Single value - use ==
+      return(rlang::call2("==", var_sym, categories[1]))
+    } else {
+      # Multiple values - use %in%
+      # Create c(...) call with all categories
+      c_call <- rlang::call2("c", !!!categories)
+      return(rlang::call2("%in%", var_sym, c_call))
+    }
+  } else {
+    cli::cli_abort(
+      "Unknown split type: {split_row$type}",
+      call = rlang::caller_env()
+    )
+  }
+}
+
+#' @rdname active_predictors
+#' @export
+active_predictors.cubist <- function(x, ...) {
+  rlang::check_installed("Cubist")
+
+  # Extract active predictors from the model object
+  # Need to compute this from splits and coefficients due to a bug in Cubist$vars$used
+
+  # Get variable names from splits (if any)
+  split_vars <- if (!is.null(x$splits)) {
+    vars <- unique(as.character(x$splits$variable))
+    # Remove quotes that may be present for factor variables
+    gsub('"', '', vars)
+  } else {
+    character(0)
+  }
+
+  # Get variable names from coefficients (those with non-NA values)
+  # Exclude non-predictor columns
+  coef_df <- x$coefficients
+  exclude_cols <- c("(Intercept)", "committee", "rule", "tmp")
+  predictor_cols <- setdiff(names(coef_df), exclude_cols)
+
+  # Find predictors with any non-NA coefficients
+  coef_vars <- predictor_cols[sapply(predictor_cols, function(col) {
+    any(!is.na(coef_df[[col]]))
+  })]
+
+  # Combine splits and coefficient variables
+  active_vars <- union(split_vars, coef_vars)
+
+  # Return using lorax constructor
+  new_active_predictors(active_vars)
+}
diff --git a/man/active_predictors.Rd b/man/active_predictors.Rd
index 56a5868..4b99af2 100644
--- a/man/active_predictors.Rd
+++ b/man/active_predictors.Rd
@@ -1,12 +1,13 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/C5.0.R, R/ObliqueForest.R, R/bart.R,
-%   R/cforest.R, R/generics.R, R/grf.R, R/lgb.Booster.R, R/party.R,
+%   R/cforest.R, R/cubist.R, R/generics.R, R/grf.R, R/lgb.Booster.R, R/party.R,
 %   R/randomForest.R, R/ranger.R, R/rpart.R, R/xgb.Booster.R
 \name{active_predictors.C5.0}
 \alias{active_predictors.C5.0}
 \alias{active_predictors.ObliqueForest}
 \alias{active_predictors.bart}
 \alias{active_predictors.cforest}
+\alias{active_predictors.cubist}
 \alias{active_predictors}
 \alias{active_predictors.grf}
 \alias{active_predictors.lgb.Booster}
@@ -25,6 +26,8 @@
 
 \method{active_predictors}{cforest}(x, tree = 1L, ...)
 
+\method{active_predictors}{cubist}(x, ...)
+
 active_predictors(x, ...)
 
 \method{active_predictors}{grf}(x, tree = 1L, ...)
diff --git a/man/extract_rules.Rd b/man/extract_rules.Rd
index a184966..4b116dd 100644
--- a/man/extract_rules.Rd
+++ b/man/extract_rules.Rd
@@ -10,7 +10,7 @@
 \alias{extract_rules.ranger}
 \title{Extract an expression that defines a path to a terminal node}
 \usage{
-\method{extract_rules}{C5.0}(x, tree = 1L, data = NULL, ...)
+\method{extract_rules}{C5.0}(x, tree = 1L, ...)
 
 \method{extract_rules}{cforest}(x, tree = 1L, ...)
 
@@ -29,10 +29,10 @@ extract_rules(x, ...)
 Default is \code{1L} for the first tree. Values must be between 1 and the
 number of trees in the forest (\code{x$num.trees}).}
 
+\item{...}{Other arguments passed to methods}
+
 \item{data}{Data.frame containing the training data. Required for ranger
 models to properly extract rules with fitted values and node summaries.}
-
-\item{...}{Other arguments passed to methods}
 }
 \value{
 A data frame with column \code{rules} (an R expression) and \code{id} (an
diff --git a/man/extract_rules.cubist.Rd b/man/extract_rules.cubist.Rd
new file mode 100644
index 0000000..7da46ac
--- /dev/null
+++ b/man/extract_rules.cubist.Rd
@@ -0,0 +1,79 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/cubist.R
+\name{extract_rules.cubist}
+\alias{extract_rules.cubist}
+\title{Extract rules from a Cubist model}
+\usage{
+\method{extract_rules}{cubist}(x, committee = 1L, ...)
+}
+\arguments{
+\item{x}{A \code{cubist} object from the Cubist package.}
+
+\item{committee}{An integer vector specifying which committee(s) to extract
+rules from. Defaults to \code{1L} (first committee). Values must be between 1
+and the total number of committees in the model.}
+
+\item{...}{Not used.}
+}
+\value{
+A tibble with columns:
+\itemize{
+\item \code{committee}: Integer committee number
+\item \code{id}: Integer rule number within the committee
+\item \code{rules}: List column containing R expressions for each rule's conditions
+}
+}
+\description{
+Extracts rule conditions from a Cubist regression model as R expressions.
+Each rule consists of conditions that define when a linear model applies.
+}
+\details{
+Cubist models use committees (similar to boosting iterations) where each
+committee contains multiple rules. Each rule has:
+\itemize{
+\item Conditions that determine when the rule applies (splits on predictors)
+\item A linear model that makes predictions when conditions are met
+}
+
+This function extracts the conditions as R expressions that can be evaluated
+on data. Rules with no conditions (applying to all data) return \code{TRUE}.
+
+The expressions use standard R operators:
+\itemize{
+\item Continuous splits: \code{>}, \code{<=}, etc.
+\item Categorical single value: \code{==}
+\item Categorical multiple values: \code{\%in\%}
+\item Missing values: \code{is.na()}
+}
+}
+\examples{
+\dontrun{
+library(Cubist)
+library(lorax)
+
+# Create sample data
+set.seed(1)
+n <- 100
+p <- 5
+X <- matrix(rnorm(n * p), n, p)
+colnames(X) <- paste0("x", 1:p)
+y <- X[, 1] + X[, 2]^2 + rnorm(n)
+
+# Fit Cubist model with multiple committees
+mod <- cubist(X, y, committees = 3)
+
+# Extract rules from first committee
+rules <- extract_rules(mod)
+rules
+
+# Extract from multiple committees
+rules_all <- extract_rules(mod, committee = 1:3)
+
+# Convert to readable text
+rule_text(rules$rules[[1]])
+}
+
+}
+\seealso{
+\code{\link[rules:tidy.cubist]{rules::tidy.cubist()}} for extracting rules as text strings
+}
diff --git a/tests/testthat/_snaps/C5.0.md b/tests/testthat/_snaps/C5.0.md
index 2e600b5..8535581 100644
--- a/tests/testthat/_snaps/C5.0.md
+++ b/tests/testthat/_snaps/C5.0.md
@@ -1,15 +1,15 @@
 # as.party.C5.0 validates tree parameter
 
     Code
-      as.party(c5_model, tree = 0, data = penguins)
+      lorax::as.party(c5_model, tree = 0, data = penguins)
     Condition
-      Error in `as.party()`:
+      Error in `lorax::as.party()`:
       ! `tree` must be >= 1, not 0.
 
 ---
 
     Code
-      as.party(c5_model, tree = 5, data = penguins)
+      lorax::as.party(c5_model, tree = 5, data = penguins)
     Condition
       Warning:
       `tree` = 5 exceeds number of trials (1). Using `tree` = 1 instead.
@@ -34,39 +34,39 @@
 ---
 
     Code
-      as.party(c5_model, tree = c(1, 2), data = penguins)
+      lorax::as.party(c5_model, tree = c(1, 2), data = penguins)
     Condition
-      Error in `as.party()`:
+      Error in `lorax::as.party()`:
       ! `tree` must be a single integer, not a double vector.
 
 ---
 
     Code
-      as.party(c5_model, tree = "1", data = penguins)
+      lorax::as.party(c5_model, tree = "1", data = penguins)
     Condition
-      Error in `as.party()`:
+      Error in `lorax::as.party()`:
       ! `tree` must be a single integer, not a string.
 
 # as.party.C5.0 requires data parameter
 
     Code
-      as.party(c5_model, tree = 1)
+      lorax::as.party(c5_model, tree = 1)
     Condition
-      Error in `as.party()`:
+      Error in `lorax::as.party()`:
       ! `data` is required for `as.party.C5.0()`.
 
 # as.party.C5.0 works with boosted models
 
     Code
-      as.party(c5_boost, tree = 0, data = penguins)
+      lorax::as.party(c5_boost, tree = 0, data = penguins)
     Condition
-      Error in `as.party()`:
+      Error in `lorax::as.party()`:
       ! `tree` must be >= 1, not 0.
 
 ---
 
     Code
-      as.party(c5_boost, tree = 11, data = penguins)
+      lorax::as.party(c5_boost, tree = 11, data = penguins)
     Condition
       Warning:
       `tree` = 11 exceeds number of trials (10). Using `tree` = 10 instead.
@@ -87,9 +87,9 @@
 # as.party.C5.0 rejects rule-based models
 
     Code
-      as.party(c5_rules, tree = 1, data = penguins)
+      lorax::as.party(c5_rules, tree = 1, data = penguins)
     Condition
-      Error in `as.party()`:
+      Error in `lorax::as.party()`:
       ! C50 rule-based models are not supported. Use tree-based models only.
 
 # extract_rules.C5.0() validates tree argument
@@ -124,22 +124,6 @@
       Error:
       ! `tree` values must be between 1 and 5.
 
-# extract_rules.C5.0() requires data parameter
-
-    Code
-      extract_rules(c5_tree, tree = 1L)
-    Condition
-      Error in `extract_rules()`:
-      ! `data` is required for `extract_rules.C5.0()`.
-
----
-
-    Code
-      extract_rules(c5_tree, tree = 1L, data = NULL)
-    Condition
-      Error in `extract_rules()`:
-      ! `data` is required for `extract_rules.C5.0()`.
-
 # extract_rules.C5.0() validates tree argument for rule models
 
     Code
@@ -159,7 +143,7 @@
 # active_predictors.C5.0() validates tree argument
 
     Code
-      active_predictors(fit, tree = "1")
+      lorax::active_predictors(fit, tree = "1")
     Condition
       Error:
       ! `tree` must be an integer vector, not a string.
@@ -167,7 +151,7 @@
 ---
 
     Code
-      active_predictors(fit, tree = 1.5)
+      lorax::active_predictors(fit, tree = 1.5)
     Condition
       Error:
       ! `tree` must be an integer vector, not a number.
@@ -175,7 +159,7 @@
 ---
 
     Code
-      active_predictors(fit, tree = 0L)
+      lorax::active_predictors(fit, tree = 0L)
     Condition
       Error:
       ! `tree` values must be between 1 and 3, not 0.
@@ -183,7 +167,7 @@
 ---
 
     Code
-      active_predictors(fit, tree = 10L)
+      lorax::active_predictors(fit, tree = 10L)
     Condition
       Error:
       ! `tree` values must be between 1 and 3, not 10.
diff --git a/tests/testthat/helper-data.R b/tests/testthat/helper-data.R
index ed988fc..7c269c1 100644
--- a/tests/testthat/helper-data.R
+++ b/tests/testthat/helper-data.R
@@ -115,3 +115,108 @@ get_single_factor_data <- function() {
   names(data)[names(data) == "island"] <- "x"
   data
 }
+
+# Cubist-specific test helpers
+
+# Get Ames housing data for Cubist regression tests
+get_ames_data <- function(n = 200) {
+  skip_if_not_installed("modeldata")
+  data("ames", package = "modeldata", envir = environment())
+  ames <- get("ames", envir = environment())
+
+  # Select relevant numeric and factor predictors
+  # Focus on variables likely to affect house price
+  ames_subset <- ames[
+    1:min(n, nrow(ames)),
+    c(
+      "Sale_Price", # outcome
+      "Gr_Liv_Area", # numeric: above grade living area
+      "Year_Built", # numeric: original construction date
+      "Overall_Cond", # ordinal: overall condition
+      "Garage_Cars", # numeric: size of garage in car capacity
+      "Total_Bsmt_SF", # numeric: total basement area
+      "Neighborhood", # factor: physical location
+      "Bldg_Type" # factor: type of dwelling
+    )
+  ]
+
+  # Convert to data frame and ensure no NAs
+  ames_subset <- na.omit(as.data.frame(ames_subset))
+
+  # Separate predictors and outcome
+  list(
+    x = ames_subset[, -1],
+    y = ames_subset[, 1]
+  )
+}
+
+# Get Sacramento housing data for Cubist regression tests
+get_sacramento_data <- function() {
+  skip_if_not_installed("modeldata")
+  data("Sacramento", package = "modeldata", envir = environment())
+  sacramento <- get("Sacramento", envir = environment())
+
+  # price is the outcome, all others are predictors
+  list(
+    x = sacramento[, c("beds", "baths", "sqft", "type", "city")],
+    y = sacramento$price
+  )
+}
+
+# Get simple regression data with clear structure for testing splits
+get_structured_regression_data <- function(n = 100) {
+  set.seed(42)
+  X <- data.frame(
+    x1 = runif(n, 0, 10),
+    x2 = runif(n, 0, 10),
+    x3 = factor(sample(c("A", "B", "C"), n, replace = TRUE))
+  )
+
+  # Create outcome with clear conditional structure
+  # This ensures Cubist will create meaningful splits
+  y <- ifelse(
+    X$x1 > 5,
+    ifelse(
+      X$x2 > 5,
+      X$x1 + X$x2, # High x1, high x2
+      X$x1 - X$x2
+    ), # High x1, low x2
+    ifelse(
+      X$x3 == "A",
+      2 * X$x2, # Low x1, category A
+      X$x1 + 3
+    )
+  ) # Low x1, category B or C
+  y <- y + rnorm(n, sd = 0.5)
+
+  list(x = X, y = y)
+}
+
+# Get data that produces simple models (for edge case testing)
+get_simple_cubist_data <- function(n = 10) {
+  set.seed(123)
+  # Very simple linear relationship, likely to produce single rule
+  x <- data.frame(x1 = seq(1, n))
+  y <- x$x1 + rnorm(n, sd = 0.01)
+  list(x = x, y = y)
+}
+
+# Get high-dimensional data for testing edge cases
+get_high_dim_cubist_data <- function(n = 20, p = 50) {
+  set.seed(456)
+  x <- matrix(rnorm(n * p), n, p)
+  colnames(x) <- paste0("x", 1:p)
+  # Only first two predictors are relevant
+  y <- x[, 1] + 0.5 * x[, 2] + rnorm(n, sd = 0.1)
+  list(x = as.data.frame(x), y = y)
+}
+
+# Get noisy data where predictors have little relationship with outcome
+# Used specifically for testing "no valid splits" scenarios
+get_noisy_cubist_data <- function(n = 100, p = 5) {
+  set.seed(789)
+  x <- matrix(rnorm(n * p), n, p)
+  colnames(x) <- paste0("x", 1:p)
+  y <- rnorm(n) # Completely independent of x
+  list(x = as.data.frame(x), y = y)
+}
diff --git a/tests/testthat/test-C5.0.R b/tests/testthat/test-C5.0.R
index d524f55..5e572dd 100644
--- a/tests/testthat/test-C5.0.R
+++ b/tests/testthat/test-C5.0.R
@@ -5,7 +5,7 @@ test_that("as.party.C5.0 returns valid party object for single tree", {
   penguins <- get_penguins_data()
 
   c5_model <- C50::C5.0(species ~ ., data = penguins)
-  p <- as.party(c5_model, tree = 1, data = penguins)
+  p <- lorax::as.party(c5_model, tree = 1, data = penguins)
 
   expect_s3_class(p, "party")
   expect_s3_class(p$node, "partynode")
@@ -18,7 +18,7 @@ test_that("as.party.C5.0 works with binary classification", {
   data <- get_binary_data()
 
   c5_model <- C50::C5.0(y ~ ., data = data)
-  p <- as.party(c5_model, tree = 1, data = data)
+  p <- lorax::as.party(c5_model, tree = 1, data = data)
 
   expect_s3_class(p, "party")
   expect_s3_class(p$node, "partynode")
@@ -30,7 +30,7 @@ test_that("as.party.C5.0 works with penguins data", {
   penguins_data <- get_penguins_data()
 
   c5_model <- C50::C5.0(species ~ ., data = penguins_data)
-  p <- as.party(c5_model, tree = 1, data = penguins_data)
+  p <- lorax::as.party(c5_model, tree = 1, data = penguins_data)
 
   expect_s3_class(p, "party")
   expect_s3_class(p$node, "partynode")
@@ -43,14 +43,20 @@ test_that("as.party.C5.0 validates tree parameter", {
   penguins <- get_penguins_data()
   c5_model <- C50::C5.0(species ~ ., data = penguins)
 
-  expect_snapshot(as.party(c5_model, tree = 0, data = penguins), error = TRUE)
+  expect_snapshot(
+    lorax::as.party(c5_model, tree = 0, data = penguins),
+    error = TRUE
+  )
   # When tree exceeds num_trials, it should warn and use max tree
-  expect_snapshot(as.party(c5_model, tree = 5, data = penguins))
+  expect_snapshot(lorax::as.party(c5_model, tree = 5, data = penguins))
+  expect_snapshot(
+    lorax::as.party(c5_model, tree = c(1, 2), data = penguins),
+    error = TRUE
+  )
   expect_snapshot(
-    as.party(c5_model, tree = c(1, 2), data = penguins),
+    lorax::as.party(c5_model, tree = "1", data = penguins),
     error = TRUE
   )
-  expect_snapshot(as.party(c5_model, tree = "1", data = penguins), error = TRUE)
 })
 
 test_that("as.party.C5.0 requires data parameter", {
@@ -60,7 +66,7 @@ test_that("as.party.C5.0 requires data parameter", {
   penguins <- get_penguins_data()
   c5_model <- C50::C5.0(species ~ ., data = penguins)
 
-  expect_snapshot(as.party(c5_model, tree = 1), error = TRUE)
+  expect_snapshot(lorax::as.party(c5_model, tree = 1), error = TRUE)
 })
 
 test_that("as.party.C5.0 works with boosted models", {
@@ -71,8 +77,8 @@ test_that("as.party.C5.0 works with boosted models", {
   c5_boost <- C50::C5.0(species ~ ., data = penguins, trials = 10)
 
   # Extract first and second trees
-  p1 <- as.party(c5_boost, tree = 1, data = penguins)
-  p2 <- as.party(c5_boost, tree = 2, data = penguins)
+  p1 <- lorax::as.party(c5_boost, tree = 1, data = penguins)
+  p2 <- lorax::as.party(c5_boost, tree = 2, data = penguins)
 
   expect_s3_class(p1, "party")
   expect_s3_class(p2, "party")
@@ -81,9 +87,12 @@ test_that("as.party.C5.0 works with boosted models", {
   expect_false(isTRUE(all.equal(p1, p2)))
 
   # Check tree number validation
-  expect_snapshot(as.party(c5_boost, tree = 0, data = penguins), error = TRUE)
+  expect_snapshot(
+    lorax::as.party(c5_boost, tree = 0, data = penguins),
+    error = TRUE
+  )
   # When tree exceeds num_trials, it should warn and use max tree
-  expect_snapshot(as.party(c5_boost, tree = 11, data = penguins))
+  expect_snapshot(lorax::as.party(c5_boost, tree = 11, data = penguins))
 })
 
 test_that("as.party.C5.0 extracts different boosted trees", {
@@ -94,9 +103,9 @@ test_that("as.party.C5.0 extracts different boosted trees", {
   c5_boost <- C50::C5.0(species ~ ., data = penguins, trials = 5)
 
   # Extract multiple trees
-  p1 <- as.party(c5_boost, tree = 1, data = penguins)
-  p2 <- as.party(c5_boost, tree = 2, data = penguins)
-  p3 <- as.party(c5_boost, tree = 5, data = penguins)
+  p1 <- lorax::as.party(c5_boost, tree = 1, data = penguins)
+  p2 <- lorax::as.party(c5_boost, tree = 2, data = penguins)
+  p3 <- lorax::as.party(c5_boost, tree = 5, data = penguins)
 
   # All should be valid
   expect_s3_class(p1, "party")
@@ -128,7 +137,10 @@ test_that("as.party.C5.0 rejects rule-based models", {
   penguins <- get_penguins_data()
   c5_rules <- C50::C5.0(species ~ ., data = penguins, rules = TRUE)
 
-  expect_snapshot(as.party(c5_rules, tree = 1, data = penguins), error = TRUE)
+  expect_snapshot(
+    lorax::as.party(c5_rules, tree = 1, data = penguins),
+    error = TRUE
+  )
 })
 
 test_that("as.party.C5.0 does not show asterisks in node summaries", {
@@ -138,7 +150,7 @@ test_that("as.party.C5.0 does not show asterisks in node summaries", {
   penguins <- get_penguins_data()
 
   c5_model <- C50::C5.0(species ~ ., data = penguins)
-  p <- as.party(c5_model, tree = 1, data = penguins)
+  p <- lorax::as.party(c5_model, tree = 1, data = penguins)
 
   output <- capture.output(print(p))
 
@@ -156,7 +168,7 @@ test_that("as.party.C5.0 properly assigns all observations to terminal nodes", {
   penguins <- get_penguins_data()
 
   c5_model <- C50::C5.0(species ~ ., data = penguins)
-  p <- as.party(c5_model, tree = 1, data = penguins)
+  p <- lorax::as.party(c5_model, tree = 1, data = penguins)
 
   # Get fitted node IDs
   fitted_ids <- p$fitted[["(fitted)"]]
@@ -185,7 +197,7 @@ test_that("as.party.C5.0 properly routes observations through categorical splits
   wa_trees <- get_wa_trees_data()
 
   c5_model <- C50::C5.0(class ~ ., data = wa_trees)
-  p <- as.party(c5_model, tree = 1, data = wa_trees)
+  p <- lorax::as.party(c5_model, tree = 1, data = wa_trees)
 
   # Get fitted node IDs
   fitted_ids <- p$fitted[["(fitted)"]]
@@ -219,7 +231,7 @@ test_that("as.party.C5.0 handles multiway categorical splits correctly", {
   wa_trees <- get_wa_trees_data()
 
   c5_model <- C50::C5.0(class ~ ., data = wa_trees)
-  p <- as.party(c5_model, tree = 1, data = wa_trees)
+  p <- lorax::as.party(c5_model, tree = 1, data = wa_trees)
 
   # Find a node with a categorical split
   # The root's right child (precip_annual >= 418) should split on county
@@ -259,9 +271,9 @@ test_that("extract_rules.C5.0() returns correct structure", {
   data <- get_factor_data()
   set.seed(327)
   c5_tree <- C50::C5.0(y ~ bill_length_mm + island + bill_depth_mm, data = data)
-  rules <- extract_rules(c5_tree, tree = 1L, data = data)
+  rules <- lorax::extract_rules(c5_tree, tree = 1L, data = data)
 
-  expect_s3_class(rules, "rule_set_party")
+  expect_s3_class(rules, "rule_set_C5.0")
   expect_s3_class(rules, "rule_set")
   expect_s3_class(rules, "tbl_df")
   expect_named(rules, c("id", "rules", "tree"))
@@ -276,7 +288,7 @@ test_that("extract_rules.C5.0() extracts from single tree", {
   data <- get_factor_data()
   set.seed(438)
   c5_tree <- C50::C5.0(y ~ bill_length_mm + island + bill_depth_mm, data = data)
-  rules <- extract_rules(c5_tree, tree = 1L, data = data)
+  rules <- lorax::extract_rules(c5_tree, tree = 1L, data = data)
 
   expect_equal(unique(rules$tree), 1L)
   expect_true(nrow(rules) > 0)
@@ -303,7 +315,7 @@ test_that("extract_rules.C5.0() extracts from multiple boosted trees", {
   trees_to_extract <- min(3L, n_trials)
   tree_nums <- seq_len(trees_to_extract)
 
-  rules <- extract_rules(c5_boost, tree = tree_nums, data = wa_trees)
+  rules <- lorax::extract_rules(c5_boost, tree = tree_nums, data = wa_trees)
 
   expect_equal(sort(unique(rules$tree)), tree_nums)
   expect_true(nrow(rules) > 0)
@@ -340,15 +352,23 @@ test_that("extract_rules.C5.0() validates tree argument", {
   )
 })
 
-test_that("extract_rules.C5.0() requires data parameter", {
+test_that("extract_rules.C5.0() works without data parameter", {
   skip_if_not_installed("C50")
 
   data <- get_factor_data()
   set.seed(762)
   c5_tree <- C50::C5.0(y ~ bill_length_mm + island + bill_depth_mm, data = data)
 
-  expect_snapshot(extract_rules(c5_tree, tree = 1L), error = TRUE)
-  expect_snapshot(extract_rules(c5_tree, tree = 1L, data = NULL), error = TRUE)
+  # Data parameter is now optional
+  rules_no_data <- lorax::extract_rules(c5_tree, tree = 1L)
+  expect_s3_class(rules_no_data, "rule_set_C5.0")
+
+  # Should also work with data for backward compatibility
+  rules_with_data <- lorax::extract_rules(c5_tree, tree = 1L, data = data)
+  expect_s3_class(rules_with_data, "rule_set_C5.0")
+
+  # Results should be identical
+  expect_equal(rules_no_data$rules, rules_with_data$rules)
 })
 
 test_that("extract_rules.C5.0() works with numeric predictors", {
@@ -358,9 +378,9 @@ test_that("extract_rules.C5.0() works with numeric predictors", {
   data$y_cat <- cut(data$y, breaks = 3, labels = c("low", "med", "high"))
   set.seed(873)
   c5_tree <- C50::C5.0(y_cat ~ predictor_01 + predictor_02, data = data)
-  rules <- extract_rules(c5_tree, tree = 1L, data = data)
+  rules <- lorax::extract_rules(c5_tree, tree = 1L, data = data)
 
-  expect_s3_class(rules, "rule_set_party")
+  expect_s3_class(rules, "rule_set_C5.0")
   expect_true(nrow(rules) > 0)
 })
 
@@ -370,9 +390,9 @@ test_that("extract_rules.C5.0() works with factor predictors", {
   data <- get_factor_data()
   set.seed(984)
   c5_tree <- C50::C5.0(y ~ island + sex, data = data)
-  rules <- extract_rules(c5_tree, tree = 1L, data = data)
+  rules <- lorax::extract_rules(c5_tree, tree = 1L, data = data)
 
-  expect_s3_class(rules, "rule_set_party")
+  expect_s3_class(rules, "rule_set_C5.0")
   expect_true(nrow(rules) > 0)
 })
 
@@ -382,9 +402,9 @@ test_that("extract_rules.C5.0() works with mixed predictors", {
   wa_trees <- get_wa_trees_data()[1:200, ]
   set.seed(195)
   c5_tree <- C50::C5.0(class ~ elevation + county, data = wa_trees)
-  rules <- extract_rules(c5_tree, tree = 1L, data = wa_trees)
+  rules <- lorax::extract_rules(c5_tree, tree = 1L, data = wa_trees)
 
-  expect_s3_class(rules, "rule_set_party")
+  expect_s3_class(rules, "rule_set_C5.0")
   expect_true(nrow(rules) > 0)
 })
 
@@ -404,7 +424,7 @@ test_that("extract_rules.C5.0() rules are sorted by tree then id", {
   trees_to_extract <- min(3L, n_trials)
   tree_nums <- rev(seq_len(trees_to_extract))
 
-  rules <- extract_rules(c5_boost, tree = tree_nums, data = wa_trees)
+  rules <- lorax::extract_rules(c5_boost, tree = tree_nums, data = wa_trees)
 
   # Check sorting
   expect_true(all(diff(rules$tree) >= 0))
@@ -426,14 +446,18 @@ test_that("extract_rules.C5.0() handles duplicate tree numbers", {
   # Get actual number of trials
   n_trials <- c5_boost$trials["Actual"]
   if (n_trials >= 2) {
-    rules <- extract_rules(c5_boost, tree = c(1L, 1L, 2L), data = wa_trees)
+    rules <- lorax::extract_rules(
+      c5_boost,
+      tree = c(1L, 1L, 2L),
+      data = wa_trees
+    )
 
     # Should have results for tree 1 twice
     tree_counts <- table(rules$tree)
     expect_equal(as.numeric(names(tree_counts)), c(1, 2))
   } else {
     # If only 1 trial, just test with duplicates of that
-    rules <- extract_rules(c5_boost, tree = c(1L, 1L), data = wa_trees)
+    rules <- lorax::extract_rules(c5_boost, tree = c(1L, 1L), data = wa_trees)
     tree_counts <- table(rules$tree)
     expect_equal(as.numeric(names(tree_counts)), 1)
   }
@@ -452,7 +476,7 @@ test_that("extract_rules.C5.0() works with all trees", {
 
   # Get actual number of trials
   n_trials <- c5_boost$trials["Actual"]
-  rules <- extract_rules(c5_boost, tree = 1:n_trials, data = wa_trees)
+  rules <- lorax::extract_rules(c5_boost, tree = 1:n_trials, data = wa_trees)
 
   expect_equal(sort(unique(rules$tree)), 1:n_trials)
   expect_true(nrow(rules) > 0)
@@ -469,9 +493,9 @@ test_that("extract_rules.C5.0() handles tree with no valid splits", {
   )
   set.seed(549)
   c5_tree <- C50::C5.0(y ~ ., data = null_data)
-  rules <- extract_rules(c5_tree, tree = 1L, data = null_data)
+  rules <- lorax::extract_rules(c5_tree, tree = 1L, data = null_data)
 
-  expect_s3_class(rules, "rule_set_party")
+  expect_s3_class(rules, "rule_set_C5.0")
   # Even with no splits, should return at least one rule (TRUE)
   expect_true(nrow(rules) >= 1)
 })
@@ -483,7 +507,7 @@ test_that("extract_rules.C5.0() works with rule-based models", {
   set.seed(123)
   c5_rules <- C50::C5.0(species ~ ., data = data, rules = TRUE)
 
-  rules <- extract_rules(c5_rules, tree = 1L, data = data)
+  rules <- lorax::extract_rules(c5_rules, tree = 1L, data = data)
 
   expect_s3_class(rules, "rule_set_C5.0")
   expect_s3_class(rules, "rule_set")
@@ -498,7 +522,7 @@ test_that("extract_rules.C5.0() handles boosted rule models", {
   set.seed(456)
   c5_rules <- C50::C5.0(species ~ ., data = data, rules = TRUE, trials = 3)
 
-  rules <- extract_rules(c5_rules, tree = c(1L, 2L), data = data)
+  rules <- lorax::extract_rules(c5_rules, tree = c(1L, 2L), data = data)
 
   expect_equal(sort(unique(rules$tree)), c(1L, 2L))
   expect_true(all(rules$tree %in% c(1L, 2L)))
@@ -511,7 +535,7 @@ test_that("extract_rules.C5.0() parses rule conditions correctly", {
   set.seed(789)
   c5_rules <- C50::C5.0(species ~ ., data = data, rules = TRUE)
 
-  rules <- extract_rules(c5_rules, tree = 1L, data = data)
+  rules <- lorax::extract_rules(c5_rules, tree = 1L, data = data)
 
   # Check that rules contain expected variable names
   rule_text_all <- sapply(rules$rules, function(r) deparse1(r))
@@ -530,7 +554,7 @@ test_that("extract_rules.C5.0() handles rule models with numeric predictors", {
   set.seed(111)
   c5_rules <- C50::C5.0(vs ~ mpg + hp + wt, data = mtcars_factor, rules = TRUE)
 
-  rules <- extract_rules(c5_rules, tree = 1L, data = mtcars_factor)
+  rules <- lorax::extract_rules(c5_rules, tree = 1L, data = mtcars_factor)
 
   expect_s3_class(rules, "rule_set_C5.0")
   expect_true(nrow(rules) > 0)
@@ -574,7 +598,7 @@ test_that("extract_rules.C5.0() handles categorical predictors in rules", {
   )
 
   c5_rules <- C50::C5.0(outcome ~ ., data = test_data, rules = TRUE)
-  rules <- extract_rules(c5_rules, tree = 1L, data = test_data)
+  rules <- lorax::extract_rules(c5_rules, tree = 1L, data = test_data)
 
   expect_s3_class(rules, "rule_set_C5.0")
   expect_true(nrow(rules) >= 0)
@@ -614,7 +638,7 @@ test_that("extract_rules.C5.0() handles rules with only categorical predictors",
   )
 
   c5_rules <- C50::C5.0(outcome ~ ., data = test_data, rules = TRUE)
-  rules <- extract_rules(c5_rules, tree = 1L, data = test_data)
+  rules <- lorax::extract_rules(c5_rules, tree = 1L, data = test_data)
 
   expect_s3_class(rules, "rule_set_C5.0")
   # C5.0 might produce 0 rules if no good splits found
@@ -643,7 +667,7 @@ test_that("extract_rules.C5.0() handles mixed numeric and categorical conditions
   set.seed(555)
   c5_rules <- C50::C5.0(class ~ ., data = wa_trees, rules = TRUE)
 
-  rules <- extract_rules(c5_rules, tree = 1L, data = wa_trees)
+  rules <- lorax::extract_rules(c5_rules, tree = 1L, data = wa_trees)
 
   expect_s3_class(rules, "rule_set_C5.0")
   expect_true(nrow(rules) > 0)
@@ -661,6 +685,104 @@ test_that("extract_rules.C5.0() handles mixed numeric and categorical conditions
   expect_true(has_numeric || has_categorical)
 })
 
+# Internal C5.0 tree parsing functions ----------------------------------------
+
+test_that("c5_parse_tree_line correctly parses tree lines", {
+  # Test parsing various line formats
+  line1 <- 'type="2" class="setosa" freq="50,50,50" att="Petal.Length" forks="2" cut="2.45"'
+  parsed1 <- lorax:::c5_parse_tree_line(line1)
+
+  expect_equal(parsed1$type, 2L)
+  expect_equal(parsed1$class, "setosa")
+  expect_equal(parsed1$freq, "50,50,50")
+  expect_equal(parsed1$att, "Petal.Length")
+  expect_equal(parsed1$forks, 2L)
+  expect_equal(parsed1$cut, 2.45)
+
+  # Test leaf node
+  line2 <- 'type="0" class="versicolor" freq="0,47,3"'
+  parsed2 <- lorax:::c5_parse_tree_line(line2)
+
+  expect_equal(parsed2$type, 0L)
+  expect_equal(parsed2$class, "versicolor")
+  expect_equal(parsed2$freq, "0,47,3")
+})
+
+test_that("c5_build_tree_structure builds correct hierarchy", {
+  # Simple tree with one split and two leaves
+  lines <- c(
+    'type="2" class="setosa" freq="50,50,50" att="Petal.Length" forks="2" cut="2.45"',
+    'type="0" class="setosa" freq="50,0,0"',
+    'type="0" class="versicolor" freq="0,50,50"'
+  )
+
+  parsed_lines <- lapply(lines, lorax:::c5_parse_tree_line)
+  tree <- lorax:::c5_build_tree_structure(parsed_lines)
+
+  # Root should have type 2 and 2 children
+  expect_equal(tree$properties$type, 2L)
+  expect_length(tree$children, 2)
+
+  # Children should be leaf nodes
+  expect_equal(tree$children[[1]]$properties$type, 0L)
+  expect_equal(tree$children[[2]]$properties$type, 0L)
+})
+
+test_that("c5_extract_tree_paths extracts correct paths", {
+  # Build a simple tree
+  lines <- c(
+    'type="2" class="setosa" freq="50,50,50" att="Petal.Length" forks="2" cut="2.45"',
+    'type="0" class="setosa" freq="50,0,0"',
+    'type="0" class="versicolor" freq="0,50,50"'
+  )
+
+  parsed_lines <- lapply(lines, lorax:::c5_parse_tree_line)
+  tree <- lorax:::c5_build_tree_structure(parsed_lines)
+  paths <- lorax:::c5_extract_tree_paths(tree)
+
+  # Should have 2 paths (one to each leaf)
+  expect_length(paths, 2)
+
+  # First path: Petal.Length <= 2.45 -> setosa
+  expect_length(paths[[1]]$conditions, 1)
+  expect_equal(paths[[1]]$conditions[[1]]$att, "Petal.Length")
+  expect_equal(paths[[1]]$conditions[[1]]$op, "<=")
+  expect_equal(paths[[1]]$class, "setosa")
+
+  # Second path: Petal.Length > 2.45 -> versicolor
+  expect_length(paths[[2]]$conditions, 1)
+  expect_equal(paths[[2]]$conditions[[1]]$att, "Petal.Length")
+  expect_equal(paths[[2]]$conditions[[1]]$op, ">")
+  expect_equal(paths[[2]]$class, "versicolor")
+})
+
+test_that("c5_extract_rules_from_tree handles different node types", {
+  # Simple tree with threshold split
+  tree_text <- 'id="See5/C5.0 2.07 GPL Edition 2024-01-01"
+entries="1"
+type="2" class="setosa" freq="50,50,50" att="Petal.Length" forks="2" cut="2.45"
+type="0" class="setosa" freq="50,0,0"
+type="0" class="versicolor" freq="0,50,50"'
+
+  # Parse the tree
+  result <- lorax:::c5_extract_rules_from_tree(tree_text, trial = 1L)
+
+  # Should have 2 rules (one for each leaf)
+  expect_equal(nrow(result), 2)
+
+  # Check first rule (Petal.Length <= 2.45)
+  rule1 <- result$rules[[1]]
+  expect_true(is.language(rule1))
+  expect_true(grepl("Petal.Length", deparse(rule1)))
+  expect_true(grepl("<=", deparse(rule1)))
+
+  # Check second rule (Petal.Length > 2.45)
+  rule2 <- result$rules[[2]]
+  expect_true(is.language(rule2))
+  expect_true(grepl("Petal.Length", deparse(rule2)))
+  expect_true(grepl(">", deparse(rule2)))
+})
+
 # active_predictors() tests ---------------------------------------------------
 
 test_that("active_predictors.C5.0() has correct structure for tree models", {
@@ -669,7 +791,7 @@ test_that("active_predictors.C5.0() has correct structure for tree models", {
   penguins <- get_penguins_data()
 
   fit <- C50::C5.0(species ~ ., data = penguins)
-  result <- active_predictors(fit)
+  result <- lorax::active_predictors(fit)
 
   expect_s3_class(result, "tbl_df")
   expect_named(result, c("active_predictors", "tree"))
@@ -683,7 +805,7 @@ test_that("active_predictors.C5.0() extracts correct variables", {
   penguins <- get_penguins_data()
 
   fit <- C50::C5.0(species ~ ., data = penguins)
-  result <- active_predictors(fit)
+  result <- lorax::active_predictors(fit)
 
   active_vars <- result$active_predictors[[1]]
   expect_true(is.character(active_vars))
@@ -699,7 +821,7 @@ test_that("active_predictors.C5.0() works with numeric predictors", {
   mtcars_factor <- mtcars
   mtcars_factor$vs <- factor(mtcars_factor$vs)
   fit <- C50::C5.0(vs ~ mpg + hp + wt, data = mtcars_factor)
-  result <- active_predictors(fit)
+  result <- lorax::active_predictors(fit)
 
   expect_s3_class(result, "tbl_df")
   active_vars <- result$active_predictors[[1]]
@@ -712,7 +834,7 @@ test_that("active_predictors.C5.0() works with factor predictors", {
   wa_trees <- get_wa_trees_data()
 
   fit <- C50::C5.0(class ~ county + roughness, data = wa_trees)
-  result <- active_predictors(fit)
+  result <- lorax::active_predictors(fit)
 
   active_vars <- result$active_predictors[[1]]
   expect_true(all(active_vars %in% c("county", "roughness")))
@@ -728,7 +850,7 @@ test_that("active_predictors.C5.0() handles tree with no splits", {
     y = factor(c("a", "a", "a"))
   )
   fit <- C50::C5.0(y ~ x, data = small_data)
-  result <- active_predictors(fit)
+  result <- lorax::active_predictors(fit)
 
   expect_s3_class(result, "tbl_df")
   # A tree with no splits will have an empty active_predictors list
@@ -742,7 +864,7 @@ test_that("active_predictors.C5.0() extracts from single boosted tree", {
   penguins <- get_penguins_data()
 
   fit <- C50::C5.0(species ~ ., data = penguins, trials = 5)
-  result <- active_predictors(fit, tree = 2L)
+  result <- lorax::active_predictors(fit, tree = 2L)
 
   expect_equal(nrow(result), 1)
   expect_equal(result$tree, 2L)
@@ -756,7 +878,7 @@ test_that("active_predictors.C5.0() extracts from multiple boosted trees", {
   penguins <- get_penguins_data()
 
   fit <- C50::C5.0(species ~ ., data = penguins, trials = 5)
-  result <- active_predictors(fit, tree = c(1L, 2L, 3L))
+  result <- lorax::active_predictors(fit, tree = c(1L, 2L, 3L))
 
   expect_equal(nrow(result), 3)
   expect_equal(result$tree, c(1L, 2L, 3L))
@@ -771,7 +893,7 @@ test_that("active_predictors.C5.0() extracts from all trees", {
 
   fit <- C50::C5.0(species ~ ., data = penguins, trials = 5)
   num_trials <- as.integer(fit$trials["Actual"])
-  result <- active_predictors(fit, tree = 1:num_trials)
+  result <- lorax::active_predictors(fit, tree = 1:num_trials)
 
   expect_equal(nrow(result), num_trials)
   expect_equal(result$tree, 1:num_trials)
@@ -784,10 +906,10 @@ test_that("active_predictors.C5.0() validates tree argument", {
 
   fit <- C50::C5.0(species ~ ., data = penguins, trials = 3)
 
-  expect_snapshot(active_predictors(fit, tree = "1"), error = TRUE)
-  expect_snapshot(active_predictors(fit, tree = 1.5), error = TRUE)
-  expect_snapshot(active_predictors(fit, tree = 0L), error = TRUE)
-  expect_snapshot(active_predictors(fit, tree = 10L), error = TRUE)
+  expect_snapshot(lorax::active_predictors(fit, tree = "1"), error = TRUE)
+  expect_snapshot(lorax::active_predictors(fit, tree = 1.5), error = TRUE)
+  expect_snapshot(lorax::active_predictors(fit, tree = 0L), error = TRUE)
+  expect_snapshot(lorax::active_predictors(fit, tree = 10L), error = TRUE)
 })
 
 test_that("active_predictors.C5.0() returns sorted unique results", {
@@ -796,7 +918,7 @@ test_that("active_predictors.C5.0() returns sorted unique results", {
   penguins <- get_penguins_data()
 
   fit <- C50::C5.0(species ~ ., data = penguins)
-  result <- active_predictors(fit)
+  result <- lorax::active_predictors(fit)
 
   active_vars <- result$active_predictors[[1]]
   expect_equal(active_vars, active_vars[order(tolower(active_vars))])
@@ -809,7 +931,7 @@ test_that("active_predictors.C5.0() automatically deduplicates tree numbers", {
   penguins <- get_penguins_data()
 
   fit <- C50::C5.0(species ~ ., data = penguins, trials = 3)
-  result <- active_predictors(fit, tree = c(1L, 1L, 2L))
+  result <- lorax::active_predictors(fit, tree = c(1L, 1L, 2L))
 
   expect_equal(nrow(result), 2)
   expect_equal(result$tree, c(1L, 2L))
@@ -821,7 +943,7 @@ test_that("active_predictors.C5.0() rule model has correct structure", {
   penguins <- get_penguins_data()
 
   fit <- C50::C5.0(species ~ ., data = penguins, rules = TRUE)
-  result <- active_predictors(fit)
+  result <- lorax::active_predictors(fit)
 
   expect_s3_class(result, "tbl_df")
   expect_named(result, "active_predictors")
@@ -836,7 +958,7 @@ test_that("active_predictors.C5.0() extracts from rule conditions", {
   penguins <- get_penguins_data()
 
   fit <- C50::C5.0(species ~ ., data = penguins, rules = TRUE)
-  result <- active_predictors(fit)
+  result <- lorax::active_predictors(fit)
 
   active_vars <- result$active_predictors[[1]]
   expect_true(is.character(active_vars))
@@ -850,7 +972,7 @@ test_that("active_predictors.C5.0() rule model with factors", {
   wa_trees <- get_wa_trees_data()
 
   fit <- C50::C5.0(class ~ county + roughness, data = wa_trees, rules = TRUE)
-  result <- active_predictors(fit)
+  result <- lorax::active_predictors(fit)
 
   active_vars <- result$active_predictors[[1]]
   expect_true(all(active_vars %in% c("county", "roughness")))
@@ -864,7 +986,7 @@ test_that("active_predictors.C5.0() rule model with numerics", {
   mtcars_factor <- mtcars
   mtcars_factor$vs <- factor(mtcars_factor$vs)
   fit <- C50::C5.0(vs ~ mpg + hp + wt, data = mtcars_factor, rules = TRUE)
-  result <- active_predictors(fit)
+  result <- lorax::active_predictors(fit)
 
   active_vars <- result$active_predictors[[1]]
   expect_true(all(active_vars %in% c("mpg", "hp", "wt")))
@@ -877,7 +999,7 @@ test_that("extract_rules.C5.0() works with single numeric predictor", {
   # C5.0 requires factor outcome - use single numeric data
   data <- get_single_numeric_data()
   model <- C50::C5.0(y ~ x, data = data, trials = 3)
-  rules <- extract_rules(model, tree = 1, data = data)
+  rules <- lorax::extract_rules(model, tree = 1, data = data)
 
   expect_s3_class(rules, "rule_set")
   expect_true(nrow(rules) > 0)
@@ -894,7 +1016,7 @@ test_that("extract_rules.C5.0() works with single factor predictor", {
 
   data <- get_single_factor_data()
   model <- C50::C5.0(y ~ x, data = data, trials = 3)
-  rules <- extract_rules(model, tree = 1, data = data)
+  rules <- lorax::extract_rules(model, tree = 1, data = data)
 
   expect_s3_class(rules, "rule_set")
   expect_true(nrow(rules) > 0)
@@ -913,7 +1035,7 @@ test_that("active_predictors.C5.0() works with single numeric predictor", {
   data <- get_single_factor_data()
   data$x_num <- rnorm(nrow(data))
   model <- C50::C5.0(y ~ x_num, data = data, trials = 3)
-  active <- active_predictors(model)
+  active <- lorax::active_predictors(model)
 
   expect_s3_class(active, "tbl_df")
   # If the model made splits, should have "x_num" as active predictor
@@ -928,7 +1050,7 @@ test_that("active_predictors.C5.0() works with single factor predictor", {
 
   data <- get_single_factor_data()
   model <- C50::C5.0(y ~ x, data = data, trials = 3)
-  active <- active_predictors(model)
+  active <- lorax::active_predictors(model)
 
   expect_s3_class(active, "tbl_df")
   # If the model made splits, should have "x" as active predictor
diff --git a/tests/testthat/test-cforest.R b/tests/testthat/test-cforest.R
index 83481a5..9bf8af8 100644
--- a/tests/testthat/test-cforest.R
+++ b/tests/testthat/test-cforest.R
@@ -1,4 +1,7 @@
 # Tests for extract_rules.cforest() --------------------------------------------
+skip_on_cran()
+skip_if(runif(1) <= 0.1)
+
 
 test_that("extract_rules.cforest() returns correct structure", {
   skip_if_not_installed("partykit")
diff --git a/tests/testthat/test-cubist.R b/tests/testthat/test-cubist.R
new file mode 100644
index 0000000..fd48510
--- /dev/null
+++ b/tests/testthat/test-cubist.R
@@ -0,0 +1,430 @@
+test_that("extract_rules.cubist returns correct structure", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Use Ames housing data - real estate pricing with meaningful predictors
+  data <- get_ames_data(n = 200)
+  mod <- cubist(data$x, data$y)
+
+  rules <- extract_rules(mod)
+
+  # Check class
+  expect_s3_class(rules, "rule_set_cubist")
+  expect_s3_class(rules, "rule_set")
+  expect_s3_class(rules, "tbl_df")
+
+  # Check columns
+  expect_named(rules, c("committee", "id", "rules"))
+
+  # Check column types
+  expect_type(rules$committee, "integer")
+  expect_type(rules$id, "integer")
+  expect_type(rules$rules, "list")
+
+  # Check that rules are expressions
+  for (rule in rules$rules) {
+    expect_true(is.language(rule) || is.logical(rule))
+  }
+})
+
+test_that("extract_rules.cubist handles committee parameter", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Use Sacramento housing data with multiple committees
+  data <- get_sacramento_data()
+  mod <- cubist(data$x, data$y, committees = 3)
+
+  # Default should be first committee
+  rules1 <- extract_rules(mod)
+  expect_equal(unique(rules1$committee), 1L)
+
+  # Extract specific committee
+  rules2 <- extract_rules(mod, committee = 2L)
+  expect_equal(unique(rules2$committee), 2L)
+
+  # Extract multiple committees
+  rules_multi <- extract_rules(mod, committee = c(1L, 3L))
+  expect_equal(sort(unique(rules_multi$committee)), c(1L, 3L))
+
+  # Check ordering
+  expect_true(all(diff(rules_multi$committee) >= 0))
+
+  # Invalid committee should error
+  expect_error(
+    extract_rules(mod, committee = 0L),
+    "must be between 1 and"
+  )
+  expect_error(
+    extract_rules(mod, committee = 4L),
+    "must be between 1 and"
+  )
+  expect_error(
+    extract_rules(mod, committee = "1"),
+    "must be an integer vector"
+  )
+})
+
+test_that("extract_rules.cubist handles continuous splits correctly", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Use structured data with clear continuous splits
+  data <- get_structured_regression_data(n = 100)
+  mod <- cubist(data$x, data$y)
+  rules <- extract_rules(mod)
+
+  # Should have at least one rule
+  expect_gt(nrow(rules), 0)
+
+  # Check that expressions can be evaluated
+  if (nrow(rules) > 0 && !isTRUE(rules$rules[[1]])) {
+    expr <- rules$rules[[1]]
+    result <- eval(expr, data$x)
+    expect_type(result, "logical")
+    expect_length(result, nrow(data$x))
+  }
+
+  # Verify that rules contain expected variables
+  # With structured data, we expect splits on x1 and/or x2
+  rule_text_all <- sapply(rules$rules, deparse)
+  expect_true(
+    any(grepl("x1", rule_text_all)) || any(grepl("x2", rule_text_all))
+  )
+})
+
+test_that("extract_rules.cubist handles categorical splits correctly", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Use Ames data which has factor predictors like Neighborhood and Bldg_Type
+  data <- get_ames_data(n = 200)
+  mod <- cubist(data$x, data$y)
+
+  rules <- extract_rules(mod)
+
+  # Should have rules
+  expect_gt(nrow(rules), 0)
+
+  # Check if we can evaluate the rules
+  for (i in seq_len(min(3, nrow(rules)))) {
+    expr <- rules$rules[[i]]
+    if (!isTRUE(expr)) {
+      result <- eval(expr, data$x)
+      expect_type(result, "logical")
+      expect_length(result, nrow(data$x))
+    }
+  }
+
+  # Check for categorical splits in the rules
+  # Ames data has factor predictors that should appear in some rules
+  rule_text_all <- paste(sapply(rules$rules, deparse), collapse = " ")
+  has_categorical <- grepl("==", rule_text_all) || grepl("%in%", rule_text_all)
+  # May or may not have categorical splits depending on the model
+  # Just verify it doesn't error
+})
+
+test_that("extract_rules.cubist handles rules with no conditions", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Use simple data that may produce rules with no conditions
+  data <- get_simple_cubist_data(n = 10)
+  mod <- cubist(data$x, data$y, control = cubistControl(rules = 1))
+  rules <- extract_rules(mod)
+
+  # Should have at least one rule
+  expect_gt(nrow(rules), 0)
+
+  # If there's a rule with no conditions, it should be TRUE
+  for (rule in rules$rules) {
+    if (is.logical(rule) && length(rule) == 1) {
+      expect_true(rule)
+    }
+  }
+})
+
+test_that("extract_rules.cubist works with multiple committees", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Use Ames data with multiple committees for ensemble modeling
+  data <- get_ames_data(n = 300)
+  mod <- cubist(data$x, data$y, committees = 5)
+
+  # Extract all committees
+  rules_all <- extract_rules(mod, committee = 1:5)
+
+  # Should have rules from all committees
+  expect_equal(sort(unique(rules_all$committee)), 1:5)
+
+  # Each committee should have at least one rule
+  for (com in 1:5) {
+    com_rules <- rules_all[rules_all$committee == com, ]
+    expect_gt(nrow(com_rules), 0)
+  }
+
+  # Rule IDs should restart for each committee
+  for (com in 1:5) {
+    com_rules <- rules_all[rules_all$committee == com, ]
+    expect_equal(min(com_rules$id), 1L)
+  }
+})
+
+test_that("active_predictors.cubist returns correct structure", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Use Sacramento data - well-understood housing predictors
+  data <- get_sacramento_data()
+  mod <- cubist(data$x, data$y)
+
+  active <- active_predictors(mod)
+
+  # Check structure
+  expect_s3_class(active, "tbl_df")
+  expect_named(active, "active_predictors")
+
+  # Check content
+  expect_type(active$active_predictors, "list")
+  expect_length(active$active_predictors, 1)
+
+  vars <- active$active_predictors[[1]]
+  expect_type(vars, "character")
+
+  # Should be subset of original predictors
+  all_vars <- colnames(data$x)
+  expect_true(all(vars %in% all_vars))
+
+  # Should be sorted alphabetically (case-insensitive)
+  expect_equal(vars, vars[order(tolower(vars))])
+})
+
+test_that("active_predictors.cubist identifies correct variables", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Use structured data where we know which predictors are important
+  data <- get_structured_regression_data(n = 200)
+  mod <- cubist(data$x, data$y)
+
+  active <- active_predictors(mod)
+  vars <- active$active_predictors[[1]]
+
+  # With structured data, x1 and x2 should definitely be included
+  # as they directly determine the outcome
+  expect_true("x1" %in% vars)
+  expect_true("x2" %in% vars)
+
+  # x3 might be included depending on the model
+  # All variables should be from the original set
+  expect_true(all(vars %in% c("x1", "x2", "x3")))
+})
+
+test_that("active_predictors.cubist works with real estate data", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Use Ames housing data
+  data <- get_ames_data(n = 200)
+  mod <- cubist(data$x, data$y)
+
+  active <- active_predictors(mod)
+  vars <- active$active_predictors[[1]]
+
+  # Should include important predictors
+  expect_gt(length(vars), 0)
+
+  # All should be valid column names
+  expect_true(all(vars %in% colnames(data$x)))
+
+  # Likely to include key housing predictors
+  # At minimum, living area and overall condition should be important
+  important_vars <- c("Gr_Liv_Area", "Overall_Cond")
+  expect_true(any(important_vars %in% vars))
+})
+
+test_that("active_predictors.cubist works with multiple committees", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Use Sacramento data with multiple committees
+  data <- get_sacramento_data()
+  mod <- cubist(data$x, data$y, committees = 5)
+
+  active <- active_predictors(mod)
+
+  # Should still return single row (aggregated across committees)
+  expect_equal(nrow(active), 1)
+
+  vars <- active$active_predictors[[1]]
+  expect_type(vars, "character")
+
+  # Variables should be from the original set
+  all_vars <- colnames(data$x)
+  expect_true(all(vars %in% all_vars))
+
+  # With real estate data, we expect most variables to be used
+  expect_gte(length(vars), 3) # At least beds, baths, sqft
+})
+
+test_that("extract_rules.cubist integrates with rule_text", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Use Ames data for realistic rules
+  data <- get_ames_data(n = 150)
+  mod <- cubist(data$x, data$y)
+
+  rules <- extract_rules(mod)
+
+  # Should be able to convert to text
+  for (i in seq_len(min(3, nrow(rules)))) {
+    rule_expr <- rules$rules[[i]]
+    if (!isTRUE(rule_expr)) {
+      text <- rule_text(rule_expr)
+      expect_type(text, "character")
+      expect_length(text, 1)
+      expect_gt(nchar(text), 0)
+
+      # Text should contain variable names from the data
+      # Check that at least one predictor appears
+      any_var_present <- any(sapply(colnames(data$x), function(v) {
+        grepl(v, text)
+      }))
+      expect_true(any_var_present)
+    }
+  }
+})
+
+test_that("extract_rules.cubist handles edge case models", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Very small dataset - edge case
+  data <- get_simple_cubist_data(n = 5)
+  mod <- cubist(data$x, data$y)
+  rules <- extract_rules(mod)
+
+  # Should still work
+  expect_s3_class(rules, "rule_set_cubist")
+  expect_gt(nrow(rules), 0)
+
+  # High dimensional data (p > n) - edge case
+  data2 <- get_high_dim_cubist_data(n = 20, p = 50)
+  mod2 <- cubist(data2$x, data2$y)
+  rules2 <- extract_rules(mod2)
+
+  expect_s3_class(rules2, "rule_set_cubist")
+  # Model should still produce some rules
+  expect_gt(nrow(rules2), 0)
+})
+
+test_that("extract_rules and active_predictors are consistent", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Use Sacramento data for consistency check
+  data <- get_sacramento_data()
+  mod <- cubist(data$x, data$y)
+
+  rules <- extract_rules(mod)
+  active <- active_predictors(mod)
+
+  active_vars <- active$active_predictors[[1]]
+
+  # Extract variables mentioned in rules
+  rule_vars <- character()
+  for (rule_expr in rules$rules) {
+    if (!isTRUE(rule_expr)) {
+      # Get all symbols from the expression
+      vars_in_rule <- all.vars(rule_expr)
+      rule_vars <- c(rule_vars, vars_in_rule)
+    }
+  }
+  rule_vars <- unique(rule_vars)
+
+  # Variables in rules should be subset of active predictors
+  # (active predictors may include variables only in linear models)
+  if (length(rule_vars) > 0) {
+    expect_true(all(rule_vars %in% active_vars))
+  }
+
+  # Active predictors should include at least the rule variables
+  if (length(rule_vars) > 0) {
+    expect_gte(length(active_vars), length(rule_vars))
+  }
+})
+
+test_that("extract_rules.cubist handles Sacramento categorical predictors", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Sacramento has 'type' (Condo, Multi_Family, Residential) and 'city'
+  data <- get_sacramento_data()
+  mod <- cubist(data$x, data$y)
+
+  rules <- extract_rules(mod)
+
+  # Should be able to evaluate rules with factor predictors
+  for (i in seq_len(min(5, nrow(rules)))) {
+    expr <- rules$rules[[i]]
+    if (!isTRUE(expr)) {
+      result <- eval(expr, data$x)
+      expect_type(result, "logical")
+      expect_length(result, nrow(data$x))
+      # Should partition the data (some TRUE, some FALSE ideally)
+      # But at minimum shouldn't error
+    }
+  }
+})
+
+test_that("extract_rules.cubist handles models with no valid splits", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Use noisy data where predictors have no relationship with outcome
+  # This is the exception case where we use random data
+  data <- get_noisy_cubist_data(n = 50, p = 3)
+
+  # Force a simple model with few rules
+  mod <- cubist(data$x, data$y, control = cubistControl(rules = 1))
+  rules <- extract_rules(mod)
+
+  # Should still produce valid output
+  expect_s3_class(rules, "rule_set_cubist")
+  expect_gt(nrow(rules), 0)
+
+  # With noisy data, might produce rules with no conditions (TRUE)
+  # or very simple conditions
+  for (rule in rules$rules) {
+    # Should be evaluable
+    if (!isTRUE(rule)) {
+      result <- eval(rule, data$x)
+      expect_type(result, "logical")
+    }
+  }
+})
+
+test_that("active_predictors.cubist handles models with no active predictors", {
+  skip_if_not_installed("Cubist")
+  library(Cubist)
+
+  # Edge case: very simple model that might use no predictors
+  # (just intercept)
+  set.seed(999)
+  x <- data.frame(x1 = rep(1, 10)) # Constant predictor
+  y <- rnorm(10) # Random outcome
+
+  mod <- cubist(x, y, control = cubistControl(rules = 1))
+  active <- active_predictors(mod)
+
+  # Should still return valid structure
+  expect_s3_class(active, "tbl_df")
+  expect_named(active, "active_predictors")
+
+  vars <- active$active_predictors[[1]]
+  expect_type(vars, "character")
+  # Might be empty or contain x1
+  expect_lte(length(vars), 1)
+})