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


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

df <- data |> 
  enframe("status", "time") |> 
  unnest(time) |> 
    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() +

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

  • $\begingroup$ 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. $\endgroup$
    – Peter H.
    Jan 19 at 17:28

1 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 

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 enter image description here

  • $\begingroup$ Ah! That explains it - good catch. $\endgroup$
    – Peter H.
    Jan 20 at 7:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.