diff --git a/R/Boost.R b/R/Boost.R index d430f6b16..61db34613 100644 --- a/R/Boost.R +++ b/R/Boost.R @@ -8,23 +8,24 @@ #' @param type \code{NULL} (default). To perform a classification of discrete integer classes from factor target variable \code{(DV.train)} with a base category of 1, set to \code{(type = "CLASS")}, else for continuous \code{(DV.train)} set to \code{(type = NULL)}. #' @param depth options: (integer, NULL, "max"); \code{(depth = NULL)}(default) Specifies the \code{order} parameter in the \link{NNS.reg} routine, assigning a number of splits in the regressors, analogous to tree depth. #' @param learner.trials integer; 100 (default) Sets the number of trials to obtain an accuracy \code{threshold} level. If the number of all possible feature combinations is less than selected value, the minimum of the two values will be used. -#' @param epochs integer; \code{2*length(DV.train)} (default) Total number of feature combinations to run. +#' @param epochs integer; \code{2*length(DV.train)} (default) Number of repeated holdout re-evaluations of the learner-trial feature subsets that pass the accuracy threshold. #' @param CV.size numeric [0, 1]; \code{NULL} (default) Sets the cross-validation size. Defaults to a random value between 0.2 and 0.33 for a random sampling of the training set. #' @param balance logical; \code{FALSE} (default) Uses both up and down sampling to balance the classes. \code{type="CLASS"} required. #' @param ts.test integer; NULL (default) Sets the length of the test set for time-series data; typically \code{2*h} parameter value from \link{NNS.ARMA} or double known periods to forecast. -#' @param threshold numeric; \code{NULL} (default) Sets the \code{obj.fn} threshold to keep feature combinations. +#' @param threshold numeric [0, 1]; \code{NULL} (default) Probability supplied to \link{LPM.VaR} over the learner-trial objective distribution to determine the objective cutoff for keeping feature combinations. Defaults to 0.80 when \code{objective = "max"} and 0.20 when \code{objective = "min"}. It is not a literal objective-score cutoff. #' @param obj.fn expression; #' \code{expression( sum((predicted - actual)^2) )} (default) Sum of squared errors is the default objective function. Any \code{expression(...)} using the specific terms \code{predicted} and \code{actual} can be used. Automatically selects an accuracy measure when \code{(type = "CLASS")}. #' @param objective options: ("min", "max") \code{"max"} (default) Select whether to minimize or maximize the objective function \code{obj.fn}. -#' @param extreme logical; \code{FALSE} (default) Uses the maximum (minimum) \code{threshold} obtained from the \code{learner.trials}, rather than the upper (lower) quintile level for maximization (minimization) \code{objective}. +#' @param extreme logical; \code{FALSE} (default) Sets the \link{LPM.VaR} probability to 1 (0) for maximization (minimization) \code{objective}, i.e. the most extreme learner-trial objective value becomes the cutoff. Overrides \code{threshold}. #' @param features.only logical; \code{FALSE} (default) Returns only the final feature loadings along with the final feature frequencies. -#' @param feature.importance logical; \code{TRUE} (default) Plots the frequency of features used in the final estimate. +#' @param feature.importance logical; \code{TRUE} (default) Draws a two-panel diagnostic: the learner-trial objective distribution with its \link{LPM.VaR} cutoff, and the frequency of features used in the final estimate. #' @param pred.int numeric [0,1]; \code{NULL} (default) Returns the associated prediction intervals for the final estimate. #' @param status logical; \code{TRUE} (default) Prints status update message in console. #' @param seed Optional integer random seed used for reproducible resampling, fold construction, and stochastic fitting steps. If `NULL`, the current random-number-generator state is used. #' @param dist options:(NULL, "NNS", "L1", "L2", "FACTOR") the method of distance calculation passed to delegated \link{NNS.reg} and \link{NNS.stack} calls. \code{dist = NULL} is the default and selects the native blended NNS distance; \code{dist = "NNS"} is an explicit alias for the default. +#' @param folds integer; 5 (default) Number of cross-validation \code{folds} passed to the final \link{NNS.stack} call. #' -#' @return Returns a vector of fitted values for the dependent variable test set \code{$results}, prediction intervals \code{$pred.int}, and the final feature loadings \code{$feature.weights}, along with final feature frequencies \code{$feature.frequency}. +#' @return Returns a vector of fitted values for the dependent variable test set \code{$results}, prediction intervals \code{$pred.int}, the final feature loadings \code{$feature.weights}, final feature frequencies \code{$feature.frequency}, and (for classification) the class labels \code{$class.levels}. Classification results are numeric: a factor or character \code{DV.train} yields integer class codes with a base category of 1 (label recoverable as \code{class.levels[results]}), and a numeric \code{DV.train} yields its original numeric class values. #' #' @note #' \itemize{ @@ -44,6 +45,9 @@ #' #' ## Test accuracy #' mean(a$results == as.numeric(iris[141:150, 5])) +#' +#' ## Recover the labels +#' a$class.levels[a$results] #' } #' #' @export @@ -68,7 +72,8 @@ NNS.boost <- function(IVs.train, pred.int = NULL, status = TRUE, seed = 123L, - dist = NULL) { + dist = NULL, + folds = 5) { dist <- .nns_reg_validate_dist(dist) # --------------------------------------------------------------------------- @@ -348,6 +353,7 @@ NNS.boost <- function(IVs.train, learner.trials <- .scalar_integer(learner.trials, "learner.trials", minimum = 1L) epochs <- .scalar_integer(epochs, "epochs", minimum = 0L, allow_null = TRUE) + folds <- .scalar_integer(folds, "folds", minimum = 1L) if (!is.null(CV.size)) { if (!is.numeric(CV.size) || @@ -599,15 +605,15 @@ NNS.boost <- function(IVs.train, balanced <- .balance_training(train_x, train_y) + # Every learner trial must remain a genuine NNS.reg base learner on the + # sampled feature subset. Do not collapse multivariate subsets to an + # equal-weight synthetic X* before scoring them. fit <- suppressWarnings( NNS.reg( balanced$x, balanced$y, point.est = test_x, - dim.red.method = if (ncol(balanced$x) > 1L) - "equal" - else - NULL, + dim.red.method = NULL, plot = FALSE, residual.plot = FALSE, order = depth, @@ -672,9 +678,7 @@ NNS.boost <- function(IVs.train, } learner_count <- length(test.features) - results <- rep(NA_real_, learner_count) - train_index <- setdiff(seq_len(n_obs), validation_index) - actual <- y[validation_index] + learner.results <- rep(NA_real_, learner_count) for (i in seq_len(learner_count)) { if (status) { @@ -685,92 +689,130 @@ NNS.boost <- function(IVs.train, appendLF = FALSE) } - predicted <- .fit_subset(test.features[[i]], train_index, validation_index) - results[i] <- .score(predicted, actual) + # Cross-sectional learner trials draw a fresh validation holdout per trial + # (the same resampling geometry as the epoch stage). Time-series trials + # keep the chronological terminal block. + trial_validation_index <- if (is.null(ts.test)) { + .random_validation_index() + } else { + validation_index + } + trial_train_index <- setdiff(seq_len(n_obs), trial_validation_index) + + predicted <- .fit_subset(test.features[[i]], + trial_train_index, + trial_validation_index) + learner.results[i] <- .score(predicted, y[trial_validation_index]) } - finite_results <- which(is.finite(results)) + finite_results <- which(is.finite(learner.results)) if (!length(finite_results)) { stop("No learner trial produced a finite objective value.", call. = FALSE) } - supplied_threshold <- !is.null(threshold) - if (!supplied_threshold) { - if (extreme) { - threshold <- if (objective == "max") { - max(results[finite_results]) - } else { - min(results[finite_results]) - } - } else { - threshold <- as.numeric(stats::quantile( - results[finite_results], - probs = if (objective == "max") - 0.75 - else - 0.25, - names = FALSE, - type = 2 - )) - } - } + # The public [threshold] argument is a probability supplied to LPM.VaR over + # the learner-trial objective distribution; the objective-score cutoff is the + # distinct value LPM.VaR returns. Neither variable overwrites the other. + threshold_info <- .nns_boost_threshold( + threshold = threshold, + objective = objective, + extreme = extreme, + learner.results = learner.results[finite_results] + ) + threshold.probability <- threshold_info$probability + learner.threshold <- threshold_info$cutoff if (status) { - message(paste0( - "\nLearner Accuracy Threshold = ", - format(threshold, digits = 4, nsmall = 2), - " " - ), - appendLF = TRUE) + message( + sprintf( + "\nLearner threshold probability = %.2f; objective cutoff = %.6f", + threshold.probability, + learner.threshold + ) + ) } - passes_learner <- is.finite(results) & if (objective == "max") { - results >= threshold + passes_learner <- is.finite(learner.results) & if (objective == "max") { + learner.results >= learner.threshold } else { - results <= threshold + learner.results <= learner.threshold } reduced.test.features <- test.features[passes_learner] if (!length(reduced.test.features)) { - if (supplied_threshold) { - if (objective == "min") { - stop("No learner subset met [threshold]; increase the threshold.", - call. = FALSE) - } else { - stop("No learner subset met [threshold]; reduce the threshold.", - call. = FALSE) - } - } - + # Defensive: LPM.VaR returns a value inside the observed range, so the + # best trial always passes; guard anyway for degenerate distributions. best_index <- if (objective == "min") { - finite_results[which.min(results[finite_results])] + finite_results[which.min(learner.results[finite_results])] } else { - finite_results[which.max(results[finite_results])] + finite_results[which.max(learner.results[finite_results])] } reduced.test.features <- list(test.features[[best_index]]) } - # Preserve repeated survivor counts. Features with no survivor count retain a very - - # small probability so every predictor can still be explored during epochs. - - feature_count <- tabulate(unlist(reduced.test.features), nbins = n_features) - feature_prob <- as.numeric(feature_count) - feature_prob <- feature_prob + max(1, sum(feature_prob)) * 1e-12 - feature_prob <- feature_prob / sum(feature_prob) - # --------------------------------------------------------------------------- - # Weighted epoch stage + # Epoch stability stage + + # The learner trials establish the objective threshold and identify the + # survivor feature combinations. Epochs then repeatedly re-test only those + # survivor combinations. For ordinary cross-sectional data, each epoch uses + # a fresh holdout so the final feature frequencies measure out-of-sample + # stability rather than repeated scoring on one fixed validation sample. + # + # Exhaustive learner trials do not disable this stage: complete enumeration + # answers which combinations passed once, while epochs answer which of those + # passing combinations continue to pass under repeated holdouts. # --------------------------------------------------------------------------- keeper.features <- list() - if (!exhaustive && epochs > 0L) { + if (epochs > 0L) { keeper.features <- vector("list", epochs) + survivor_count <- length(reduced.test.features) + + # Give every survivor approximately equal re-test exposure. Randomizing the + # balanced schedule avoids ordering effects without allowing exposure count + # alone to masquerade as feature importance. + epoch_survivor_id <- rep(seq_len(survivor_count), length.out = epochs) + if (length(epoch_survivor_id) > 1L) { + epoch_survivor_id <- sample(epoch_survivor_id, + length(epoch_survivor_id), + replace = FALSE) + } + + # For time-series data, construct expanding-window chronological holdouts of + # the same size as ts.test. This makes repeated epochs genuine stability + # checks rather than repetitions of the same terminal block. + chronological_splits <- NULL + epoch_split_id <- NULL + if (!is.null(ts.test)) { + possible_blocks <- floor((n_obs - 1L) / ts.test) + block_starts <- n_obs - seq_len(possible_blocks) * ts.test + 1L + chronological_splits <- lapply(block_starts, function(start) { + validation <- seq.int(start, start + ts.test - 1L) + training <- seq_len(start - 1L) + list(train = training, validation = validation) + }) + valid_split <- vapply(chronological_splits, function(split) { + length(split$train) >= 3L && .has_all_classes(y[split$train]) + }, logical(1L)) + chronological_splits <- chronological_splits[valid_split] + if (!length(chronological_splits)) { + stop("No chronological epoch split retained enough training observations and every response class.", + call. = FALSE) + } + epoch_split_id <- rep(seq_along(chronological_splits), + length.out = epochs) + if (length(epoch_split_id) > 1L) { + epoch_split_id <- sample(epoch_split_id, + length(epoch_split_id), + replace = FALSE) + } + } for (j in seq_len(epochs)) { if (status) { @@ -781,21 +823,29 @@ NNS.boost <- function(IVs.train, appendLF = FALSE) } - k <- sample.int(n_features, 1L) - features_j <- sort(sample.int( - n_features, - size = k, - replace = FALSE, - prob = feature_prob - )) + features_j <- reduced.test.features[[epoch_survivor_id[j]]] - predicted <- .fit_subset(features_j, train_index, validation_index) - new_result <- .score(predicted, actual) + # Cross-sectional epochs draw a fresh random holdout. Time-series epochs + # cycle through expanding-window chronological holdouts. + if (is.null(ts.test)) { + epoch_validation_index <- .random_validation_index() + epoch_train_index <- setdiff(seq_len(n_obs), epoch_validation_index) + } else { + epoch_split <- chronological_splits[[epoch_split_id[j]]] + epoch_train_index <- epoch_split$train + epoch_validation_index <- epoch_split$validation + } + epoch_actual <- y[epoch_validation_index] + + predicted <- .fit_subset(features_j, + epoch_train_index, + epoch_validation_index) + new_result <- .score(predicted, epoch_actual) passes <- is.finite(new_result) && if (objective == "max") { - new_result >= threshold + new_result >= learner.threshold } else { - new_result <= threshold + new_result <= learner.threshold } keeper.features[[j]] <- if (passes) @@ -810,20 +860,14 @@ NNS.boost <- function(IVs.train, } if (!length(keeper.features)) { - if (supplied_threshold) { - if (objective == "min") { - stop("No epoch subset met [threshold]; increase the threshold.", - call. = FALSE) - } else { - stop("No epoch subset met [threshold]; reduce the threshold.", - call. = FALSE) - } - } - + warning( + "No feature combination re-passed the objective cutoff during epochs; using the best learner-trial combination. Consider a lower [threshold] probability.", + call. = FALSE + ) best_index <- if (objective == "min") { - finite_results[which.min(results[finite_results])] + finite_results[which.min(learner.results[finite_results])] } else { - finite_results[which.max(results[finite_results])] + finite_results[which.max(learner.results[finite_results])] } keeper.features <- list(test.features[[best_index]]) } @@ -832,6 +876,18 @@ NNS.boost <- function(IVs.train, names(plot.table) <- names(x) plot.table <- sort(plot.table[plot.table > 0L], decreasing = TRUE) + # Both diagnostic panels are drawn once the final feature frequencies exist, + # so they are complete before an early features.only return and cannot be + # disturbed by the final NNS.stack call (whose NNS.reg fits never plot). + if (feature.importance) { + .nns_boost_plot_diagnostics( + learner.results = learner.results[finite_results], + threshold.probability = threshold.probability, + learner.threshold = learner.threshold, + feature.frequency = plot.table + ) + } + if (features.only) { return(.NNS.out( list( @@ -846,156 +902,67 @@ NNS.boost <- function(IVs.train, # --------------------------------------------------------------------------- - # Build a training-fitted numeric design and frequency-weighted X*. - - # Categorical features are one-hot encoded with training levels only. Each active - - # dummy receives its original feature's weight, so a categorical feature's total - - # per-row contribution remains equal to that feature weight. + # Final estimate: replicate the original keeper predictors by their relative + # epoch frequencies and fit a genuine multivariate Method 1 NNS.stack. + # + # The historical scaling rule converts positive keeper-feature counts into + # integer replication factors, so a more stable feature contributes + # proportionally more columns to the final design. No synthetic scalar X* is + # constructed and the frequencies are never passed through dim.red.method: + # method = 1 with stack = FALSE keeps the final estimator a multivariate + # NNS.reg whose n.best is selected by NNS.stack's cross-validation. # --------------------------------------------------------------------------- - .numeric_design <- function(train, test) { - train_blocks <- list() - test_blocks <- list() - source_feature <- character() - - for (j in seq_along(train)) { - nm <- names(train)[j] - tr <- train[[j]] - te <- test[[j]] - - if (is.factor(tr)) { - lev <- levels(tr) - if (!length(lev)) { - stop(sprintf("Predictor [%s] has no observed levels.", nm), - call. = FALSE) - } - tr_block <- vapply(lev, function(level) - as.numeric(as.character(tr) == level), numeric(length(tr))) - te_block <- vapply(lev, function(level) - as.numeric(as.character(te) == level), numeric(length(te))) - if (is.null(dim(tr_block))) - tr_block <- matrix(tr_block, ncol = 1L) - if (is.null(dim(te_block))) - te_block <- matrix(te_block, ncol = 1L) - colnames(tr_block) <- paste0(nm, "__", make.names(lev, unique = TRUE)) - colnames(te_block) <- colnames(tr_block) - source_feature <- c(source_feature, rep(nm, length(lev))) - } else { - tr_block <- matrix(as.numeric(tr), - ncol = 1L, - dimnames = list(NULL, nm)) - te_block <- matrix(as.numeric(te), - ncol = 1L, - dimnames = list(NULL, nm)) - source_feature <- c(source_feature, nm) - } - - train_blocks[[j]] <- tr_block - test_blocks[[j]] <- te_block - } - - train_matrix <- do.call(cbind, train_blocks) - test_matrix <- do.call(cbind, test_blocks) - storage.mode(train_matrix) <- "double" - storage.mode(test_matrix) <- "double" - - list(train = train_matrix, - test = test_matrix, - source = source_feature) - } - - design <- .numeric_design(x, z) - - train_min <- apply(design$train, 2L, min) - train_max <- apply(design$train, 2L, max) - train_range <- train_max - train_min - train_range[!is.finite(train_range) | train_range == 0] <- 1 - - train_norm <- sweep(design$train, 2L, train_min, "-") - train_norm <- sweep(train_norm, 2L, train_range, "/") - test_norm <- sweep(design$test, 2L, train_min, "-") - test_norm <- sweep(test_norm, 2L, train_range, "/") - - feature_weights <- as.numeric(plot.table / sum(plot.table)) - names(feature_weights) <- names(plot.table) - coef_design <- feature_weights[design$source] - coef_design[is.na(coef_design)] <- 0 - - if (!any(coef_design > 0)) { - stop("No positive feature weights were available for the final estimate.", - call. = FALSE) - } - - xstar_train <- as.numeric(train_norm %*% coef_design) - xstar_test <- as.numeric(test_norm %*% coef_design) - - if (any(!is.finite(xstar_train)) || any(!is.finite(xstar_test))) { - stop("The final synthetic predictor contains non-finite values.", - call. = FALSE) - } - - xstar_frame <- function(v) { - data.frame(xstar = v, - xstar2 = v, - check.names = FALSE) - } - - # --------------------------------------------------------------------------- - - # Final estimate via NNS.stack Method 1 on the duplicated synthetic X*. - # - # X* is a univariate synthetic predictor. Duplicating it invokes the - # multivariate NNS.reg path used by NNS.stack so that n.best can be selected - # on the regression-point matrix without maintaining a second local search - # implementation inside NNS.boost. - # - # method = 1 uses only the NNS.reg component because NNS.boost has already - # performed feature screening and constructed the frequency-weighted X*. - # stack = FALSE prevents another dimension-reduction stage. - # optimize.threshold = FALSE preserves NNS.boost's existing 0.5 class-rounding - # convention. Passing CV.size = cv_fraction reproduces the repeated-holdout - # geometry used above. For ts.test, folds = 1 preserves the single terminal - # validation block used by the former local implementation. - # --------------------------------------------------------------------------- + feature.frequency <- plot.table + relative.frequency <- as.numeric(feature.frequency) / min(as.numeric(feature.frequency)) + replication.count <- pmax(1L, as.integer(round(relative.frequency))) + replicated_names <- rep(names(feature.frequency), times = replication.count) + replicated.train <- x[, replicated_names, drop = FALSE] + replicated.test <- z[, replicated_names, drop = FALSE] + names(replicated.train) <- make.unique(names(replicated.train), sep = ".") + names(replicated.test) <- names(replicated.train) - final_stack <- suppressWarnings( + final.stack <- suppressWarnings( NNS.stack( - IVs.train = xstar_frame(xstar_train), + IVs.train = replicated.train, DV.train = y, - IVs.test = xstar_frame(xstar_test), + IVs.test = replicated.test, type = type, obj.fn = obj.fn, objective = objective, optimize.threshold = FALSE, dist = dist, - CV.size = cv_fraction, - balance = balance, + CV.size = CV.size, + balance = FALSE, ts.test = ts.test, - folds = if (is.null(ts.test)) - 5L - else - 1L, + folds = folds, order = depth, method = 1L, stack = FALSE, pred.int = pred.int, status = status, + ncores = 1, seed = seed ) ) - estimates_code <- .sanitize_predictions(final_stack$reg, y) - pred_int_out <- final_stack$reg.pred.int + results <- final.stack$reg + pred.int.output <- final.stack$reg.pred.int + + estimates_code <- .sanitize_predictions(results, y) + pred_int_out <- pred.int.output if (is_class) { estimates_code <- pmin(pmax(estimates_code, 1L), length(class_values)) estimates_code <- as.integer(round(estimates_code)) + # Classification results keep the historical numeric coding: integer + # class codes with a base category of 1 (numeric responses recover their + # original numeric class values). $class.levels supplies the label for + # each code, so labels are always recoverable via class.levels[results]. if (response_was_numeric) { estimates <- class_values[estimates_code] } else { @@ -1016,29 +983,141 @@ NNS.boost <- function(IVs.train, estimates <- estimates_code } - if (feature.importance) { - old_par <- graphics::par(no.readonly = TRUE) - on.exit(graphics::par(old_par), add = TRUE) - - top <- tail(sort(plot.table, decreasing = FALSE), min(length(plot.table), 10L)) - label_margin <- max(graphics::strwidth(names(top), "inch") + 0.4, na.rm = TRUE) - graphics::par(mai = c(1.0, label_margin, 0.8, 0.5)) - graphics::barplot( - top, - horiz = TRUE, - col = "steelblue", - main = "Feature Frequency in Final Estimate", - xlab = "Frequency", - las = 1 - ) - } - .NNS.out( list( results = estimates, pred.int = pred_int_out, feature.weights = plot.table / sum(plot.table), - feature.frequency = plot.table + feature.frequency = plot.table, + class.levels = if (is_class) class_values else NULL + ) + ) +} + + +#' Learner threshold from the trial objective distribution +#' +#' Maps the public \code{threshold} probability (with \code{objective} and +#' \code{extreme} defaults) onto the objective-score cutoff via +#' \link{LPM.VaR}. The probability and the cutoff remain distinct values. +#' +#' @keywords internal +#' @noRd +.nns_boost_threshold <- function(threshold, objective, extreme, learner.results) { + threshold.probability <- threshold + + if (is.null(threshold.probability)) { + threshold.probability <- if (objective == "max") 0.80 else 0.20 + } + + if (extreme) { + threshold.probability <- if (objective == "max") 1 else 0 + } + + learner.threshold <- as.numeric( + LPM.VaR( + percentile = threshold.probability, + degree = 1, + x = learner.results ) ) + + list(probability = threshold.probability, cutoff = learner.threshold) +} + + +#' Panel 1: learner-trial objective distribution with its LPM.VaR cutoff +#' +#' @keywords internal +#' @noRd +.nns_boost_plot_learner_distribution <- function(learner.results, + threshold.probability, + learner.threshold) { + graphics::hist( + learner.results, + main = "Distribution of Learner Trials Objective Function", + xlab = "Objective Function", + col = "steelblue" + ) + + graphics::abline( + v = learner.threshold, + col = "red", + lty = 2, + lwd = 2 + ) + + graphics::mtext( + sprintf( + "LPM.VaR(p = %.2f) = %.4f", + threshold.probability, + learner.threshold + ), + side = 3, + col = "red" + ) + + invisible(NULL) +} + + +#' Panel 2: horizontal feature-frequency bar plot +#' +#' @keywords internal +#' @noRd +.nns_boost_plot_feature_frequency <- function(feature.frequency) { + sorted <- sort(feature.frequency, decreasing = FALSE) + + # Widen this panel's left margin so horizontal las = 1 labels are not + # clipped by the default margin. The change is scoped to this helper (and + # margins are per-panel, so the surrounding mfrow layout is untouched). + old_mai <- graphics::par("mai") + on.exit(graphics::par(mai = old_mai), add = TRUE) + label_margin <- max(graphics::strwidth(names(sorted), units = "inches") + 0.4, + na.rm = TRUE) + graphics::par(mai = c(old_mai[1L], + max(old_mai[2L], label_margin), + old_mai[3L], + old_mai[4L])) + + graphics::barplot( + sorted, + horiz = TRUE, + col = "steelblue", + main = "Feature Frequency in Final Estimate", + xlab = "Frequency", + las = 1 + ) + + invisible(NULL) +} + + +#' Two-panel NNS.boost diagnostic on one graphics device +#' +#' Saves the caller's graphics parameters, draws the learner-trial +#' distribution and the feature-frequency panels under +#' \code{par(mfrow = c(2, 1))}, and restores the original parameters only +#' after both panels are complete. +#' +#' @keywords internal +#' @noRd +.nns_boost_plot_diagnostics <- function(learner.results, + threshold.probability, + learner.threshold, + feature.frequency) { + old_par <- graphics::par(no.readonly = TRUE) + on.exit(graphics::par(old_par), add = TRUE) + + graphics::par(mfrow = c(2, 1)) + + .nns_boost_plot_learner_distribution( + learner.results = learner.results, + threshold.probability = threshold.probability, + learner.threshold = learner.threshold + ) + + .nns_boost_plot_feature_frequency(feature.frequency = feature.frequency) + + invisible(NULL) } diff --git a/man/NNS.boost.Rd b/man/NNS.boost.Rd index 3304a3280..92904e75f 100644 --- a/man/NNS.boost.Rd +++ b/man/NNS.boost.Rd @@ -24,7 +24,8 @@ NNS.boost( pred.int = NULL, status = TRUE, seed = 123L, - dist = NULL + dist = NULL, + folds = 5 ) } \arguments{ @@ -40,7 +41,7 @@ NNS.boost( \item{learner.trials}{integer; 100 (default) Sets the number of trials to obtain an accuracy \code{threshold} level. If the number of all possible feature combinations is less than selected value, the minimum of the two values will be used.} -\item{epochs}{integer; \code{2*length(DV.train)} (default) Total number of feature combinations to run.} +\item{epochs}{integer; \code{2*length(DV.train)} (default) Number of repeated holdout re-evaluations of the learner-trial feature subsets that pass the accuracy threshold.} \item{CV.size}{numeric [0, 1]; \code{NULL} (default) Sets the cross-validation size. Defaults to a random value between 0.2 and 0.33 for a random sampling of the training set.} @@ -48,18 +49,18 @@ NNS.boost( \item{ts.test}{integer; NULL (default) Sets the length of the test set for time-series data; typically \code{2*h} parameter value from \link{NNS.ARMA} or double known periods to forecast.} -\item{threshold}{numeric; \code{NULL} (default) Sets the \code{obj.fn} threshold to keep feature combinations.} +\item{threshold}{numeric [0, 1]; \code{NULL} (default) Probability supplied to \link{LPM.VaR} over the learner-trial objective distribution to determine the objective cutoff for keeping feature combinations. Defaults to 0.80 when \code{objective = "max"} and 0.20 when \code{objective = "min"}. It is not a literal objective-score cutoff.} \item{obj.fn}{expression; \code{expression( sum((predicted - actual)^2) )} (default) Sum of squared errors is the default objective function. Any \code{expression(...)} using the specific terms \code{predicted} and \code{actual} can be used. Automatically selects an accuracy measure when \code{(type = "CLASS")}.} \item{objective}{options: ("min", "max") \code{"max"} (default) Select whether to minimize or maximize the objective function \code{obj.fn}.} -\item{extreme}{logical; \code{FALSE} (default) Uses the maximum (minimum) \code{threshold} obtained from the \code{learner.trials}, rather than the upper (lower) quintile level for maximization (minimization) \code{objective}.} +\item{extreme}{logical; \code{FALSE} (default) Sets the \link{LPM.VaR} probability to 1 (0) for maximization (minimization) \code{objective}, i.e. the most extreme learner-trial objective value becomes the cutoff. Overrides \code{threshold}.} \item{features.only}{logical; \code{FALSE} (default) Returns only the final feature loadings along with the final feature frequencies.} -\item{feature.importance}{logical; \code{TRUE} (default) Plots the frequency of features used in the final estimate.} +\item{feature.importance}{logical; \code{TRUE} (default) Draws a two-panel diagnostic: the learner-trial objective distribution with its \link{LPM.VaR} cutoff, and the frequency of features used in the final estimate.} \item{pred.int}{numeric [0,1]; \code{NULL} (default) Returns the associated prediction intervals for the final estimate.} @@ -68,9 +69,11 @@ NNS.boost( \item{seed}{Optional integer random seed used for reproducible resampling, fold construction, and stochastic fitting steps. If `NULL`, the current random-number-generator state is used.} \item{dist}{options:(NULL, "NNS", "L1", "L2", "FACTOR") the method of distance calculation passed to delegated \link{NNS.reg} and \link{NNS.stack} calls. \code{dist = NULL} is the default and selects the native blended NNS distance; \code{dist = "NNS"} is an explicit alias for the default.} + +\item{folds}{integer; 5 (default) Number of cross-validation \code{folds} passed to the final \link{NNS.stack} call.} } \value{ -Returns a vector of fitted values for the dependent variable test set \code{$results}, prediction intervals \code{$pred.int}, and the final feature loadings \code{$feature.weights}, along with final feature frequencies \code{$feature.frequency}. +Returns a vector of fitted values for the dependent variable test set \code{$results}, prediction intervals \code{$pred.int}, the final feature loadings \code{$feature.weights}, final feature frequencies \code{$feature.frequency}, and (for classification) the class labels \code{$class.levels}. Classification results are numeric: a factor or character \code{DV.train} yields integer class codes with a base category of 1 (label recoverable as \code{class.levels[results]}), and a numeric \code{DV.train} yields its original numeric class values. } \description{ Ensemble method for classification using the NNS multivariate regression \link{NNS.reg} as the base learner instead of trees. @@ -92,6 +95,9 @@ Ensemble method for classification using the NNS multivariate regression \link{N ## Test accuracy mean(a$results == as.numeric(iris[141:150, 5])) + + ## Recover the labels + a$class.levels[a$results] } } diff --git a/tests/testthat/test-boost-iris-integration.R b/tests/testthat/test-boost-iris-integration.R new file mode 100644 index 000000000..c5de6ff94 --- /dev/null +++ b/tests/testthat/test-boost-iris-integration.R @@ -0,0 +1,130 @@ +# Full iris integration run mirroring the documented NNS.boost example. Heavy: +# 100 learner trials (exhaustive: 15 subsets), 100 epochs, final NNS.stack. +test_that("NNS.boost iris run has probability threshold, epochs, and a replicated multivariate final stack", { + skip_on_cran() + + stack_args <- new.env(parent = emptyenv()) + suppressMessages(trace( + NNS.stack, + where = asNamespace("NNS"), + print = FALSE, + tracer = bquote({ + assign("captured", list( + method = method, + folds = folds, + stack = stack, + n_predictors = ncol(IVs.train), + predictor_names = names(IVs.train) + ), envir = .(stack_args)) + }) + )) + on.exit(suppressMessages(untrace(NNS.stack, where = asNamespace("NNS"))), + add = TRUE) + + pdf_file <- withr::local_tempfile(fileext = ".pdf") + grDevices::pdf(pdf_file) + + set.seed(123) + msgs <- capture_messages( + a <- NNS.boost( + iris[1:140, 1:4], + iris[1:140, 5], + IVs.test = iris[141:150, 1:4], + epochs = 100, + learner.trials = 100, + type = "CLASS", + depth = NULL, + balance = TRUE, + feature.importance = TRUE, + status = TRUE + ) + ) + + grDevices::dev.off() + + # 1-3: the console states the probability and the distinct LPM.VaR cutoff, + # and never prints the obsolete status string. + expect_true(any(grepl("Learner threshold probability = 0.80", msgs, + fixed = TRUE))) + cutoff_line <- grep("objective cutoff = ", msgs, fixed = TRUE, value = TRUE) + expect_length(cutoff_line, 1L) + cutoff <- as.numeric(sub(".*objective cutoff = ([0-9.]+).*", "\\1", + cutoff_line)) + expect_true(is.finite(cutoff)) + expect_false(isTRUE(all.equal(cutoff, 0.80))) + expect_false(any(grepl("Learner Accuracy Threshold", msgs, fixed = TRUE))) + + # 4-5: both diagnostic panels were drawn (nonempty device output). + expect_true(file.exists(pdf_file)) + expect_gt(file.size(pdf_file), 0) + + # 6: iris has 15 unique feature subsets < 100 trials (exhaustive), yet + # epochs still ran. + expect_true(any(grepl("% of epochs", msgs, fixed = TRUE))) + + # 7-9: results and feature outputs. Classification results keep the + # historical numeric coding (integer codes, base category 1) and + # $class.levels recovers the label for each code. + expect_length(a$results, 10L) + expect_true(is.numeric(a$results)) + expect_true(all(a$results %in% seq_along(levels(iris$Species)))) + expect_identical(a$class.levels, levels(iris$Species)) + expect_gte(mean(a$results == as.numeric(iris[141:150, 5])), 0.8) + expect_gte(mean(a$class.levels[a$results] == + as.character(iris[141:150, 5])), 0.8) + expect_equal(sum(a$feature.weights), 1) + expect_true(all(names(a$feature.frequency) %in% names(iris)[1:4])) + expect_true(length(a$feature.frequency) >= 1L) + + # 10: the final call used replicated multivariate predictors through + # NNS.stack(method = 1, folds = 5, stack = FALSE). + captured <- get("captured", envir = stack_args) + expect_equal(as.numeric(captured$method), 1) + expect_equal(as.numeric(captured$folds), 5) + expect_false(captured$stack) + expect_gt(captured$n_predictors, 1L) + base_names <- sub("[.][0-9]+$", "", captured$predictor_names) + expect_true(all(base_names %in% names(iris)[1:4])) + expect_false(any(grepl("xstar", captured$predictor_names, fixed = TRUE))) + + # Replication counts follow the historical scaling rule. + expected_replication <- pmax(1L, as.integer(round( + as.numeric(a$feature.frequency) / min(as.numeric(a$feature.frequency)) + ))) + expect_identical(captured$n_predictors, sum(expected_replication)) +}) + +test_that("NNS.boost passes a supplied threshold to LPM.VaR as a probability", { + skip_on_cran() + + set.seed(123) + msgs <- capture_messages( + b <- NNS.boost( + iris[1:140, 1:4], + iris[1:140, 5], + IVs.test = iris[141:150, 1:4], + epochs = 100, + learner.trials = 100, + type = "CLASS", + depth = NULL, + balance = TRUE, + threshold = 0.65, + feature.importance = FALSE, + status = TRUE + ) + ) + + expect_true(any(grepl("Learner threshold probability = 0.65", msgs, + fixed = TRUE))) + + # The supplied probability produces a distinct objective cutoff: 0.65 is + # not applied as a literal accuracy cutoff. + cutoff_line <- grep("objective cutoff = ", msgs, fixed = TRUE, value = TRUE) + expect_length(cutoff_line, 1L) + cutoff <- as.numeric(sub(".*objective cutoff = ([0-9.]+).*", "\\1", + cutoff_line)) + expect_true(is.finite(cutoff)) + expect_false(isTRUE(all.equal(cutoff, 0.65))) + + expect_length(b$results, 10L) +}) diff --git a/tests/testthat/test-boost-plots.R b/tests/testthat/test-boost-plots.R new file mode 100644 index 000000000..9efb5e178 --- /dev/null +++ b/tests/testthat/test-boost-plots.R @@ -0,0 +1,131 @@ +.boost_plot_fixture <- function() { + set.seed(42) + list( + learner.results = runif(40, 0.5, 1), + threshold.probability = 0.80, + learner.threshold = 0.9, + feature.frequency = c(a = 5L, b = 3L, c = 1L) + ) +} + +test_that("diagnostics draws both panels in one two-panel layout", { + fx <- .boost_plot_fixture() + seen <- new.env(parent = emptyenv()) + + local_mocked_bindings( + .nns_boost_plot_learner_distribution = function(learner.results, + threshold.probability, + learner.threshold) { + seen$dist_layout <- graphics::par("mfrow") + seen$dist_threshold <- learner.threshold + invisible(NULL) + }, + .nns_boost_plot_feature_frequency = function(feature.frequency) { + seen$freq_layout <- graphics::par("mfrow") + invisible(NULL) + }, + .package = "NNS" + ) + + pdf_file <- withr::local_tempfile(fileext = ".pdf") + grDevices::pdf(pdf_file) + on.exit(grDevices::dev.off(), add = TRUE) + + .nns_boost_plot_diagnostics( + learner.results = fx$learner.results, + threshold.probability = fx$threshold.probability, + learner.threshold = fx$learner.threshold, + feature.frequency = fx$feature.frequency + ) + + # Both helpers were invoked, under the same undisturbed two-panel layout. + expect_identical(seen$dist_layout, c(2L, 1L)) + expect_identical(seen$freq_layout, c(2L, 1L)) + + # The threshold line receives learner.threshold, not the probability. + expect_identical(seen$dist_threshold, fx$learner.threshold) + + # The original graphics parameters are restored after both panels. + expect_identical(graphics::par("mfrow"), c(1L, 1L)) +}) + +test_that("the learner-distribution panel draws the cutoff line and label", { + src <- paste(deparse(body(.nns_boost_plot_learner_distribution)), + collapse = "\n") + expect_true(grepl("abline", src, fixed = TRUE)) + expect_true(grepl("v = learner.threshold", src, fixed = TRUE)) + expect_true(grepl("LPM.VaR(p = %.2f)", src, fixed = TRUE)) + expect_true(grepl("hist", src, fixed = TRUE)) + + freq_src <- paste(deparse(body(.nns_boost_plot_feature_frequency)), + collapse = "\n") + expect_true(grepl("barplot", freq_src, fixed = TRUE)) + expect_true(grepl("horiz = TRUE", freq_src, fixed = TRUE)) +}) + +test_that("NNS.boost invokes both plot helpers only when feature.importance = TRUE", { + calls <- new.env(parent = emptyenv()) + calls$dist <- 0L + calls$freq <- 0L + + local_mocked_bindings( + .nns_boost_plot_learner_distribution = function(...) { + calls$dist <- calls$dist + 1L + invisible(NULL) + }, + .nns_boost_plot_feature_frequency = function(...) { + calls$freq <- calls$freq + 1L + invisible(NULL) + }, + .package = "NNS" + ) + + pdf_file <- withr::local_tempfile(fileext = ".pdf") + grDevices::pdf(pdf_file) + on.exit(grDevices::dev.off(), add = TRUE) + + run_boost <- function(feature.importance) { + NNS.boost( + iris[1:60, 1:4], + iris[1:60, 5], + IVs.test = iris[61:70, 1:4], + epochs = 3, + learner.trials = 5, + type = "CLASS", + feature.importance = feature.importance, + status = FALSE, + seed = 123L + ) + } + + invisible(run_boost(feature.importance = TRUE)) + expect_identical(calls$dist, 1L) + expect_identical(calls$freq, 1L) + + invisible(run_boost(feature.importance = FALSE)) + expect_identical(calls$dist, 1L) + expect_identical(calls$freq, 1L) +}) + +test_that("NNS.boost completes its diagnostics on a PDF device without errors", { + pdf_file <- withr::local_tempfile(fileext = ".pdf") + grDevices::pdf(pdf_file) + + expect_no_error( + invisible(NNS.boost( + iris[1:60, 1:4], + iris[1:60, 5], + IVs.test = iris[61:70, 1:4], + epochs = 3, + learner.trials = 5, + type = "CLASS", + feature.importance = TRUE, + status = FALSE, + seed = 123L + )) + ) + + grDevices::dev.off() + expect_true(file.exists(pdf_file)) + expect_gt(file.size(pdf_file), 0) +}) diff --git a/tests/testthat/test-boost-threshold.R b/tests/testthat/test-boost-threshold.R new file mode 100644 index 000000000..1b6eaf16e --- /dev/null +++ b/tests/testthat/test-boost-threshold.R @@ -0,0 +1,50 @@ +test_that("learner threshold maps probabilities through LPM.VaR", { + scores <- c(0.50, 0.60, 0.70, 0.80, 0.90, 1.00) + + default_max <- .nns_boost_threshold(NULL, "max", FALSE, scores) + expect_identical(default_max$probability, 0.80) + expect_equal(default_max$cutoff, as.numeric(LPM.VaR(0.80, 1, scores))) + + default_min <- .nns_boost_threshold(NULL, "min", FALSE, scores) + expect_identical(default_min$probability, 0.20) + expect_equal(default_min$cutoff, as.numeric(LPM.VaR(0.20, 1, scores))) + + extreme_max <- .nns_boost_threshold(NULL, "max", TRUE, scores) + expect_identical(extreme_max$probability, 1) + expect_equal(extreme_max$cutoff, as.numeric(LPM.VaR(1, 1, scores))) + + extreme_min <- .nns_boost_threshold(NULL, "min", TRUE, scores) + expect_identical(extreme_min$probability, 0) + expect_equal(extreme_min$cutoff, as.numeric(LPM.VaR(0, 1, scores))) +}) + +test_that("a supplied threshold is used as the LPM.VaR probability", { + scores <- c(0.50, 0.60, 0.70, 0.80, 0.90, 1.00) + + supplied <- .nns_boost_threshold(0.65, "max", FALSE, scores) + expect_identical(supplied$probability, 0.65) + expect_equal(supplied$cutoff, as.numeric(LPM.VaR(0.65, 1, scores))) + + # The probability is not treated as a literal objective-score cutoff. + expect_false(isTRUE(all.equal(supplied$cutoff, 0.65))) + + # extreme overrides a supplied probability. + overridden <- .nns_boost_threshold(0.65, "max", TRUE, scores) + expect_identical(overridden$probability, 1) +}) + +test_that("threshold selection does not use fivenum or empirical quantiles", { + boost_source <- paste(deparse(body(NNS.boost)), collapse = "\n") + helper_source <- paste(deparse(body(.nns_boost_threshold)), collapse = "\n") + + for (src in c(boost_source, helper_source)) { + expect_false(grepl("fivenum", src, fixed = TRUE)) + expect_false(grepl("quantile(", src, fixed = TRUE)) + expect_false(grepl("stats::quantile", src, fixed = TRUE)) + } + + # The cutoff comes from LPM.VaR, and epochs reuse the same cutoff variable. + expect_true(grepl("LPM.VaR", helper_source, fixed = TRUE)) + expect_true(grepl("learner.threshold", boost_source, fixed = TRUE)) + expect_false(grepl("Learner Accuracy Threshold", boost_source, fixed = TRUE)) +}) diff --git a/tests/testthat/test-native-nns-distance-default.R b/tests/testthat/test-native-nns-distance-default.R index 282f97531..2f3743dcf 100644 --- a/tests/testthat/test-native-nns-distance-default.R +++ b/tests/testthat/test-native-nns-distance-default.R @@ -66,7 +66,7 @@ test_that("NNS.boost NULL distance matches explicit NNS and propagates L2", { b = c(0, 0.2, 1, 0.9, 0.7, 0.1, 0.6, 0.4)) y <- c(0, 2, 4, 8, 5, 3, 6, 7) args <- list(IVs.train = x, DV.train = y, IVs.test = x[1:2, ], learner.trials = 2, - epochs = 1, status = FALSE, seed = 42) + epochs = 0, status = FALSE, seed = 42) a <- do.call(NNS.boost, c(args, list(dist = NULL))) b <- do.call(NNS.boost, c(args, list(dist = "NNS"))) expect_equal(a$results, b$results, tolerance = 1e-12) @@ -128,11 +128,11 @@ test_that("NNS.boost preserves old positional objective argument order", { test <- x[1:2, ] obj <- expression(sum(abs(predicted - actual))) - positional <- NNS.boost(x, y, test, NULL, NULL, 2, 1, 0.25, FALSE, + positional <- NNS.boost(x, y, test, NULL, NULL, 2, 0, 0.25, FALSE, NULL, NULL, obj, "min", FALSE, FALSE, FALSE, NULL, FALSE, 42) named <- NNS.boost(IVs.train = x, DV.train = y, IVs.test = test, - learner.trials = 2, epochs = 1, CV.size = 0.25, + learner.trials = 2, epochs = 0, CV.size = 0.25, balance = FALSE, ts.test = NULL, threshold = NULL, obj.fn = obj, objective = "min", extreme = FALSE, features.only = FALSE, feature.importance = FALSE,