Draft Surv to classif pipeline #194

DESCRIPTION

@@ -35,6 +35,10 @@ Authors@R:
              email = "[email protected]",
              role = "ctb",
              comment = c(ORCID = "0000-0001-7528-3795")),
+      person(given = "Philip",
+             family = "Studener",
+             role = "aut",
+             email = "[email protected]"),
       person(given = "Maximilian",
              family = "Muecke",
              email = "[email protected]",
@@ -79,7 +83,9 @@ Suggests:
     vdiffr,
     abind,
     Ecdat,
-    coxed
+    coxed,
+    mlr3learners,
+    pammtools
 LinkingTo:
     Rcpp
 Remotes:
@@ -133,13 +139,15 @@ Collate:
     'PipeOpBreslow.R'
     'PipeOpCrankCompositor.R'
     'PipeOpDistrCompositor.R'
+    'PipeOpPredClassifSurv.R'
     'PipeOpTransformer.R'
     'PipeOpPredTransformer.R'
     'PipeOpPredRegrSurv.R'
     'PipeOpPredSurvRegr.R'
     'PipeOpProbregrCompositor.R'
     'PipeOpSurvAvg.R'
     'PipeOpTaskRegrSurv.R'
+    'PipeOpTaskSurvClassif.R'
     'PipeOpTaskSurvRegr.R'
     'PipeOpTaskTransformer.R'
     'PredictionDataDens.R'

NAMESPACE

@@ -72,12 +72,14 @@ export(MeasureSurvXuR2)
 export(PipeOpBreslow)
 export(PipeOpCrankCompositor)
 export(PipeOpDistrCompositor)
+export(PipeOpPredClassifSurv)
 export(PipeOpPredRegrSurv)
 export(PipeOpPredSurvRegr)
 export(PipeOpPredTransformer)
 export(PipeOpProbregr)
 export(PipeOpSurvAvg)
 export(PipeOpTaskRegrSurv)
+export(PipeOpTaskSurvClassif)
 export(PipeOpTaskSurvRegr)
 export(PipeOpTaskTransformer)
 export(PipeOpTransformer)
@@ -95,6 +97,7 @@ export(as_task_surv)
 export(assert_surv)
 export(breslow)
 export(pecs)
+export(pipeline_survtoclassif)
 export(pipeline_survtoregr)
 export(plot_probregr)
 import(checkmate)

NEWS.md

@@ -1,7 +1,8 @@
 # mlr3proba 0.6.4
 
+* Add `PipeTaskSurvClassif`, `PipeOpPredClassifSurv` and `pipeline_survtoclassif` to transform a survival task into a classification task by discretizing the status.
 * Add useR! 2024 tutorial
-* Lots of refactoring, improve code quality (thanks to @m-muecke)
+* Lots of refactoring, improving code quality, migration to testthat v3, etc. (thanks to @m-muecke)
 
 # mlr3proba 0.6.3
 

R/LearnerSurvCoxPH.R

@@ -42,30 +42,20 @@ LearnerSurvCoxPH = R6Class("LearnerSurvCoxPH",
         pv$weights = task$weights$weight
       }
 
-      invoke(survival::coxph, formula = task$formula(), data = task$data(), .args = pv, x = TRUE)
+      invoke(survival::coxph, formula = task$formula(), data = task$data(),
+             .args = pv, x = TRUE)
     },
 
     .predict = function(task) {
-
-      newdata = task$data(cols = task$feature_names)
-
-      # We move the missingness checks here manually as if any NAs are made in predictions then the
-      # distribution object cannot be create (initialization of distr6 objects does not handle NAs)
-      if (anyMissing(newdata)) {
-        stopf(
-          "Learner %s on task %s failed to predict: Missing values in new data (line(s) %s)\n",
-          self$id, task$id,
-          toString(which(!complete.cases(newdata)))
-        )
-      }
-
+      newdata = ordered_features(task, self)
       pv = self$param_set$get_values(tags = "predict")
 
-      # Get predicted values
+      # Get survival predictions via `survfit`
       fit = invoke(survival::survfit, formula = self$model, newdata = newdata,
-        se.fit = FALSE, .args = pv)
+                   se.fit = FALSE, .args = pv)
 
-      lp = predict(self$model, type = "lp", newdata = newdata)
+      # Get linear predictors
+      lp = invoke(predict, self$model, type = "lp", newdata = newdata)
 
       .surv_return(times = fit$time, surv = t(fit$surv), lp = lp)
     }

R/LearnerSurvRpart.R

@@ -79,14 +79,15 @@ LearnerSurvRpart = R6Class("LearnerSurvRpart",
         pv = insert_named(pv, list(weights = task$weights$weight))
       }
 
-      invoke(rpart::rpart,
-        formula = task$formula(), data = task$data(),
-        method = "exp", .args = pv)
+      invoke(rpart::rpart, formula = task$formula(), data = task$data(),
+             method = "exp", .args = pv)
     },
 
     .predict = function(task) {
-      preds = invoke(predict, object = self$model, newdata = task$data(cols = task$feature_names))
-      list(crank = preds)
+      newdata = ordered_features(task, self)
+      p = invoke(predict, object = self$model, newdata = newdata)
+
+      list(crank = p)
     }
   )
 )

R/PipeOpPredClassifSurv.R

@@ -0,0 +1,112 @@
+#' @title PipeOpPredClassifSurv
+#' @name mlr_pipeops_trafopred_classifsurv
+#'
+#' @description
+#' Transform [PredictionClassif] to [PredictionSurv] by converting
+#' event probabilities of a pseudo status variable (discrete time hazards)
+#' to survival probabilities.
+#'
+#' @section Input and Output Channels:
+#' Input and output channels are inherited from [PipeOp][mlr3pipelines::PipeOp].
+#'
+#' The output is the input [PredictionClassif] transformed to a [PredictionSurv].
+#'
+#' @examples
+#' \dontrun{
+#' if (requireNamespace("mlr3pipelines", quietly = TRUE)) {
+#'   library(mlr3)
+#'   library(mlr3pipelines)
+#'
+#'   task = tsk("rats")
+#'
+#'   if (requireNamespace("mlr3learners", quietly = TRUE)) {
+#'     library(mlr3learners)
+#'     po_tasktoclassif = po("trafotask_survclassif")
+#'     po_tasktoclassif$train(list(task))
+#'     task_classif = po_tasktoclassif$predict(list(task))[[1]]
+#'     trafo_data = po_tasktoclassif$predict(list(task))[[2]]
+#'
+#'     learner = lrn("classif.log_reg", predict_type = "prob")
+#'     learner$train(task_classif)
+#'     pred = learner$predict(task_classif)
+#'
+#'     po_predtosurv = po("trafopred_classifsurv")
+#'     po_predtosurv$train(list(pred, trafo_data))
+#'     po_predtosurv$predict(list(pred, trafo_data))
+#'   }
+#' }
+#' }
+#' @family PipeOps
+#' @family Transformation PipeOps
+#' @export
+PipeOpPredClassifSurv = R6Class(
+  "PipeOpPredClassifSurv",
+  inherit = mlr3pipelines::PipeOp,
+
+  public = list(
+    #' @description
+    #' Creates a new instance of this [R6][R6::R6Class] class.
+    #' @param id (character(1))\cr
+    #' Identifier of the resulting object.
+    initialize = function(id = "trafopred_classifsurv") {
+      super$initialize(
+        id = id,
+        input = data.table::data.table(
+          name = c("input", "transformed_data"),
+          train = c("NULL", "data.table"),
+          predict = c("PredictionClassif", "data.table")
+        ),
+        output = data.table::data.table(
+          name = "output",
+          train = "NULL",
+          predict = "PredictionSurv"
+        )
+      )
+    }
+  ),
+
+  private = list(
+    .predict = function(input) {
+      pred = input[[1]]
+      data = input[[2]]
+      assert_true(!is.null(pred$prob))
+      data = cbind(data, pred = pred$prob[, "0"])
+
+      ## convert hazards to surv as prod(1 - h(t))
+      rows_per_id = nrow(data)/length(unique(data$id))
+      surv = t(vapply(unique(data$id), function(unique_id) {
+        x = cumprod((data[data$id == unique_id, ][["pred"]]))
+        x
+      }, numeric(rows_per_id)))
+
+      pred_list = list()
+      unique_end_times = sort(unique(data$tend))
+      ## coerce to distribution and crank
+      pred_list = .surv_return(times = unique_end_times, surv = surv)
+
+      # select the real tend values by only selecting the last row of each id
+      # basically a slightly more complex unique()
+      real_tend = data$time2[seq_len(nrow(data)) %% rows_per_id == 0]
+
+      # select last row for every id
+      data = as.data.table(data)
+      id = ped_status = NULL # to fix note
+      data = data[, .SD[.N, list(ped_status)], by = id]
+
+      ## create prediction object
+      p = PredictionSurv$new(
+        row_ids = seq_row(data),
+        crank = pred_list$crank, distr = pred_list$distr,
+        truth = Surv(real_tend, as.integer(as.character(data$ped_status))))
+
+      list(p)
+    },
+
+    .train = function(input) {
+      self$state = list()
+      list(input)
+    }
+  )
+)
+
+register_pipeop("trafopred_classifsurv", PipeOpPredClassifSurv)

R/PipeOpTaskSurvClassif.R

@@ -0,0 +1,150 @@
+#' @title PipeOpTaskSurvClassif
+#' @name mlr_pipeops_trafotask_survclassif
+#' @template param_pipelines
+#'
+#' @description
+#' Transform [TaskSurv] to [TaskClassif][mlr3::TaskClassif] by creating multiple
+#' interval observations for each subject based on `cut`, with a `ped_status` variable
+#' indicating whether an event occurred in each interval.
+#'
+#' @section Input and Output Channels:
+#' Input and output channels are inherited from [PipeOp][mlr3pipelines::PipeOp].
+#'
+#' The output is the input [TaskSurv] transformed to a [TaskClassif][mlr3::TaskClassif]
+#' as well as the transformed data during prediction.
+#'
+#' @section State:
+#' The `$state` contains information about the `cut` parameter used
+#' as well as `time_var` and `event_var`, the names of the two target
+#' columns of the survival task.
+#'
+#' @section Parameters:
+#' The parameters are
+#'
+#' * `cut :: numeric()`\cr
+#' Split points, used to partition the data into intervals.
+#' If unspecified, all unique event times will be used.
+#' If `cut` is a single integer, it will be interpreted as the number of equidistant
+#' intervals from 0 until the maximum event time.
+#' * `max_time :: numeric(1)`\cr
+#' If cut is unspecified, this will be the last possible event time.
+#' All event times after max_time will be administratively censored at max_time.
+#' Needs to be greater than the minimum event time.
+#'
+#' @examples
+#' \dontrun{
+#' if (requireNamespace("mlr3pipelines", quietly = TRUE)) {
+#'   library(mlr3)
+#'   library(mlr3pipelines)
+#'
+#'   task = tsk("lung")
+#'   po = po("trafotask_survclassif")
+#'   po$train(list(task))
+#'   po$predict(list(task))[[1]]
+#' }
+#' }
+#'
+#' @references
+#' `r format_bib("tutz_2016")`
+#'
+#' @family PipeOps
+#' @family Transformation PipeOps
+#' @export
+PipeOpTaskSurvClassif = R6Class("PipeOpTaskSurvClassif",
+  inherit = mlr3pipelines::PipeOp,
+
+  public = list(
+    #' @description
+    #' Creates a new instance of this [R6][R6::R6Class] class.
+    initialize = function(id = "trafotask_survclassif") {
+      param_set = ps(
+        cut = p_uty(default = NULL),
+        max_time = p_dbl(default = NULL, special_vals = list(NULL))
+      )
+      super$initialize(
+        id = id,
+        param_set = param_set,
+        input = data.table::data.table(
+          name    = "input",
+          train   = "TaskSurv",
+          predict = "TaskSurv"
+        ),
+        output = data.table::data.table(
+          name    = c("output", "transformed_data"),
+          train   = c("TaskClassif", "data.table"),
+          predict = c("TaskClassif", "data.table")
+        )
+      )
+    }
+  ),
+
+  private = list(
+    .train = function(input) {
+      task = input[[1]]
+      data = task$data()
+      assert_true(task$censtype == "right")
+
+      cut = assert_numeric(self$param_set$values$cut, null.ok = TRUE, lower = 0)
+      max_time = self$param_set$values$max_time
+
+      time_var = task$target_names[1]
+      event_var = task$target_names[2]
+      if (testInt(cut, lower = 1)) {
+        cut = seq(0, data[get(event_var) == 1, max(get(time_var))], length.out = cut + 1)
+      }
+      if (!is.null(max_time)) {
+        assert(max_time > data[get(event_var) == 1, min(get(time_var))],
+               "max_time must be greater than the minimum event time.")
+      }
+
+      form = mlr3misc::formulate(sprintf("Surv(%s, %s)", time_var, event_var), ".")
+
+      # TODO: do without pammtools
+      long_data = pammtools::as_ped(data = data, formula = form, cut = cut, max_time = max_time)
+      self$state$cut = attributes(long_data)$trafo_args$cut
+      self$state$event_var = event_var
+      self$state$time_var = time_var
+      long_data = as.data.table(long_data)
+      long_data$ped_status = factor(long_data$ped_status, levels = c("0", "1"))
+
+      # remove offset, tstart, interval for dataframe long_data
+      long_data[, c("offset", "tstart", "interval") := NULL]
+
+      task = TaskClassif$new(paste0(task$id, "_disc"), long_data, target = "ped_status", positive = "1")
+      task$set_col_roles("id", roles = "name")
+
+      list(task, data.table())
+    },
+
+    .predict = function(input) {
+      task = input[[1]]
+      data = task$data()
+
+      # extract required data from `state`
+      cut = self$state$cut
+      time_var = self$state$time_var
+      event_var = self$state$event_var
+
+      max_time = max(cut)
+      time = data[[time_var]]
+      data[[time_var]] = max_time
+
+      # update form
+      form = mlr3misc::formulate(sprintf("Surv(%s, %s)", time_var, event_var), ".")
+
+      new_data = pammtools::as_ped(data, formula = form, cut = cut)
+      new_data = as.data.table(new_data)
+      new_data$ped_status = factor(new_data$ped_status, levels = c("0", "1"))
+
+      # remove offset, tstart, interval for dataframe long_data
+      new_data[, c("offset", "tstart", "interval") := NULL]
+      task = TaskClassif$new(paste0(task$id, "_disc"), new_data, target = "ped_status", positive = "1")
+      task$set_col_roles("id", roles = "name")
+
+      new_data$time2 = rep(time, each = sum(new_data$id == 1))
+      list(task, new_data)
+    }
+  )
+)
+
+register_pipeop("trafotask_survclassif", PipeOpTaskSurvClassif)