further clean up to the point explain_forest works for ranger

Author: Yue-Jiang

R/explain_forest.R

@@ -22,15 +22,29 @@
 #'
 #' @export
 explain_forest <- function(forest, interactions = FALSE, data = NULL, vars = NULL, no_of_pred_plots = 3, pred_grid = 100,
-                             measures = if(forest$type == "classification")
-                               c("mean_min_depth", "accuracy_decrease", "gini_decrease", "no_of_nodes", "times_a_root") else
-                                 c("mean_min_depth", "mse_increase", "node_purity_increase", "no_of_nodes", "times_a_root")){
-  if(any(c("accuracy_decrease", "mse_increase") %in% measures) & dim(forest$importance)[2] == 1) {
+                           measures = NULL){
+  if(is.null(measures)){
+    if("randomForest" %in% class(forest)){
+      if(forest$type == "classification"){
+        measures <- c("mean_min_depth", "accuracy_decrease", "gini_decrease", "no_of_nodes", "times_a_root")
+      } else{
+        measures <- c("mean_min_depth", "mse_increase", "node_purity_increase", "no_of_nodes", "times_a_root")
+      }
+    } else if("ranger" %in% class(forest)){
+      measures <- c("mean_min_depth", forest$importance.mode, "no_of_nodes", "times_a_root")
+    }
+  }
+  if("randomForest" %in% class(forest) && dim(forest$importance)[2] == 1){
     stop(paste("Your forest does not contain information on local importance so",
-               paste(intersect(c("accuracy_decrease", "mse_increase"), measures), sep=", "),
+               ifelse(forest$type == "classification", "accuracy_decrease", "mse_increase"),
                "measure cannot be extracted.",
                "To add it regrow the forest with the option localImp = TRUE and run this function again."))
   }
+  if("ranger" %in% class(forest) && forest$importance.mode == "none"){
+   stop(paste("Your forest does not contain importance information so",
+              "importance cannot be extracted.",
+              "To add it regrow the forest with the option importance other than 'none' and run this function again."))
+  }
   environment <- new.env()
   environment$forest <- forest
   environment$data <- data
@@ -42,7 +56,7 @@ explain_forest <- function(forest, interactions = FALSE, data = NULL, vars = NUL
   directory <- getwd()
   path_to_templates <- file.path(path.package("randomForestExplainer"), "templates")
   template_name <- grep('explain_forest_template.rmd', list.files(path_to_templates),
-                   ignore.case = TRUE, value = TRUE)
+                        ignore.case = TRUE, value = TRUE)
 
   rmarkdown::render(file.path(path_to_templates, template_name),
                     "html_document", output_file = paste0(directory, "/Your_forest_explained.html"),

R/measure_importance.R

@@ -60,7 +60,7 @@ measure_vimp <- function(forest, only_nonlocal = FALSE){
 # Extract randomForest variable importance measures
 # ranger
 measure_vimp_ranger <- function(forest){
-  if (forest$importance.mode == "none") {
+  if (forest$importance.mode == "none"){
     stop("No variable importance available, regenerate forest by ranger(..., importance='impurity').")
   }
   frame <- data.frame(importance=forest$variable.importance,
@@ -95,8 +95,8 @@ measure_times_a_root <- function(min_depth_frame){
 measure_p_value <- function(importance_frame){
   total_no_of_nodes <- sum(importance_frame$no_of_nodes)
   p_value <- unlist(lapply(importance_frame$no_of_nodes,
-                  function(x) stats::binom.test(x, total_no_of_nodes, 1/nrow(importance_frame),
-                                         alternative = "greater")$p.value))
+                           function(x) stats::binom.test(x, total_no_of_nodes, 1/nrow(importance_frame),
+                                                         alternative = "greater")$p.value))
   return(p_value)
 }
 
@@ -132,7 +132,7 @@ measure_importance.randomForest <- function(forest, mean_sample = "top_trees", m
     if(forest$type == "classification"){
       measures <- c("mean_min_depth", "no_of_nodes", "accuracy_decrease",
                     "gini_decrease", "no_of_trees", "times_a_root", "p_value")
-    } else if(forest$type =="regression") {
+    } else if(forest$type =="regression"){
       measures <- c("mean_min_depth", "no_of_nodes", "mse_increase", "node_purity_increase",
                     "no_of_trees", "times_a_root", "p_value")
     }
@@ -161,7 +161,7 @@ measure_importance.randomForest <- function(forest, mean_sample = "top_trees", m
   }
   if(forest$type == "classification"){
     vimp <- c("accuracy_decrease", "gini_decrease")
-  } else if(forest$type =="regression") {
+  } else if(forest$type =="regression"){
     vimp <- c("mse_increase", "node_purity_increase")
   }
   if(all(vimp %in% measures)){
@@ -193,7 +193,7 @@ measure_importance.randomForest <- function(forest, mean_sample = "top_trees", m
 measure_importance.ranger <- function(forest, mean_sample = "top_trees", measures = NULL){
   tree <- NULL; splitvarName <- NULL; depth <- NULL
   if(is.null(measures)){
-      measures <- c("mean_min_depth", "no_of_nodes", forest$importance.mode, "no_of_trees", "times_a_root", "p_value")
+    measures <- c("mean_min_depth", "no_of_nodes", forest$importance.mode, "no_of_trees", "times_a_root", "p_value")
   }
   if(("p_value" %in% measures) && !("no_of_nodes" %in% measures)){
     measures <- c(measures, "no_of_nodes")
@@ -261,7 +261,7 @@ important_variables <- function(importance_frame, k = 15,
     if("predicted" %in% measures){
       measures <- names(importance_frame)[2:5]
     }
-  } else if ("ranger" %in% class(importance_frame)) {
+  } else if ("ranger" %in% class(importance_frame)){
     importance_frame <- measure_importance(importance_frame)
   }
   rankings <- data.frame(variable = importance_frame$variable, mean_min_depth =
@@ -376,14 +376,16 @@ plot_multi_way_importance <- function(importance_frame, x_measure = "mean_min_de
 #' plot_importance_ggpairs(frame, measures = c("mean_min_depth", "times_a_root"))
 #'
 #' @export
-plot_importance_ggpairs <- function(importance_frame, measures =
-                                      names(importance_frame)[c(2, 4, 5, 3, 7)],
+plot_importance_ggpairs <- function(importance_frame, measures = NULL,
                                     main = "Relations between measures of importance"){
-  if("randomForest" %in% class(importance_frame)){
+  if(any(c("randomForest", "ranger") %in% class(importance_frame))){
     importance_frame <- measure_importance(importance_frame)
-    if("predicted" %in% measures){
-      names(importance_frame)[c(2, 4, 5, 3, 7)]
-    }
+  }
+  if (is.null(measures)){
+    default_measures <- c("gini_decrease", "node_purity_increase", # randomForest
+                          "impurity", "impurity_corrected", "permutation", # ranger
+                          "mean_min_depth", "no_of_trees", "no_of_nodes", "p_value")
+    measures <- intersect(default_measures, colnames(importance_frame))
   }
   plot <- ggpairs(importance_frame[, measures]) + theme_bw()
   if(!is.null(main)){
@@ -397,7 +399,7 @@ plot_importance_ggpairs <- function(importance_frame, measures =
 #' Plot against each other rankings of variables according to various measures of importance
 #'
 #' @param importance_frame A result of using the function measure_importance() to a random forest or a randomForest object
-#' @param measures A character vector specifying the measures of importance to be used
+#' @param measures A character vector specifying the measures of importance to be used.
 #' @param main A string to be used as title of the plot
 #'
 #' @return A ggplot object
@@ -411,22 +413,29 @@ plot_importance_ggpairs <- function(importance_frame, measures =
 #' plot_importance_ggpairs(frame, measures = c("mean_min_depth", "times_a_root"))
 #'
 #' @export
-plot_importance_rankings <- function(importance_frame, measures =
-                                       names(importance_frame)[c(2, 4, 5, 3, 7)],
+plot_importance_rankings <- function(importance_frame, measures = NULL,
                                      main = "Relations between rankings according to different measures"){
-  if("randomForest" %in% class(importance_frame)){
+  if(any(c("randomForest", "ranger") %in% class(importance_frame))){
     importance_frame <- measure_importance(importance_frame)
-    if("predicted" %in% measures){
-      names(importance_frame)[c(2, 4, 5, 3, 7)]
-    }
   }
-  rankings <- data.frame(variable = importance_frame$variable, mean_min_depth =
-                           frankv(importance_frame$mean_min_depth, ties.method = "dense"),
-                         apply(importance_frame[, -c(1, 2)], 2,
+  if (is.null(measures)){
+    default_measures <- c("gini_decrease", "node_purity_increase", # randomForest
+                          "impurity", "impurity_corrected", "permutation", # ranger
+                          "mean_min_depth", "no_of_trees", "no_of_nodes", "p_value")
+    measures <- intersect(default_measures, colnames(importance_frame))
+  }
+  rankings <- data.frame(variable = importance_frame$variable,
+                         apply(importance_frame[, !colnames(importance_frame) %in% c("variable", "mean_min_depth", "p_value")], 2,
                                function(x) frankv(x, order = -1, ties.method = "dense")))
+  if ("mean_min_depth" %in% measures){
+    rankings$mean_min_depth = frankv(importance_frame$mean_min_depth, ties.method = "dense")
+  }
+  if ("p_value" %in% measures){
+    rankings$p_value = frankv(importance_frame$p_value, ties.method = "dense")
+  }
   plot <- ggpairs(rankings[, measures], lower = list(continuous = function(data, mapping){
     ggplot(data = data, mapping = mapping) + geom_point() +  geom_smooth(method = "loess")
-    }))+ theme_bw()
+  })) + theme_bw()
   if(!is.null(main)){
     plot <- plot + ggtitle(main)
   }

R/min_depth_interactions.R

@@ -39,11 +39,11 @@ conditional_depth_ranger <- function(frame, vars){
         df <- frame[begin:nrow(frame), setdiff(names(frame), setdiff(vars, j))]
         df[[j]][1] <- 0
         for(k in 2:nrow(df)){
-          if(length(df[df[, "leftChild"] == as.numeric(df[k, "number"]) |
-                       df[, "rightChild"] == as.numeric(df[k, "number"]), j]) != 0){
+          if(length(df[(!is.na(df[, "leftChild"]) & df[, "leftChild"] == as.numeric(df[k, "number"])) |
+                       (!is.na(df[, "rightChild"]) & df[, "rightChild"] == as.numeric(df[k, "number"])), j]) != 0){
             df[k, j] <-
-              df[df[, "leftChild"] == as.numeric(df[k, "number"]) |
-                   df[, "rightChild"] == as.numeric(df[k, "number"]), j] + 1
+              df[(!is.na(df[, "leftChild"]) & df[, "leftChild"] == as.numeric(df[k, "number"])) |
+                   (!is.na(df[, "rightChild"]) & df[, "rightChild"] == as.numeric(df[k, "number"])), j] + 1
           }
         }
         frame[begin:nrow(frame), setdiff(names(frame), setdiff(vars, j))] <- df
@@ -84,8 +84,8 @@ min_depth_interactions_values <- function(forest, vars){
 min_depth_interactions_values_ranger <- function(forest, vars){
   `.` <- NULL; .SD <- NULL; tree <- NULL; splitvarName <- NULL
   interactions_frame <-
-    lapply(1:forest$num.trees, function(i) ranger::treeInfo(forest, k = i, labelVar = T) %>%
-             calculate_tree_depth() %>% cbind(., tree = i, number = 1:nrow(.))) %>%
+    lapply(1:forest$num.trees, function(i) ranger::treeInfo(forest, tree = i) %>%
+             calculate_tree_depth_ranger() %>% cbind(., tree = i, number = 1:nrow(.))) %>%
     data.table::rbindlist() %>% as.data.frame()
   interactions_frame[vars] <- as.numeric(NA)
   interactions_frame <-
@@ -95,10 +95,10 @@ min_depth_interactions_values_ranger <- function(forest, vars){
   mean_tree_depth[mean_tree_depth == -Inf] <- NA
   mean_tree_depth <- colMeans(mean_tree_depth[, vars], na.rm = TRUE)
   min_depth_interactions_frame <-
-    interactions_frame %>% dplyr::group_by(tree, `split var`) %>%
+    interactions_frame %>% dplyr::group_by(tree, splitvarName) %>%
     dplyr::summarize_at(vars, funs(min(., na.rm = TRUE))) %>% as.data.frame()
   min_depth_interactions_frame[min_depth_interactions_frame == Inf] <- NA
-  min_depth_interactions_frame <- min_depth_interactions_frame[!is.na(min_depth_interactions_frame$`split var`), ]
+  min_depth_interactions_frame <- min_depth_interactions_frame[!is.na(min_depth_interactions_frame$splitvarName), ]
   colnames(min_depth_interactions_frame)[2] <- "variable"
   min_depth_interactions_frame[, -c(1:2)] <- min_depth_interactions_frame[, -c(1:2)] - 1
   return(list(min_depth_interactions_frame, mean_tree_depth))
@@ -120,7 +120,8 @@ min_depth_interactions_values_ranger <- function(forest, vars){
 #' min_depth_interactions(forest, c("Petal.Width", "Petal.Length"))
 #'
 #' @export
-min_depth_interactions <- function(){
+min_depth_interactions <- function(forest, vars = important_variables(measure_importance(forest)),
+                                   mean_sample = "top_trees", uncond_mean_sample = mean_sample){
   UseMethod("min_depth_interactions")
 }
 
@@ -217,7 +218,7 @@ min_depth_interactions.ranger <- function(forest, vars = important_variables(mea
     dplyr::summarize(min(depth))
   colnames(min_depth_frame) <- c("tree", "variable", "minimal_depth")
   min_depth_frame <- as.data.frame(min_depth_frame[!is.na(min_depth_frame$variable),])
-  importance_frame <- get_min_depth_means_ranger(min_depth_frame, min_depth_count(min_depth_frame), uncond_mean_sample)
+  importance_frame <- get_min_depth_means(min_depth_frame, min_depth_count(min_depth_frame), uncond_mean_sample)
   colnames(importance_frame)[2] <- "uncond_mean_min_depth"
   interactions_frame <- merge(interactions_frame, importance_frame)
 }
@@ -273,6 +274,9 @@ plot_min_depth_interactions <- function(interactions_frame, k = 30,
 #' @param variable2 A character string with the name a numerical predictor that will on Y-axis
 #' @param grid The number of points on the one-dimensional grid on x and y-axis
 #' @param main A string to be used as title of the plot
+#' @param time A numeric value specifying the time at which to predict survival probability, only
+#'  applies to survival forests. If not specified, the time closest to predicted median survival
+#'  time is used
 #'
 #' @return A ggplot2 object
 #'
@@ -285,7 +289,8 @@ plot_min_depth_interactions <- function(interactions_frame, k = 30,
 #' @export
 plot_predict_interaction <- function(forest, data, variable1, variable2, grid = 100,
                                      main = paste0("Prediction of the forest for different values of ",
-                                                   paste0(variable1, paste0(" and ", variable2)))){
+                                                   paste0(variable1, paste0(" and ", variable2))),
+                                     time = NULL){
   UseMethod("plot_predict_interaction")
 }
 
@@ -296,7 +301,8 @@ plot_predict_interaction <- function(forest, data, variable1, variable2, grid =
 #' @export
 plot_predict_interaction.randomForest <- function(forest, data, variable1, variable2, grid = 100,
                                                   main = paste0("Prediction of the forest for different values of ",
-                                                                paste0(variable1, paste0(" and ", variable2)))){
+                                                                paste0(variable1, paste0(" and ", variable2))),
+                                                  time = NULL){
   newdata <- expand.grid(seq(min(data[[variable1]]), max(data[[variable1]]), length.out = grid),
                          seq(min(data[[variable2]]), max(data[[variable2]]), length.out = grid))
   colnames(newdata) <- c(variable1, variable2)
@@ -343,7 +349,8 @@ plot_predict_interaction.randomForest <- function(forest, data, variable1, varia
 #' @export
 plot_predict_interaction.ranger <- function(forest, data, variable1, variable2, grid = 100,
                                             main = paste0("Prediction of the forest for different values of ",
-                                                          paste0(variable1, paste0(" and ", variable2)))){
+                                                          paste0(variable1, paste0(" and ", variable2))),
+                                            time = NULL){
   newdata <- expand.grid(seq(min(data[[variable1]]), max(data[[variable1]]), length.out = grid),
                          seq(min(data[[variable2]]), max(data[[variable2]]), length.out = grid))
   colnames(newdata) <- c(variable1, variable2)
@@ -358,7 +365,7 @@ plot_predict_interaction.ranger <- function(forest, data, variable1, variable2,
     newdata[[i]] <- data[[i]][sample(1:n, nrow(newdata), replace = TRUE)]
   }
   if(forest$treetype == "Regression"){
-    newdata$prediction <- predict(forest, newdata, type = "response")
+    newdata$prediction <- predict(forest, newdata, type = "response")$predictions
     plot <- ggplot(newdata, aes_string(x = variable1, y = variable2, fill = "prediction")) +
       geom_raster() + theme_bw() +
       scale_fill_gradient2(midpoint = min(newdata$prediction) + 0.5 * (max(newdata$prediction) - min(newdata$prediction)),
@@ -377,8 +384,23 @@ plot_predict_interaction.ranger <- function(forest, data, variable1, variable2,
       geom_raster() + theme_bw() + facet_wrap(~ variable) +
       scale_fill_gradient2(midpoint = min(newdata$prediction) + 0.5 * (max(newdata$prediction) - min(newdata$prediction)),
                            low = "blue", high = "red")
-  } else if(forest$treetype == "Classification") {
+  } else if(forest$treetype == "Classification"){
     stop("Ranger forest for classification needs to be generated by ranger(..., probability = TRUE).")
+  } else if(forest$treetype == "Survival"){
+    pred <- predict(forest, newdata, type = "response")
+    if (is.null(time)) {
+      time <- pred$unique.death.times[which.min(abs(colMeans(pred$survival, na.rm = TRUE) - 0.5))]
+      message(sprintf("Using unique death time %s which is the closest to predicted median survival time.", time))
+    } else if (!time %in% pred$unique.death.times) {
+      new_time <- pred$unique.death.times[which.min(abs(pred$unique.death.times - time))]
+      message(sprintf("Using closest unique death time %s instead of %s.", new_time, time))
+      time <- new_time
+    }
+    newdata$prediction <- pred$survival[, pred$unique.death.times == time, drop = TRUE]
+    plot <- ggplot(newdata, aes_string(x = variable1, y = variable2, fill = "prediction")) +
+      geom_raster() + theme_bw() +
+      scale_fill_gradient2(midpoint = min(newdata$prediction) + 0.5 * (max(newdata$prediction) - min(newdata$prediction)),
+                           low = "blue", high = "red")
   } else {
     stop(sprintf("Ranger forest type '%s' is currently not supported.", forest$treetype))
   }