Develop surv estimating functions

R-package/targeted/DESCRIPTION

@@ -52,7 +52,9 @@ Suggests:
     SuperLearner (>= 2.0-28),
     knitr,
     xgboost,
-    viridisLite
+    viridisLite,
+    ranger,
+    randomForestSRC
 BugReports: https://github.com/kkholst/targeted/issues
 URL: https://kkholst.github.io/targeted/
 License: Apache License (== 2.0)

R-package/targeted/NAMESPACE

@@ -49,6 +49,7 @@ export(calibration)
 export(cate)
 export(cate_link)
 export(crr)
+export(cumhaz)
 export(cv)
 export(design)
 export(expand.list)
@@ -84,6 +85,7 @@ import(Rcpp)
 import(methods)
 importFrom(R6,R6Class)
 importFrom(data.table,data.table)
+importFrom(data.table,dcast)
 importFrom(data.table,is.data.table)
 importFrom(grDevices,nclass.Sturges)
 importFrom(graphics,abline)

R-package/targeted/R/H_constructor.R

@@ -0,0 +1,131 @@
+## H(u|X_i, A_i) = E[I\{T_i \leq \tau\} | T_i > u, X_i, A_i] = I\{u \leq \tau \} \frac{S(u|X_i, A_i) - S(\tau|X_i, A_i)}{S(u|X_i, A_i)}
+H_constructor_risk <- function(T_model, tau, individual_time, ...) {
+  force(T_model)
+  force(tau)
+  force(individual_time)
+
+  H <- function(u, data) {
+    S <- cumhaz(
+      T_model,
+      newdata = data,
+      times = u,
+      individual.time = individual_time
+    )$surv
+
+    S_tau <- cumhaz(
+      T_model,
+      newdata = data,
+      times = tau,
+      individual_time = individual_time
+    )$surv
+    S_tau <- as.vector(S_tau)
+
+
+    if (individual_time == FALSE) {
+      res <- apply(S, 2, function(x) x - S_tau)
+      res <- res / S
+      indicator <- (u <= tau)
+      res <- apply(res, 1, function(x) x * indicator)
+      res <- t(res)
+    } else {
+      res <- (S - S_tau) / S * (u <= tau)
+      res <- as.vector(res)
+    }
+    return(res)
+  }
+  return(H)
+}
+
+## Hu: H(u|X_i, A_i) = E[I\{T_i > \tau\} | T_i \geq u, X_i, A_i] = I\{u \leq \tau \} \frac{S(\tau|X_i, A_i)}{S(u|X_i, A_i)} +  I\{u > \tau \}
+H_constructor_surv <- function(T_model, tau, individual_time, ...) {
+  force(T_model)
+  force(tau)
+  force(individual_time)
+
+  H <- function(u, data) {
+    S <- cumhaz(
+      T_model,
+      newdata = data,
+      times = u,
+      individual.time = individual_time
+    )$surv
+
+    S_tau <- cumhaz(
+      T_model,
+      newdata = data,
+      times = tau,
+      individual_time = individual_time
+    )$surv
+    S_tau <- as.vector(S_tau)
+
+    if (individual_time == FALSE) {
+      res <- 1 / S
+      res <- apply(res, 2, function(x) x * S_tau)
+      indicator_1 <- (u <= tau)
+      indicator_2 <- (u > tau)
+      res <- apply(res, 1, function(x) x * indicator_1 + indicator_2)
+      res <- t(res)
+    } else {
+      res <- S_tau / S * (u <= tau) + (u > tau)
+      res <- as.vector(res)
+    }
+
+    return(res)
+  }
+  return(H)
+}
+
+## Hu: H_\tau(u|X_i, A_i) = E[\min(T, \tau) | T_i \geq u, X_i, A_i] = u + \frac{1}{S(u|X,A)} \int_u^\tau S(t|X,A) dt
+H_constructor_rmst <- function(T_model, time, event, tau, individual_time, sample = 0) {
+  force(T_model)
+  force(tau)
+  force(individual_time)
+  force(time)
+  force(event)
+  force(sample)
+
+  H <- function(u, data) {
+    S <- cumhaz(
+      T_model,
+      newdata = data,
+      times = u,
+      individual.time = individual_time
+    )$surv
+
+    tt <- time[event == 1]
+    if (sample > 0) {
+      tt <- subjumps(tt, size = sample, tau = tau)
+    }
+
+    S_T <- cumhaz(
+      T_model,
+      newdata = data,
+      times = tt,
+      individual.time = FALSE
+    )$surv
+    if (individual_time == FALSE) {
+      int_S <- apply(
+        S_T,
+        1,
+        function(x) {
+          int_surv(times = tt, surv = x, start = u, stop = tau, extend = FALSE)
+        },
+        simplify = FALSE
+      )
+      int_S <- do.call(what = "rbind", int_S)
+      res <- 1 / S
+      res <- res * int_S
+      res <- apply(res, 1, function(x) x + pmin(u, tau))
+      res <- t(res)
+    } else {
+      int_S <- numeric(length = length(u))
+      for (k in seq_along(u)) {
+        int_S[k] <- int_surv(times = tt, surv = S_T[k, ], start = u[k], stop = tau, extend = FALSE)
+      }
+      res <- pmin(u, tau) + 1 / S * int_S
+    }
+    return(res)
+  }
+
+  return(H)
+}

R-package/targeted/R/RATE.R

@@ -328,53 +328,17 @@ RATE.surv <- function(response, post.treatment, treatment, censoring,
   return(out)
 }
 
-cumhaz <- function(object, newdata, times=NULL, ...) {
-  if (inherits(object, "phreg")) {
-    if (is.null(times)) times <- object$times
-    pp <- predict(object, newdata=newdata,
-                  times=times,
-                  individual.times=FALSE, ...)
-    chf <- t(pp$cumhaz)
-    tt <- pp$times
-  } else if (inherits(object, "rfsrc")) {
-    pp <- predict(object, newdata=newdata, oob=TRUE, ...)
-    chf <- t(rbind(pp$chf))
-    tt <- pp$time.interest
-    if (!is.null(times)) {
-      idx <- mets::fast.approx(tt, times)
-      chf <- chf[idx,,drop=FALSE]
-      tt <- times
-    }
-  } else if (inherits(object, "ranger")) {
-    num.threads <- object$call$num.threads
-    pp <- predict(object, type="response", data=newdata, num.threads = num.threads, ...)
-    chf <- t(rbind(pp$chf))
-    tt <- pp$unique.death.times
-    if (!is.null(times)) {
-      idx <- mets::fast.approx(tt, times)
-      chf <- chf[idx,,drop=FALSE]
-      tt <- times
-    }
-  } else if (inherits(object, "coxph")) {
-    pp <- survfit(object, newdata=newdata)
-    pp <- summary(pp, time=times)
-    chf <- rbind(pp$cumhaz)
-    tt <- pp$time
-  }
-  list(time=tt, chf=chf, surv=exp(-chf), dchf=diff(rbind(0,chf)))
-}
-
 F.tau <- function(T.est, D.est, data, tau, a, treatment, post.treatment){
 
   data[lava::getoutcome(treatment)] <- a
   pred.D <- predict(D.est, type = "response", data)
 
   data[lava::getoutcome(post.treatment)] <- 1
-  surv.T.D1 <- cumhaz(T.est, newdata = data, times = tau)$surv[1,]
+  surv.T.D1 <- cumhaz(T.est, newdata = data, times = tau)$surv[ , 1]
 
   data[lava::getoutcome(post.treatment)] <- 0
-  surv.T.D0 <- cumhaz(T.est, newdata = data, times = tau)$surv[1,]
-
+  surv.T.D0 <- cumhaz(T.est, newdata = data, times = tau)$surv[ , 1]
+g
   surv <- pred.D * surv.T.D1 + (1 - pred.D) * surv.T.D0
 
   f.tau <- 1 - surv
@@ -392,7 +356,7 @@ HMc.tau <- function(T.est, C.est, data, time, event, tau){
   time.C <- time[event == 0]
 
   S <- diag(cumhaz(T.est, newdata = data.C, times = time.C)$surv)
-  S.tau <- cumhaz(T.est, newdata = data.C, times = tau)$surv[1,]
+  S.tau <- cumhaz(T.est, newdata = data.C, times = tau)$surv[ , 1]
   Sc <- diag(cumhaz(C.est, newdata = data.C, times = time.C)$surv)
   stopifnot(all(S * Sc> 0))