Bug in `tidy.survfit()` coming in the next release of the {survival} pkg

[ddsjoberg]

Hello, The authors of the {survival} package are prepping a release to CRAN and plan to submit in two weeks. The tidy.survfit() function does not error with the new version; however, the results are not correct for multi-state models (the most common type is competing risks).

I haven't looked into the details too much, but it seems that in the new version the number of censored subjects is stored as a matrix, where previously it was a vector.

I prepared an example illustrating the differences in the result from the current release of CRAN to the dev version.
https://rpubs.com/ddsjoberg/1174489

Below is a table showing the difference in the reported number of censored subjects from the current CRAN release (n.censor.cran) compared to dev (n.censor.dev). The dev version of the package is here www.github.com/therneau/survival

## # A tibble: 58 × 4
##    state         time n.censor.cran n.censor.dev
##    <chr>        <dbl>         <int>        <dbl>
##  1 cancer death 0.130             1            0
##  2 cancer death 0.153             1            0
##  3 cancer death 0.192             1            0
##  4 cancer death 0.192             1            0
##  5 cancer death 0.202             1            0
##  6 cancer death 0.211             1            0
##  7 cancer death 0.269             1            0
##  8 cancer death 0.288             1            0
##  9 cancer death 0.292             1            0
## 10 cancer death 0.343             1            0
## # ℹ 48 more rows

I would submit a PR myself, but I can't this week and I am traveling next week. I hope someone from the community has the bandwidth to address this important issue for an incredibly common package and method.

Thanks,
Daniel

[simonpcouch]

Thank you for the issue, @ddsjoberg! Just wanted to affirm that we've seen this issue and we're on it. :)

[ddsjoberg]

Thank you for the issue, @ddsjoberg! Just wanted to affirm that we've seen this issue and we're on it. :)

Awesome, thank you @simonpcouch !! 🎸

[simonpcouch]

Okay, I believe that this is an artifact of an "inconsistency" in survival that is now fixed. It looks like broom had surfaced this issue some years ago when transitioning to tibble but reintroduced the issue by reverting to data.frame()'s recycling behavior. The new survival release will fix the issue—as in, broom output has propagated the "inconsistent" result from survival for some time, and will now be correct. Once the new survival release is out, I'll switch back to tibble under the hood here so we're more vigilant about vector recycling there.

Some more context:

library(broom)

path <- fs::path_temp()
path_cran <- fs::path(path, "cran")
path_dev <- fs::path(path, "dev")
fs::dir_create(path_cran)
fs::dir_create(path_dev)

install.packages("survival", lib = path_cran)
#> 
#> The downloaded binary packages are in
#>  /var/folders/6c/w21prsj167b_x82q4_s45t340000gn/T//RtmpKmSXN4/downloaded_packages
pak::pak("therneau/survival", lib = path_dev)
#> ℹ Loading metadata database
#> ✔ Loading metadata database ... done
#> 
#> 
#> → Will update 1 package.
#> → The package (0 B) is cached.
#> + survival 3.5-8 → 3.6-1 👷🏿‍♂️🔧 (GitHub: f6016a7)
#> ℹ No downloads are needed, 1 pkg is cached
#> ✔ Installed survival 3.6-1 (github::therneau/survival@f6016a7) (26ms)
#> ✔ 1 pkg + 2 deps: upd 1 [3.7s]

library("survival", lib.loc = path_cran)
packageVersion("survival")
#> [1] '3.5.8'

set.seed(1)
n <- 100
d <- 
  dplyr::tibble(
    death_cr = 
      sample.int(n = 3, size = n, replace = TRUE) |> 
      factor(labels = c("censor", "cancer death", "other cause death")),
    ttdeath = rgamma(n = n, shape = 1)
  )

set.seed(1)
cran_fit <- survfit(Surv(ttdeath, death_cr) ~ 1, data = d)
cran_tidy <- cran_fit |> tidy()

detach("package:survival", unload = TRUE)
library("survival", lib.loc = path_dev)
packageVersion("survival")
#> [1] '3.6.1'

set.seed(1)
dev_fit <- survfit(Surv(ttdeath, death_cr) ~ 1, data = d)
dev_tidy <- dev_fit |> tidy()

dim(cran_fit$n.censor)
#> NULL
length(cran_fit$n.censor)
#> [1] 100
dim(dev_fit$n.censor)
#> [1] 100   3
length(dev_fit$n.censor)
#> [1] 300

# note that the rest of the recycled components are 
# in cran survival are actually arrays
dim(cran_fit$n.risk)
#> [1] 100   3
dim(cran_fit$n.event)
#> [1] 100   3

^{Created on 2024-04-17 with reprex v2.1.0}

The survfit tidiers currently just rely on data.frame()'s vector recycling to gloss over that inconsistency. With tibble() switched out under the hood, we'd have seen:

cran_tidy <- cran_fit |> tidy()
#> Error in `tibble::tibble()`:
#> ! Tibble columns must have compatible sizes.
#> • Size 300: Existing data.
#> • Size 100: Column `n.censor`.
#> ℹ Only values of size one are recycled.

[ddsjoberg]

If it helps, the author of {survival} also added, "For a multi-state, the survfit routine now returns n.censor as a matrix, with one column per state." Competing risks falls under the umbrella of multi-state models.

[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.