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

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

• Please let me know if this is off-topic for Cross-Validated and instead should be posted to stackoverflow. My impression was that more people here would be familiar with the packages in question. Jan 19 at 17:28

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

• Ah! That explains it - good catch. Jan 20 at 7:44