# Survival estimation when death/censoring is probabilistic

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?

• 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. – tristan Apr 8 '15 at 22:33
• @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. – Piotr Migdal Apr 9 '15 at 7:12
• 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. – tristan Apr 9 '15 at 14:37
• @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 – Cam.Davidson.Pilon Apr 17 '15 at 12:30
• @Cam.Davidson.Pilon Thanks! Yesterday I learnt about R Buy 'Til You Die, but I prefer Python (and I am already using Lifelines). – Piotr Migdal Apr 17 '15 at 17:19

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.

EDIT

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).

The code for it is:

library(ggplot2)

# 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)

# PRODUCT LIMIT
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")

• It won't work that way. I mean, the tricky part with KM is that the product is only of death events, not censoring. – Piotr Migdal Apr 9 '15 at 7:14
• 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. – tristan Apr 9 '15 at 8:57
• @PiotrMigdal I have updated my answer – tristan Apr 9 '15 at 13:38
• Now I see it. I.e. the hazard function estimator should be just $-\log[(n_i - p_i)/n_i]$, giving this result. – Piotr Migdal Apr 13 '15 at 21:24