From 819dab496e1b76eb1ed5417d8d89a1822477038c Mon Sep 17 00:00:00 2001 From: Emily de la Rua <59304861+edelarua@users.noreply.github.com> Date: Thu, 23 May 2024 18:39:06 -0400 Subject: [PATCH] Update `ard_survival_survfit` (#152) **What changes are proposed in this pull request?** * Updated `ard_survival_survfit` to return `std.error` and `n.risk` statistics. (#139 ) @ddsjoberg as I currently have it, `std.error` gets transformed along with `estimate`, `conf.low`, and `conf.high` (with options `"survival"`, `"cumhaz"`, and `"risk"`), but if you think it's more informative to always return the untransformed `std.error` that can be updated. Closes #139 -------------------------------------------------------------------------------- Pre-review Checklist (if item does not apply, mark is as complete) - [x] **All** GitHub Action workflows pass with a :white_check_mark: - [x] PR branch has pulled the most recent updates from master branch: `usethis::pr_merge_main()` - [x] If a bug was fixed, a unit test was added. - [x] If a new `ard_*()` function was added, it passes the ARD structural checks from `cards::check_ard_structure()`. - [x] If a new `ard_*()` function was added, `set_cli_abort_call()` has been set. - [x] If a new `ard_*()` function was added and it depends on another package (such as, `broom`), `is_pkg_installed("broom", reference_pkg = "cardx")` has been set in the function call and the following added to the roxygen comments: `@examplesIf do.call(asNamespace("cardx")$is_pkg_installed, list(pkg = "broom"", reference_pkg = "cardx"))` - [x] Code coverage is suitable for any new functions/features (generally, 100% coverage for new code): `devtools::test_coverage()` Reviewer Checklist (if item does not apply, mark is as complete) - [ ] If a bug was fixed, a unit test was added. - [ ] Code coverage is suitable for any new functions/features: `devtools::test_coverage()` When the branch is ready to be merged: - [ ] Update `NEWS.md` with the changes from this pull request under the heading "`# cardx (development version)`". If there is an issue associated with the pull request, reference it in parentheses at the end update (see `NEWS.md` for examples). - [ ] **All** GitHub Action workflows pass with a :white_check_mark: - [ ] Approve Pull Request - [ ] Merge the PR. Please use "Squash and merge" or "Rebase and merge". --- R/ard_survival_survfit.R | 74 +++---- man/dot-process_survfit_time.Rd | 5 +- tests/testthat/_snaps/ard_survival_survfit.md | 194 +++++++++++------- 3 files changed, 154 insertions(+), 119 deletions(-) diff --git a/R/ard_survival_survfit.R b/R/ard_survival_survfit.R index 9a30adf5f..e430eedc0 100644 --- a/R/ard_survival_survfit.R +++ b/R/ard_survival_survfit.R @@ -111,6 +111,8 @@ ard_survival_survfit <- function(x, times = NULL, probs = NULL, type = NULL) { #' #' @inheritParams cards::tidy_as_ard #' @inheritParams ard_survival_survfit +#' @param start.time (`numeric`)\cr +#' default starting time. See [survival::survfit0()] for more details. #' #' @return a `tibble` #' @@ -119,42 +121,39 @@ ard_survival_survfit <- function(x, times = NULL, probs = NULL, type = NULL) { #' cardx:::.process_survfit_time(times = c(60, 180), type = "risk") #' #' @keywords internal -.process_survfit_time <- function(x, times, type) { - # tidy survfit results - tidy_x <- broom::tidy(x) +.process_survfit_time <- function(x, times, type, start.time = NULL) { + # add start time + min_time <- min(x$time) + if (is.null(start.time) && min_time < 0) { + cli::cli_inform(paste( + "The {.arg start.time} argument has not been set and negative times have been observed. Please set start", + "time via the {.arg start.time} argument, otherwise the minimum observed time will be used by default." + )) + start.time <- min_time + } else if (is.null(start.time)) { + start.time <- 0 + } + x <- survival::survfit0(x, start.time) %>% + summary(times) # process competing risks/multi-state models - multi_state <- inherits(x, "survfitms") + multi_state <- inherits(x, "summary.survfitms") - if (multi_state == TRUE) { + if (multi_state) { # selecting state to show - state <- setdiff(unique(tidy_x$state), "(s0)")[[1]] + state <- setdiff(unique(x$states), "(s0)")[[1]] cli::cli_inform("Multi-state model detected. Showing probabilities into state '{state}'.") - tidy_x <- dplyr::filter(tidy_x, .data$state == .env$state) + x$n.risk <- x$n.risk[, 1] + ms_cols <- c("pstate", "std.err", "upper", "lower") + state_col <- which(colnames(x$pstate) == state) + x[ms_cols] <- lapply(x[ms_cols], function(m) m[, state_col]) + x$surv <- x$pstate } - # adding time 0 to data frame - tidy_x <- tidy_x %>% - # make strata a fct to preserve ordering - dplyr::mutate(dplyr::across(dplyr::any_of("strata"), ~ factor(., levels = unique(.)))) %>% - # if CI is missing and SE is 0, use estimate as the CI - dplyr::mutate_at( - dplyr::vars("conf.high", "conf.low"), - ~ ifelse(is.na(.) & .data$std.error == 0, .data$estimate, .) - ) %>% - dplyr::select(dplyr::any_of(c("time", "estimate", "conf.high", "conf.low", "strata"))) %>% - # add data for time 0 - dplyr::bind_rows( - dplyr::group_by_at(., dplyr::vars(dplyr::any_of("strata"))) %>% - dplyr::slice(1) %>% - dplyr::mutate( - time = 0, - estimate = ifelse(multi_state, 0, 1), - conf.low = ifelse(multi_state, 0, 1), - conf.high = ifelse(multi_state, 0, 1) - ) - ) %>% - dplyr::ungroup() + # tidy survfit results + x_cols <- intersect(names(x), c("time", "n.risk", "surv", "std.err", "upper", "lower", "strata")) + tidy_x <- data.frame(x[x_cols]) %>% + dplyr::rename(estimate = "surv", std.error = "std.err", conf.high = "upper", conf.low = "lower") strat <- "strata" %in% names(tidy_x) @@ -182,16 +181,7 @@ ard_survival_survfit <- function(x, times = NULL, probs = NULL, type = NULL) { } df_stat <- df_stat %>% - # if user-specifed time is unobserved, fill estimate with previous value dplyr::arrange(.data$time) %>% - dplyr::group_by_at(dplyr::vars(dplyr::any_of("strata"))) %>% - tidyr::fill( - "estimate", "conf.high", "conf.low", "time_max", - .direction = "down" - ) %>% - dplyr::ungroup() %>% - # keep only user-specified times - dplyr::filter(!is.na(.data$col_name)) %>% # if user-specified time is after max time, make estimate NA dplyr::mutate_at( dplyr::vars("estimate", "conf.high", "conf.low"), @@ -236,7 +226,7 @@ ard_survival_survfit <- function(x, times = NULL, probs = NULL, type = NULL) { as.data.frame() %>% set_names(c("estimate", "conf.low", "conf.high")) %>% dplyr::mutate(strata = row.names(.)) %>% - dplyr::select(dplyr::any_of(c("strata", "estimate", "conf.low", "conf.high"))) %>% + dplyr::select(dplyr::any_of(c("n.risk", "strata", "estimate", "std.error", "conf.low", "conf.high"))) %>% dplyr::mutate(prob = .x) ) %>% dplyr::bind_rows() %>% @@ -293,10 +283,10 @@ extract_multi_strata <- function(x, df_stat) { ret <- tidy_survfit %>% dplyr::mutate(dplyr::across( - dplyr::any_of(c("estimate", "conf.high", "conf.low", "time", "prob")), ~ as.list(.) + dplyr::any_of(c("n.risk", "estimate", "std.error", "conf.high", "conf.low", "time", "prob")), ~ as.list(.) )) %>% tidyr::pivot_longer( - cols = dplyr::any_of(c("estimate", "conf.high", "conf.low")), + cols = dplyr::any_of(c("n.risk", "estimate", "std.error", "conf.high", "conf.low")), names_to = "stat_name", values_to = "stat" ) %>% @@ -342,7 +332,9 @@ extract_multi_strata <- function(x, df_stat) { .df_survfit_stat_labels <- function() { dplyr::tribble( ~stat_name, ~stat_label, + "n.risk", "Number of Subjects at Risk", "estimate", "Survival Probability", + "std.error", "Standard Error (untransformed)", "conf.low", "CI Lower Bound", "conf.high", "CI Upper Bound", "conf.level", "CI Confidence Level", diff --git a/man/dot-process_survfit_time.Rd b/man/dot-process_survfit_time.Rd index c6f9bee25..cc3b6682c 100644 --- a/man/dot-process_survfit_time.Rd +++ b/man/dot-process_survfit_time.Rd @@ -4,7 +4,7 @@ \alias{.process_survfit_time} \title{Process Survival Fit For Time Estimates} \usage{ -.process_survfit_time(x, times, type) +.process_survfit_time(x, times, type, start.time = NULL) } \arguments{ \item{x}{(\code{\link[survival:survfit]{survival::survfit()}})\cr @@ -22,6 +22,9 @@ Must be one of the following:\tabular{ll}{ \code{"risk"} \tab \code{1 - x} \cr \code{"cumhaz"} \tab \code{-log(x)} \cr }} + +\item{start.time}{(\code{numeric})\cr +default starting time. See \code{\link[survival:survfit0]{survival::survfit0()}} for more details.} } \value{ a \code{tibble} diff --git a/tests/testthat/_snaps/ard_survival_survfit.md b/tests/testthat/_snaps/ard_survival_survfit.md index 6697507dd..7f2a2529d 100644 --- a/tests/testthat/_snaps/ard_survival_survfit.md +++ b/tests/testthat/_snaps/ard_survival_survfit.md @@ -5,27 +5,39 @@ CNSR) ~ TRTA, cards::ADTTE), times = c(60, 180)), stat = lapply(stat, function(x) ifelse(is.numeric(x), cards::round5(x, 3), x))), n = Inf) Message - {cards} data frame: 18 x 11 + {cards} data frame: 30 x 11 Output group1 group1_level variable variable_level stat_name stat_label stat - 1 TRTA Placebo time 60 estimate Survival… 0.893 - 2 TRTA Placebo time 60 conf.high CI Upper… 0.966 - 3 TRTA Placebo time 60 conf.low CI Lower… 0.825 - 4 TRTA Placebo time 180 estimate Survival… 0.651 - 5 TRTA Placebo time 180 conf.high CI Upper… 0.783 - 6 TRTA Placebo time 180 conf.low CI Lower… 0.541 - 7 TRTA Xanomeli… time 60 estimate Survival… 0.694 - 8 TRTA Xanomeli… time 60 conf.high CI Upper… 0.849 - 9 TRTA Xanomeli… time 60 conf.low CI Lower… 0.568 - 10 TRTA Xanomeli… time 180 estimate Survival… 0.262 - 11 TRTA Xanomeli… time 180 conf.high CI Upper… 0.749 - 12 TRTA Xanomeli… time 180 conf.low CI Lower… 0.092 - 13 TRTA Xanomeli… time 60 estimate Survival… 0.732 - 14 TRTA Xanomeli… time 60 conf.high CI Upper… 0.878 - 15 TRTA Xanomeli… time 60 conf.low CI Lower… 0.61 - 16 TRTA Xanomeli… time 180 estimate Survival… 0.381 - 17 TRTA Xanomeli… time 180 conf.high CI Upper… 0.743 - 18 TRTA Xanomeli… time 180 conf.low CI Lower… 0.195 + 1 TRTA Placebo time 60 n.risk Number o… 59 + 2 TRTA Placebo time 60 estimate Survival… 0.893 + 3 TRTA Placebo time 60 std.error Standard… 0.036 + 4 TRTA Placebo time 60 conf.high CI Upper… 0.966 + 5 TRTA Placebo time 60 conf.low CI Lower… 0.825 + 6 TRTA Placebo time 180 n.risk Number o… 35 + 7 TRTA Placebo time 180 estimate Survival… 0.651 + 8 TRTA Placebo time 180 std.error Standard… 0.061 + 9 TRTA Placebo time 180 conf.high CI Upper… 0.783 + 10 TRTA Placebo time 180 conf.low CI Lower… 0.541 + 11 TRTA Xanomeli… time 60 n.risk Number o… 14 + 12 TRTA Xanomeli… time 60 estimate Survival… 0.694 + 13 TRTA Xanomeli… time 60 std.error Standard… 0.071 + 14 TRTA Xanomeli… time 60 conf.high CI Upper… 0.849 + 15 TRTA Xanomeli… time 60 conf.low CI Lower… 0.568 + 16 TRTA Xanomeli… time 180 n.risk Number o… 3 + 17 TRTA Xanomeli… time 180 estimate Survival… 0.262 + 18 TRTA Xanomeli… time 180 std.error Standard… 0.14 + 19 TRTA Xanomeli… time 180 conf.high CI Upper… 0.749 + 20 TRTA Xanomeli… time 180 conf.low CI Lower… 0.092 + 21 TRTA Xanomeli… time 60 n.risk Number o… 20 + 22 TRTA Xanomeli… time 60 estimate Survival… 0.732 + 23 TRTA Xanomeli… time 60 std.error Standard… 0.068 + 24 TRTA Xanomeli… time 60 conf.high CI Upper… 0.878 + 25 TRTA Xanomeli… time 60 conf.low CI Lower… 0.61 + 26 TRTA Xanomeli… time 180 n.risk Number o… 5 + 27 TRTA Xanomeli… time 180 estimate Survival… 0.381 + 28 TRTA Xanomeli… time 180 std.error Standard… 0.13 + 29 TRTA Xanomeli… time 180 conf.high CI Upper… 0.743 + 30 TRTA Xanomeli… time 180 conf.low CI Lower… 0.195 Message i 4 more variables: context, fmt_fn, warning, error @@ -36,27 +48,39 @@ CNSR) ~ TRTA, cards::ADTTE), times = c(60, 180), type = "risk"), stat = lapply( stat, function(x) ifelse(is.numeric(x), cards::round5(x, 3), x))), n = Inf) Message - {cards} data frame: 18 x 11 + {cards} data frame: 30 x 11 Output group1 group1_level variable variable_level stat_name stat_label stat - 1 TRTA Placebo time 60 estimate Survival… 0.107 - 2 TRTA Placebo time 60 conf.high CI Upper… 0.175 - 3 TRTA Placebo time 60 conf.low CI Lower… 0.034 - 4 TRTA Placebo time 180 estimate Survival… 0.349 - 5 TRTA Placebo time 180 conf.high CI Upper… 0.459 - 6 TRTA Placebo time 180 conf.low CI Lower… 0.217 - 7 TRTA Xanomeli… time 60 estimate Survival… 0.306 - 8 TRTA Xanomeli… time 60 conf.high CI Upper… 0.432 - 9 TRTA Xanomeli… time 60 conf.low CI Lower… 0.151 - 10 TRTA Xanomeli… time 180 estimate Survival… 0.738 - 11 TRTA Xanomeli… time 180 conf.high CI Upper… 0.908 - 12 TRTA Xanomeli… time 180 conf.low CI Lower… 0.251 - 13 TRTA Xanomeli… time 60 estimate Survival… 0.268 - 14 TRTA Xanomeli… time 60 conf.high CI Upper… 0.39 - 15 TRTA Xanomeli… time 60 conf.low CI Lower… 0.122 - 16 TRTA Xanomeli… time 180 estimate Survival… 0.619 - 17 TRTA Xanomeli… time 180 conf.high CI Upper… 0.805 - 18 TRTA Xanomeli… time 180 conf.low CI Lower… 0.257 + 1 TRTA Placebo time 60 n.risk Number o… 59 + 2 TRTA Placebo time 60 estimate Survival… 0.107 + 3 TRTA Placebo time 60 std.error Standard… 0.036 + 4 TRTA Placebo time 60 conf.high CI Upper… 0.175 + 5 TRTA Placebo time 60 conf.low CI Lower… 0.034 + 6 TRTA Placebo time 180 n.risk Number o… 35 + 7 TRTA Placebo time 180 estimate Survival… 0.349 + 8 TRTA Placebo time 180 std.error Standard… 0.061 + 9 TRTA Placebo time 180 conf.high CI Upper… 0.459 + 10 TRTA Placebo time 180 conf.low CI Lower… 0.217 + 11 TRTA Xanomeli… time 60 n.risk Number o… 14 + 12 TRTA Xanomeli… time 60 estimate Survival… 0.306 + 13 TRTA Xanomeli… time 60 std.error Standard… 0.071 + 14 TRTA Xanomeli… time 60 conf.high CI Upper… 0.432 + 15 TRTA Xanomeli… time 60 conf.low CI Lower… 0.151 + 16 TRTA Xanomeli… time 180 n.risk Number o… 3 + 17 TRTA Xanomeli… time 180 estimate Survival… 0.738 + 18 TRTA Xanomeli… time 180 std.error Standard… 0.14 + 19 TRTA Xanomeli… time 180 conf.high CI Upper… 0.908 + 20 TRTA Xanomeli… time 180 conf.low CI Lower… 0.251 + 21 TRTA Xanomeli… time 60 n.risk Number o… 20 + 22 TRTA Xanomeli… time 60 estimate Survival… 0.268 + 23 TRTA Xanomeli… time 60 std.error Standard… 0.068 + 24 TRTA Xanomeli… time 60 conf.high CI Upper… 0.39 + 25 TRTA Xanomeli… time 60 conf.low CI Lower… 0.122 + 26 TRTA Xanomeli… time 180 n.risk Number o… 5 + 27 TRTA Xanomeli… time 180 estimate Survival… 0.619 + 28 TRTA Xanomeli… time 180 std.error Standard… 0.13 + 29 TRTA Xanomeli… time 180 conf.high CI Upper… 0.805 + 30 TRTA Xanomeli… time 180 conf.low CI Lower… 0.257 Message i 4 more variables: context, fmt_fn, warning, error @@ -98,15 +122,19 @@ status) ~ 1, data = survival::lung), times = c(60, 180)), stat = lapply(stat, function(x) ifelse(is.numeric(x), cards::round5(x, 3), x))), n = Inf) Message - {cards} data frame: 6 x 9 + {cards} data frame: 10 x 9 Output - variable variable_level context stat_name stat_label stat - 1 time 60 survival estimate Survival… 0.925 - 2 time 60 survival conf.high CI Upper… 0.96 - 3 time 60 survival conf.low CI Lower… 0.892 - 4 time 180 survival estimate Survival… 0.722 - 5 time 180 survival conf.high CI Upper… 0.783 - 6 time 180 survival conf.low CI Lower… 0.666 + variable variable_level context stat_name stat_label stat + 1 time 60 survival n.risk Number o… 213 + 2 time 60 survival estimate Survival… 0.925 + 3 time 60 survival std.error Standard… 0.017 + 4 time 60 survival conf.high CI Upper… 0.96 + 5 time 60 survival conf.low CI Lower… 0.892 + 6 time 180 survival n.risk Number o… 160 + 7 time 180 survival estimate Survival… 0.722 + 8 time 180 survival std.error Standard… 0.03 + 9 time 180 survival conf.high CI Upper… 0.783 + 10 time 180 survival conf.low CI Lower… 0.666 Message i 3 more variables: fmt_fn, warning, error @@ -146,20 +174,20 @@ 4 sex 1 ph.ecog 0 5 sex 1 ph.ecog 0 6 sex 1 ph.ecog 0 - 7 sex 1 ph.ecog 1 - 8 sex 1 ph.ecog 1 - 9 sex 1 ph.ecog 1 - 10 sex 1 ph.ecog 1 + 7 sex 1 ph.ecog 0 + 8 sex 1 ph.ecog 0 + 9 sex 1 ph.ecog 0 + 10 sex 1 ph.ecog 0 11 sex 1 ph.ecog 1 12 sex 1 ph.ecog 1 - 13 sex 1 ph.ecog 2 - 14 sex 1 ph.ecog 2 - 15 sex 1 ph.ecog 2 - 16 sex 1 ph.ecog 2 - 17 sex 1 ph.ecog 2 - 18 sex 1 ph.ecog 2 - 19 sex 1 ph.ecog 3 - 20 sex 1 ph.ecog 3 + 13 sex 1 ph.ecog 1 + 14 sex 1 ph.ecog 1 + 15 sex 1 ph.ecog 1 + 16 sex 1 ph.ecog 1 + 17 sex 1 ph.ecog 1 + 18 sex 1 ph.ecog 1 + 19 sex 1 ph.ecog 1 + 20 sex 1 ph.ecog 1 --- @@ -201,27 +229,39 @@ ifelse(is.numeric(x), cards::round5(x, 3), x))), n = Inf) Message Multi-state model detected. Showing probabilities into state 'death from cancer'. - {cards} data frame: 18 x 11 + {cards} data frame: 30 x 11 Output group1 group1_level variable variable_level stat_name stat_label stat - 1 TRTA Placebo time 60 estimate Survival… 0.054 - 2 TRTA Placebo time 60 conf.high CI Upper… 0.14 - 3 TRTA Placebo time 60 conf.low CI Lower… 0.021 - 4 TRTA Placebo time 180 estimate Survival… 0.226 - 5 TRTA Placebo time 180 conf.high CI Upper… 0.361 - 6 TRTA Placebo time 180 conf.low CI Lower… 0.142 - 7 TRTA Xanomeli… time 60 estimate Survival… 0.137 - 8 TRTA Xanomeli… time 60 conf.high CI Upper… 0.311 - 9 TRTA Xanomeli… time 60 conf.low CI Lower… 0.06 - 10 TRTA Xanomeli… time 180 estimate Survival… 0.51 - 11 TRTA Xanomeli… time 180 conf.high CI Upper… 0.892 - 12 TRTA Xanomeli… time 180 conf.low CI Lower… 0.292 - 13 TRTA Xanomeli… time 60 estimate Survival… 0.162 - 14 TRTA Xanomeli… time 60 conf.high CI Upper… 0.33 - 15 TRTA Xanomeli… time 60 conf.low CI Lower… 0.08 - 16 TRTA Xanomeli… time 180 estimate Survival… 0.244 - 17 TRTA Xanomeli… time 180 conf.high CI Upper… 0.516 - 18 TRTA Xanomeli… time 180 conf.low CI Lower… 0.115 + 1 TRTA Placebo time 60 n.risk Number o… 59 + 2 TRTA Placebo time 60 estimate Survival… 0.054 + 3 TRTA Placebo time 60 std.error Standard… 0.026 + 4 TRTA Placebo time 60 conf.high CI Upper… 0.14 + 5 TRTA Placebo time 60 conf.low CI Lower… 0.021 + 6 TRTA Placebo time 180 n.risk Number o… 35 + 7 TRTA Placebo time 180 estimate Survival… 0.226 + 8 TRTA Placebo time 180 std.error Standard… 0.054 + 9 TRTA Placebo time 180 conf.high CI Upper… 0.361 + 10 TRTA Placebo time 180 conf.low CI Lower… 0.142 + 11 TRTA Xanomeli… time 60 n.risk Number o… 14 + 12 TRTA Xanomeli… time 60 estimate Survival… 0.137 + 13 TRTA Xanomeli… time 60 std.error Standard… 0.057 + 14 TRTA Xanomeli… time 60 conf.high CI Upper… 0.311 + 15 TRTA Xanomeli… time 60 conf.low CI Lower… 0.06 + 16 TRTA Xanomeli… time 180 n.risk Number o… 3 + 17 TRTA Xanomeli… time 180 estimate Survival… 0.51 + 18 TRTA Xanomeli… time 180 std.error Standard… 0.145 + 19 TRTA Xanomeli… time 180 conf.high CI Upper… 0.892 + 20 TRTA Xanomeli… time 180 conf.low CI Lower… 0.292 + 21 TRTA Xanomeli… time 60 n.risk Number o… 20 + 22 TRTA Xanomeli… time 60 estimate Survival… 0.162 + 23 TRTA Xanomeli… time 60 std.error Standard… 0.059 + 24 TRTA Xanomeli… time 60 conf.high CI Upper… 0.33 + 25 TRTA Xanomeli… time 60 conf.low CI Lower… 0.08 + 26 TRTA Xanomeli… time 180 n.risk Number o… 5 + 27 TRTA Xanomeli… time 180 estimate Survival… 0.244 + 28 TRTA Xanomeli… time 180 std.error Standard… 0.093 + 29 TRTA Xanomeli… time 180 conf.high CI Upper… 0.516 + 30 TRTA Xanomeli… time 180 conf.low CI Lower… 0.115 Message i 4 more variables: context, fmt_fn, warning, error