Question: Get predicted survival curve for estimated ATE on a new population

[krz]

Hi,

grf is a great package and I'm following the excellent documentation ATE transport example here: example https://grf-labs.github.io/grf/articles/ate_transport.html

I adapted it to a use case where I use a causal survival forest instead of the causal forest in the example.

Is it possible to get predicted survival curve from this setup? Ideally two survival curves with treatment=1 vs treatment=0 on the target population.

Or produce a RMST plot like this?

Thanks!

[erikcs]

Hi @krz, causal survival forest estimates a single point on the survival curve difference, but you could re-estimate the forest at many horizon time points to trace out a survival curve.

In order to get estimates for a single arm, you could follow the same logic as in this vignette. The implied point on the treatment counterfactual curve would be P[Y(1) > h | X] = m(X) - e(X) tau^h(X), where m(X) is the usual Y.hat argument in CSF.

For the RMST, it could maybe be useful to recall that it is equal to the area under the survival curve (or in case of CSF, it estimates the area between the counterfactual survival curves up to horizon).

[krz]

Hi @erikcs , thanks for the explanation. I know that the RMST is equal to the area between the counterfactual survival curves, so I hoped these survival curves are already computed and can be extracted somehow.

Using your advice, I tried the following to estimate the curves at different horizons, but the resulting curves look very weird.
I'm probably doing something wrong.

# data generation
n <- 700
p <- 5
X <- matrix(runif(n * p), n, p)
W <- rbinom(n, 1, 0.5)
horizon <- 10
event.time <- pmin(rexp(n, rate = 0.1 + W * X[, 1]), horizon)
censor.time <- runif(n, 0, 12)
Y <- round(pmin(event.time, censor.time), 1)
D <- as.integer(event.time <= censor.time)

# random test data
X.test <- matrix(runif(n * p), n, p)
n.rct <- length(Y)
n.obs <- nrow(X.test)

# estimate on a grid
horizon_max <- 10
grid <- seq(0.2, horizon_max, by = 0.2)

S1_curve <- numeric(length(grid))
S0_curve <- numeric(length(grid))

for (k in seq_along(grid)) {
  h <- grid[k]
  forest <- causal_survival_forest(X, Y, W, D, W.hat=0.5, horizon=h, 
                        target="survival.probability")
  S.forest <- regression_forest(X = rbind(forest$X.orig, X.test),
                                  Y = c(rep(1, n.rct), rep(0, n.obs)))
  p.hat <- predict(S.forest)$predictions[1:n.rct]

  tau.hat.test <- predict(forest, X.test)$predictions
  tau.hat.train <- predict(forest)$predictions

  Y.hat <- forest$Y.hat
  Y.orig <- forest$Y.orig
  W.hat <- forest$W.hat
  W.orig <- forest$W.orig

  debiasing.weights <- (1 - p.hat) / p.hat * (W.orig - W.hat) / (W.hat * (1 - W.hat))
  Y.residual <- Y.orig - (Y.hat + tau.hat.train * (W.orig - W.hat))

  tau.hat <- mean(tau.hat.test) + sum(debiasing.weights * Y.residual) / n.obs

  mu.hat.0 <- Y.hat - W.hat * tau.hat
  mu.hat.1 <- Y.hat + (1 - W.hat) * tau.hat

  S0_curve[k] <- mean(mu.hat.0)
  S1_curve[k] <- mean(mu.hat.1)

}

but the plot looks very weird

plot(grid, S0_curve, type = "l", col = "steelblue", lwd = 2,
     ylim=c(-1.5,1.5),
     xlab = "Time (horizon)", ylab = "Survival probability")
lines(grid, S1_curve, col = "tomato", lwd = 2)