tagteam / prodlim

Product limit estimation for survival analysis
7 stars 9 forks source link

Different cumulative incidence estimates obtained with `prodlim` and `survfit` for a simple competing risks dataset #3

Closed peterchristofferholm closed 1 year ago

peterchristofferholm commented 1 year ago

Thanks for a very useful package.

I have a question regarding the cumulative incidence curves obtained using prodlim, which seems to differ quite a lot compared to the ones that I can obtain from the survfit function in the survival package. In a relatively basic competing risks setup, I would expect the two implementations to give more similar results.

As a reprex, we can use some example data from “Kleinbaum, David G., and Mitchel Klein. Survival analysis: a self-learning text. Vol. 3. New York: Springer, 2012.” which looks like this:

data <- list(
  cens = c(3.2, 7.6, 10, 11, 15, 24.4),
  cause_1 = c(0.7, 3, 4.9, 6, 6, 6.9, 10, 10.8, 17.1, 20.3),
  cause_2 = c(1.5, 2.8, 3.8, 4.7, 7, 10, 10, 11.2)
)

For reference, I'm using the following libraries

library(tidyverse)
library(survival)
library(prodlim)

First I wrangle the data into a tabular format that can be used by prodlim and survfit

df <- data |> 
  enframe("status", "time") |> 
  unnest(time) |> 
  mutate(
    status = factor(status, levels = c("cens", "cause_1", "cause_2"))
  )

And then, fitting and tidying the results from survfit

sfit <- survfit(Surv(time, status) ~ 1, data = df)
pstate <- summary(sfit)$pstate
colnames(pstate) <- c("cens", "cause_1", "cause_2")
sfit_tidy <- as_tibble(pstate) |> 
  add_column(time = summary(sfit)$time) |> 
  pivot_longer(-time, names_to = "cause", values_to = "cuminc")

... and similarly for prodlim

pfit <- prodlim(Hist(time, status) ~ 1, data = df)
pfit_tidy <- map_dfr(summary(pfit)$table, as_tibble, .id = "cause") |> 
  select(cause, time, cuminc)

If we compare the cumulative incidence curves they differ quite alot. As far as I understand, the two methods both use the Aalen-Johansen estimator and I have a hard time understanding why the shouldn’t be identical.

bind_rows(sfit = sfit_tidy, pfit = pfit_tidy, .id = "method") |> 
  filter(cause != "cens") |> 
  ggplot(aes(time, cuminc, color = method)) +
  geom_step() +
  facet_wrap(~cause)

Created on 2023-01-19 with reprex v2.0.2

I can also add, that the cumulative incidence from survfit corresponds exactly to the ones that are presented in the Kleinbaum & Klein textbook where they "manually" go through the calculation of cumulative incidence curves.

Kind regards, Peter

peterchristofferholm commented 1 year ago

Turns out i have to explicitly state which level corresponds to censoring, and then the two approaches agree.

https://stats.stackexchange.com/questions/602547/different-cumulative-incidence-estimates-obtained-with-prodlim-and-survfit-f/602572#602572