diff --git a/R/ranger.R b/R/ranger.R index bebe3ea21..69f801323 100644 --- a/R/ranger.R +++ b/R/ranger.R @@ -580,6 +580,8 @@ ranger <- function(formula = NULL, data = NULL, num.trees = 500, mtry = NULL, } else { importance.mode <- 3 } + } else if (importance == 'sobolMDA'){ + importance.mode <- 7 } else { stop("Error: Unknown importance mode.") } diff --git a/man/ranger.Rd b/man/ranger.Rd index 410cd7fe5..b6f9ed9a1 100644 --- a/man/ranger.Rd +++ b/man/ranger.Rd @@ -55,7 +55,7 @@ ranger( \item{mtry}{Number of variables to possibly split at in each node. Default is the (rounded down) square root of the number variables. Alternatively, a single argument function returning an integer, given the number of independent variables.} -\item{importance}{Variable importance mode, one of 'none', 'impurity', 'impurity_corrected', 'permutation'. The 'impurity' measure is the Gini index for classification, the variance of the responses for regression and the sum of test statistics (see \code{splitrule}) for survival.} +\item{importance}{Variable importance mode, one of 'none', 'sobolMDA', 'impurity', 'impurity_corrected', 'permutation'. The sobolMDA provides the proportion of explained output variance lost by removing a given variable. The 'impurity' measure is the Gini index for classification, the variance of the responses for regression and the sum of test statistics (see \code{splitrule}) for survival.} \item{write.forest}{Save \code{ranger.forest} object, required for prediction. Set to \code{FALSE} to reduce memory usage if no prediction intended.} @@ -250,6 +250,7 @@ ranger(trait ~ ., data = dat.gwaa) } \references{ \itemize{ + \item Benard, C., Da Veiga, S., & Scornet, E. (2022). Mean decrease accuracy for random forests: inconsistency, and a practical solution via the Sobol-MDA. Biometrika. \doi{https://doi.org/10.1093/biomet/asac017}. \item Wright, M. N. & Ziegler, A. (2017). ranger: A fast implementation of random forests for high dimensional data in C++ and R. J Stat Softw 77:1-17. \doi{10.18637/jss.v077.i01}. \item Schmid, M., Wright, M. N. & Ziegler, A. (2016). On the use of Harrell's C for clinical risk prediction via random survival forests. Expert Syst Appl 63:450-459. \doi{10.1016/j.eswa.2016.07.018}. \item Wright, M. N., Dankowski, T. & Ziegler, A. (2017). Unbiased split variable selection for random survival forests using maximally selected rank statistics. Stat Med 36:1272-1284. \doi{10.1002/sim.7212}. diff --git a/src/Forest.cpp b/src/Forest.cpp index f914d44ee..d067139b5 100644 --- a/src/Forest.cpp +++ b/src/Forest.cpp @@ -760,6 +760,80 @@ void Forest::computePermutationImportance() { } } +// Compute the Sobol-MDA importance measures. +std::vector Forest::computeSobolMDA(){ + // Compute output variance. + double varY = computeVarianceY(); + // Compute the output unexplained variances of the projected forests. + std::vector SMDA = computeProjectedVariance(); + // Subtract the unexplained variance of the original forest, and normalize with the output variance. + for (size_t j = 0; j < num_independent_variables; j++){ + SMDA[j] = (SMDA[j] - overall_prediction_error)/varY; + } + return(SMDA); +} + +// Compute output variance. +double Forest::computeVarianceY(){ + std::vector y; + for (size_t i = 0; i < num_samples; ++i) { + y.push_back(data->get_y(i, 0)); + } + double meanY = std::accumulate(y.begin(), y.end(), 0.0); + meanY = meanY/num_samples; + double varY = 0; + for (auto & z : y){ + varY += (z - meanY)*(z - meanY); + } + varY = varY/(num_samples-1); + return(varY); +} + +// Compute the output unexplained variances of the projected forests. +std::vector Forest::computeProjectedVariance(){ + // Compute projected forest predictions. + std::vector> forest_predictions = predictProjectedForest(); + // Compute the output unexplained variances associated with these projected predictions. + std::vector projected_variances(num_independent_variables, 0); + for (size_t i = 0; i < num_samples; ++i) { + double y = data->get_y(i, 0); + for (size_t j = 0; j < num_independent_variables; ++j){ + projected_variances[j] += (y - forest_predictions[j][i])*(y - forest_predictions[j][i]); + } + } + for (size_t j = 0; j < num_independent_variables; ++j){ + projected_variances[j] = projected_variances[j]/num_samples; + } + return(projected_variances); +} + +// Compute projected forest predictions. +std::vector> Forest::predictProjectedForest(){ + std::vector> forest_predictions_oob(num_independent_variables, std::vector (num_samples, 0)); + std::vector num_oob(num_samples, 0); + for (size_t i = 0; i < num_trees; ++i) { + // Compute projected tree predictions for all input variables. + std::vector> tree_predictions = trees[i] -> predictProjectedTree(); + for (size_t j = 0; j < num_independent_variables; ++j){ + for (size_t i = 0; i < num_samples; ++i) { + forest_predictions_oob[j][i] += tree_predictions[j][i]; // Sum predictions of all trees. + } + } + // Compute the number of trees used for each oob observation prediction. + std::vector oob_sampleIDs = trees[i] -> getOobSampleIDs(); + for (auto & i : oob_sampleIDs){ + num_oob[i] += 1; + } + } + // Recover average tree predictions by normalizing with the number of trees involved in each case. + for (size_t j = 0; j < num_independent_variables; ++j){ + for (size_t i = 0; i < num_samples; ++i) { + forest_predictions_oob[j][i] = forest_predictions_oob[j][i]/num_oob[i]; + } + } + return(forest_predictions_oob); +} + #ifndef OLD_WIN_R_BUILD void Forest::growTreesInThread(uint thread_idx, std::vector* variable_importance) { if (thread_ranges.size() > thread_idx + 1) { diff --git a/src/Forest.h b/src/Forest.h index a1275b43e..c72f6a25c 100644 --- a/src/Forest.h +++ b/src/Forest.h @@ -107,6 +107,7 @@ class Forest { const std::vector& getVariableImportanceCasewise() const { return variable_importance_casewise; } + std::vector computeSobolMDA(); double getOverallPredictionError() const { return overall_prediction_error; } @@ -155,6 +156,11 @@ class Forest { virtual void computePredictionErrorInternal() = 0; void computePermutationImportance(); + + // SobolMDA methods + double computeVarianceY(); + std::vector computeProjectedVariance(); + std::vector> predictProjectedForest(); // Multithreading methods for growing/prediction/importance, called by each thread void growTreesInThread(uint thread_idx, std::vector* variable_importance); diff --git a/src/Tree.cpp b/src/Tree.cpp index 4305ab15f..94941a9d3 100644 --- a/src/Tree.cpp +++ b/src/Tree.cpp @@ -138,9 +138,12 @@ void Tree::grow(std::vector* variable_importance) { } // Delete sampleID vector to save memory - sampleIDs.clear(); - sampleIDs.shrink_to_fit(); + if (importance_mode != IMP_SOBOL_MDA){ + sampleIDs.clear(); + sampleIDs.shrink_to_fit(); + } cleanUpInternal(); + } void Tree::predict(const Data* prediction_data, bool oob_prediction) { @@ -252,6 +255,211 @@ void Tree::computePermutationImportance(std::vector& forest_importance, } } +// Compute predictions of projected trees for all input variables. +std::vector> Tree::predictProjectedTree() { + + // Initialize variables. + size_t p = data->getNumCols(); // Set p as the number of input variables. + // The vectors nodeIDs_samplesIDs_inb and nodeIDs_samplesIDs_oob are of size p, and respectively store the in-bag and out-of-bag sampleIDs of the projected partitions. + // The j-th component of these two vectors stores the projected partitions of the tree at all levels when the splits on variable j are ignored (i.e., observations are sent to both children nodes). + // Each projected partition is a map: the key is a list of nodeIDs, and the value is the list of the sampleIDs falling simultaneously in all nodes of the key. + std::vector, std::vector>> nodeIDs_samplesIDs_inb(p); + std::vector, std::vector>> nodeIDs_samplesIDs_oob(p); + std::map> leaves_variables; // Lists of variables used in the sequence of splits to reach each terminal leave. + std::map> leaves_sampleIDs; // Lists of sampleIDs falling in each terminal leave. + + // Deduplication of in-bag sample to speed up computations. + std::vector sampleIDs_unique; // Deduplicated sampleIDs. + std::map sampleIDs_count; // Number of repetitions of each observation of sampleIDs_unique in sampleIDs. + for (auto & i : sampleIDs){ + std::pair::iterator, bool> itIDs = sampleIDs_count.insert(std::pair (i, 1)); + if (!itIDs.second){ + itIDs.first->second += 1; + }else{ + sampleIDs_unique.push_back(i); + } + } + + // Drop observations in the projected partitions. First consider in-bag sample (s = 0), and then the out-of-bag sample (s = 1). + for (size_t s = 0; s <= 1; ++s){ + std::vector sampleIDs_temp; + if (s == 0){ + sampleIDs_temp = sampleIDs_unique; // Deduplicated in-bag sample. + } else { + sampleIDs_temp = oob_sampleIDs; // Out-of-bag sample. + } + // Loop through all observations of the considered sample. + for (auto & i : sampleIDs_temp) { + std::vector split_variables; // Track the list of variables involved in the splits of the previous tree levels. + bool is_terminal = false; // Track whether terminal leave is reached. + size_t nodeIDorigin = 0; // nodeIDorigin is the nodeID where observation i falls at each level of the original tree. Initialized here at the root node. + // Loop through the tree levels, and stop when observation i has reached the terminal leave of the original tree. + while (!is_terminal){ + if (child_nodeIDs[0][nodeIDorigin] == 0 && child_nodeIDs[1][nodeIDorigin] == 0) { + is_terminal = true; // If nodeIDorigin is a terminal leave, break the loop through the levels of the original tree. + }else{ + std::vector nodeIDtemp; // Store the list of nodeIDs at each level when observation i is dropped down the projected tree. + nodeIDtemp.push_back(nodeIDorigin); + size_t varIDorigin = split_varIDs[nodeIDorigin]; // varIDorigin is the splitting variable at nodeIDorigin. + nodeIDorigin = getChildID(i, nodeIDorigin, varIDorigin); // Update nodeIDorigin as the child node where observation i falls at the next level of the original tree. + // If varIDorigin is NOT already involved in the nodes hit by observation i at the previous tree levels, + // then compute the projected cell where observation i falls, by ignoring all splits using varIDorigin down the tree, otherwise send observation i to the next level of the original tree. + if (find(split_variables.begin(), split_variables.end(), varIDorigin) == split_variables.end()){ + split_variables.push_back(varIDorigin); // Add varIDorigin to the list of variables already met by observation i at the previous tree levels. + std::vector nodeIDnext; // List of inner nodeIDs of the next tree level where observation i falls. + std::vector nodeIDpred; // List of terminal nodeIDs where observation i falls. + std::vector nodeIDpredtemp = nodeIDtemp; // Store projected cell of the previous level of the projected tree. + bool next_level = true; + // Loop through the remaining tree levels, ignoring all splits using varIDorigin, i.e., + // observation i is sent to both children nodes when a split on varIDorigin is hit, and then, observation i ends up in multiple terminal leaves. + while (next_level) { + for (auto & nodeID : nodeIDtemp){ + if (child_nodeIDs[0][nodeID] == 0 && child_nodeIDs[1][nodeID] == 0) { + nodeIDpred.push_back(nodeID); + }else{ + size_t split_varID = split_varIDs[nodeID]; + if (split_varID != varIDorigin){ + // Move to child node if splitting variable is NOT varIDorigin. + nodeIDnext.push_back(getChildID(i, nodeID, split_varID)); + }else{ + // Move to both left and right children nodes if splitting variable is varIDorigin. + nodeIDnext.push_back(child_nodeIDs[0][nodeID]); + nodeIDnext.push_back(child_nodeIDs[1][nodeID]); + } + } + } + // Store observation i in the projected partition of the current tree level. + if (s == 0){ + // In-bag observation case. + if (nodeIDnext.size() > 0 || nodeIDtemp == std::vector {0}){ + std::vector nodeIDlevel = nodeIDpred; // nodeIDlevel is the list of nodeIDs defining the projected cell at the current tree level. + nodeIDlevel.insert(nodeIDlevel.end(), nodeIDnext.begin(), nodeIDnext.end()); // Assemble nodeIDs of terminal leaves and inner nodes in nodeIDlevel. + std::pair, std::vector>::iterator, bool> it_bool; + // Insert the list of nodeIDs defining the projected cell at the current tree level (nodeIDlevel) and the associated sampleIDs. + it_bool = nodeIDs_samplesIDs_inb[varIDorigin].insert(std::pair, std::vector> (nodeIDlevel, {i})); + if (!it_bool.second){ + it_bool.first->second.push_back(i); + } + } + nodeIDtemp = nodeIDnext; + nodeIDnext = {}; + next_level = nodeIDtemp.size() > 0; + }else{ + // Out-of-bag observation case. + std::vector nodeIDlevel = nodeIDpred; + nodeIDlevel.insert(nodeIDlevel.end(), nodeIDnext.begin(), nodeIDnext.end()); // Assemble nodeIDs of terminal leaves and inner nodes. + // Check whether the projected cell exists in the pre-computed in-bag projected partition. + std::map, std::vector>::iterator it = nodeIDs_samplesIDs_inb[varIDorigin].find(nodeIDlevel); + if (it == nodeIDs_samplesIDs_inb[varIDorigin].end() || nodeIDnext.size() == 0){ + // If the projected cell is NOT in the in-bag projected partition or the projected cell is empty, + // set the projected cell of the previous tree level (nodeIDpredtemp) as the terminal projected leave. + next_level = false; + std::pair, std::vector>::iterator, bool> it_bool; + it_bool = nodeIDs_samplesIDs_oob[varIDorigin].insert(std::pair, std::vector> (nodeIDpredtemp, {i})); + if (!it_bool.second){ + it_bool.first->second.push_back(i); + } + }else{ + // If the projected cell is in the in-bag projected partition and not empty, keep splitting. + nodeIDpredtemp = nodeIDlevel; + nodeIDtemp = nodeIDnext; + nodeIDnext = {}; + } + } + } + } + } + } + // For out-of-bag observations, store terminal leave information. + if (s == 1){ + // Store list of variables used in the sequence of splits. + leaves_variables.insert(std::pair> (nodeIDorigin, split_variables)); + std::pair>::iterator, bool> it_leaves; + // Store oob sampleIDs in terminal leaves. + it_leaves = leaves_sampleIDs.insert(std::pair> (nodeIDorigin, {i})); + if (!it_leaves.second){ + it_leaves.first->second.push_back(i); + } + } + } + } + + // Generate out-of-bag predictions from the obtained projected partitions. + int nsample = data->getNumRows(); // Sample size. + std::vector> predictions_oob; // Out-of-bag projected tree predictions for all input variables. + std::map, double> nodeIDs_response; + std::pair, double>::iterator, bool> it_response; + for (size_t j = 0; j < p; ++j){ + std::vector predictions_varID(nsample, 0); // Out-of-bag projected tree predictions when variable varID is ignored. + for (std::map, std::vector>::iterator it = nodeIDs_samplesIDs_oob[j].begin(); it != nodeIDs_samplesIDs_oob[j].end(); ++it){ + double response; + std::vector nodeIDs = it->first; + it_response = nodeIDs_response.insert(std::pair, double> (nodeIDs, 0)); + // Compute projected cell output at its first occurence, using in-bag data. + if (it_response.second){ + std::vector sampleIDs_inb = nodeIDs_samplesIDs_inb[j][nodeIDs]; + response = 0; + size_t inb_size = 0; + for (auto & i : sampleIDs_inb){ + size_t size_temp = sampleIDs_count[i]; + response += data->get_y(i, 0)*size_temp; + inb_size += size_temp; + } + response = response/inb_size; + it_response.first->second = response; + }else{ + response = it_response.first->second; + } + // Assign projected cell output to all oob observations falling in this projected cell. + for (auto & i : it->second){ + predictions_varID[i] = response; + } + } + predictions_oob.push_back(predictions_varID); + } + // For all variables not involved in the splits of the original tree path of a given observation, the projected prediction is given by the original tree. + for (std::map>::iterator it = leaves_variables.begin(); it != leaves_variables.end(); ++it){ + double response = split_values[it->first]; + for (size_t j = 0; j < p; ++j){ + if (find(it->second.begin(), it->second.end(), j) == it->second.end()){ + for (auto & i : leaves_sampleIDs[it->first]){ + predictions_oob[j][i] = response; + } + } + } + } + + return(predictions_oob); + +} + +// Return child nodeID. +size_t Tree::getChildID(size_t i, size_t nodeID, size_t split_varID){ + size_t childID; + double value = data->get_x(i, split_varID); + if (data->isOrderedVariable(split_varID)) { + if (value <= split_values[nodeID]) { + // Move to left child. + childID = child_nodeIDs[0][nodeID]; + } else { + // Move to right child. + childID = child_nodeIDs[1][nodeID]; + } + } else { + size_t factorID = floor(value) - 1; + size_t splitID = floor(split_values[nodeID]); + // Left if 0 found at position factorID. + if (!(splitID & (1ULL << factorID))) { + // Move to left child. + childID = child_nodeIDs[0][nodeID]; + } else { + // Move to right child. + childID = child_nodeIDs[1][nodeID]; + } + } + return(childID); +} + // #nocov start void Tree::appendToFile(std::ofstream& file) { diff --git a/src/Tree.h b/src/Tree.h index cb409b39c..f97faa440 100644 --- a/src/Tree.h +++ b/src/Tree.h @@ -16,6 +16,7 @@ #include #include #include +#include #include "globals.h" #include "Data.h" @@ -51,6 +52,10 @@ class Tree { void computePermutationImportance(std::vector& forest_importance, std::vector& forest_variance, std::vector& forest_importance_casewise); + + std::vector> predictProjectedTree(); + + size_t getChildID(size_t i, size_t nodeID, size_t split_varID); void appendToFile(std::ofstream& file); virtual void appendToFileInternal(std::ofstream& file) = 0; diff --git a/src/globals.h b/src/globals.h index b794273be..38cbb00a8 100644 --- a/src/globals.h +++ b/src/globals.h @@ -57,7 +57,8 @@ enum ImportanceMode { IMP_PERM_LIAW = 4, IMP_PERM_RAW = 3, IMP_GINI_CORRECTED = 5, - IMP_PERM_CASEWISE = 6 + IMP_PERM_CASEWISE = 6, + IMP_SOBOL_MDA = 7 }; const uint MAX_IMP_MODE = 6; diff --git a/src/rangerCpp.cpp b/src/rangerCpp.cpp index 991fc9609..39b0e4afc 100644 --- a/src/rangerCpp.cpp +++ b/src/rangerCpp.cpp @@ -227,9 +227,13 @@ Rcpp::List rangerCpp(uint treetype, Rcpp::NumericMatrix& input_x, Rcpp::NumericM result.push_back(forest->getMtry(), "mtry"); result.push_back(forest->getMinNodeSize(), "min.node.size"); if (importance_mode != IMP_NONE) { - result.push_back(forest->getVariableImportance(), "variable.importance"); - if (importance_mode == IMP_PERM_CASEWISE) { - result.push_back(forest->getVariableImportanceCasewise(), "variable.importance.local"); + if (importance_mode == IMP_SOBOL_MDA && treetype == TREE_REGRESSION){ + result.push_back(forest->computeSobolMDA(), "variable.importance"); + }else{ + result.push_back(forest->getVariableImportance(), "variable.importance"); + if (importance_mode == IMP_PERM_CASEWISE) { + result.push_back(forest->getVariableImportanceCasewise(), "variable.importance.local"); + } } } result.push_back(forest->getOverallPredictionError(), "prediction.error");