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 far as I understand, they are both implementing the Aalen-Johansen estimator. My question is basically
- Am I wrong to expect that the two implementations should be identical?
- Do you know if the two implementations differ in a methodological perspective?
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
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