tlverse / sl3
Showing 12 of 198 files from the diff.
Other files ignored by Codecov
man/Lrnr_glm.Rd has changed.
man/Stack.Rd has changed.
man/Lrnr_hts.Rd has changed.
NAMESPACE has changed.
man/Lrnr_pca.Rd has changed.
man/Lrnr_gam.Rd has changed.
man/Lrnr_nnls.Rd has changed.
docs/index.html has changed.
docs/pkgdown.yml has changed.
man/Lrnr_solnp.Rd has changed.
man/Lrnr_tsDyn.Rd has changed.
man/importance.Rd has changed.
man/Lrnr_arima.Rd has changed.
man/Lrnr_svm.Rd has changed.
man/Lrnr_gts.Rd has changed.
man/Lrnr_base.Rd has changed.
man/Lrnr_gbm.Rd has changed.
man/Lrnr_rpart.Rd has changed.
docs/sitemap.xml has changed.
man/Lrnr_grf.Rd has changed.
man/Lrnr_nnet.Rd has changed.
man/Lrnr_caret.Rd has changed.
man/Lrnr_earth.Rd has changed.
man/CV_lrnr_sl.Rd has changed.
man/cv_risk.Rd has changed.
man/Lrnr_sl.Rd has changed.
README.md has changed.
man/Lrnr_optim.Rd has changed.
NEWS.md has changed.
man/Lrnr_cv.Rd has changed.
man/Lrnr_mean.Rd has changed.
README.Rmd has changed.
DESCRIPTION has changed.
man/Pipeline.Rd has changed.

@@ -4,11 +4,10 @@
Loading
4 4
#'
5 5
#' @docType class
6 6
#'
7 -
#' @importFrom R6 R6Class
8 -
#' @importFrom assertthat assert_that is.count is.flag
9 -
#'
10 7
#' @export
11 8
#'
9 +
#' @importFrom origami make_folds
10 +
#'
12 11
#' @keywords data
13 12
#'
14 13
#' @return Learner object with methods for training and prediction. See
@@ -31,14 +30,8 @@
Loading
31 30
  class = TRUE,
32 31
  public = list(
33 32
    initialize = function(categorical_learner = NULL, type = "equal_mass",
34 -
                          n_bins = 20, ...) {
35 -
      if (is.null(categorical_learner)) {
36 -
        categorical_learner <- make_learner(Lrnr_glmnet)
37 -
      }
38 -
      params <- list(
39 -
        type = type, n_bins = n_bins,
40 -
        categorical_learner = categorical_learner, ...
41 -
      )
33 +
                          n_bins = 20, breaks = NULL, ...) {
34 +
      params <- args_to_list()
42 35
      super$initialize(params = params, ...)
43 36
    }
44 37
  ),
@@ -48,21 +41,34 @@
Loading
48 41
      discretized <- discretize_variable(task$Y,
49 42
        type = self$params$type,
50 43
        n_bins = self$params$n_bins,
51 -
        breaks = self$params_breaks
44 +
        breaks = self$params$breaks
52 45
      )
53 46
54 47
      # make discretized task
55 -
      new_columns <-
56 -
        task$add_columns(data.table(
57 -
          discrete_Y =
58 -
            factor(discretized$x_discrete)
59 -
        ))
48 +
      folds <- tryCatch(
49 +
        {
50 +
          origami::make_folds(strata_ids = factor(discretized$x_discrete))
51 +
        },
52 +
        warning = function(c) {
53 +
          message("Cannot construct stratified CV due to insufficient sample size for at least one level of discretized outcome, using folds from original task")
54 +
          task$folds
55 +
        }
56 +
      )
57 +
      new_columns <- task$add_columns(
58 +
        data.table(discrete_Y = factor(discretized$x_discrete))
59 +
      )
60 60
      discrete_task <- task$next_in_chain(
61 61
        outcome = "discrete_Y",
62 -
        column_names = new_columns
62 +
        column_names = new_columns,
63 +
        folds = folds
63 64
      )
64 65
      # fit categorical learner to discretized task
65 -
      categorical_fit <- self$params$categorical_learner$train(discrete_task)
66 +
      categorical_learner <- self$params$categorical_learner
67 +
      if (is.null(categorical_learner)) {
68 +
        message("Categorical learner NULL, fitting with LASSO regression")
69 +
        categorical_learner <- make_learner(Lrnr_glmnet)
70 +
      }
71 +
      categorical_fit <- categorical_learner$train(discrete_task)
66 72
67 73
      fit_object <- list(
68 74
        categorical_fit = categorical_fit,
@@ -75,14 +81,22 @@
Loading
75 81
      discretized <- discretize_variable(task$Y,
76 82
        breaks = self$fit_object$breaks
77 83
      )
78 -
      new_columns <-
79 -
        task$add_columns(data.table(
80 -
          discrete_Y =
81 -
            factor(discretized$x_discrete)
82 -
        ))
84 +
      folds <- tryCatch(
85 +
        {
86 +
          origami::make_folds(strata_ids = factor(discretized$x_discrete))
87 +
        },
88 +
        warning = function(c) {
89 +
          message("Cannot construct stratified CV due to insufficient sample size for at least one level of discretized outcome, using folds from original task")
90 +
          task$folds
91 +
        }
92 +
      )
93 +
      new_columns <- task$add_columns(
94 +
        data.table(discrete_Y = factor(discretized$x_discrete))
95 +
      )
83 96
      discrete_task <- task$next_in_chain(
84 97
        outcome = "discrete_Y",
85 -
        column_names = new_columns
98 +
        column_names = new_columns,
99 +
        folds = folds
86 100
      )
87 101
88 102
      # predict categorical learner on discretized task

@@ -1,10 +1,14 @@
Loading
1 1
#' Nonlinear Optimization via Augmented Lagrange
2 2
#'
3 -
#' This meta-learner provides fitting procedures for any pairing of loss
3 +
#' This meta-learner provides fitting procedures for any pairing of loss or risk
4 4
#' function and metalearner function, subject to constraints. The optimization
5 5
#' problem is solved by making use of \code{\link[Rsolnp]{solnp}}, using
6 -
#' Lagrange multipliers. For further details, consult the documentation of the
7 -
#' \code{Rsolnp} package.
6 +
#' Lagrange multipliers. An important note from the \code{\link[Rsolnp]{solnp}}
7 +
#' documentation states that the control parameters \code{tol} and \code{delta}
8 +
#' are key in getting any possibility of successful convergence, therefore it
9 +
#' is suggested that the user change these appropriately to reflect their
10 +
#' problem specification. For further details, consult the documentation of the
11 +
#' \pkg{Rsolnp} package.
8 12
#'
9 13
#' @docType class
10 14
#'
@@ -14,55 +18,80 @@
Loading
14 18
#'
15 19
#' @keywords data
16 20
#'
17 -
#' @return Learner object with methods for training and prediction. See
18 -
#'  \code{\link{Lrnr_base}} for documentation on learners.
21 +
#' @return A learner object inheriting from \code{\link{Lrnr_base}} with
22 +
#'  methods for training and prediction. For a full list of learner
23 +
#'  functionality, see the complete documentation of \code{\link{Lrnr_base}}.
19 24
#'
20 -
#' @format \code{\link{R6Class}} object.
25 +
#' @format An \code{\link[R6]{R6Class}} object inheriting from
26 +
#'  \code{\link{Lrnr_base}}.
21 27
#'
22 28
#' @family Learners
23 29
#'
24 30
#' @section Parameters:
25 -
#' \describe{
26 -
#'   \item{\code{learner_function=metalearner_linear}}{A function(alpha, X) that
27 -
#'     takes a vector of covariates and a matrix of data and combines them into
28 -
#'     a vector of predictions. See \link{metalearners} for options.}
29 -
#'   \item{\code{loss_function=loss_squared_error}}{A function(pred, truth) that
30 -
#'     takes prediction and truth vectors and returns a loss vector. See
31 -
#'     \link{loss_functions} for options.}
32 -
#'   \item{\code{make_sparse=TRUE}}{If TRUE, zeros out small alpha values.}
33 -
#'   \item{\code{convex_combination=TRUE}}{If \code{TRUE}, constrain alpha to
34 -
#'     sum to 1.}
35 -
#'   \item{\code{init_0=FALSE}}{If TRUE, alpha is initialized to all 0's, useful
36 -
#'     for TMLE. Otherwise, it is initialized to equal weights summing to 1,
37 -
#'     useful for SuperLearner.}
38 -
#'   \item{\code{...}}{Not currently used.}
39 -
#' }
31 +
#'   - \code{learner_function = metalearner_linear}: A function(alpha, X) that
32 +
#'       takes a vector of covariates and a matrix of data and combines them
33 +
#'       into a vector of predictions. See \code{\link{metalearners}} for
34 +
#'       options.
35 +
#'   - \code{eval_function = loss_squared_error}: A function(pred, truth) that
36 +
#'       takes prediction and truth vectors and returns a loss vector or a risk
37 +
#'       scalar. See \code{\link{loss_functions}} and
38 +
#'       \code{\link{risk_functions}} for options and more detail.
39 +
#'   - \code{make_sparse = TRUE}: If \code{TRUE}, zeros out small alpha values.
40 +
#'   - \code{convex_combination = TRUE}: If \code{TRUE}, constrain alpha to sum
41 +
#'       to 1.
42 +
#'   - \code{init_0 = FALSE}: If \code{TRUE}, alpha is initialized to all 0's,
43 +
#'       useful for TMLE. Otherwise, it is initialized to equal weights summing
44 +
#'       to 1, useful for Super Learner.
45 +
#'   - \code{outer.iter = 400}: Maximum number of major (outer) iterations.
46 +
#'   - \code{inner.iter = 800}: Maximum number of minor (inner) iterations.
47 +
#'   - \code{delta = 1e-7}:Relative step size in forward difference evaluation.
48 +
#'   - \code{tol = 1e-8}: Relative tolerance on feasibility and optimality.
49 +
#'   - \code{trace = FALSE}: The value of the objective function and the
50 +
#'       parameters are printed at every major iteration.
51 +
#'   - \code{...}: Additional arguments defined in \code{\link{Lrnr_base}},
52 +
#'     such as \code{params} (like \code{formula}) and \code{name}.
40 53
#'
41 -
#' @template common_parameters
42 -
#
54 +
#' @examples
55 +
#' # define ML task
56 +
#' data(cpp_imputed)
57 +
#' covs <- c("apgar1", "apgar5", "parity", "gagebrth", "mage", "meducyrs")
58 +
#' task <- sl3_Task$new(cpp_imputed, covariates = covs, outcome = "haz")
59 +
#'
60 +
#' # build relatively fast learner library (not recommended for real analysis)
61 +
#' lasso_lrnr <- Lrnr_glmnet$new()
62 +
#' glm_lrnr <- Lrnr_glm$new()
63 +
#' ranger_lrnr <- Lrnr_ranger$new()
64 +
#' lrnrs <- c(lasso_lrnr, glm_lrnr, ranger_lrnr)
65 +
#' names(lrnrs) <- c("lasso", "glm", "ranger")
66 +
#' lrnr_stack <- make_learner(Stack, lrnrs)
67 +
#'
68 +
#' # instantiate SL with GA metalearner
69 +
#' solnp_meta <- Lrnr_solnp$new()
70 +
#' sl <- Lrnr_sl$new(lrnr_stack, solnp_meta)
71 +
#' sl_fit <- sl$train(task)
43 72
Lrnr_solnp <- R6Class(
44 73
  classname = "Lrnr_solnp",
45 74
  inherit = Lrnr_base, portable = TRUE,
46 75
  class = TRUE,
47 76
  public = list(
48 77
    initialize = function(learner_function = metalearner_linear,
49 -
                          loss_function = loss_squared_error,
78 +
                          eval_function = loss_squared_error,
50 79
                          make_sparse = TRUE, convex_combination = TRUE,
51 -
                          init_0 = FALSE, tol = 1e-5, ...) {
80 +
                          init_0 = FALSE, outer.iter = 400, inner.iter = 800,
81 +
                          delta = 1e-7, tol = 1e-8, trace = FALSE, ...) {
52 82
      params <- args_to_list()
53 -
      super$initialize(params = params, ...)
83 +
      super$initialize(params = params)
54 84
    }
55 85
  ),
56 86
  private = list(
57 87
    .properties = c(
58 -
      "continuous", "binomial", "categorical", "weights",
59 -
      "offset"
88 +
      "continuous", "binomial", "categorical", "weights", "offset"
60 89
    ),
61 90
    .train = function(task) {
62 91
      verbose <- getOption("sl3.verbose")
63 92
      params <- self$params
64 93
      learner_function <- params$learner_function
65 -
      loss_function <- params$loss_function
94 +
      eval_function <- params$eval_function
66 95
      outcome_type <- self$get_outcome_type(task)
67 96
68 97
      # specify data
@@ -82,8 +111,18 @@
Loading
82 111
        } else {
83 112
          preds <- learner_function(alphas, X)
84 113
        }
85 -
        losses <- loss_function(preds, Y)
86 -
        risk <- weighted.mean(losses, weights)
114 +
        eval_result <- eval_function(preds, Y)
115 +
116 +
        if (!is.null(attr(eval_result, "loss")) && !attr(eval_result, "loss")) {
117 +
          risk <- eval_result
118 +
        } else {
119 +
          loss <- eval_result
120 +
          risk <- weighted.mean(loss, weights)
121 +
        }
122 +
        if (!is.null(attr(eval_result, "optimize")) &&
123 +
          attr(eval_result, "optimize") == "maximize") {
124 +
          risk <- risk * -1
125 +
        }
87 126
        return(risk)
88 127
      }
89 128
      if (params$convex_combination) {
@@ -108,7 +147,10 @@
Loading
108 147
        init_alphas, risk,
109 148
        eqfun = eq_fun, eqB = eqB,
110 149
        LB = LB,
111 -
        control = list(trace = 0, tol = params$tol)
150 +
        control = list(
151 +
          outer.iter = params$outer.iter, inner.iter = params$inner.iter,
152 +
          delta = params$delta, tol = params$tol, trace = params$trace
153 +
        )
112 154
      )
113 155
      coefs <- fit_object$pars
114 156
      names(coefs) <- colnames(task$X)
@@ -117,7 +159,10 @@
Loading
117 159
        max_coef <- max(coefs)
118 160
        threshold <- max_coef / 1000
119 161
        coefs[coefs < threshold] <- 0
120 -
        coefs <- coefs / sum(coefs)
162 +
        if (params$convex_combination) {
163 +
          # renormalize so coefficients sum to 1
164 +
          coefs <- coefs / sum(coefs)
165 +
        }
121 166
      }
122 167
      fit_object$coefficients <- coefs
123 168
      fit_object$training_offset <- task$has_node("offset")

@@ -34,7 +34,7 @@
Loading
34 34
      ylab("Predicted") +
35 35
      theme_bw() +
36 36
      geom_smooth(se = FALSE)
37 -
  } else {
37 +
  } else if (outcome_type$type == "categorical") {
38 38
    unpacked <- unpack_predictions(predictions)
39 39
    unpacked <- as.data.table(unpacked)
40 40
    setnames(unpacked, outcome_type$levels)
@@ -57,6 +57,17 @@
Loading
57 57
      ylab("True Positive Rate") +
58 58
      scale_color_discrete("Observed") +
59 59
      coord_equal()
60 +
  } else if (outcome_type$type == "binomial") {
61 +
    pred_data <- data.table(pred = predictions, obs = observed)
62 +
    pred_plot <- ggplot(pred_data, aes_(x = ~obs, y = ~pred, fill = ~obs)) +
63 +
      geom_violin(trim = FALSE) +
64 +
      geom_boxplot(width = 0.1, fill = "white") +
65 +
      xlab("Observed") +
66 +
      ylab("Predicted") +
67 +
      theme_bw() +
68 +
      theme(legend.position = "none")
69 +
  } else {
70 +
    stop(sprintf("Outcome type of %s is not supported", outcome_type$type))
60 71
  }
61 72
62 73
  return(pred_plot)

@@ -80,8 +80,13 @@
Loading
80 80
        } else {
81 81
          preds <- learner_function(alphas, X)
82 82
        }
83 -
        losses <- loss_function(preds, Y)
84 -
        risk <- weighted.mean(losses, weights)
83 +
        eval_result <- loss_function(preds, Y)
84 +
        if (!is.null(attr(eval_result, "loss")) && !attr(eval_result, "loss")) {
85 +
          risk <- eval_result
86 +
        } else {
87 +
          losses <- eval_result
88 +
          risk <- weighted.mean(losses, weights)
89 +
        }
85 90
        return(risk)
86 91
      }
87 92
      p <- ncol(X)

@@ -0,0 +1,196 @@
Loading
1 +
#' Nonlinear Optimization via Genetic Algorithm (GA)
2 +
#'
3 +
#' This metalearner provides fitting procedures for any pairing of loss or risk
4 +
#' function and metalearner function, subject to constraints. The optimization
5 +
#' problem is solved by making use of the \code{\link[GA]{ga}} function in the
6 +
#' \pkg{GA} R package. For further consult the documentation of this package.
7 +
#'
8 +
#' @docType class
9 +
#'
10 +
#' @importFrom R6 R6Class
11 +
#'
12 +
#' @export
13 +
#'
14 +
#' @keywords data
15 +
#'
16 +
#' @return A learner object inheriting from \code{\link{Lrnr_base}} with
17 +
#'  methods for training and prediction. For a full list of learner
18 +
#'  functionality, see the complete documentation of \code{\link{Lrnr_base}}.
19 +
#'
20 +
#' @format An \code{\link[R6]{R6Class}} object inheriting from
21 +
#'  \code{\link{Lrnr_base}}.
22 +
#'
23 +
#' @family Learners
24 +
#'
25 +
#' @section Parameters:
26 +
#'   - \code{learner_function = metalearner_linear}: A function(alpha, X) that
27 +
#'       takes a vector of covariates and a matrix of data and combines them
28 +
#'       into a vector of predictions. See \code{\link{metalearners}} for
29 +
#'       options.
30 +
#'   - \code{eval_function = loss_squared_error}: A function(pred, truth) that
31 +
#'       takes prediction and truth vectors and returns a loss vector or a risk
32 +
#'       scalar. See \code{\link{loss_functions}} and
33 +
#'       \code{\link{risk_functions}} for options and more detail.
34 +
#'   - \code{make_sparse = TRUE}: If \code{TRUE}, zeros out small alpha values.
35 +
#'   - \code{convex_combination = TRUE}: If \code{TRUE}, constrain alpha to sum
36 +
#'       to 1.
37 +
#'   - \code{maxiter = 100}: The maximum number of iterations to run before the
38 +
#'       GA search is halted.
39 +
#'   - \code{run = 10}: The number of consecutive generations without any
40 +
#'       improvement in the best fitness value before the GA is stopped.
41 +
#'   - \code{optim = TRUE}: A logical determining whether or not a local search
42 +
#'       using general-purpose optimization algorithms should be used. Argument
43 +
#'       \code{optimArgs} of \code{\link[GA]{ga}} provides further details and
44 +
#'       finer control.
45 +
#'   - \code{...}: Additional arguments to \code{\link[GA]{ga}} and/or
46 +
#'       \code{\link{Lrnr_base}}.
47 +
#'
48 +
#' @examples
49 +
#' # define ML task
50 +
#' data(cpp_imputed)
51 +
#' covs <- c("apgar1", "apgar5", "parity", "gagebrth", "mage", "meducyrs")
52 +
#' task <- sl3_Task$new(cpp_imputed, covariates = covs, outcome = "haz")
53 +
#'
54 +
#' # build relatively fast learner library (not recommended for real analysis)
55 +
#' lasso_lrnr <- Lrnr_glmnet$new()
56 +
#' glm_lrnr <- Lrnr_glm$new()
57 +
#' ranger_lrnr <- Lrnr_ranger$new()
58 +
#' lrnrs <- c(lasso_lrnr, glm_lrnr, ranger_lrnr)
59 +
#' names(lrnrs) <- c("lasso", "glm", "ranger")
60 +
#' lrnr_stack <- make_learner(Stack, lrnrs)
61 +
#'
62 +
#' # instantiate SL with GA metalearner
63 +
#' ga <- Lrnr_ga$new()
64 +
#' sl <- Lrnr_sl$new(lrnr_stack, ga)
65 +
#' sl_fit <- sl$train(task)
66 +
Lrnr_ga <- R6Class(
67 +
  classname = "Lrnr_ga",
68 +
  inherit = Lrnr_base, portable = TRUE,
69 +
  class = TRUE,
70 +
  public = list(
71 +
    initialize = function(learner_function = metalearner_linear,
72 +
                          eval_function = loss_squared_error,
73 +
                          make_sparse = TRUE, convex_combination = TRUE,
74 +
                          maxiter = 100, run = 10, optim = TRUE, ...) {
75 +
      params <- args_to_list()
76 +
      super$initialize(params = params)
77 +
    }
78 +
  ),
79 +
  private = list(
80 +
    .properties = c(
81 +
      "continuous", "binomial", "categorical", "weights", "offset"
82 +
    ),
83 +
    .train = function(task) {
84 +
      verbose <- getOption("sl3.verbose")
85 +
      params <- self$params
86 +
      learner_function <- params$learner_function
87 +
      eval_function <- params$eval_function
88 +
      outcome_type <- self$get_outcome_type(task)
89 +
90 +
      # specify data
91 +
      X <- as.matrix(task$X)
92 +
      Y <- outcome_type$format(task$Y)
93 +
94 +
      if (task$has_node("offset")) {
95 +
        offset <- task$offset
96 +
      } else {
97 +
        offset <- NULL
98 +
      }
99 +
100 +
      weights <- task$weights
101 +
102 +
      # Borrow the risk code from Lrnr_solnp
103 +
      # NB: We enforce convex combination by rescaling inside risk calculation
104 +
      # Probably better to use lagrange multipliers
105 +
      risk <- function(alphas) {
106 +
        if (sum(alphas) == 0) {
107 +
          return(NA)
108 +
        }
109 +
        alphas <- alphas / sum(alphas)
110 +
        if (!is.null(offset)) {
111 +
          preds <- learner_function(alphas, X, offset)
112 +
        } else {
113 +
          preds <- learner_function(alphas, X)
114 +
        }
115 +
        eval_result <- eval_function(preds, Y)
116 +
117 +
        if (!is.null(attr(eval_result, "loss")) && !attr(eval_result, "loss")) {
118 +
          risk <- eval_result
119 +
        } else {
120 +
          loss <- eval_result
121 +
          risk <- weighted.mean(loss, weights)
122 +
        }
123 +
        if (!is.null(attr(eval_result, "optimize")) &&
124 +
          attr(eval_result, "optimize") == "maximize") {
125 +
          risk <- risk * -1
126 +
        }
127 +
        return(risk)
128 +
      }
129 +
130 +
      # build a matrix of suggestions
131 +
      p <- ncol(X)
132 +
      discrete <- diag(p) # first all the discrete SL solutions
133 +
      equal <- rep(1 / p, p) # then equal weights
134 +
135 +
      # maybe a nnls for good measure
136 +
      nnls_coef <- tryCatch(
137 +
        {
138 +
          nnls_fit <- nnls(X, Y)
139 +
          coef(nnls_fit)
140 +
        },
141 +
        error = function(error) {
142 +
          return(equal)
143 +
        }
144 +
      )
145 +
146 +
      suggestions <- rbind(discrete, equal, nnls_coef)
147 +
148 +
      # note we flip back to fitness because GA is a maximizer
149 +
      args <- c(list(
150 +
        type = "real-valued", fitness = function(x) {
151 +
          -1 * risk(x)
152 +
        },
153 +
        lower = rep(0, p), upper = rep(1, p), suggestions = suggestions,
154 +
        popSize = 10 * p, keepBest = TRUE
155 +
      ), params)
156 +
      GA1 <- call_with_args(
157 +
        GA::ga, args,
158 +
        ignore = c(
159 +
          "learner_function", "eval_function", "make_sparse",
160 +
          "convex_combination"
161 +
        )
162 +
      )
163 +
164 +
      coefs <- as.vector(GA1@bestSol[[1]])
165 +
      names(coefs) <- colnames(task$X)
166 +
167 +
      fit_object <- list(ga_fit <- GA1)
168 +
      if (params$make_sparse) {
169 +
        max_coef <- max(coefs)
170 +
        threshold <- max_coef / 1000
171 +
        coefs[coefs < threshold] <- 0
172 +
      }
173 +
      if (params$convex_combination) {
174 +
        # renormalize so coefficients sum to 1
175 +
        coefs <- coefs / sum(coefs)
176 +
      }
177 +
      fit_object$coefficients <- coefs
178 +
      fit_object$training_offset <- task$has_node("offset")
179 +
      fit_object$name <- "ga"
180 +
      return(fit_object)
181 +
    },
182 +
    .predict = function(task = NULL) {
183 +
      verbose <- getOption("sl3.verbose")
184 +
      X <- as.matrix(task$X)
185 +
      alphas <- self$fit_object$coefficients
186 +
187 +
      if (self$fit_object$training_offset) {
188 +
        predictions <- self$params$learner_function(alphas, X, task$offset)
189 +
      } else {
190 +
        predictions <- self$params$learner_function(alphas, X)
191 +
      }
192 +
      return(predictions)
193 +
    },
194 +
    .required_packages = c("GA")
195 +
  )
196 +
)

@@ -91,8 +91,8 @@
Loading
91 91
      params <- list(learner = learner, folds = folds, full_fit = full_fit, ...)
92 92
      super$initialize(params = params, ...)
93 93
    },
94 -
    cv_risk = function(loss_fun) {
95 -
      return(cv_risk(self, loss_fun))
94 +
    cv_risk = function(eval_fun) {
95 +
      return(cv_risk(self, eval_fun))
96 96
    },
97 97
    print = function() {
98 98
      print("Lrnr_cv")

@@ -16,6 +16,7 @@
Loading
16 16
#' @export
17 17
loss_squared_error <- function(pred, observed) {
18 18
  out <- (pred - observed)^2
19 +
  attributes(out)$name <- "MSE"
19 20
  return(out)
20 21
}
21 22
@@ -28,6 +29,7 @@
Loading
28 29
#' @export
29 30
loss_loglik_true_cat <- function(pred, observed) {
30 31
  out <- -log(bound(pred))
32 +
  attributes(out)$name <- "NLL"
31 33
  return(out)
32 34
}
33 35
@@ -38,6 +40,7 @@
Loading
38 40
#' @export
39 41
loss_loglik_binomial <- function(pred, observed) {
40 42
  out <- -1 * ifelse(observed == 1, log(bound(pred)), log(1 - bound(pred)))
43 +
  attributes(out)$name <- "NLL"
41 44
  return(out)
42 45
}
43 46
@@ -53,7 +56,9 @@
Loading
53 56
  index_mat <- cbind(seq_along(observed), observed)
54 57
  unpacked <- unpack_predictions(pred)
55 58
  class_liks <- log(bound(unpacked[index_mat]))
56 -
  return(-1 * class_liks)
59 +
  out <- -1 * class_liks
60 +
  attributes(out)$name <- "NLL"
61 +
  return(out)
57 62
}
58 63
59 64
#' SQUARED ERROR LOSS FOR MULTIVARIATE (LOSS AVERAGED ACROSS OUTCOMES)
@@ -66,6 +71,7 @@
Loading
66 71
loss_squared_error_multivariate <- function(pred, observed) {
67 72
  unpacked <- unpack_predictions(pred)
68 73
  losses <- rowMeans((unpacked - observed)^2)
74 +
  attributes(losses)$name <- "MSE"
69 75
  return(losses)
70 76
}
71 77
@@ -82,28 +88,110 @@
Loading
82 88
#'
83 89
#' @export
84 90
risk <- function(pred, observed, loss = loss_squared_error, weights = NULL) {
91 +
  # get losses
92 +
  losses <- loss(pred, observed)
93 +
94 +
  # get weights
85 95
  if (is.null(weights)) {
86 96
    weights <- rep(1, length(observed))
87 97
  }
88 -
  risk <- weighted.mean(loss(observed, pred), weights)
98 +
99 +
  # calculate risk
100 +
  risk <- weighted.mean(losses, weights)
89 101
  return(risk)
90 102
}
91 103
104 +
#' FACTORY RISK FUNCTION FOR ROCR PERFORMANCE MEASURES WITH BINARY OUTCOMES
105 +
#'
106 +
#' Factory function for estimating an ROCR-based risk for a given ROCR measure,
107 +
#' and the risk is defined as one minus the performance measure.
108 +
#'
109 +
#' @note This risk does not take into account weights. In order to use this
110 +
#' risk, it must first be instantiated with respect to the \pkg{ROCR}
111 +
#' performance measure of interest, and then the user-defined function can be
112 +
#' used.
113 +
#'
114 +
#' @name risk_functions
115 +
#'
116 +
#' @param pred A vector of predicted values.
117 +
#' @param observed A vector of binary observed values.
118 +
#' @param measure A character indicating which \code{ROCR} performance measure
119 +
#' to use for evaluation. The \code{measure} must be either cutoff-dependent
120 +
#' so a single value can be selected (e.g., "tpr"), or it's value is a scalar
121 +
#' (e.g., "aucpr"). For more information, see \code{\link[ROCR]{performance}}.
122 +
#' @param cutoff A numeric value specifying the cutoff for choosing a single
123 +
#' performance measure from the returned set. Only used for performance measures
124 +
#' that are cutoff-dependent and default is 0.5. See
125 +
#' \code{\link[ROCR]{performance}} for more detail.
126 +
#' @param name An optional character string for user to supply their desired
127 +
#' name for the performance measure, which will be used for naming subsequent
128 +
#' risk-related tables and metrics (e.g., \code{cv_risk} column names). When
129 +
#' \code{name} is not supplied, the \code{measure} will be used for naming.
130 +
#' @param ... Optional arguments to specific \code{ROCR} performance
131 +
#' measures. See \code{\link[ROCR]{performance}} for more detail.
132 +
#'
133 +
#' @rdname risk_functions
134 +
#'
135 +
#' @importFrom ROCR prediction performance
136 +
#'
137 +
#' @export
138 +
#'
139 +
custom_ROCR_risk <- function(measure, cutoff = 0.5, name = NULL, ...) {
140 +
  function(pred, observed) {
141 +
142 +
    # remove NA, NaN, Inf values
143 +
    if (any(is.na(pred)) | any(is.nan(pred)) | any(is.infinite(pred))) {
144 +
      to_rm <- unique(which(is.na(pred) | is.nan(pred) | is.infinite(pred)))
145 +
      pred <- pred[-to_rm]
146 +
      observed <- observed[-to_rm]
147 +
    }
148 +
149 +
    # get ROCR performance object
150 +
    args <- list(...)
151 +
    args$measure <- measure
152 +
    args$prediction.obj <- ROCR::prediction(pred, observed)
153 +
    performance_object <- do.call(ROCR::performance, args)
154 +
155 +
    # get single performance value from the performance object
156 +
    performance <- unlist(performance_object@y.values)
157 +
    if (length(performance) != 1) {
158 +
      if (performance_object@x.name != "Cutoff") {
159 +
        stop(
160 +
          "Multiple performance values returned, but the measure is not",
161 +
          "cutoff-depedent. Check that the supplied measure is valid."
162 +
        )
163 +
      } else {
164 +
        # select the performance closest to the supplied cutoff
165 +
        performance <- performance[
166 +
          which.min(abs(unlist(performance_object@x.values) - cutoff))
167 +
        ]
168 +
      }
169 +
    }
170 +
    attributes(performance)$loss <- FALSE
171 +
    attributes(performance)$optimize <- "maximize"
172 +
    attributes(performance)$name <- ifelse(is.null(name), measure, name)
173 +
    return(performance)
174 +
  }
175 +
}
176 +
92 177
utils::globalVariables(c("id", "loss", "obs", "pred", "wts"))
93 178
#' Cross-validated Risk Estimation
94 179
#'
95 -
#' Estimates the cross-validated risk for a given learner and loss function.
180 +
#' Estimates the cross-validated risk for a given learner and evaluation
181 +
#' function, which can be either a loss or a risk function.
96 182
#'
97 183
#' @param learner A trained learner object.
98 -
#' @param loss_fun A valid loss function. See \code{\link{loss_functions}}.
184 +
#' @param eval_fun A valid loss or risk function. See
185 +
#' \code{\link{loss_functions}} and \code{\link{risk_functions}}.
99 186
#' @param coefs A \code{numeric} vector of coefficients.
100 187
#'
101 188
#' @importFrom assertthat assert_that
102 189
#' @importFrom data.table data.table ":=" set setnames setorderv
103 190
#' @importFrom origami cross_validate validation fold_index
104 191
#' @importFrom stats sd
105 -
cv_risk <- function(learner, loss_fun, coefs = NULL) {
106 -
  assertthat::assert_that("cv" %in% learner$properties,
192 +
cv_risk <- function(learner, eval_fun = NULL, coefs = NULL) {
193 +
  assertthat::assert_that(
194 +
    "cv" %in% learner$properties,
107 195
    msg = "learner is not cv-aware"
108 196
  )
109 197
@@ -115,7 +203,7 @@
Loading
115 203
  }
116 204
117 205
  get_obsdata <- function(fold, task) {
118 -
    list(loss_dt = data.table::data.table(
206 +
    list(dt = data.table::data.table(
119 207
      fold_index = origami::fold_index(),
120 208
      index = origami::validation(),
121 209
      obs = origami::validation(task$Y),
@@ -123,44 +211,67 @@
Loading
123 211
      wts = origami::validation(task$weights)
124 212
    ))
125 213
  }
126 -
127 -
  loss_dt <- origami::cross_validate(get_obsdata, task$folds, task)$loss_dt
128 -
  data.table::setorderv(loss_dt, c("index", "fold_index"))
129 -
  loss_dt <- cbind(loss_dt, preds)
130 -
  pred_cols <- colnames(preds)
131 -
  loss_long <- melt(loss_dt,
132 -
    measure.vars = pred_cols,
214 +
  dt <- origami::cross_validate(get_obsdata, task$folds, task)$dt
215 +
  data.table::setorderv(dt, c("index", "fold_index"))
216 +
  dt <- cbind(dt, preds)
217 +
  dt_long <- melt(
218 +
    dt,
219 +
    measure.vars = colnames(preds),
133 220
    variable.name = "learner",
134 221
    value.name = "pred"
135 222
  )
136 -
  loss_long[, `:=`(loss = wts * loss_fun(pred, obs))]
137 223
138 -
  # average loss in id-fold cluster
139 -
  loss_by_id <- loss_long[, list(loss = mean(loss, na.rm = TRUE)),
140 -
    by = list(learner, id, fold_index)
141 -
  ]
224 +
  eval_result <- eval_fun(dt_long[["pred"]], dt_long[["obs"]])
225 +
226 +
  if (!is.null(attr(eval_result, "loss")) && !attr(eval_result, "loss")) {
227 +
    # try stratifying by id
228 +
    eval_by_id <- tryCatch(
229 +
      {
230 +
        dt_long[, list(
231 +
          eval_result = eval_fun(pred, obs)
232 +
        ), by = list(learner, id, fold_index)]
233 +
      },
234 +
      error = function(c) {
235 +
        dt_long[, list(
236 +
          eval_result = eval_fun(pred, obs)
237 +
        ), by = list(learner, fold_index)]
238 +
      }
239 +
    )
240 +
  } else {
241 +
    dt_long <- cbind(dt_long, eval_result = dt_long[["wts"]] * eval_result)
242 +
243 +
    eval_by_id <- dt_long[, list(
244 +
      eval_result = mean(eval_result, na.rm = TRUE)
245 +
    ), by = list(learner, id, fold_index)]
246 +
  }
142 247
143 -
  # get learner level loss statistics
144 -
  loss_stats <- loss_by_id[, list(
248 +
  # get learner-level evaluation statistics
249 +
  eval_stats <- eval_by_id[, list(
145 250
    coefficients = NA_real_,
146 -
    risk = mean(loss, na.rm = TRUE),
147 -
    se = (1 / sqrt(.N)) * stats::sd(loss)
251 +
    risk = mean(eval_result, na.rm = TRUE),
252 +
    se = (1 / sqrt(.N)) * stats::sd(eval_result)
148 253
  ), by = list(learner)]
149 254
150 -
  # get fold-learner level loss statistics
151 -
  loss_fold_stats <- loss_by_id[, list(risk = mean(loss, na.rm = TRUE)),
152 -
    by = list(learner, fold_index)
153 -
  ]
154 -
  loss_stats_fold <- loss_fold_stats[, list(
255 +
  # get learner-level evaluation statistics by fold
256 +
  eval_fold_stats <- eval_by_id[, list(
257 +
    risk = mean(eval_result, na.rm = TRUE)
258 +
  ), by = list(learner, fold_index)]
259 +
  eval_stats_fold <- eval_fold_stats[, list(
155 260
    fold_sd = stats::sd(risk, na.rm = TRUE),
156 261
    fold_min_risk = min(risk, na.rm = TRUE),
157 262
    fold_max_risk = max(risk, na.rm = TRUE)
158 263
  ), by = list(learner)]
159 264
160 -
  risk_dt <- loss_stats <- merge(loss_stats, loss_stats_fold, by = "learner")
265 +
  risk_dt <- eval_stats <- merge(eval_stats, eval_stats_fold, by = "learner")
161 266
162 267
  if (!is.null(coefs)) {
163 268
    data.table::set(risk_dt, , "coefficients", coefs)
164 269
  }
270 +
271 +
  if (!is.null(attr(eval_result, "name"))) {
272 +
    colnames(risk_dt) <- gsub(
273 +
      "risk", attr(eval_result, "name"), colnames(risk_dt)
274 +
    )
275 +
  }
165 276
  return(risk_dt)
166 277
}

@@ -1,166 +1,305 @@
Loading
1 +
utils::globalVariables(c("score"))
1 2
#' Variable Importance
2 3
#'
3 4
#' Function that takes a cross-validated fit (i.e., cross-validated learner that
4 5
#' has already been trained on a task), which could be a cross-validated single
5 -
#' learner or super learner, to generate a loss-based variable importance
6 -
#' measure for each predictor, where the predictors are the covariates in the
7 -
#' trained task. This function generates a \code{data.table} in which each row
8 -
#' corresponds to the risk difference or risk ratio between the following
9 -
#' two risks: the risk when a predictor is permuted or removed, and the original
10 -
#' risk (i.e., when all predictors are included). A higher risk ratio/difference
11 -
#' corresponds to a more important predictor. A plot can be generated from the
12 -
#' returned \code{data.table} by calling companion function
13 -
#' \code{plot_importance}.
14 -
#'
15 -
#' @param fit A trained cross-validated learner (e.g., cv stack, super learner),
16 -
#'  from which cross-validated predictions can be generated.
17 -
#' @param loss The loss function for evaluating the risk. Defaults according to
18 -
#'  outcome type: squared error loss for continuous outcomes, and negative
19 -
#'  log-likelihood loss for discrete outcomes. See \code{\link{loss_functions}}.
20 -
#' @param fold_number The fold number to use for obtaining the predictions
21 -
#'  from the fit. Either a positive integer for obtaining predictions from a
22 -
#'  specific fold's fit; \code{"full"} for obtaining predictions from a fit on
23 -
#'  all of the data, or \code{"validation"} (default) for obtaining
24 -
#'  cross-validated predictions, where the data used for training and prediction
25 -
#'  never overlaps across the folds. Note that if a positive integer or
26 -
#'  \code{"full"} is supplied, then there will be overlap between the  data
27 -
#'  used for training and prediction.
28 -
#' @param type Which method should be used to obscure the relationship between
29 -
#'  the covariate and the outcome. When type is \code{"remove"} (default), each
30 -
#'  covariate is removed from the task and the cross-validated learner is refit
31 -
#'  to this modified task and then predictions are obtained from this refit.
32 -
#'  When type is \code{"permute"}, each covariate is permuted (sampled without
33 -
#'  replacement), and then predictions are obtained from this modified data with
34 -
#'  one permuted covariate.
35 -
#' @param importance_metric Either \code{"ratio"} (default), which returns
36 -
#'  the risk with the permuted/removed X divided by observed risk with all X; or
37 -
#'  \code{"difference"}, which returns the difference between the risk with the
38 -
#'  permuted/removed X and the observed risk.
6 +
#' learner or super learner, and generates a risk-based variable importance
7 +
#' score for either each covariate or each group of covariates in the task.
8 +
#' This function outputs a \code{data.table}, where each row corresponds to the
9 +
#' risk difference or the risk ratio between the following two risks: the risk
10 +
#' when a covariate (or group of covariates) is permuted or removed, and the
11 +
#' original risk (i.e., when all covariates are included as they were in the
12 +
#' observed data). A higher risk ratio/difference corresponds to a more
13 +
#' important covariate/group. A plot can be generated from the returned
14 +
#' \code{data.table} by calling companion function \code{\link{importance_plot}}.
39 15
#'
40 -
#' @importFrom data.table data.table
41 -
#'
42 -
#' @return A \code{data.table} of variable importance for each covariate.
16 +
#' @export
43 17
#'
44 18
#' @name importance
19 +
#' @rdname importance
45 20
#' @keywords variable importance
46 21
#'
47 -
#' @export
48 -
importance <- function(fit, loss = NULL, fold_number = "validation",
22 +
#' @importFrom data.table data.table setorder setnames
23 +
#'
24 +
#' @return A \code{data.table} of variable importance for each covariate.
25 +
#'
26 +
#' @section Parameters:
27 +
#'   - \code{fit}: A trained cross-validated (CV) learner (such as a CV stack
28 +
#'       or super learner), from which cross-validated predictions can be
29 +
#'       generated.
30 +
#'   - \code{eval_fun = NULL}: The evaluation function (risk or loss function)
31 +
#'       for evaluating the risk. Defaults vary based on the outcome type,
32 +
#'       matching defaults in \code{\link{default_metalearner}}. See
33 +
#'       \code{\link{loss_functions}} and \code{\link{risk_functions}} for
34 +
#'       options.
35 +
#'   - \code{fold_number}: The fold number to use for obtaining the predictions
36 +
#'       from the fit. Either a positive integer for obtaining predictions from
37 +
#'       a specific fold's fit; \code{"full"} for obtaining predictions from a
38 +
#'       fit on all of the data, or \code{"validation"} (default) for obtaining
39 +
#'       cross-validated predictions, where the data used for training and
40 +
#'       prediction never overlaps across the folds. Note that if a positive
41 +
#'       integer or \code{"full"} is supplied here then there will be overlap
42 +
#'       between the data used for training and validation, so
43 +
#'       \code{fold_number ="validation"} is recommended.
44 +
#'   - \code{type}: Which method should be used to obscure the relationship
45 +
#'       between each covariate / covariate group and the outcome? When
46 +
#'       \code{type} is \code{"remove"} (default), each covariate / covariate
47 +
#'       group is removed one at a time from the task; the cross-validated
48 +
#'       learner is refit to this modified task; and finally, predictions are
49 +
#'       obtained from this refit. When \code{type} is \code{"permute"}, each
50 +
#'       covariate / covariate group is permuted (sampled without replacement)
51 +
#'       one at a time, and then predictions are obtained from this modified
52 +
#'       data.
53 +
#'   - \code{importance_metric}: Either \code{"ratio"} or \code{"difference"}
54 +
#'       (default). For each covariate / covariate group, \code{"ratio"}
55 +
#'       returns the risk of the permuted/removed covariate / covariate group
56 +
#'       divided by observed/original risk (i.e., the risk with all covariates
57 +
#'       as they existed in the sample) and \code{"difference"} returns the
58 +
#'       difference between the risk with the permuted/removed covariate /
59 +
#'       covariate group and the observed risk.
60 +
#'   - \code{covariate_groups}: Optional named list covariate groups which will
61 +
#'       invoke variable importance evaluation at the group-level, by
62 +
#'       removing/permuting all covariates in the same group together. If
63 +
#'       covariates in the task are not specified in the list of groups, then
64 +
#'       those covariates will be added as additional single-covariate groups.
65 +
#'
66 +
#' @examples
67 +
#' # define ML task
68 +
#' data(cpp_imputed)
69 +
#' covs <- c("apgar1", "apgar5", "parity", "gagebrth", "mage", "meducyrs")
70 +
#' task <- sl3_Task$new(cpp_imputed, covariates = covs, outcome = "haz")
71 +
#'
72 +
#' # build relatively fast learner library (not recommended for real analysis)
73 +
#' lasso_lrnr <- Lrnr_glmnet$new()
74 +
#' glm_lrnr <- Lrnr_glm$new()
75 +
#' ranger_lrnr <- Lrnr_ranger$new()
76 +
#' lrnrs <- c(lasso_lrnr, glm_lrnr, ranger_lrnr)
77 +
#' names(lrnrs) <- c("lasso", "glm", "ranger")
78 +
#' lrnr_stack <- make_learner(Stack, lrnrs)
79 +
#'
80 +
#' # instantiate SL with default metalearner
81 +
#' sl <- Lrnr_sl$new(lrnr_stack)
82 +
#' sl_fit <- sl$train(task)
83 +
#'
84 +
#' importance_result <- importance(sl_fit)
85 +
#' importance_result
86 +
#'
87 +
#' # importance with groups of covariates
88 +
#' groups <- list(
89 +
#'   scores = c("apgar1", "apgar5"),
90 +
#'   maternal = c("parity", "mage", "meducyrs")
91 +
#' )
92 +
#' importance_result_groups <- importance(sl_fit, covariate_groups = groups)
93 +
#' importance_result_groups
94 +
importance <- function(fit, eval_fun = NULL,
95 +
                       fold_number = "validation",
49 96
                       type = c("remove", "permute"),
50 -
                       importance_metric = c("ratio", "difference")) {
97 +
                       importance_metric = c("difference", "ratio"),
98 +
                       covariate_groups = NULL) {
51 99
52 -
  # check arguments
100 +
  # check fit is trained
53 101
  if (!fit$is_trained) {
54 102
    stop("Fit is not trained.")
55 103
  }
56 104
57 -
  # set defaults
105 +
  ################################ set defaults ################################
58 106
  importance_metric <- match.arg(importance_metric)
59 107
  type <- match.arg(type)
60 108
61 -
  # extract task and data
62 -
  task <- fit$training_task
63 -
  d <- task$data
64 -
  X <- task$nodes$covariates
65 -
  Y <- task$Y
66 -
67 -
  if (is.null(loss)) {
68 -
    outcome_type <- task$outcome_type$type
109 +
  if (is.null(eval_fun)) {
110 +
    outcome_type <- fit$training_task$outcome_type$type
69 111
    if (outcome_type %in% c("constant", "binomial")) {
70 -
      loss <- loss_loglik_binomial
112 +
      eval_fun <- loss_squared_error
71 113
    } else if (outcome_type == "categorical") {
72 -
      loss <- loss_loglik_multinomial
114 +
      eval_fun <- loss_loglik_multinomial
73 115
    } else if (outcome_type == "continuous") {
74 -
      loss <- loss_squared_error
116 +
      eval_fun <- loss_squared_error
75 117
    } else if (outcome_type == "multivariate") {
76 -
      loss <- loss_squared_error_multivariate
118 +
      eval_fun <- loss_squared_error_multivariate
77 119
    } else {
78 120
      stop(paste0(
79 -
        "No default loss for outcome type ", outcome_type,
121 +
        "No default eval_fun for outcome type ", outcome_type,
80 122
        ". Please specify your own."
81 123
      ))
82 124
    }
83 125
  }
84 126
85 -
  # get predictions and risk
86 -
  pred <- fit$predict_fold(task, fold_number = fold_number)
87 -
  original_risk <- mean(loss(pred, Y))
127 +
  ########################### extract nodes and data ###########################
128 +
  task <- fit$training_task
129 +
  d <- task$data
130 +
  X <- task$nodes$covariates
131 +
  Y <- task$Y
132 +
  weights <- task$weights
88 133
89 -
  # X-length list of importance scores
134 +
  ############################ check covariate groups ##########################
135 +
  if (!is.null(covariate_groups)) {
136 +
    if (!is.list(covariate_groups)) {
137 +
      stop("Covariate groups must be a list.")
138 +
    }
139 +
140 +
    # check that all covariates in the groups are also in the task's covariates
141 +
    if (!all(unlist(covariate_groups) %in% X)) {
142 +
      stop("Groups contain covariates that are not in the task's covariates.")
143 +
    }
144 +
145 +
    # check if all covariates in task are in groups, and if not add them
146 +
    if (!all(X %in% unlist(covariate_groups))) {
147 +
      missingX <- as.list(X[which(!X %in% unlist(covariate_groups))])
148 +
      names(missingX) <- X[which(!X %in% unlist(covariate_groups))]
149 +
      covariate_groups <- c(covariate_groups, missingX)
150 +
    }
151 +
    # check that groups with more than one covariate are named, and add name to
152 +
    # unnamed groups but only if unnamed group is a single covariate
153 +
    if (any(is.null(names(covariate_groups))) |
154 +
      any(names(covariate_groups) == "") |
155 +
      any(is.na(names(covariate_groups)))) {
156 +
      no_name <- unique(which(
157 +
        is.null(names(covariate_groups)) | names(covariate_groups) == "" |
158 +
          is.na(names(covariate_groups))
159 +
      ))
160 +
      if (any(sapply(covariate_groups[unique(no_name)], length)) != 1) {
161 +
        stop("Covariate groups with more than one covariate must be named.")
162 +
      } else if (all(sapply(covariate_groups[unique(no_name)], length)) == 1) {
163 +
        names(covariate_groups[no_name]) <- unlist(covariate_groups[no_name])
164 +
      }
165 +
    }
166 +
    X <- covariate_groups # for convenience in remaining parts of function
167 +
  } else {
168 +
    names(X) <- X
169 +
  }
170 +
171 +
  ########################## risk under observed data ##########################
172 +
  original_pred <- fit$predict_fold(task, fold_number = fold_number)
173 +
  original_eval <- eval_fun(original_pred, Y)
174 +
  if (!is.null(attr(original_eval, "loss")) && !attr(original_eval, "loss")) {
175 +
    original_risk <- original_eval
176 +
  } else {
177 +
    original_losses <- original_eval
178 +
    original_risk <- weighted.mean(original_losses, weights)
179 +
  }
180 +
181 +
  ######################## list of importance results ##########################
90 182
  res_list <- lapply(X, function(x) {
183 +
    # get the risk when the covariate/group is permuted/removed
91 184
    if (type == "permute") {
92 -
      # permute covariate x, and give it the same name as the original x
93 -
      x_permuted <- data.table(sample(unlist(d[, x, with = FALSE]), nrow(d)))
94 -
      names(x_permuted) <- x
185 +
      # get the permutation rows
186 +
      perm <- sample(1:nrow(d), nrow(d))
187 +
      # permute x (the covariate/group), and name it as the original x
188 +
      x_perm <- d[perm, x, with = FALSE]
189 +
      data.table::setnames(x_perm, x)
95 190
      # replace original x with permuted x, and update task with permuted x
96 -
      x_permuted_name <- task$add_columns(x_permuted)
97 -
      task_x_permuted <- task$next_in_chain(column_names = x_permuted_name)
191 +
      x_perm_name <- task$add_columns(x_perm)
192 +
      task_x_perm <- task$next_in_chain(column_names = x_perm_name)
98 193
      # obtain predictions & risk on the new task with permuted x
99 -
      x_permuted_pred <- fit$predict_fold(task_x_permuted, fold_number)
100 -
      no_x_risk <- mean(loss(x_permuted_pred, Y))
194 +
      x_perm_pred <- fit$predict_fold(task_x_perm, fold_number = fold_number)
195 +
      x_perm_eval <- eval_fun(x_perm_pred, Y)
196 +
      if (!is.null(attr(original_eval, "loss")) && !attr(original_eval, "loss")) {
197 +
        no_x_risk <- x_perm_eval
198 +
      } else {
199 +
        no_x_losses <- x_perm_eval
200 +
        no_x_risk <- weighted.mean(no_x_losses, weights)
201 +
      }
101 202
    } else if (type == "remove") {
102 -
      # modify learner to not include covariate x
103 -
      x_removed_lrnr <- fit$reparameterize(list(covariates = setdiff(X, x)))
104 -
      x_removed_fit <- x_removed_lrnr$train(task)
105 -
      x_removed_pred <- x_removed_fit$predict_fold(task, fold_number)
106 -
      no_x_risk <- mean(loss(x_removed_pred, Y))
203 +
      # modify learner to not include x (covariate/group)
204 +
      if (!is.null(covariate_groups)) {
205 +
        x_rm_covars <- setdiff(unlist(X), x)
206 +
      } else {
207 +
        x_rm_covars <- setdiff(X, x)
208 +
      }
209 +
      x_rm_lrnr <- fit$reparameterize(list(covariates = x_rm_covars))
210 +
      x_rm_fit <- x_rm_lrnr$train(task)
211 +
      x_rm_pred <- x_rm_fit$predict_fold(task, fold_number = fold_number)
212 +
      x_rm_eval <- eval_fun(x_rm_pred, Y)
213 +
      if (!is.null(attr(original_eval, "loss")) && !attr(original_eval, "loss")) {
214 +
        no_x_risk <- x_rm_eval
215 +
      } else {
216 +
        no_x_losses <- x_rm_eval
217 +
        no_x_risk <- weighted.mean(no_x_losses, weights)
218 +
      }
107 219
    }
220 +
108 221
    # evaluate importance
109 222
    if (importance_metric == "ratio") {
110 -
      return(no_x_risk / original_risk)
223 +
      result <- no_x_risk / original_risk
111 224
    } else if (importance_metric == "difference") {
112 -
      return(no_x_risk - original_risk)
225 +
      result <- no_x_risk - original_risk
113 226
    }
227 +
    return(result)
114 228
  })
115 -
  names(res_list) <- X
116 229
230 +
  ############################## prep output ###################################
117 231
  # importance results ordered by decreasing importance
118 -
  if (importance_metric == "ratio") {
119 -
    res <- data.table(X = names(res_list), risk_ratio = unlist(res_list))
120 -
    return(res[order(-res$risk_ratio)])
121 -
  } else if (importance_metric == "difference") {
122 -
    res <- data.table(X = names(res_list), risk_difference = unlist(res_list))
123 -
    return(res[order(-res$risk_difference)])
232 +
  result <- data.table::data.table(covariate = names(X), metric = unlist(res_list))
233 +
  if (!is.null(covariate_groups)) {
234 +
    colnames(result)[1] <- "covariate_group"
124 235
  }
236 +
  data.table::setorder(result, -metric)
237 +
238 +
  # name the importance metric appropriately
239 +
  metric_name <- paste0("risk_", importance_metric)
240 +
  if (!is.null(attr(original_eval, "name"))) {
241 +
    metric_name <- gsub("risk", attr(original_eval, "name"), metric_name)
242 +
  }
243 +
  colnames(result)[2] <- metric_name
244 +
  return(result)
125 245
}
126 246
127 247
#' Variable Importance Plot
128 248
#'
129 -
#' @param x The 2-column \code{data.table} returned by the call to
130 -
#'  \code{importance}, where the first column is the predictors and the second
131 -
#'  col
132 -
#' @param nvar The maximum number of predictors to be plotted. Defaults to 30.
133 -
#' @param ... Other parameters passed to \code{\link[graphics]{dotchart}}.
249 +
#' @section Parameters:
250 +
#'   - \code{x}: The 2-column \code{data.table} returned by
251 +
#'       \code{\link{importance}}, where the first column is the
252 +
#'       covariate/groups and the second column is the importance score.
253 +
#'   - \code{nvar}: The maximum number of predictors to be plotted. Defaults to
254 +
#'       the minimum between 30 and the number of rows in \code{x}.
134 255
#'
135 -
#' @return A \code{\link[graphics]{dotchart}} of variable importance.
256 +
#' @return A \code{\link[ggplot2]{ggplot}} of variable importance.
136 257
#'
137 -
#' @importFrom graphics dotchart
258 +
#' @importFrom ggplot2 ggplot geom_point coord_flip scale_x_discrete labs
259 +
#' @importFrom data.table data.table
138 260
#'
261 +
#' @rdname importance_plot
139 262
#' @name importance_plot
140 263
#' @keywords variable importance
141 264
#'
142 265
#' @export
143 -
importance_plot <- function(x, nvar = min(30, nrow(x)), ...) {
266 +
#'
267 +
#' @examples
268 +
#' # define ML task
269 +
#' data(cpp_imputed)
270 +
#' covs <- c("apgar1", "apgar5", "parity", "gagebrth", "mage", "meducyrs")
271 +
#' task <- sl3_Task$new(cpp_imputed, covariates = covs, outcome = "haz")
272 +
#'
273 +
#' # build relatively fast learner library (not recommended for real analysis)
274 +
#' lasso_lrnr <- Lrnr_glmnet$new()
275 +
#' glm_lrnr <- Lrnr_glm$new()
276 +
#' ranger_lrnr <- Lrnr_ranger$new()
277 +
#' lrnrs <- c(lasso_lrnr, glm_lrnr, ranger_lrnr)
278 +
#' names(lrnrs) <- c("lasso", "glm", "ranger")
279 +
#' lrnr_stack <- make_learner(Stack, lrnrs)
280 +
#'
281 +
#' # instantiate SL with default metalearner
282 +
#' sl <- Lrnr_sl$new(lrnr_stack)
283 +
#' sl_fit <- sl$train(task)
284 +
#' importance_result <- importance(sl_fit)
285 +
#' importance_plot(importance_result)
286 +
importance_plot <- function(x, nvar = min(30, nrow(x))) {
144 287
145 288
  # get the importance metric
146 -
  if (grepl("ratio", colnames(x)[2])) {
147 -
    xlab <- "Risk Ratio"
148 -
  } else if (grepl("difference", colnames(x)[2])) {
149 -
    xlab <- "Risk Difference"
150 -
  }
289 +
  xlab <- colnames(x)[2]
151 290
152 -
  # sort by decreasing importance
153 -
  x_sorted <- x[order(-x[, 2])]
154 -
  # subset to include most at most nvar
155 -
  x_sorted <- x_sorted[1:(min(nvar, nrow(x_sorted))), ]
156 291
  # sort by increasing importance
157 -
  x_sorted <- x_sorted[order(x_sorted[, 2])]
158 -
159 -
  # make dotchart with most important variables on top
160 -
  # x_sorted[[2]] is importance scores & x_sorted[[1]] is covariate names
161 -
  dotchart(
162 -
    x = x_sorted[[2]], labels = x_sorted[[1]],
163 -
    xlab = xlab, ylab = "",
164 -
    xlim = c(min(x_sorted[[2]]), max(x_sorted[[2]])), ...
165 -
  )
292 +
  x <- x[order(-x[, 2]), ]
293 +
  # subset to include most at most nvar
294 +
  x <- x[1:(min(nvar, nrow(x))), ]
295 +
296 +
  # format for ggplot
297 +
  d <- data.table::data.table(vars = factor(x[[1]], levels = x[[1]]), score = x[[2]])
298 +
299 +
  ggplot2::ggplot(d, aes(x = vars, y = score)) +
300 +
    ggplot2::geom_point() +
301 +
    ggplot2::ylim(c(min(d$score), max(d$score))) +
302 +
    ggplot2::labs(x = "", y = xlab, title = "sl3 Variable Importance Plot") +
303 +
    ggplot2::scale_x_discrete(limits = rev(levels(d$covs))) +
304 +
    ggplot2::coord_flip()
166 305
}

@@ -474,24 +474,31 @@
Loading
474 474
    folds = function(new_folds) {
475 475
      if (!missing(new_folds)) {
476 476
        private$.folds <- new_folds
477 -
      } else if (is.numeric(private$.folds)) {
478 -
        # if an integer, create new_folds object but pass integer to V argument
479 -
        if (self$has_node("id")) {
480 -
          new_folds <- origami::make_folds(
481 -
            cluster_ids = self$id,
482 -
            V = private$.folds
483 -
          )
484 -
        } else {
485 -
          new_folds <- origami::make_folds(n = self$nrow, V = private$.folds)
477 +
      } else if (is.numeric(private$.folds) | is.null(private$.folds)) {
478 +
        args <- list()
479 +
        args$n <- self$nrow
480 +
        if (is.numeric(private$.folds)) {
481 +
          # pass integer to V argument when it's supplied
482 +
          args$V <- private$.folds
486 483
        }
487 -
        private$.folds <- new_folds
488 -
      } else if (is.null(private$.folds)) {
489 -
        # generate folds now if never specified
490 484
        if (self$has_node("id")) {
491 -
          new_folds <- origami::make_folds(cluster_ids = self$id)
492 -
        } else {
493 -
          new_folds <- origami::make_folds(n = self$nrow)
485 +
          # clustered cross-validation for clustered data
486 +
          args$cluster_ids <- self$id
487 +
        }
488 +
        if (self$outcome_type$type %in% c("binomial", "categorical")) {
489 +
          # stratified cross-validation folds for discrete outcomes
490 +
          args$strata_ids <- self$Y
491 +
          if (self$has_node("id")) {
492 +
            # don't use stratified CV if clusters are not nested in strata
493 +
            is_nested <- all(
494 +
              rowSums(table(args$cluster_ids, args$strata_ids) > 0) == 1
495 +
            )
496 +
            if (!is_nested) {
497 +
              args <- args[!(names(args) == "strata_ids")]
498 +
            }
499 +
          }
494 500
        }
501 +
        new_folds <- do.call(origami::make_folds, args)
495 502
        private$.folds <- new_folds
496 503
      }
497 504
      return(private$.folds)

@@ -21,10 +21,10 @@
Loading
21 21
#' @family Learners
22 22
#'
23 23
#' @section Parameters:
24 -
#'   - \code{loss_function}: A function that takes a vector of predictions as
24 +
#'   - \code{eval_function}: A function that takes a vector of predictions as
25 25
#'     it's first argument, and a vector of truths/observations as it's second
26 -
#'     argument, and then returns a vector of losses. See \link{loss_functions}
27 -
#'     for options.
26 +
#'     argument, and then returns a vector of losses or a numeric risk. See
27 +
#'     \link{loss_functions} and \link{risk_functions} for options.
28 28
#'
29 29
#' @examples
30 30
#' data(cpp_imputed)
@@ -50,7 +50,7 @@
Loading
50 50
  classname = "Lrnr_cv_selector",
51 51
  inherit = Lrnr_base, portable = TRUE, class = TRUE,
52 52
  public = list(
53 -
    initialize = function(loss_function = loss_squared_error) {
53 +
    initialize = function(eval_function = loss_squared_error) {
54 54
      params <- args_to_list()
55 55
      super$initialize(params = params)
56 56
    }
@@ -60,7 +60,7 @@
Loading
60 60
      "continuous", "binomial", "categorical", "weights", "wrapper"
61 61
    ),
62 62
    .train = function(task) {
63 -
      loss_function <- self$params$loss_function
63 +
      eval_function <- self$params$eval_function
64 64
65 65
      # specify data
66 66
      outcome_type <- self$get_outcome_type(task)
@@ -70,8 +70,13 @@
Loading
70 70
      weights <- task$weights
71 71
72 72
      risk <- function(preds) {
73 -
        loss <- loss_function(preds, Y)
74 -
        risk <- weighted.mean(loss, weights)
73 +
        eval_result <- eval_function(preds, Y)
74 +
        if (!is.null(attr(eval_result, "loss")) && !attr(eval_result, "loss")) {
75 +
          risk <- eval_result
76 +
        } else {
77 +
          loss <- eval_result
78 +
          risk <- weighted.mean(loss, weights)
79 +
        }
75 80
        return(risk)
76 81
      }
77 82
      risks <- apply(X, 2, risk)

@@ -73,15 +73,15 @@
Loading
73 73
          if (is.null(private$.cv_risk)) {
74 74
            tryCatch(
75 75
              {
76 -
                # try using loss function based on outcome type
77 -
                loss_fun <- private$.params$metalearner$params$loss_function
78 -
                private$.cv_risk <- self$cv_risk(loss_fun)
76 +
                # try using eval function based on outcome type
77 +
                eval_fun <- private$.params$metalearner$params$eval_function
78 +
                private$.cv_risk <- self$cv_risk(eval_fun)
79 79
              },
80 80
              error = function(c) {
81 81
                # check training outcome type explicitly
82 82
                metalearner <- default_metalearner(self$training_outcome_type)
83 -
                loss_fun <- metalearner$params$loss_function
84 -
                private$.cv_risk <- self$cv_risk(loss_fun)
83 +
                eval_fun <- metalearner$params$eval_function
84 +
                private$.cv_risk <- self$cv_risk(eval_fun)
85 85
              }
86 86
            )
87 87
          }
@@ -97,10 +97,10 @@
Loading
97 97
      self$assert_trained()
98 98
      return(private$.fit_object$cv_meta_fit$fit_object)
99 99
    },
100 -
    cv_risk = function(loss_fun) {
100 +
    cv_risk = function(eval_fun) {
101 101
      # get risks for cv learners (nested cv)
102 102
      cv_stack_fit <- self$fit_object$cv_fit
103 -
      stack_risks <- cv_stack_fit$cv_risk(loss_fun)
103 +
      stack_risks <- cv_stack_fit$cv_risk(eval_fun)
104 104
105 105
      coefs <- self$coefficients
106 106
      if (!is.null(coefs)) {
@@ -119,7 +119,7 @@
Loading
119 119
      )
120 120
121 121
      # get risks for super learner ("revere" CV)
122 -
      sl_risk <- cv_risk(self, loss_fun)
122 +
      sl_risk <- cv_risk(self, eval_fun)
123 123
      set(sl_risk, , "learner", "SuperLearner")
124 124
125 125
      # combine and return

@@ -1,35 +1,21 @@
Loading
1 -
#' Estimate Cross-Validated Risk of Super Learner
2 -
#' @param lrnr_sl a Lrnr_sl object specifying the Super Learner
3 -
#' @param task the task used for training and performance assessment
4 -
#' @param loss_fun the loss function used to evaluate Super Learner
1 +
#' Estimates cross-validated risk of the Super Learner
2 +
#' @param lrnr_sl a \code{Lrnr_sl} object specifying the Super Learner.
3 +
#' @param task the task used for training and performance assessment.
4 +
#' @param eval_fun the evaluation function, either a loss or risk function, for
5 +
#'  evaluating the Super Learner's predictions.
5 6
#' @export
6 -
CV_lrnr_sl <- function(lrnr_sl, task, loss_fun) {
7 +
CV_lrnr_sl <- function(lrnr_sl, task, eval_fun) {
8 +
  # check arguments
7 9
  if (!inherits(lrnr_sl, "Lrnr_sl")) {
8 10
    stop("lrnr_sl must be a Lrnr_sl object")
9 11
  }
12 +
  # cross-validate the SL
10 13
  cv_sl <- make_learner(Lrnr_cv, lrnr_sl, full_fit = TRUE)
11 14
  cv_sl_fit <- cv_sl$train(task)
12 -
  #
13 -
  # # to avoid refitting the stack to the full data,
14 -
  # # extract the full data stack fits and combine
15 -
  # fold_fits <- cv_sl_fit$fit_object$fold_fits
16 -
  #
17 -
  # combined_fold_fits <- lapply(fold_fits, function(fold_fit){
18 -
  #   full_fit <- fold_fit$fit_object$full_fit
19 -
  #   full_stack <- fold_fit$fit_object$cv_fit$fit_object$full_fit
20 -
  #   stack_combined <- make_learner(Stack, full_stack, full_fit)
21 -
  # })
22 -
  # combined_fit_object <- cv_sl_fit$fit_object
23 -
  # combined_fit_object$fold_fits <- combined_fold_fits
24 -
  # cv_combined_fit <- copy(cv_sl)
25 -
  # cv_combined_fit$set_train(combined_fit_object, task)
26 -
  #
27 15
  full_fit <- cv_sl_fit$fit_object$full_fit
28 -
29 16
  # TODO: extract loss function from lrnr_sl where possible
30 -
  full_risk <- full_fit$cv_risk(loss_fun)
31 -
  sl_risk <- cv_sl_fit$cv_risk(loss_fun)
32 -
17 +
  full_risk <- full_fit$cv_risk(eval_fun)
18 +
  sl_risk <- cv_sl_fit$cv_risk(eval_fun)
33 19
  # replace revere cv sl risk with nested cv sl risk
34 20
  stack_risks <- full_risk[full_risk$learner != "SuperLearner"]
35 21
  set(sl_risk, , "learner", "SuperLearner")
Files Coverage
R 79.77%
Project Totals (91 files) 79.77%
1
codecov:
2
  token: 0e172078-76a1-4b94-bd7f-2e10d0647b3d
3

4
comment: false
5

6
ignore:
7
  - "R/utils.R"
8

9
coverage:
10
  status:
11
    project:
12
      default:
13
        against: parent
14
        target: auto
15
        threshold: 1%
16
        branches:
17
          - master
18
        if_no_uploads: error
19
        if_not_found: success
20
        if_ci_failed: error
21
        only_pulls: false
22
        flags:
23
          - integration
24
        paths:
25
          - folder
26

27
    patch:
28
      default:
29
        against: parent
30
        target: 80%
31
        branches: null
32
        if_no_uploads: success
33
        if_not_found: success
34
        if_ci_failed: error
35
        only_pulls: false
36
        flags:
37
          - integration
38
        paths:
39
          - folder
40

41
    changes:
42
      default:
43
        against: parent
44
        branches: null
45
        if_no_uploads: error
46
        if_not_found: success
47
        if_ci_failed: error
48
        only_pulls: false
49
        flags:
50
          - integration
51
        paths:
52
          - folder
53

54
  flags:
55
    integration:
56
      branches:
57
        - master
58
      ignore:
59
        - app/ui
60

61
  ignore: # files and folders for processing
62
    - tests/*
63

64
  fixes:
65
    - "old_path::new_path"
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading