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

Question

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 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}

Answer 1

This could easily have been a CrossValidated question, but it turns out to be more a StackOverflow one.

If you look at pfit

> pfit

Call: prodlim(formula = Hist(time, status) ~ 1, data = df)


No covariates

Uncensored response of a competing.risks model

No.Observations: 24 

Pattern:
         
Cause     event
  cens        6
  cause_1    10
  cause_2     8
  unknown     0

So the problem is with the censoring (what I noticed is that the curves start to separate after the first censoring). In survfit() the first level of a factor is automatically taken as censoring, but in Hist() there's a cens.code argument to set the censoring level. If you do

pfit <- prodlim(Hist(time, status,cens.code="cens") ~ 1, data = df)

you get