I am trying to estimate survival function, but in case where each event is censoring with probability $p_i$. (That is, I am never sure if the event is right-censoring or death, but I can estimate the probability of each.)

Typical Kaplan-Meier works only for binary right-censoring. Is there any straightforward generalization?

  • $\begingroup$ Out of interest, do you have an application in mind? I'm struggling to think of a situation where you would know $t_i$ but not whether the event is right censoring or death. $\endgroup$
    – tristan
    Apr 8, 2015 at 22:33
  • $\begingroup$ @tristan I am looking at user churn rate. But what I have are points in time (e.g. if transactions), rather than anything continuous. So I need to estimate, whether the time between the last interaction and end of measurement corresponds to a typical time between transactions, or did user resigned. $\endgroup$ Apr 9, 2015 at 7:12
  • $\begingroup$ So how have you come up with the $p_i$ in this case? It sounds to me that all of the users are right-censored in fact. $\endgroup$
    – tristan
    Apr 9, 2015 at 14:37
  • 1
    $\begingroup$ @PiotrMigdal have you looked at using customer CLV to measure this? See a good intro in "“Counting Your Customers” the Easy Way: An Alternative to the Pareto/NBD Model" by Fader et al. An implementation is available at github.com/CamDavidsonPilon/lifetimes $\endgroup$ Apr 17, 2015 at 12:30
  • $\begingroup$ @Cam.Davidson.Pilon Thanks! Yesterday I learnt about R Buy 'Til You Die, but I prefer Python (and I am already using Lifelines). $\endgroup$ Apr 17, 2015 at 17:19

1 Answer 1


If we can assume that for each item/individual you have $ x_i $ (covariates), $ t_i $ (last observation) and $ p_i $ (probability of death rather than censoring), and that the $ p_i $ are conditionally independent given $ x_i $ and $ t_i $: I don't see a problem with a simple modification to the Kaplan-Meier product limit formula.

If the standard formula is:

$$ \hat S (t) = \prod_{t_i<t}{\frac{n_i-d_i}{n_i}} $$

Where $ n_i $ is the number of individuals at risk at time $ t_i $ and $ d_i $ is the number dying at time $ t_i $.

Then you should be able to use

$$ \hat S (t) = \prod_{t_i <t}{\frac{n_i-\sum_{t_j=t_i}{p_j}}{n_i}} $$

Where $ n_i $ is now the expected number of individuals at risk. (Note that I have actually used $ p_i $ as probability of death rather than censoring for notational convenience.)

For other analyses (e.g., Cox proportional hazards) you may find it more difficult to construct the likelihood function in a computationally convenient or tractable form. What you should be able to do if nothing else is sample $ D_i \sim Bernoulli (p_i) $ and perform your analyses multiple times in a Monte Carlo manner.


I have now simulated this in R and I am happy with it. Briefly I created the following data generating process:

\begin{align} T_{death} &\sim Weibull(\alpha_{death},\beta_{death}) \\ T_{cens} &\sim Weibull(\alpha_{cens}, \beta_{cens} ) \\ T &= \min\{T_{death},T_{cens}\} \\ \lambda_{death}(t) &= \frac{\alpha_{death}}{\beta_{death}}\left(\frac{t}{\beta_{death}}\right)^{\alpha_{death}-1} \\ \lambda_{cens}(t) &= \frac{\alpha_{cens}}{\beta_{cens}}\left(\frac{t}{\beta_{cens}}\right)^{\alpha_{cens}-1} \\ p(t) &= \frac{\lambda_{death}(t)}{\lambda_{death}(t)+\lambda_{cens}(t)} \end{align}

Using the estimator as I suggested above gives the following good match to the survival curve for $T_{death}$ (i.e., the censoring process has not interfered with estimation). Comparison of product-limit estimator and actual survival curve

The code for it is:


# Data generating function parameters
a_death = 1.5
b_death = 10.0
a_cens = 1.0
b_cens = 7.5
n = 1000

# Generate data
t <- c(rweibull(n*mix_deathcens, a_death, b_death), rweibull(n*(1-mix_deathcens), a_cens, b_cens))
t_death <- rweibull(n, a_death, b_death)
t_cens  <- rweibull(n, a_cens,  b_cens )
t <- pmin(t_death, t_cens)
g <- (t_death < t_cens) * 1

# Assign probabilities based on hazard rates
hr_death <- (a_death/b_death) * (t/b_death)^(a_death-1)
hr_cens  <- (a_cens /b_cens ) * (t/b_cens )^(a_cens -1)
p <- hr_death/(hr_death+hr_cens)

ord <- order(t)
df <- data.frame(t=t[ord], p=p[ord], S=rep(NA,n))
S_i <- 1
n_i <- n
for (i in 1:length(df$t)) {
  n_i <- n_i - (1 - df$p[i])
  df$S[i] <- S_i * (n_i - df$p[i])/n_i
  S_i <- df$S[i]
  n_i <- n_i - df$p[i]

# Calculate desired survival curve
df$S_death <- pweibull(df$t, a_death, b_death, lower.tail = FALSE)

# Plot
ggplot(df) + geom_step(aes(x = t, y = S)) + geom_line(aes(x = t, y = S_death)) + xlab("Time") + ylab("Survival")
  • $\begingroup$ It won't work that way. I mean, the tricky part with KM is that the product is only of death events, not censoring. $\endgroup$ Apr 9, 2015 at 7:14
  • $\begingroup$ So in this case it's the product of all events. I have simulated a dataset and analysed it in R and I'm happy with the result - will add to my answer. $\endgroup$
    – tristan
    Apr 9, 2015 at 8:57
  • $\begingroup$ @PiotrMigdal I have updated my answer $\endgroup$
    – tristan
    Apr 9, 2015 at 13:38
  • $\begingroup$ Now I see it. I.e. the hazard function estimator should be just $-\log[(n_i - p_i)/n_i]$, giving this result. $\endgroup$ Apr 13, 2015 at 21:24

Your Answer

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

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