bgreenwell / investr

Inverse estimation in R
22 stars 6 forks source link

Add functions for survival models #30

Open bgreenwell opened 7 years ago

bgreenwell commented 7 years ago

Here's a good start for single Kaplan-Meier curves:

#' @keywords internal
minmin <- function(y, x, prob = 0.5) {
  # adapted from function inside survival:::survmean
  tolerance <- .Machine$double.eps ^ 0.5
  keep <- (!is.na(y) & y < (prob + tolerance))
  if (!any(keep)) {
    NA
  } else {
    x <- x[keep]
    y <- y[keep]
    if (abs(y[1L] - prob) < tolerance && any(y < y[1L])) 
      (x[1L] + x[min(which(y < y[1L]))]) / 2
    else x[1L]
  }
}

invest.Surv <- function(object, prob = 0.5, level = 0.95, ...) {
  km <- survival::survfit(object ~ 1, conf.int = level, ...)
  est <- minmin(km$surv, km$time, prob = prob)
  upr <- minmin(km$upper, km$time, prob = prob)
  lwr <- minmin(km$lower, km$time, prob = prob)
  plot(km, conf.int = TRUE)
  abline(h = prob, v = c(lwr, est, upr))
  c("est" = est, "lwr" = lwr, "upr" = upr)
}

# Example
library(survival)
invest.Surv(Surv(aml$time, aml$status), prob = 0.5)
bgreenwell commented 7 years ago

And for "coxph" objects:

invest.coxph <- function(object, prob = 0.5, level = 0.95, ...) {
  sf <- survfit(object, conf.int = level, ...)
  est <- minmin(sf$surv, sf$time, prob = prob)
  upr <- minmin(sf$upper, sf$time, prob = prob)
  lwr <- minmin(sf$lower, sf$time, prob = prob)
  plot(sf, conf.int = TRUE)
  abline(h = prob, v = c(lwr, est, upr))
  c("est" = est, "lwr" = lwr, "upr" = upr)
}

# Example
fit <- coxph(Surv(futime, fustat) ~ age, data = ovarian) 
invest.coxph(fit, prob = 0.9)
invest.coxph(fit, prob = 0.9, newdata = data.frame(age = 60))