6. Time-to-event Models: Gastric Cancer*

import arviz as az
import numpy as np
import pymc as pm
from pymc.math import exp

rng = np.random.default_rng(6420)

6. Time-to-event Models: Gastric Cancer*#

Adapted from code for Unit 8: gastric.odc.

Data can be found here.

Problem statement#

Stablein, Carter, and Novak [1981] provide data on 90 patients affected by locally advanced, nonresectable gastric carcinoma. The patients are randomized to two treatments: chemotherapy alone (coded as 0) and chemotherapy plus radiation (coded as 1). Survival time is reported in days. Recorded times are censored either if the patient’s death was not observed during the study.

Data#

Columns are, from left to right:

  • type: Treatment type, chemotherapy (0) or chemotherapy + radiation (1)

  • censored: If censored, meaning the patient survived the observation period, the time in days appears here rather than in the times column. 0 if not censored.

  • times: Recorded days until death. NaN if censored.

Model changes#

PyMC really did not like the noninformative exponential prior on v (α in this model). For some reason, the equivalent Gamma distribution is more stable. I found passing an initial value also helps avoid divergences here.

Method 1: pm.Censored#

The way PyMC censoring works is described in some detail in this notebook (Vincent [2023]). For right-censoring, try this: pm.Censored("name", dist, lower=None, upper=censored, observed=y). The censored values can be an array of the same shape as the y values.

If the y value equals the right-censored value, pm.Censored returns the complement to the CDF evaluated at the censored value. If the y value is greater than the censored value, it returns -np.inf. Otherwise, the distribution you passed to the dist parameter works as normal. What I’ve been doing is setting the values in the censored array to np.inf if the corresponding y value is not censored, and equal to the y value if it should be censored.

Note

I’ve noticed that this method is unstable with some distributions. Try using the imputed censoring model (below) if this one isn’t working.

data = np.loadtxt("../data/gastric.txt")
data.shape

(90, 3)

x = data[:, 0].copy()
censored = data[:, 1].copy()
y = data[:, 2].copy()
# for pymc, right-censored values must be greater than or equal to than the "upper" value
y[np.isnan(y)] = censored[np.isnan(y)]
censored[censored == 0] = np.inf

Warning

PyMC and BUGS do not specify the Weibull distribution in the same way! The below model is meant to be a direct translation of the lecture model, so we converted the parameters:

\[ \alpha = \nu \]
\[ \beta = \lambda ** (-1 / \alpha) \]

Warning

We have noticed that using the tau (precision) parameter to define the coefficients is sometimes unstable. We aren’t sure why. In order to stay as close as possible to the lecture model, I’ve kept it here, but when using pm.Censored, it’s better to define the normal priors for your coefficients and intercept with sigma (standard deviation).

This:

beta0 = pm.Normal("beta0", 0, sigma=10)
beta1 = pm.Normal("beta1", 0, sigma=5)

Not this:

beta0 = pm.Normal("beta0", 0, tau=0.01)
beta1 = pm.Normal("beta1", 0, tau=0.1)

log2 = np.log(2)

with pm.Model() as m_censored:
    beta0 = pm.Normal("beta0", 0, tau=0.01)
    beta1 = pm.Normal("beta1", 0, tau=0.1)
    alpha = pm.Gamma("alpha", 1, 0.001)

    lamb = exp(beta0 + beta1 * x)
    beta = lamb ** (-1 / alpha)

    obs_latent = pm.Weibull.dist(alpha=alpha, beta=beta)
    likelihood = pm.Censored(
        "likelihood",
        obs_latent,
        lower=None,
        upper=censored,
        observed=y,
    )

    median0 = pm.Deterministic("median0", (log2 * exp(-beta0)) ** (1 / alpha))
    median1 = pm.Deterministic(
        "median1", (log2 * exp(-beta0 - beta1)) ** (1 / alpha)
    )

    S = pm.Deterministic("survival", exp(-lamb * (likelihood**alpha)))
    pdf = lamb * alpha * (likelihood ** (alpha - 1)) * S
    pm.Deterministic("hazard", pdf / S)

    trace = pm.sample(
        10000,
        tune=2000,
        initvals={"alpha": 0.25},
        target_accept=0.9,
        random_seed=rng,
    )

Hide code cell output

 
NUTS[nutpie]: [beta0, beta1, alpha]













az.summary(trace, var_names=["~survival", "~hazard"], ci_prob=0.9)










  
    

      
      mean
      sd
      eti90_lb
      eti90_ub
      ess_bulk
      ess_tail
      r_hat
      mcse_mean
      mcse_sd
    

  
  
    

      beta0
      -6.77
      0.66
      -7.9
      -5.7
      6429
      8007
      1.00
      0.0083
      0.0058
    

    

      beta1
      0.263
      0.234
      -0.12
      0.65
      18125
      17275
      1.00
      0.0017
      0.0012
    

    

      alpha
      1.025
      0.097
      0.87
      1.2
      6465
      7917
      1.00
      0.0012
      0.00085
    

    

      median0
      523
      89
      390
      680
      18942
      21425
      1.00
      0.64
      0.52
    

    

      median1
      404.4
      69.8
      300
      530
      21259
      24026
      1.00
      0.48
      0.38
    

  







Note


The above model was written to reproduce the original BUGS lecture as closely as possible, so that’s why we need to include that conversion from the rate. If we were building this model directly in PyMC using the Weibull scale parameterization, we might do something like the model below.


These are similar but not identical models and the coefficients have different interpretations. Ask yourself, what do the signs of the coefficients indicate for each model?








log2 = np.log(2)

with pm.Model() as m_censored:
    beta0 = pm.Normal("beta0", 0, sigma=10)
    beta1 = pm.Normal("beta1", 0, sigma=3.2)
    alpha = pm.Gamma("alpha", 1, 0.001)

    beta = exp(beta0 + beta1 * x)

    obs_latent = pm.Weibull.dist(alpha=alpha, beta=beta)
    likelihood = pm.Censored(
        "likelihood",
        obs_latent,
        lower=None,
        upper=censored,
        observed=y,
    )

    median0 = pm.Deterministic(
        "median0", pm.math.exp(beta0) * log2 ** (1 / alpha)
    )
    median1 = pm.Deterministic(
        "median1", pm.math.exp(beta0 + beta1) * log2 ** (1 / alpha)
    )

    hazard_ratio = pm.Deterministic("hazard", pm.math.exp(-alpha * beta1))

    trace = pm.sample(
        10000,
        tune=2000,
        initvals={"alpha": 0.25},
        target_accept=0.9,
        random_seed=rng,
    )

az.summary(trace, var_names=["~hazard"], ci_prob=0.9)









NUTS[nutpie]: [beta0, beta1, alpha]








  
    

      
      mean
      sd
      eti90_lb
      eti90_ub
      ess_bulk
      ess_tail
      r_hat
      mcse_mean
      mcse_sd
    

  
  
    

      beta0
      6.604
      0.166
      6.3
      6.9
      17704
      19918
      1.00
      0.0013
      0.00093
    

    

      beta1
      -0.255
      0.235
      -0.65
      0.13
      17568
      19679
      1.00
      0.0018
      0.0013
    

    

      alpha
      1.012
      0.0972
      0.86
      1.2
      22362
      24544
      1.00
      0.00065
      0.00046
    

    

      median0
      520
      89
      390
      680
      18073
      21110
      1.00
      0.67
      0.58
    

    

      median1
      402.8
      70.5
      300
      530
      35315
      29558
      1.00
      0.37
      0.31
    

  









Method 2: pm.Potential#


This method uses pm.Potential to achieve the same thing as above by evaluating the censored datapoints differently. It came from this notebook(Junpeng Lao [2023]).






x = data[:, 0].copy()
censored_vals = data[:, 1].copy()
y = data[:, 2].copy()

# we need to separate the observed values and the censored values
observed_mask = censored_vals == 0

y_censored = censored_vals[~observed_mask]
y_uncensored = y[observed_mask]
x_censored = x[~observed_mask]
x_uncensored = x[observed_mask]

n_right_censored = int(x_censored.shape[0])
n_observed = int(x_uncensored.shape[0])













# see https://www.pymc.io/projects/examples/en/latest/survival_analysis/weibull_aft.html
def weibull_lccdf(x, alpha, beta):
    """Log complementary cdf of Weibull distribution."""
    return -((x / beta) ** alpha)













log2 = np.log(2)

with pm.Model() as m_potential:
    beta0 = pm.Normal("beta0", 0, tau=0.01)
    beta1 = pm.Normal("beta1", 0, tau=0.1)
    alpha = pm.Gamma("alpha", 1, 0.001, initval=0.25)

    lamb_censored = exp(beta0 + beta1 * x_censored)
    beta_censored = lamb_censored ** (-1 / alpha)

    lamb_uncensored = exp(beta0 + beta1 * x_uncensored)
    beta_uncensored = lamb_uncensored ** (-1 / alpha)

    pm.Weibull(
        "observed",
        alpha=alpha,
        beta=beta_uncensored,
        observed=y_uncensored,
        shape=n_observed,
    )
    pm.Potential("censored", weibull_lccdf(y_censored, alpha, beta_censored))

    median0 = pm.Deterministic("median0", (log2 * exp(-beta0)) ** (1 / alpha))
    median1 = pm.Deterministic(
        "median1", (log2 * exp(-beta0 - beta1)) ** (1 / alpha)
    )

    trace = pm.sample(10000, tune=2000, target_accept=0.9, random_seed=rng)









NUTS[nutpie]: [beta0, beta1, alpha]















az.summary(trace, ci_prob=0.9)










  
    

      
      mean
      sd
      eti90_lb
      eti90_ub
      ess_bulk
      ess_tail
      r_hat
      mcse_mean
      mcse_sd
    

  
  
    

      beta0
      -6.79
      0.67
      -7.9
      -5.7
      6321
      8716
      1.00
      0.0084
      0.0061
    

    

      beta1
      0.261
      0.233
      -0.12
      0.65
      18758
      17735
      1.00
      0.0017
      0.0012
    

    

      alpha
      1.027
      0.098
      0.87
      1.2
      6321
      8694
      1.00
      0.0012
      0.00088
    

    

      median0
      523
      88
      390
      680
      19390
      23540
      1.00
      0.63
      0.52
    

    

      median1
      405.5
      69.8
      300
      530
      21377
      25175
      1.00
      0.47
      0.38
    

  









%load_ext watermark
%watermark -n -u -v -iv -p pytensor









Last updated: Mon, 19 Oct 2026

Python implementation: CPython
Python version       : 3.14.7
IPython version      : 9.16.1

pytensor: 3.3.0

arviz: 1.3.0
numpy: 2.4.6
pymc : 6.3.2