Make sure all models can handle out of domain survival probability predictions

[tidymodels]

Make sure all models can handle out of domain survival probability predictions

[mattwarkentin]

I see you using the summary.flexsurvreg method for predictions (both here and in the package source code). You may wish to use predict.flexsurvreg which I've implemented to be very tidy-friendly.

censored/R/survival_reg.R

Lines 53 to 61 in 14acd77

	flexsurv_probs <- function(object, new_data, time, type = "survival") {
	type <- rlang::arg_match(type, c("survival", "hazard"))
	res <- summary(object, newdata = new_data, type = type, t = time, ci = FALSE)
	res <- unname(res)
	col_name <- rlang::sym(paste0(".pred_", type))
	res <- purrr::map(res, ~ dplyr::select(.x, time, est))
	res <- purrr::map(res, ~ setNames(.x, c(".time", col_name)))
	tibble::tibble(.pred = res)
	}

[topepo]

Here's some code to get probability estimates from the survival package (verified against flexsurv):

library(survival)
library(flexsurv) # <- test against this

l_norm <- survreg(Surv(futime, fustat) ~ ecog.ps + rx, ovarian, dist = "lognormal")
l_norm_2 <- flexsurvreg(Surv(futime, fustat) ~ ecog.ps + rx, data = ovarian, dist = "lnorm")

para_survival <- function(object, time, location, scale = object$scale, ...) {
  distr <- object$dist
  tibble::tibble(
    .time = time, 
    .prob_survival = 1 - psurvreg(time, location, distribution = distr, scale, ...)
  )
}

para_hazard <- function(object, time, location, scale = object$scale, ...) {
  distr <- object$dist
  prob <- 
    dsurvreg(time, location, scale, distribution = distr, ...) /
    (1 - psurvreg(time, location, distribution = distr, scale, ...))
  tibble::tibble(
    .time = time, 
    .prob_hazard = prob
  )
}

surv_via_survival <- para_survival(l_norm, time = 1:300, predict(l_norm, head(ovarian, 1), type = "lp"))
surv_via_flexsurv <- summary(l_norm_2, newdata = head(ovarian, 1), type = "survival", t = 1:300)[[1]]
all.equal(surv_via_flexsurv$est, surv_via_survival$.prob_survival, tol = .0001)

haz_via_survival <- para_hazard(l_norm, time = 1:300, predict(l_norm, head(ovarian, 1), type = "lp"))
haz_via_flexsurv <- summary(l_norm_2, newdata = head(ovarian, 1), type = "hazard", t = 1:300)[[1]]
all.equal(haz_via_flexsurv$est, haz_via_survival$.prob_hazard, tol = .0001)

I'll put this in a PR soon.

[hfrick]

thanks @mattwarkentin ! 🙌

[topepo]

That's great. Would you like to do a PR for that?

(We ask them for this a while back I think)

[mattwarkentin]

I'm not sure what you mean. A PR for this package or a PR for {flexsurv}?

The predict.flexsurvreg method has been available in the dev version since 2020-06-10 and on CRAN since flexsurv version 2.0 (2021-02-22), I believe.

[github-actions]

This issue has been automatically locked. If you believe you have found a related problem, please file a new issue (with a reprex: https://reprex.tidyverse.org) and link to this issue.