Use of pm.Potential in Pymc3 Survival example

Question

I am working through the pymc3 Bayesian Survival Analysis example at the link below and I'm struggling to understand their use of pm.Potential:

https://docs.pymc.io/en/stable/pymc-examples/examples/survival_analysis/bayes_param_survival_pymc3.html

The heart of the example comes down to the following two snippets of code:

y_obs = pm.Gumbel("y_obs", η[~cens_], s, observed=y_std[~cens])
y_cens = pm.Potential("y_cens", gumbel_sf(y_std[cens], η[cens_], s))

Basically what is being done is they are modeling observed survival as a Gumbel distribution and censoring as a Gumbel survival function. In general that sort of makes sense to me. I do think I understand how observed deaths and observed censoring mathematically relate to these distributions. I just don't understand what pm.Potential is doing. In fact, if I take it out it doesn't appear to affect anything. Is it there to just calculate y_cens? If so, how do I access that value since I don't see it in the trace output. At this moment, I'm thinking the censored values are contributing nothing to this model and the fitted survival data is entirely based on a gumbel distribution fit of observed deaths. However, I also suspect, if I could figure out how to plot y_cens, it would also be able to give me a very similar estimate of survival trends.

Can anyone educate me on the function of pm.Potential in the above example?

Answer 1

Basically what is being done is they are modeling observed survival as a Gumbel distribution and censoring as a Gumbel survival function... I'm thinking the censored values are contributing nothing to this model and the fitted survival data is entirely based on a gumbel distribution fit of observed deaths.

That's not correct. A survival model has to incorporate the contributions of censored observations to the likelihood. Otherwise there is substantial risk of bias.

Location-scale survival modeling in this example is based on the form:

$$\log T = \eta + \sigma W ,$$

where $T$ is survival time, $\eta$ is a location parameter (the linear predictor as a function of covariate values), and $\sigma$ is a scale parameter (represented as s in your code, I think).

$W$ has a probability distribution corresponding to the underlying parametric model. For the Weibull model in your example, $W$ is standard minimum-extreme-value (one type of Gumbel).

The 2 lines of code that you rightly note to be the heart of the example describe the contributions to the likelihood separately for uncensored (first line) and right-censored (second line) survival times. This page summarizes contributions to likelihood for various types of censoring and truncation in survival analysis.

An uncensored observation contributes a factor proportional to the density of the lifetime distribution at its observation time. That's what your first line of code describes.

Provided that censoring isn't informative, a right-censored observation contributes a factor proportional to the survival function of the lifetime distribution at the (censored) observation time. (A right-censored observation provides information that an individual survived at least that long, which is the survival function of the distribution.) That's what your second line of code represents. In this case there is no modeling of censoring times, just an incorporation of likelihood contributions from censored observations into the model.

The Potential object provides a way to incorporate such likelihoods into the model without having to define an entire DensityDist. This page provides further information.

In this way both uncensored and right-censored observation times provide their contributions to the likelihood. As there is no modeling of the censoring times per se, your sampling from the posterior is simply based on the probability distribution of uncensored lifetimes. There is no need to sample censored observation times, which in this case are only used to construct the likelihood.