Return correct standard errors in `tidy.survfit()`

[mccarthy-m-g]

tidy.survfit() currently returns standard errors for the cumulative hazard instead of the survival probability. This PR fixes this so that standard errors for the survival probability are returned. This makes the output consistent with summary.survfit.

Here's a blog post covering the sneaky behaviour of standard errors in survival. The short version is that fit$std.err (which tidy.survfit() currently uses) returns the standard errors for the cumulative hazard, whereas summary(fit)$std.err returns the standard errors for the survival probabilities.

library(survival)
library(broom)

fit <- survfit(Surv(time, status) ~ x, data = aml)

# Cumulative hazard
fit$std.err
#>  [1] 0.09534626 0.14213381 0.19508758 0.24873417 0.24873417 0.33446777
#>  [7] 0.44181673 0.44181673 0.83378775 0.83378775 0.12909944 0.20412415
#> [13] 0.24397502 0.24397502 0.30472470 0.37796447 0.47559487 0.62678317
#> [19] 0.94491118        Inf

# Survival probability
summary(fit)$std.err
#>  [1] 0.08667842 0.11629130 0.13966497 0.15263233 0.16419327 0.16266889
#>  [7] 0.15349275 0.10758287 0.13608276 0.14231876 0.14813006 0.14698618
#> [13] 0.13871517 0.12187451 0.09186636        NaN

^{Created on 2023-06-03 with reprex v2.0.2}

Edit

After using this on some examples, I think it's better to calculate manually instead of using summary(); otherwise you get an error when using survfit objects that have been reformulated with survfit0() due to incompatible numbers of rows.

I've updated the PR to use the calculations done in the survival package instead of summary() like in the reprex above.

For reference, here are the relevant line in survival's source code.

survfit:

https://github.com/therneau/survival/blob/e8de7b732daa9a82c1ceaad93fd68ee7545c7736/R/survfitms.R#L178

survfitms:

https://github.com/therneau/survival/blob/e8de7b732daa9a82c1ceaad93fd68ee7545c7736/R/survfitms.R#L375

[simonpcouch]

Thanks for the PR! I'm seeing

── Error ('test-survival-survfit.R:30:3'): tidy.survfit ────────────────────────
Error in `data.frame(time = x$time, n.risk = c(x$n.risk), n.event = c(x$n.event), 
    n.censor = c(x$n.censor), estimate = c(x$pstate), std.error = c(x$std.err * 
        x$surv), conf.high = c(x$upper), conf.low = c(x$lower), 
    state = rep(x$states, each = nrow(x$pstate)))`: arguments imply differing number of rows: 237, 711, 0

in tests. Looks like this may be related to 03a4662—could you address that failure?

[mccarthy-m-g]

Welcome!

I addressed the failure, and also made this a bit more robust by using the same conditional statement as survival does to check if the standard errors need to be fixed. Should be good now.

[simonpcouch]

Alright, got it. Thanks again for putting this together!

After sitting with this for a bit, I think our best move here is to better document the current output. While the current output takes a moment of pause to interface with (and is certainly documented incorrectly), it is true to survival's implementation, and transitioning to a mix-and-match of slots of both the original object and its summary feels like a lateral move in terms of clarity/safety. At another point in the package's lifecycle, I may have considered adding a tidy.summary.survfit() method here, though that change is at odds with the notes in release 0.7.0 and the maintenance guidelines.

Would you be up for updating that documentation?

[mccarthy-m-g]

I'm up for updating the documentation. You're right that that would be more consistent with the information stored in survfit.object. For reference, here's how ?survfit.object defines std.err:

For a survival curve this contains standard error of the cumulative hazard or -log(survival), for a multi-state curve it contains the standard error of prev. This difference is a reflection of the fact that each is the natural calculation for that case.

However, I think this will lead to real errors in practice, as people generally assume the standard error next to an estimate to correspond to that estimate, not a different estimate that isn't shown.

As an alternative, would it be possible to add some sort of type argument to the function to return estimates, standard errors (and intervals) for the different curves? This would also be true to survival's implementation, and would be clearer and safer than only updating documentation or taking the mix-and-match approach.

[simonpcouch]

A type argument that proxies whether output mirrors the object structure of the model object or its summary feels to me like a workaround to not just newly adding a tidier method for the summary object, and thus again a lateral move in terms of clarity. Let's stick with just improving documentation here.👍

[mccarthy-m-g]

Hi @simonpcouch, small update:

As of survival version 3.7-0, survfit objects now return survival estimates (surv) and their standard errors (std.err), as well as cumulative hazard estimates (cumhaz) and their standard errors (std.chaz).

library(survival)
lung_survfit_1 <- survfit(Surv(time, status) ~ 1, data = lung)
str(lung_survfit_1)

## List of 17
##  $ n        : int 228
##  $ time     : num [1:186] 5 11 12 13 15 26 30 31 53 54 ...
##  $ n.risk   : num [1:186] 228 227 224 223 221 220 219 218 217 215 ...
##  $ n.event  : num [1:186] 1 3 1 2 1 1 1 1 2 1 ...
##  $ n.censor : num [1:186] 0 0 0 0 0 0 0 0 0 0 ...
##  $ surv     : num [1:186] 0.996 0.982 0.978 0.969 0.965 ...
##  $ std.err  : num [1:186] 0.0044 0.00885 0.00992 0.01179 0.01263 ...
##  $ cumhaz   : num [1:186] 0.00439 0.0176 0.02207 0.03103 0.03556 ...
##  $ std.chaz : num [1:186] 0.00439 0.0088 0.00987 0.01173 0.01257 ...
##  $ type     : chr "right"
##  $ logse    : logi TRUE
##  $ conf.int : num 0.95
##  $ conf.type: chr "log"
##  $ lower    : num [1:186] 0.987 0.966 0.959 0.947 0.941 ...
##  $ upper    : num [1:186] 1 1 0.997 0.992 0.989 ...
##  $ t0       : num 0
##  $ call     : language survfit(formula = Surv(time, status) ~ 1, data = lung)
##  - attr(*, "class")= chr "survfit"

Would you be open to this change now that a type argument would stay true to the model object’s structure?

New proposal

Here’s my proposed update to tidy.survfit(), which adds a type argument to choose between the different curves. This update would be in line with suggestions from the "Create your own broom tidier methods" article:

Sometimes a model will have several different types of components. For example, in mixed models, there is different information associated with fixed effects and random effects. Since this information doesn’t have the same interpretation, it doesn’t make sense to summarize the fixed and random effects in the same table. In cases like this you should add an argument that allows the user to specify which type of information they want.

tidy.survfit <- function(x, type = c("surv", "cumhaz"), ...) {
  type <- rlang::arg_match(type)
  if (inherits(x, "survfitms")) names(x)[names(x) == "pstate"] <- "surv"
  # For multi-state models, c() coerces matrices to vectors.
  ret <- data.frame(
    time = x$time,
    n.risk = c(x$n.risk),
    n.event = c(x$n.event),
    n.censor = c(x$n.censor),
    estimate = switch(type, surv = c(x$surv), cumhaz = c(x$cumhaz)),
    std.error = switch(type, surv = c(x$std.err), cumhaz = c(x$std.chaz)),
    # Confidence intervals for cumulative hazard estimates are not included in
    # survfit objects, but can be estimated as -log(survival) of the opposite
    # interval.
    conf.low = switch(type, surv = c(x$lower), cumhaz = -log(c(x$upper))),
    conf.high = switch(type, surv = c(x$upper), cumhaz = -log(c(x$lower)))
  )
  if (inherits(x, "survfitms")) {
    ret$state <- rep(x$states, each = nrow(x$surv))
    ret <- ret[ret$state != "", ]
  }
  # For multi-state models, strata are automatically recycled.
  if (!is.null(x$strata)) {
    ret$strata <- rep(names(x$strata), times = x$strata)
  }
  tibble::as_tibble(ret)
}

Results for different models

Output from a basic model:

tidy.survfit(lung_survfit_1)

## # A tibble: 186 × 8
##     time n.risk n.event n.censor estimate std.error conf.low conf.high
##    <dbl>  <dbl>   <dbl>    <dbl>    <dbl>     <dbl>    <dbl>     <dbl>
##  1     5    228       1        0    0.996   0.00440    0.987     1    
##  2    11    227       3        0    0.982   0.00885    0.966     1.00 
##  3    12    224       1        0    0.978   0.00992    0.959     0.997
##  4    13    223       2        0    0.969   0.0118     0.947     0.992
##  5    15    221       1        0    0.965   0.0126     0.941     0.989
##  6    26    220       1        0    0.961   0.0134     0.936     0.986
##  7    30    219       1        0    0.956   0.0142     0.930     0.983
##  8    31    218       1        0    0.952   0.0149     0.924     0.980
##  9    53    217       2        0    0.943   0.0163     0.913     0.974
## 10    54    215       1        0    0.939   0.0169     0.908     0.970
## # ℹ 176 more rows

tidy.survfit(lung_survfit_1, type = "cumhaz")

## # A tibble: 186 × 8
##     time n.risk n.event n.censor estimate std.error conf.low conf.high
##    <dbl>  <dbl>   <dbl>    <dbl>    <dbl>     <dbl>    <dbl>     <dbl>
##  1     5    228       1        0  0.00439   0.00439 0           0.0130
##  2    11    227       3        0  0.0176    0.00880 0.000354    0.0350
##  3    12    224       1        0  0.0221    0.00987 0.00274     0.0416
##  4    13    223       2        0  0.0310    0.0117  0.00808     0.0543
##  5    15    221       1        0  0.0356    0.0126  0.0110      0.0605
##  6    26    220       1        0  0.0401    0.0134  0.0140      0.0666
##  7    30    219       1        0  0.0447    0.0141  0.0171      0.0727
##  8    31    218       1        0  0.0493    0.0149  0.0202      0.0787
##  9    53    217       2        0  0.0585    0.0162  0.0268      0.0906
## 10    54    215       1        0  0.0631    0.0169  0.0302      0.0966
## # ℹ 176 more rows

A model with strata:

lung_survfit_2 <- update(lung_survfit_1, . ~ sex)
tidy.survfit(lung_survfit_2)

## # A tibble: 206 × 9
##     time n.risk n.event n.censor estimate std.error conf.low conf.high strata
##    <dbl>  <dbl>   <dbl>    <dbl>    <dbl>     <dbl>    <dbl>     <dbl> <chr> 
##  1    11    138       3        0    0.978    0.0127    0.954     1     sex=1 
##  2    12    135       1        0    0.971    0.0147    0.943     0.999 sex=1 
##  3    13    134       2        0    0.957    0.0181    0.923     0.991 sex=1 
##  4    15    132       1        0    0.949    0.0197    0.913     0.987 sex=1 
##  5    26    131       1        0    0.942    0.0211    0.904     0.982 sex=1 
##  6    30    130       1        0    0.935    0.0225    0.894     0.977 sex=1 
##  7    31    129       1        0    0.928    0.0238    0.885     0.972 sex=1 
##  8    53    128       2        0    0.913    0.0263    0.867     0.961 sex=1 
##  9    54    126       1        0    0.906    0.0275    0.858     0.956 sex=1 
## 10    59    125       1        0    0.899    0.0286    0.850     0.950 sex=1 
## # ℹ 196 more rows

tidy.survfit(lung_survfit_2, type = "cumhaz")

## # A tibble: 206 × 9
##     time n.risk n.event n.censor estimate std.error conf.low conf.high strata
##    <dbl>  <dbl>   <dbl>    <dbl>    <dbl>     <dbl>    <dbl>     <dbl> <chr> 
##  1    11    138       3        0   0.0217    0.0126 0           0.0469 sex=1 
##  2    12    135       1        0   0.0291    0.0146 0.000588    0.0582 sex=1 
##  3    13    134       2        0   0.0441    0.0180 0.00888     0.0800 sex=1 
##  4    15    132       1        0   0.0516    0.0195 0.0135      0.0906 sex=1 
##  5    26    131       1        0   0.0593    0.0210 0.0183      0.101  sex=1 
##  6    30    130       1        0   0.0670    0.0223 0.0234      0.112  sex=1 
##  7    31    129       1        0   0.0747    0.0236 0.0286      0.122  sex=1 
##  8    53    128       2        0   0.0904    0.0261 0.0395      0.142  sex=1 
##  9    54    126       1        0   0.0983    0.0273 0.0451      0.153  sex=1 
## 10    59    125       1        0   0.106     0.0284 0.0509      0.163  sex=1 
## # ℹ 196 more rows

A Cox model with strata:

lung_coxph <- coxph(Surv(time, status) ~ strata(sex), data = lung)
lung_coxsurv <- survfit(lung_coxph)
tidy.survfit(lung_coxsurv)

## # A tibble: 206 × 9
##     time n.risk n.event n.censor estimate std.error conf.low conf.high strata
##    <dbl>  <dbl>   <dbl>    <dbl>    <dbl>     <dbl>    <dbl>     <dbl> <chr> 
##  1    11    138       3        0    0.978    0.0126    0.954     1     sex=1 
##  2    12    135       1        0    0.971    0.0147    0.944     0.999 sex=1 
##  3    13    134       2        0    0.957    0.0181    0.923     0.991 sex=1 
##  4    15    132       1        0    0.949    0.0196    0.914     0.987 sex=1 
##  5    26    131       1        0    0.942    0.0210    0.904     0.982 sex=1 
##  6    30    130       1        0    0.935    0.0224    0.895     0.977 sex=1 
##  7    31    129       1        0    0.928    0.0237    0.886     0.972 sex=1 
##  8    53    128       2        0    0.913    0.0262    0.868     0.961 sex=1 
##  9    54    126       1        0    0.906    0.0273    0.859     0.956 sex=1 
## 10    59    125       1        0    0.899    0.0285    0.850     0.951 sex=1 
## # ℹ 196 more rows

tidy.survfit(lung_coxsurv, type = "cumhaz")

## # A tibble: 206 × 9
##     time n.risk n.event n.censor estimate std.error conf.low conf.high strata
##    <dbl>  <dbl>   <dbl>    <dbl>    <dbl>     <dbl>    <dbl>     <dbl> <chr> 
##  1    11    138       3        0   0.0219    0.0126 0           0.0467 sex=1 
##  2    12    135       1        0   0.0293    0.0147 0.000586    0.0580 sex=1 
##  3    13    134       2        0   0.0443    0.0181 0.00885     0.0797 sex=1 
##  4    15    132       1        0   0.0519    0.0196 0.0134      0.0903 sex=1 
##  5    26    131       1        0   0.0595    0.0210 0.0183      0.101  sex=1 
##  6    30    130       1        0   0.0672    0.0224 0.0233      0.111  sex=1 
##  7    31    129       1        0   0.0749    0.0237 0.0285      0.121  sex=1 
##  8    53    128       2        0   0.0906    0.0262 0.0393      0.142  sex=1 
##  9    54    126       1        0   0.0986    0.0273 0.0450      0.152  sex=1 
## 10    59    125       1        0   0.107     0.0285 0.0507      0.162  sex=1 
## # ℹ 196 more rows

A multi-state model:

mgus_survfitms_1 <- survfit(Surv(start, stop, event) ~ 1, id = id, data = mgus1)
tidy.survfit(mgus_survfitms_1)

## # A tibble: 900 × 9
##     time n.risk n.event n.censor estimate std.error conf.low conf.high state
##    <dbl>  <dbl>   <dbl>    <dbl>    <dbl>     <dbl>    <dbl>     <dbl> <chr>
##  1     6    241       0        0    0.996   0.00414    0.988     1     (s0) 
##  2     7    240       0        0    0.992   0.00584    0.980     1     (s0) 
##  3    31    239       0        0    0.988   0.00714    0.974     1     (s0) 
##  4    32    238       0        0    0.983   0.00823    0.967     1.00  (s0) 
##  5    39    237       0        0    0.979   0.00918    0.961     0.997 (s0) 
##  6    60    236       0        0    0.975   0.0100     0.956     0.995 (s0) 
##  7    61    235       0        0    0.967   0.0115     0.944     0.990 (s0) 
##  8   152    233       0        0    0.963   0.0122     0.939     0.987 (s0) 
##  9   153    232       0        0    0.959   0.0128     0.934     0.984 (s0) 
## 10   174    231       0        0    0.954   0.0134     0.928     0.981 (s0) 
## # ℹ 890 more rows

tidy.survfit(mgus_survfitms_1, type = "cumhaz")

## # A tibble: 900 × 9
##     time n.risk n.event n.censor estimate std.error conf.low conf.high state
##    <dbl>  <dbl>   <dbl>    <dbl>    <dbl>     <dbl>    <dbl>     <dbl> <chr>
##  1     6    241       0        0        0         0 0           0.0123 (s0) 
##  2     7    240       0        0        0         0 0           0.0199 (s0) 
##  3    31    239       0        0        0         0 0           0.0267 (s0) 
##  4    32    238       0        0        0         0 0.000335    0.0331 (s0) 
##  5    39    237       0        0        0         0 0.00259     0.0393 (s0) 
##  6    60    236       0        0        0         0 0.00504     0.0454 (s0) 
##  7    61    235       0        0        0         0 0.0104      0.0572 (s0) 
##  8   152    233       0        0        0         0 0.0132      0.0629 (s0) 
##  9   153    232       0        0        0         0 0.0161      0.0686 (s0) 
## 10   174    231       0        0        0         0 0.0191      0.0743 (s0) 
## # ℹ 890 more rows

A multi-state model with strata:

mgus_survfitms_2 <- update(mgus_survfitms_1, . ~ sex)
tidy.survfit(mgus_survfitms_2)

## # A tibble: 903 × 10
##     time n.risk n.event n.censor estimate std.error conf.low conf.high state strata    
##    <dbl>  <dbl>   <dbl>    <dbl>    <dbl>     <dbl>    <dbl>     <dbl> <chr> <chr>     
##  1    61    104       0        0    0.990   0.00957    0.972     1     (s0)  sex=female
##  2   152    103       0        0    0.981   0.0135     0.955     1     (s0)  sex=female
##  3   370    102       0        0    0.971   0.0164     0.940     1     (s0)  sex=female
##  4   499    101       0        0    0.962   0.0189     0.925     0.999 (s0)  sex=female
##  5   518    100       0        0    0.952   0.0210     0.912     0.994 (s0)  sex=female
##  6   533     99       0        0    0.942   0.0229     0.899     0.988 (s0)  sex=female
##  7   652     98       0        0    0.933   0.0246     0.886     0.982 (s0)  sex=female
##  8   700     97       0        0    0.923   0.0261     0.873     0.976 (s0)  sex=female
##  9   748     96       0        0    0.913   0.0276     0.861     0.969 (s0)  sex=female
## 10   954     95       0        0    0.904   0.0289     0.849     0.962 (s0)  sex=female
## # ℹ 893 more rows

tidy.survfit(mgus_survfitms_2, type = "cumhaz")

## # A tibble: 903 × 10
##     time n.risk n.event n.censor estimate std.error conf.low conf.high state strata    
##    <dbl>  <dbl>   <dbl>    <dbl>    <dbl>     <dbl>    <dbl>     <dbl> <chr> <chr>     
##  1    61    104       0        0   0         0      0           0.0286 (s0)  sex=female
##  2   152    103       0        0   0         0      0           0.0463 (s0)  sex=female
##  3   370    102       0        0   0         0      0           0.0624 (s0)  sex=female
##  4   499    101       0        0   0         0      0.000783    0.0777 (s0)  sex=female
##  5   518    100       0        0   0         0      0.00608     0.0925 (s0)  sex=female
##  6   533     99       0        0   0         0      0.0119      0.107  (s0)  sex=female
##  7   652     98       0        0   0         0      0.0181      0.121  (s0)  sex=female
##  8   700     97       0        0   0.0103    0.0103 0.0246      0.136  (s0)  sex=female
##  9   748     96       0        0   0.0103    0.0103 0.0314      0.150  (s0)  sex=female
## 10   954     95       0        0   0.0208    0.0147 0.0384      0.164  (s0)  sex=female
## # ℹ 893 more rows