adding ard_survival_survfit_diff() function (#136)

**What changes are proposed in this pull request?**
* Added the `ard_survival_survfit_diff()` function for calculating
differences in survival estimates. (#126, @<username>)


**Reference GitHub issue associated with pull request.** _e.g., 'closes
#<issue number>'_
closes #126


--------------------------------------------------------------------------------

Pre-review Checklist (if item does not apply, mark is as complete)
- [x] **All** GitHub Action workflows pass with a ✅
- [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)

- [x] If a bug was fixed, a unit test was added.
- [x] 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 ✅
- [ ] Approve Pull Request
- [ ] Merge the PR. Please use "Squash and merge" or "Rebase and merge".

---------

Signed-off-by: Daniel Sjoberg <[email protected]>
Co-authored-by: Emily de la Rua <[email protected]>

NAMESPACE

@@ -48,6 +48,7 @@ export(ard_survey_svyranktest)
 export(ard_survey_svyttest)
 export(ard_survival_survdiff)
 export(ard_survival_survfit)
+export(ard_survival_survfit_diff)
 export(bt)
 export(bt_strip)
 export(construct_model)

R/ard_stats_t_test_onesample.R

@@ -64,7 +64,7 @@ ard_stats_t_test_onesample <- function(data, variables, by = dplyr::group_vars(d
     ) |>
     dplyr::mutate(
       stat_label = dplyr::coalesce(.data$stat_label, .data$stat_name),
-      context = "ard_stats_t_test_onesample",
+      context = "stats_t_test_onesample",
     ) |>
     cards::tidy_ard_row_order() |>
     cards::tidy_ard_column_order()

R/ard_stats_wilcox_test_onesample.R

@@ -65,7 +65,7 @@ ard_stats_wilcox_test_onesample <- function(data, variables, by = dplyr::group_v
     ) |>
     dplyr::mutate(
       stat_label = dplyr::coalesce(.data$stat_label, .data$stat_name),
-      context = "ard_stats_wilcox_test_onesample",
+      context = "stats_wilcox_test_onesample",
     ) |>
     cards::tidy_ard_row_order() |>
     cards::tidy_ard_column_order()

R/ard_survival_survfit_diff.R

@@ -0,0 +1,121 @@
+#' ARD Survival Differences
+#'
+#' Calculate differences in the Kaplan-Meier estimator of survival using the
+#' results from [`survival::survfit()`].
+#'
+#' @param x (`survift`)\cr
+#'   object of class `'survfit'` typically created with [`survival::survfit()`]
+#' @param conf.level (scalar `numeric`)\cr
+#'   confidence level for confidence interval. Default is `0.95`.
+#' @inheritParams ard_survival_survfit
+#'
+#' @return an ARD data frame of class 'card'
+#' @export
+#'
+#' @examplesIf do.call(asNamespace("cardx")$is_pkg_installed, list(pkg = c("survival", "ggsurvfit"), reference_pkg = "cardx"))
+#' library(ggsurvfit)
+#' library(survival)
+#'
+#' survfit(Surv_CNSR() ~ TRTA, data = cards::ADTTE) |>
+#'   ard_survival_survfit_diff(times = c(25, 50))
+ard_survival_survfit_diff <- function(x, times, conf.level = 0.95) {
+  set_cli_abort_call()
+
+  # check installed packages ---------------------------------------------------
+  check_pkg_installed(c("survival", "broom"), reference_pkg = "cardx")
+  check_not_missing(x)
+  check_not_missing(times)
+  check_class(x, "survfit")
+
+  if (inherits(x, c("survfitms", "survfitcox"))) {
+    cli::cli_abort(
+      "Argument {.arg x} cannot be class {.cls {c('survfitms', 'survfitcox')}}.",
+      call = get_cli_abort_call()
+    )
+  }
+  check_scalar_range(conf.level, range = c(0, 1))
+  check_length(
+    as.list(x$call)[["formula"]] |> stats::as.formula() |> stats::terms() |> attr("term.labels"),
+    length = 1L,
+    message = "The {.cls survfit} object passed in argument {.arg x} must be stratified by a single variable."
+  )
+  if (length(x$strata) < 2) {
+    cli::cli_abort(
+      "The {.cls survfit} object's stratifying variable must have 2 or more levels.",
+      call = get_cli_abort_call()
+    )
+  }
+
+  # calculate the survival at the specified times
+  ard_survival_survfit <-
+    ard_survival_survfit(x = x, times = times) |>
+    dplyr::filter(.data$stat_name %in% c("estimate", "std.error")) |>
+    dplyr::select(-c("stat_label", "context", "fmt_fn"))
+
+  # transform the survival ARD into a cards object with the survival difference
+  card <-
+    ard_survival_survfit %>%
+    {dplyr::left_join( # styler: off
+        # remove the first group from the data frame (this is our reference group)
+        dplyr::filter(., .by = cards::all_ard_groups(), dplyr::cur_group_id() > 1L) |>
+          dplyr::rename(stat1 = "stat"),
+        # merge the reference group data
+        dplyr::filter(., .by = cards::all_ard_groups(), dplyr::cur_group_id() == 1L) |>
+          dplyr::select(stat0 = "stat", everything(), -c("group1_level", "error", "warning")),
+        by = c("group1", "variable", "variable_level", "stat_name")
+    )} |> # styler: off
+    # reshape to put the stats that need to be combined on the same row
+    tidyr::pivot_wider(
+      id_cols = c("group1", "group1_level", "variable", "variable_level"),
+      names_from = "stat_name",
+      values_from = c("stat0", "stat1"),
+      values_fn = unlist
+    ) |>
+    # calcualte the primary statistics to return
+    dplyr::mutate(
+      # reference level
+      reference_level = ard_survival_survfit[["group1_level"]][1],
+      # short description of method
+      method = "Survival Difference (Z-test)",
+      # survival difference
+      estimate = .data$stat0_estimate - .data$stat1_estimate,
+      # survival difference standard error
+      std.error = sqrt(.data$stat0_std.error^2 + .data$stat1_std.error^2),
+      # Z test statistic
+      statistic = .data$estimate / .data$std.error,
+      # confidence limits of the survival difference
+      conf.low = .data$estimate - .data$std.error * stats::qnorm(1 - (1 - .env$conf.level) / 2),
+      conf.high = .data$estimate + .data$std.error * stats::qnorm(1 - (1 - .env$conf.level) / 2),
+      # p-value for test where H0: no difference
+      p.value = 2 * (1 - stats::pnorm(abs(.data$statistic))),
+      across(c("reference_level", "estimate", "std.error", "statistic", "conf.low", "conf.high", "p.value", "method"), as.list)
+    ) |>
+    # reshape into the cards structure
+    dplyr::select(-starts_with("stat0_"), -starts_with("stat1_")) |>
+    tidyr::pivot_longer(
+      cols = -c(cards::all_ard_groups(), cards::all_ard_variables()),
+      names_to = "stat_name",
+      values_to = "stat"
+    )
+
+  # final prepping of the cards object -----------------------------------------
+  card |>
+    dplyr::mutate(
+      warning = ard_survival_survfit[["warning"]][1],
+      error = ard_survival_survfit[["error"]][1],
+      fmt_fn = list(1L),
+      stat_label =
+        dplyr::case_when(
+          .data$stat_name %in% "estimate" ~ "Survival Difference",
+          .data$stat_name %in% "std.error" ~ "Survival Difference Standard Error",
+          .data$stat_name %in% "conf.low" ~ "CI Lower Bound",
+          .data$stat_name %in% "conf.high" ~ "CI Upper Bound",
+          .data$stat_name %in% "statistic" ~ "z statistic",
+          .data$stat_name %in% "p.value" ~ "p-value",
+          .default = .data$stat_name
+        ),
+      context = "survival_survfit_diff",
+    ) |>
+    cards::tidy_ard_column_order() %>%
+    structure(., class = c("card", class(.)))
+}

R/construction_helpers.R

@@ -25,6 +25,11 @@
 #'   must be specified in the `package` argument.
 #' @param method.args (named `list`)\cr
 #'   named list of arguments that will be passed to `fn`.
+#'
+#'   Note that this list may contain non-standard evaluation components.
+#'   If you are wrapping this function in other functions, the argument
+#'   must be passed in a way that does not evaluate the list, e.g.
+#'   using rlang's  embrace operator `{{ . }}`.
 #' @param package (`string`)\cr
 #'   string of package name that will be temporarily loaded when function
 #'   specified in `method` is executed.

_pkgdown.yml

@@ -71,6 +71,7 @@ reference:
   - subtitle: "{survival} package"
   - contents:
     - ard_survival_survfit
+    - ard_survival_survfit_diff
     - ard_survival_survdiff
 
   - subtitle: "Other ARD functions"

inst/WORDLIST

@@ -34,6 +34,7 @@ mmrm
 pearson
 pre
 quosures
+rlang's
 sd
 strat
 vif

tests/testthat/_snaps/ard_attributes.survey.design.md

@@ -0,0 +1,13 @@
+# ard_attributes.survey.design() works
+
+    Code
+      attr(dclus1$variables$sname, "label") <- "School Name"
+      as.data.frame(ard_attributes(dclus1, variables = c(sname, dname), label = list(
+        dname = "District Name")))
+    Output
+        variable    context stat_name     stat_label          stat
+      1    sname attributes     label Variable Label   School Name
+      2    sname attributes     class Variable Class     character
+      3    dname attributes     label Variable Label District Name
+      4    dname attributes     class Variable Class     character
+

tests/testthat/_snaps/ard_survival_survfit_diff.md

@@ -0,0 +1,27 @@
+# ard_survival_survfit_diff() messaging
+
+    Code
+      ard_survival_survfit_diff(survfit(Surv(AVAL, 1 - CNSR) ~ SEX + TRTA, cards::ADTTE),
+      times = c(25, 50))
+    Condition
+      Error in `ard_survival_survfit_diff()`:
+      ! The <survfit> object passed in argument `x` must be stratified by a single variable.
+
+---
+
+    Code
+      ard_survival_survfit_diff(survfit(Surv(AVAL, 1 - CNSR) ~ constant, dplyr::mutate(
+        cards::ADTTE, constant = 1L)), times = c(25, 50))
+    Condition
+      Error in `ard_survival_survfit_diff()`:
+      ! The <survfit> object's stratifying variable must have 2 or more levels.
+
+---
+
+    Code
+      ard_survival_survfit_diff(survfit(coxph(Surv(AVAL, CNSR) ~ SEX + strata(TRTA),
+      cards::ADTTE)), times = c(25, 50))
+    Condition
+      Error in `ard_survival_survfit_diff()`:
+      ! Argument `x` cannot be class <survfitms/survfitcox>.
+

tests/testthat/test-ard_survival_survfit_diff.R

@@ -0,0 +1,63 @@
+skip_if_not(is_pkg_installed(c("survival", "broom"), reference_pkg = "cardx"))
+
+test_that("ard_survival_survfit_diff() works", {
+  withr::local_package("survival")
+  sf <- survfit(Surv(AVAL, 1 - CNSR) ~ SEX, cards::ADTTE)
+  expect_silent(
+    ard1 <- ard_survival_survfit_diff(sf, times = c(25, 50))
+  )
+
+  # check the survival differences are accurate
+  expect_equal(
+    ard1 |>
+      dplyr::filter(variable_level == 25, stat_name == "estimate") |>
+      dplyr::pull(stat) |>
+      unlist(),
+    summary(sf, times = 25) |>
+      getElement("surv") |>
+      reduce(`-`)
+  )
+  expect_equal(
+    ard1 |>
+      dplyr::filter(variable_level == 50, stat_name == "estimate") |>
+      dplyr::pull(stat) |>
+      unlist(),
+    summary(sf, times = 50) |>
+      getElement("surv") |>
+      reduce(`-`)
+  )
+
+  # check the structure of the ARD object
+  expect_silent(
+    cards::check_ard_structure(ard1)
+  )
+})
+
+test_that("ard_survival_survfit_diff() messaging", {
+  withr::local_package("survival")
+
+  # we can only do one stratifying variable at a time
+  expect_snapshot(
+    error = TRUE,
+    survfit(Surv(AVAL, 1 - CNSR) ~ SEX + TRTA, cards::ADTTE) |>
+      ard_survival_survfit_diff(times = c(25, 50))
+  )
+
+  # the stratifying variable must have 2 or more levels
+  expect_snapshot(
+    error = TRUE,
+    survfit(
+      Surv(AVAL, 1 - CNSR) ~ constant,
+      cards::ADTTE |> dplyr::mutate(constant = 1L)
+    ) |>
+      ard_survival_survfit_diff(times = c(25, 50))
+  )
+
+  # cannot pass a multi-state model or stratified Cox
+  expect_snapshot(
+    error = TRUE,
+    coxph(Surv(AVAL, CNSR) ~ SEX + strata(TRTA), cards::ADTTE) |>
+      survfit() |>
+      ard_survival_survfit_diff(times = c(25, 50))
+  )
+})