version 0.0-2

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: tramicp
 Type: Package
 Title: Model-Based Causal Feature Selection for General Response Types
-Version: 0.0-1
+Version: 0.0-2
 Authors@R: c(person("Lucas", "Kook", role = c("aut", "cre"),
     comment = c("ORCID" = "0000-0002-7546-7356"),
     email = "lucasheinrich.kook@gmail.com"), person("Sorawit", "Saengkyongam",
@@ -15,19 +15,19 @@ Description: Extends invariant causal prediction (Peters et al., 2016,
     <doi:10.48550/arXiv.2309.12833>).
 License: GPL-3
 Encoding: UTF-8
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
 Depends: R (>= 4.0.0)
 Imports: tram, mlt, coin, multcomp, survival, variables, basefun, MASS,
         cotram, dHSIC, ranger, sandwich
 Suggests: testthat (>= 3.0.0)
 Config/testthat/edition: 3
 NeedsCompilation: no
-Packaged: 2023-09-29 13:26:28 UTC; stigler
+Packaged: 2024-03-14 15:08:22 UTC; lkook
 Author: Lucas Kook [aut, cre] (<https://orcid.org/0000-0002-7546-7356>),
   Sorawit Saengkyongam [ctb],
   Anton Rask Lundborg [ctb],
   Torsten Hothorn [ctb],
   Jonas Peters [ctb]
 Maintainer: Lucas Kook <lucasheinrich.kook@gmail.com>
 Repository: CRAN
-Date/Publication: 2023-09-29 13:52:43 UTC
+Date/Publication: 2024-03-14 15:30:02 UTC

MD5

@@ -1,23 +1,23 @@
-7f0c353ee5ccbc7ea096092524e828a9 *DESCRIPTION
-01f51cfad7e0eb0be21c69208701e24a *NAMESPACE
-f0d4e7659a4c1c35174b8932ed2f41f3 *R/alias.R
+5d525c4141fb24f81bffadd660d221e9 *DESCRIPTION
+87bbb2b2e92e87c75bbd1369e8fd3f98 *NAMESPACE
+b94231a2cf30f6cbc1d7da77202e22e6 *R/alias.R
 b8bf164fce1e55b435fd60ae6b67da5f *R/boostrap-stability.R
-4268f00196c8c6634a004876b642a9b4 *R/controls.R
+b3923230b96b4aa0a498dc734ad9abc0 *R/controls.R
 96acd7e17489733996341194190c968d *R/dgps.R
-f7b53792e9080725e6c33196463f0e0f *R/helpers.R
-a79be13ee4dfd8252be7fd786c1e2038 *R/invariance-types.R
+35c690076af486c4988e86ffaf530ccb *R/helpers.R
+7918d523dbcb8d104221b181b4c05d95 *R/invariance-types.R
 d8827edbb242c3ec907b26ee14cf44a6 *R/methods.R
 678d5bfb15106bb0b5b445e8acd02e30 *R/tests.R
 67722902b9592c98c4e3300994540e0d *R/tram-sim.R
-d41c89281b25b7dee35e1ee4ef5aaece *R/tramicp.R
-6dcc673e53b5d5b2ab8b2f78399d60ba *build/partial.rdb
+08426fc98a54520516be978318f6c980 *R/tramicp.R
+20c8d3f97d93e2e318a22508275b5a02 *build/partial.rdb
 20ea6e052a5b798fd68a35a71c73d0c5 *man/bootstrap_stability.Rd
 3115c939b6e6b2532e9c0983e6c8e2a9 *man/dgp_dicp.Rd
 d7eb6f6aa1634d8e536f8e50ce80500d *man/dicp.Rd
 cfad8a947f0b967dafd83a8b7283f2c6 *man/dicp_controls.Rd
 be0d0d1fb2b74a941540ca99e3f9836f *man/invariant_sets.Rd
 bae6d437fe79a6c61849cfe7bda312b9 *man/pvalues.Rd
-824f229f685b2572af64445c7930c5f9 *man/tramicp-alias.Rd
+d1f16b22025d8d13d7644604f6fcf574 *man/tramicp-alias.Rd
 15fc7658313e6bf54c386a2fada38b53 *tests/testthat.R
 89e508aee768d0a419c6baa334a2b03b *tests/testthat/test-dgps.R
 57b9b77a72c6d1c8138334b8e7a7457f *tests/testthat/test-greedy.R

NAMESPACE

@@ -1,9 +1,12 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(predict,qrf)
 S3method(print,dICP)
 S3method(print,dICPtest)
 S3method(residuals,polr)
+S3method(residuals,qrf)
 S3method(residuals,ranger)
+S3method(residuals,survforest)
 S3method(summary,dICP)
 export(BoxCoxICP)
 export(ColrICP)
@@ -23,6 +26,8 @@ export(invariant_sets)
 export(lmICP)
 export(polrICP)
 export(pvalues)
+export(qrfICP)
 export(rangerICP)
+export(survforestICP)
 export(survregICP)
 import(tram)

R/alias.R

@@ -6,7 +6,9 @@
 #' @description ICP for Box-Cox-type transformed normal regression, parametric
 #'     and semiparametric survival models, continuous outcome logistic
 #'     regression, linear regression, cumulative ordered regression, generalized
-#'     linear models; and nonparametric ICP via ranger.
+#'     linear models; and nonparametric ICP via ranger. While TRAMICP based on
+#'     quantile and survival random forests is also supported, for these methods
+#'     it comes without theoretical guarantees as of yet.
 #' @rdname tramicp-alias
 #'
 #' @inheritParams dicp
@@ -45,6 +47,8 @@ BoxCoxICP <- function(formula, data, env, verbose = TRUE, type = "residual",
 #' d <- dgp_dicp(mod = "weibull", n = 300)
 #' SurvregICP(Y ~ X1 + X2 + X3, data = d, env = ~ E)
 #' ### or
+#' library("survival")
+#' d$Y <- Surv(d$Y)
 #' survregICP(Y ~ X1 + X2 + X3, data = d, env = ~ E)
 #' CoxphICP(Y ~ X2, data = d, env = ~ E)
 #' coxphICP(Y ~ X2, data = d, env = ~ E)
@@ -315,14 +319,10 @@ glmICP <- function(formula, data, env, verbose = TRUE, type = "residual",
   call <- match.call()
   if (is.character(test))
     test <- match.arg(test, .implemented_tests())
-  if (!is.null(fam <- list(...)$family) && (identical(fam, stats::binomial) ||
-                                            fam == "binomial")) {
-    resid <- "residuals.binglm"
-  }
   if (is.null(controls))
     controls <- dicp_controls(type, test, alpha = alpha,
                               baseline_fixed = baseline_fixed,
-                              residuals = resid)
+                              residuals = "residuals.binglm")
   ret <- dicp(formula = formula, data = data, env = env, modFUN = stats::glm,
               verbose = verbose, type = type, test = test, controls = controls,
               alpha = alpha, baseline_fixed = baseline_fixed, greedy = greedy,
@@ -384,3 +384,60 @@ rangerICP <- function(formula, data, env, verbose = TRUE, type = "residual",
   ret$call <- call
   ret
 }
+
+#' nonparametric ICP for right-censored observations with survival forest GCM
+#' @rdname tramicp-alias
+#'
+#' @inheritParams dicp
+#'
+#' @export
+#'
+#' @examples
+#' \donttest{
+#' set.seed(12)
+#' d <- dgp_dicp(mod = "coxph", n = 3e2)
+#' d$Y <- survival::Surv(d$Y, sample(0:1, 3e2, TRUE, prob = c(0.1, 0.9)))
+#' survforestICP(Y ~ X1 + X2 + X3, data = d, env = ~ E)
+#' }
+#'
+survforestICP <- function(formula, data, env, verbose = TRUE, type = "residual",
+                          test = "gcm.test", controls = NULL, alpha = 0.05,
+                          baseline_fixed = TRUE, greedy = FALSE, max_size = NULL,
+                          mandatory = NULL, ...) {
+  call <- match.call()
+  message("Note: `survforestICP()` does not come with theoretical guarantees.")
+  ret <- dicp(formula = formula, data = data, env = env, modFUN = survforest,
+              verbose = verbose, type = type, test = test, controls = controls,
+              alpha = alpha, baseline_fixed = baseline_fixed, greedy = greedy,
+              max_size = max_size, mandatory = mandatory, ... = ...)
+  ret$call <- call
+  ret
+}
+
+#' nonparametric ICP with quantile forest GCM
+#' @rdname tramicp-alias
+#'
+#' @inheritParams dicp
+#'
+#' @export
+#'
+#' @examples
+#' \donttest{
+#' set.seed(12)
+#' d <- dgp_dicp(mod = "boxcox", n = 3e2)
+#' qrfICP(Y ~ X1 + X2 + X3, data = d, env = ~ E)
+#' }
+#'
+qrfICP <- function(formula, data, env, verbose = TRUE, type = "residual",
+                   test = "gcm.test", controls = NULL, alpha = 0.05,
+                   baseline_fixed = TRUE, greedy = FALSE, max_size = NULL,
+                   mandatory = NULL, ...) {
+  call <- match.call()
+  message("Note: `qrfICP()` does not come with theoretical guarantees.")
+  ret <- dicp(formula = formula, data = data, env = env, modFUN = qrf,
+              verbose = verbose, type = type, test = test, controls = controls,
+              alpha = alpha, baseline_fixed = baseline_fixed, greedy = greedy,
+              max_size = max_size, mandatory = mandatory, ... = ...)
+  ret$call <- call
+  ret
+}

R/controls.R

@@ -79,7 +79,7 @@ dicp_controls <- function(
 
 .test_fun <- function(type, test, ctest) {
   if (is.function(test))
-    return(list(test = "custom", test_fun = test_fun, test_name = ctest))
+    return(list(test = "custom", test_fun = test, test_name = ctest))
 
   if (type %in% c("wald", "partial")) {
     ctest <- "wald"

R/helpers.R

@@ -129,10 +129,7 @@ intersect_intervals <- function(...) {
 
 #' @method residuals binglm
 residuals.binglm <- function(object, ...) {
-  resp <- stats::model.response(stats::model.frame(object))
-  success <- if (is.factor(resp)) levels(resp)[2] else sort(unique(resp))[2]
-    y <- c(0, 1)[1 + as.numeric(resp == success)]
-  y - stats::predict(object, type = "response")
+  stats::residuals.glm(object, type = "response")
 }
 
 .check_depth <- function(x) {
@@ -184,7 +181,7 @@ residuals.binglm <- function(object, ...) {
   res <- lapply(terms, \(term) suppressWarnings(
     max(pvals[!grepl(term, names(pvals), fixed = TRUE)], na.rm = TRUE)))
   ret <- structure(unlist(res), names = terms)
-  if (all(ret < alpha))
+  if (all(pvals < alpha, na.rm = TRUE) & length(ret) > 1)
     ret[] <- 1
   ret
 }
@@ -309,3 +306,71 @@ residuals.ranger <- function(object, newdata = NULL, newy = NULL, ...) {
     return(character(0))
   ret
 }
+
+# Survival and quantile RF
+
+survforest <- function(formula, data, ...) {
+  tms <- .get_terms(formula)
+  if (identical(tms$me, character(0))) {
+    return(survival::coxph(formula, data))
+  }
+  rf <- ranger::ranger(formula, data, ...)
+  class(rf) <- c("survforest", class(rf))
+  rf$y <- stats::model.response(stats::model.frame(formula, data))
+  rf$data <- data
+  rf
+}
+
+#' @exportS3Method residuals survforest
+residuals.survforest <- function(object, ...) {
+  times <- object$y[, 1]
+  status <- object$y[, 2]
+  pred <- stats::predict(object, data = object$data)
+  idx <- sapply(times, \(x) which.min(abs(x - pred$unique.death.times))[1])
+  preds <- pred$survival
+  ipreds <- sapply(seq_len(nrow(preds)), \(smpl) {
+    -log(preds[smpl, idx[smpl]])
+  })
+  status - ipreds
+}
+
+qrf <- \(formula, data, ...) {
+  rY <- stats::model.response(stats::model.frame(formula, data))
+  tms <- .get_terms(formula)
+  if (identical(tms$me, character(0))) {
+    ret <- list(m = stats::ecdf(rY), data = data, response = rY,
+                unconditional = TRUE)
+    class(ret) <- c("qrf")
+    return(ret)
+  }
+  rf <- ranger::ranger(formula, data, ...)
+  rf$response <- rY
+  rf$data <- data
+  rf$unconditional <- FALSE
+  class(rf) <- c("qrf", class(rf))
+  rf
+}
+
+#' @exportS3Method predict qrf
+predict.qrf <- \(object, data, ...) {
+  if (object$unconditional)
+    return(object$m(object$response))
+  class(object) <- class(object)[-1]
+  tn <- stats::predict(object, data = data, type = "terminalNodes")$predictions
+  K <- matrix(0, nrow = N <- nrow(tn), ncol = N)
+  for (tree in seq_len(B <- object$num.trees)) {
+    K <- K + sapply(seq_len(nrow(tn)), \(obs) {
+      as.numeric(tn[obs, tree] == tn[, tree])
+    }, simplify = "matrix")
+  }
+  K <- K / B
+  diag(K) <- 0
+  K <- K / pmax(colSums(K), .Machine$double.eps)
+  pred <- \(y) mean(K %*% as.numeric(object$response <= y))
+  sapply(object$response, pred)
+}
+
+#' @exportS3Method residuals qrf
+residuals.qrf <- \(object, ...) {
+  2 * predict.qrf(object, object$data) - 1
+}

R/invariance-types.R

@@ -27,8 +27,8 @@
 
   ### Return
   if (set == 1) tset <- "Empty"
-  structure(list(set = tset, test = tst, coef = stats::coef(m), tram = m$tram),
-            class = "dICPtest")
+  structure(list(set = tset, test = tst, coef = stats::coef(m), tram = m$tram,
+                 rYX = r, rEX = e), class = "dICPtest")
 
 }
 

R/tramicp.R

@@ -211,8 +211,8 @@ invariant_sets <- function(object, with_pvalues = FALSE) {
         modFUN = modFUN, data = data, controls = controls,
         mandatory = mandatory, ... = ...
       )
-
-      if (.get_pvalue(ret$test) > controls$alpha) {
+      tpv <- .get_pvalue(ret$test)
+      if (!is.nan(tpv) && tpv > controls$alpha) {
         MI <- c(MI, lps[[set]])
       }