require(survival)
prob_predict.coxph <- function(object, newdata, times) {
fit.detail <- suppressWarnings(basehaz(object))
cum.haz <- fit.detail$hazard[sapply(times, function(x) max(which(fit.detail$time <= x)))]
predX <- predict(object = object, newdata = newdata, type = "risk")
res <- matrix(NA, ncol = length(times), nrow = length(predX))
for (ti in seq_len(length(times))) {
res[, ti] <- exp(-predX * cum.haz[ti])
}
res
}
set.seed(68)
n <- 500
Z <- rnorm(n)
X <- rbinom(n, 1, prob = 0.5)
T <- rexp(n, rate = exp(X + Z + X * Z)) # survival time
C <- rexp(n, rate = exp(X + Z + X * Z)) # censoring time
U <- pmin(T, C) # time at risk
D <- as.numeric(T < C) # event indicator
dd <- data.frame(Z, X, U, D)
x <- standardize_level(
fitter_list = list("coxph", "coxph"),
arguments = list(
list(
formula = Surv(U, D) ~ X + Z + X * Z,
method = "breslow",
x = TRUE,
y = TRUE
),
list(
formula = Surv(U, D) ~ X,
method = "breslow",
x = TRUE,
y = TRUE
)
),
predict_fun_list = list(prob_predict.coxph, prob_predict.coxph),
data = dd,
times = seq(1, 5, 0.1),
values = list(X = c(0, 1)),
B = 100,
reference = 0,
contrasts = "difference"
)
print(x)
Run the code above in your browser using DataLab