development version

DESCRIPTION

@@ -1,10 +1,11 @@
 Package: survex
 Title: Explainable Machine Learning in Survival Analysis
-Version: 1.0.0.9000
+Version: 1.0.0.9102
 Authors@R: 
     c(
         person("Mikołaj", "Spytek", email = "[email protected]", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-7111-2286")),
         person("Mateusz", "Krzyziński", role = c("aut"), comment = c(ORCID = "0000-0001-6143-488X")),
+        person("Sophie", "Langbein", role = c("aut")), 
         person("Hubert", "Baniecki", role = c("aut"), comment = c(ORCID = "0000-0001-6661-5364")),
         person("Przemyslaw", "Biecek", role = c("aut"), comment = c(ORCID = "0000-0001-8423-1823"))
     )
@@ -18,11 +19,11 @@ Description: Survival analysis models are commonly used in medicine and other ar
 License: GPL (>= 3)
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3
 Depends: R (>= 3.5.0)
 Imports:
     DALEX (>= 2.2.1),
-    ggplot2,
+    ggplot2 (>= 3.4.0),
     kernelshap,
     pec,
     survival,
@@ -34,13 +35,14 @@ Suggests:
     generics,
     glmnet,
     ingredients,
-    knitr,
+    knitr (>= 1.42),
     mboost,
     parsnip,
     progressr,
     randomForestSRC,
     ranger,
     rmarkdown,
+    rms,
     testthat (>= 3.0.0),
     withr,
     xgboost

NAMESPACE

@@ -12,9 +12,12 @@ S3method(model_performance,default)
 S3method(model_performance,surv_explainer)
 S3method(model_profile,default)
 S3method(model_profile,surv_explainer)
-S3method(plot,feature_importance_explainer)
+S3method(model_profile_2d,surv_explainer)
+S3method(model_survshap,surv_explainer)
+S3method(plot,aggregated_surv_shap)
 S3method(plot,model_parts_survival)
 S3method(plot,model_performance_survival)
+S3method(plot,model_profile_2d_survival)
 S3method(plot,model_profile_survival)
 S3method(plot,predict_parts_survival)
 S3method(plot,predict_profile_survival)
@@ -29,6 +32,7 @@ S3method(predict_parts,default)
 S3method(predict_parts,surv_explainer)
 S3method(predict_profile,default)
 S3method(predict_profile,surv_explainer)
+S3method(print,model_profile_2d_survival)
 S3method(print,model_profile_survival)
 S3method(print,surv_ceteris_paribus)
 S3method(print,surv_feature_importance)
@@ -59,6 +63,8 @@ export(loss_one_minus_integrated_cd_auc)
 export(model_parts)
 export(model_performance)
 export(model_profile)
+export(model_profile_2d)
+export(model_survshap)
 export(predict_parts)
 export(predict_profile)
 export(risk_from_chf)
@@ -78,6 +84,9 @@ importFrom(DALEX,theme_ema)
 importFrom(DALEX,theme_ema_vertical)
 importFrom(stats,aggregate)
 importFrom(stats,as.formula)
+importFrom(stats,ave)
+importFrom(stats,cmdscale)
+importFrom(stats,ecdf)
 importFrom(stats,median)
 importFrom(stats,model.frame)
 importFrom(stats,model.matrix)
@@ -86,8 +95,10 @@ importFrom(stats,optim)
 importFrom(stats,predict)
 importFrom(stats,quantile)
 importFrom(stats,reorder)
+importFrom(stats,reshape)
 importFrom(stats,rnorm)
 importFrom(stats,stepfun)
+importFrom(stats,xtabs)
 importFrom(utils,head)
 importFrom(utils,stack)
 importFrom(utils,tail)

NEWS.md

@@ -1,6 +1,14 @@
 # survex (development)
 
-* Fix not being able to plot or print SurvLIME results for the cph model sometimes. ([#72](https://github.com/ModelOriented/survex/issues/72))
+* fixed not being able to plot or print SurvLIME results for the cph model sometimes. ([#72](https://github.com/ModelOriented/survex/issues/72))
+* added global explanations via the SurvSHAP(t) method (see `model_survshap()` function)
+* added plots for global SurvSHAP(t) explanations (see `plot.aggregated_surv_shap()`)
+* added Accumulated Local Effects (ALE) explanations (see `model_profile(..., type = "accumulated")`)
+* added 2-dimensional PDP and ALE plots (see `model_profile_2d()` function)
+* added `plot2()` function for plotting PDP and ALE explanations without the time dimension
+* made improvements on the vignettes for the package (see `vignette("pdp")` and `vignette("global-survshap")`)
+* increased the test coverage of the pacakge
+* reduced the number of expensive `requireNamespace()` calls ([#83](https://github.com/ModelOriented/survex/issues/83))
 
 # survex 1.0.0
 

R/explain.R

@@ -126,19 +126,15 @@ explain_survival <-
 
         if (is.null(predict_survival_function) &&
             !is.null(predict_cumulative_hazard_function)) {
-            predict_survival_function <- function(model, newdata, times) {
-                cumulative_hazard_to_survival(predict_cumulative_hazard_function(model, newdata, times))
-            }
+            predict_survival_function <- function(model, newdata, times) cumulative_hazard_to_survival(predict_cumulative_hazard_function(model, newdata, times))
             attr(predict_survival_function, "verbose_info") <- "exp(-predict_cumulative_hazard_function) will be used"
             attr(predict_survival_function, "is.default") <- TRUE
         }
 
         if (is.null(predict_cumulative_hazard_function) &&
             !is.null(predict_survival_function)) {
             predict_cumulative_hazard_function <-
-                function(model, newdata, times) {
-                    survival_to_cumulative_hazard(predict_survival_function(model, newdata, times))
-                }
+                function(model, newdata, times) survival_to_cumulative_hazard(predict_survival_function(model, newdata, times))
             attr(predict_cumulative_hazard_function, "verbose_info") <- "-log(predict_survival_function) will be used"
             attr(predict_cumulative_hazard_function, "is.default") <- TRUE
         }
@@ -250,9 +246,7 @@ explain_survival <-
         # verbose predict function
         if (is.null(predict_function)) {
             if (!is.null(predict_cumulative_hazard_function)) {
-                predict_function <- function(model, newdata) {
-                    risk_from_chf(predict_cumulative_hazard_function(model, newdata, times = times))
-                }
+                predict_function <- function(model, newdata) risk_from_chf(predict_cumulative_hazard_function(model, newdata, times = times))
                 verbose_cat("  -> predict function  : ", "sum over the predict_cumulative_hazard_function will be used",  is.default = TRUE, verbose = verbose)
             } else {
                 verbose_cat("  -> predict function   :  not specified! (", color_codes$red_start, "WARNING", color_codes$red_end, ")", verbose = verbose)
@@ -726,9 +720,7 @@ explain.model_fit <-
 
         if (is.null(predict_cumulative_hazard_function)) {
             predict_cumulative_hazard_function <-
-                function(object, newdata, times) {
-                    survival_to_cumulative_hazard(predict_survival_function(object, newdata, times))
-                }
+                function(object, newdata, times) survival_to_cumulative_hazard(predict_survival_function(object, newdata, times))
             attr(predict_cumulative_hazard_function, "verbose_info") <- "-log(predict_survival_function) will be used"
             attr(predict_cumulative_hazard_function, "is.default") <- TRUE
         } else {
@@ -737,14 +729,10 @@ explain.model_fit <-
 
         if (is.null(predict_function)) {
             if (model$spec$engine %in% c("mboost", "survival", "glmnet", "flexsurv")){
-                predict_function <- function(model, newdata, times) {
-                    predict(model, new_data = newdata, type = "linear_pred")$.pred_linear_pred
-                }
+                predict_function <- function(model, newdata, times) predict(model, new_data = newdata, type = "linear_pred")$.pred_linear_pred
                 attr(predict_function, "verbose_info") <- "predict.model_fit with type = 'linear_pred' will be used"
             } else {
-                predict_function <- function(model, newdata, times) {
-                    rowSums(predict_cumulative_hazard_function(model, newdata, times = times))
-                }
+                predict_function <- function(model, newdata, times) rowSums(predict_cumulative_hazard_function(model, newdata, times = times))
                 attr(predict_function, "verbose_info") <- "sum over the predict_cumulative_hazard_function will be used"
             }
             attr(predict_function, "use.times") <- TRUE

R/metrics.R

@@ -27,38 +27,15 @@ loss_integrate <- function(loss_function, ..., normalization = NULL , max_quanti
         times <- times[quantile_mask]
         surv <- surv[,quantile_mask]
 
-        loss_values <- loss_function(y_true = y_true, risk = risk, surv = surv,times = times)
+        loss_values <- loss_function(y_true = y_true, risk = risk, surv = surv, times = times)
 
         na_mask <- (!is.na(loss_values))
 
         times <- times[na_mask]
         loss_values <- loss_values[na_mask]
         surv <- surv[na_mask]
 
-        n <- length(loss_values)
-        # integral using trapezoid method
-
-        if (is.null(normalization)){
-            tmp <- (loss_values[1:(n - 1)] + loss_values[2:n]) * diff(times) / 2
-            integrated_metric <- sum(tmp) / (max(times) - min(times))
-            return(integrated_metric)
-        }
-        else if (normalization == "t_max") {
-            tmp <- (loss_values[1:(n - 1)] + loss_values[2:n]) * diff(times) / 2
-            integrated_metric <- sum(tmp)
-            return(integrated_metric/max(times))
-        } else if (normalization == "survival"){
-
-            km <- survival::survfit(y_true ~ 1)
-            estimator <- stepfun(km$time, c(1, km$surv))
-
-            dwt <- 1 - estimator(times)
-
-            tmp <- (loss_values[1:(n - 1)] + loss_values[2:n]) * diff(dwt) / 2
-            integrated_metric <- sum(tmp)
-            return(integrated_metric/(1 - estimator(max(times))))
-        }
-
+        calculate_integral(loss_values, times, normalization, y_true=y_true)
     }
 
     attr(integrated_loss_function, "loss_type") <- "integrated"

R/model_profile.R

@@ -10,10 +10,11 @@
 #' @param ... other parameters passed to `DALEX::model_profile` if `output_type == "risk"`, otherwise ignored
 #' @param categorical_variables character, a vector of names of additional variables which should be treated as categorical (factors are automatically treated as categorical variables). If it contains variable names not present in the `variables` argument, they will be added at the end.
 #' @param grid_points maximum number of points for profile calculations. Note that the final number of points may be lower than grid_points. Will be passed to internal function. By default `51`.
+#' @param variable_splits_type character, decides how variable grids should be calculated. Use `"quantiles"` for percentiles or `"uniform"` (default) to get uniform grid of points.
 #' @param groups if `output_type == "risk"` a variable name that will be used for grouping. By default `NULL`, so no groups are calculated. If `output_type == "survival"` then ignored
 #' @param k passed to `DALEX::model_profile` if `output_type == "risk"`, otherwise ignored
 #' @param center logical, should profiles be centered before clustering
-#' @param type the type of variable profile. If `output_type == "survival"` then only `"partial"` is implemented, otherwise passed to `DALEX::model_profile`.
+#' @param type the type of variable profile, `"partial"` for Partial Dependence, `"accumulated"` for Accumulated Local Effects, or `"conditional"` (available only for `output_type == "risk"`)
 #' @param output_type either `"survival"` or `"risk"` the type of survival model output that should be considered for explanations. If `"survival"` the explanations are based on the survival function. Otherwise the scalar risk predictions are used by the `DALEX::model_profile` function.
 #'
 #' @return An object of class `model_profile_survival`. It is a list with the element `result` containing the results of the calculation.
@@ -39,7 +40,8 @@
 #' plot(cph_model_profile)
 #'
 #' rsf_model_profile <- model_profile(rsf_src_exp, output_type = "survival",
-#'                                    variables = c("age", "celltype"))
+#'                                    variables = c("age", "celltype"),
+#'                                    type = "accumulated")
 #'
 #' head(rsf_model_profile$result)
 #'
@@ -54,7 +56,7 @@ model_profile <- function(explainer,
                           ...,
                           groups = NULL,
                           k = NULL,
-                          center = TRUE,
+                          center = FALSE,
                           type = "partial",
                           output_type = "survival") UseMethod("model_profile", explainer)
 
@@ -66,6 +68,7 @@ model_profile.surv_explainer <- function(explainer,
                                          ...,
                                          categorical_variables = NULL,
                                          grid_points = 51,
+                                         variable_splits_type = "uniform",
                                          groups = NULL,
                                          k = NULL,
                                          center = TRUE,
@@ -84,28 +87,40 @@ model_profile.surv_explainer <- function(explainer,
                                                    type = type),
             "survival" = {
                 test_explainer(explainer, "model_profile", has_data = TRUE, has_survival = TRUE)
-
                 data <- explainer$data
                 if (!is.null(N) && N < nrow(data)) {
-                    ndata <- data[sample(1:nrow(data), N), ]
+                    ndata <- data[sample(1:nrow(data), N), , drop = FALSE]
                 } else {
-                    ndata <- data
+                    ndata <- data[1:nrow(data), , drop = FALSE]
                 }
 
-                cp_profiles <- surv_ceteris_paribus(explainer,
+                if (type == "partial"){
+                    cp_profiles <- surv_ceteris_paribus(explainer,
                                                         new_observation = ndata,
                                                         variables = variables,
                                                         categorical_variables = categorical_variables,
                                                         grid_points = grid_points,
+                                                        variable_splits_type = variable_splits_type,
                                                         ...)
+                    result <- surv_aggregate_profiles(cp_profiles, ...,
+                                                            variables = variables,
+                                                            center = center)
+                } else if (type == "accumulated"){
+                    cp_profiles <- list(variable_values = data.frame(ndata))
+                    result <- surv_ale(explainer,
+                                       data = ndata,
+                                       variables = variables,
+                                       categorical_variables = categorical_variables,
+                                       grid_points = grid_points,
+                                       ...)
+                } else {
+                    stop("Currently only `partial` and `accumulated` types are implemented")
+                }
 
-                agr_profiles <- surv_aggregate_profiles(cp_profiles, ...,
-                                                             groups = groups,
-                                                             type = type,
-                                                             variables = variables,
-                                                             center = center)
-
-                ret <- list(eval_times = unique(agr_profiles$`_times_`), cp_profiles = cp_profiles, result = agr_profiles)
+                ret <- list(eval_times = unique(result$`_times_`),
+                            cp_profiles = cp_profiles,
+                            result = result,
+                            type = type)
                 class(ret) <- c("model_profile_survival", "list")
                 ret$event_times <- explainer$y[explainer$y[, 1] <= max(explainer$times), 1]
                 ret$event_statuses <- explainer$y[explainer$y[, 1] <= max(explainer$times), 2]