First commit of isotonic constraints

NAMESPACE

@@ -21,6 +21,7 @@ export(lmtp_tmle)
 export(static_binary_off)
 export(static_binary_on)
 export(tidy)
+importFrom(OrdMonReg,BoundedIsoMean)
 importFrom(R6,R6Class)
 importFrom(checkmate,assert_character)
 importFrom(checkmate,assert_function)

R/curve.R

@@ -77,10 +77,17 @@ estimate_curve_sdr <- function(task, fold, ratios, learners, control, pb) {
       fits[[t]] <- extract_sl_weights(fit)
     }
 
-    mnt <- update_m(mnt, fit, natural$train, t, task$tau)
-    mst <- update_m(mst, fit, ust, t, task$tau)
-    mnv <- update_m(mnv, fit, natural$valid, t, task$tau)
-    msv <- update_m(msv, fit, usv, t, task$tau)
+    lambda <- function(x) x
+    if(control$.isotonic_constraint == TRUE && task$outcome_type == "binomial") {
+      pred <- predict(fit, natural$train[at_risk & observed & time, vars])
+      out <- natural$train[at_risk & observed & time, "..i..Y_1"]
+      lambda <- isotonic_constraint(pred, out)
+    }
+
+    mnt <- update_m(mnt, fit, natural$train, t, task$tau, lambda)
+    mst <- update_m(mst, fit, ust, t, task$tau, lambda)
+    mnv <- update_m(mnv, fit, natural$valid, t, task$tau, lambda)
+    msv <- update_m(msv, fit, usv, t, task$tau, lambda)
 
     if ((t + 1) <= task$tau) {
       psuedo <- unlist(lapply((t + 1):task$tau, function(x) {
@@ -113,14 +120,14 @@ predict_long <- function(fit, newdata, t, tau) {
   ans[, 1]
 }
 
-update_m <- function(m, fit, newdata, t, tau) {
+update_m <- function(m, fit, newdata, t, tau, lambda) {
   pred <- predict_long(fit, newdata, t, tau)
   time <- as.numeric(newdata$time) >= t
 
   for (l in seq_along(t:tau)) {
     x <- (t:tau)[l]
     j <- as.numeric(newdata$time[time]) == x
-    m[[x]][, x - t + 1] <- pred[j]
+    m[[x]][, x - t + 1] <- lambda(pred[j])
   }
   m
 }

R/isotonic_constraint.R

@@ -0,0 +1,21 @@
+#' @importFrom OrdMonReg BoundedIsoMean
+isotonic_constraint <- function(x, y) {
+  na <- is.na(x) | is.na(y)
+  x <- x[!na]
+  y <- y[!na]
+
+  o <- order(x)
+  x <- x[o]
+  y <- y[o]
+
+  a <- rep(max(0, min(x)), length(x))
+  b <- rep(max(0, min(1, max(x))), length(x))
+
+  w <- rep(1, length(x))
+
+  iso <- BoundedIsoMean(y, w, a, b)
+
+  change_points <- x[iso[2:length(iso)] != lag(iso)[2:length(iso)]]
+
+  approxfun(x, iso, method = "constant", rule = 2)
+}

R/lmtp_control.R

@@ -13,6 +13,8 @@
 #'  The number of cross-validation folds for \code{learners_trt}.
 #' @param .return_full_fits \[\code{logical(1)}\]\cr
 #'  Return full 'mlr3superlearner' fits? Default is \code{FALSE}, return only 'mlr3superlearner' weights.
+#' @param .iosotonic_constraint \[\code{logical(1)}]\cr
+#'  Constrain pseudo-outcomes using isotonic regression in curve algorithm
 #' @param .discrete \[\code{logical(1)}\]\cr
 #'  Use discrete or ensemble super learner?
 #' @param .info \[\code{logical(1)}\]\cr
@@ -28,6 +30,7 @@ lmtp_control <- function(.bound = 1e5,
                          .learners_outcome_folds = 10,
                          .learners_trt_folds = 10,
                          .return_full_fits = FALSE,
+                         .isotonic_constraint = TRUE,
                          .discrete = TRUE,
                          .info = FALSE) {
 
@@ -36,6 +39,7 @@ lmtp_control <- function(.bound = 1e5,
   assert_number(.bound)
   assert_number(.trim, upper = 1)
   assert_logical(.return_full_fits, len = 1)
+  assert_logical(.isotonic_constraint, len = 1)
   assert_logical(.discrete, len = 1)
   assert_logical(.info, len = 1)
 
@@ -44,6 +48,7 @@ lmtp_control <- function(.bound = 1e5,
        .learners_outcome_folds = .learners_outcome_folds,
        .learners_trt_folds = .learners_trt_folds,
        .return_full_fits = .return_full_fits,
+       .isotonic_constraint = .isotonic_constraint,
        .discrete = .discrete,
        .info = .info)
 }

tests/testthat/test-time_varying_treatment.R

@@ -25,10 +25,12 @@ truth <- 0.305
 
 tmle <- sw(lmtp_tmle(tmp, a, "Y", time_vary = time_varying, shift = d, mtp = T, folds = 1))
 sdr <- sw(lmtp_sdr(tmp, a, "Y", time_vary = time_varying, shift = d, mtp = T, folds = 1))
-curve <- sw(lmtp_curve(tmp, a, "Y", time_vary = time_varying, shift = d, mtp = T, folds = 1, k = 0))
+curve <- sw(lmtp_curve(tmp, a, "Y", time_vary = time_varying, shift = d, mtp = T, folds = 1, k = 0, control = lmtp_control(.isotonic_constraint = FALSE)))
+curve_isotonic <- sw(lmtp_curve(tmp, a, "Y", time_vary = time_varying, shift = d, mtp = T, folds = 1, k = 0))
 
 test_that("time varying treatment fidelity, t = 4", {
   expect_equal(truth, tmle$estimate@x, tolerance = 0.01)
   expect_equal(truth, sdr$estimate@x, tolerance = 0.01)
   expect_equal(truth, curve$estimates[[4]]@x, tolerance = 0.01)
+  expect_equal(truth, curve_isotonic$estimates[[4]]@x, tolerance = 0.01)
 })