Survival models with multiple censoring methods

Modeling
1

Hello,

I’m having issues specifying a stan model for a problem with multiple forms of censoring. When running our reliability experiment we add stress to the system in discrete pulses. We can assess whether something has experienced the survival event only after a pulse has ended. Every point y \sim Weibull(\alpha, \sigma whose event is "observed’ is interval censored between [floor(y), ceil(y)]. Additionally, some subjects survive the test and are right censored at test end. Here is my stan code.

functions {
    real my_weibull_lcdf(real t, real alpha, real sigma) {
        //required to avoid zero point during accumulation, see
        //(https://discourse.mc-stan.org/t/
        //interval-censored-data-fails-with-weibull-but-not-gamma/28780/3)
        return log1m_exp(-pow(t / sigma, alpha));
    }
}

data {
    int<lower=0> N_intervals; //number of survival interval bins
    array[N_intervals] int<lower=0> observed; //censoring codes, 
                                                                       //0 for interval, 1 for right
    vector<lower=0>[N_intervals] low_bound; //lower bound of a bin
    vector<lower=0>[N_intervals] high_bound;//upper bound of a bin
    vector<lower=0>[N_intervals] N_cens; //num censors in a bin
    int<lower=0> N_draws; //control posterior predictive draw count
}

parameters {
    real<lower=0> alpha;
    real<lower=0> sigma;
}

model {
    //weak priors for testing
    alpha ~ lognormal(0,1);
    sigma ~ lognormal(0, 1);

    //likelihood
    for (n in 1:N_intervals) {
        if (observed[n] == 0) { //interval censoring
            target += N_cens[n] * log_diff_exp(
                my_weibull_lcdf(high_bound[n] | alpha, sigma),
                my_weibull_lcdf(low_bound[n] | alpha, sigma)
            );
        }
        else if (observed[n] == 1) { //right censoring
            target += N_cens[n] * weibull_lccdf(low_bound[n] | alpha, sigma);
       }
    }
}

generated quantities {
   vector<lower=0>[N_draws] y_tilde;
   for (n in 1:N_draws) {
        y_tilde[n] = weibull_rng(alpha, sigma);
   }
}

This works as I would hope; here is an example fitting some simulated data (here a test that ends at 10):

image
image1920×1446 100 KB

Where things go off the rails is when I have multiple right censoring times. An example of how this would be achieved in practice is a second measurement on the experimental probe; if a certain signal/noise threshold is breached, we no longer trust the results of our experiment. I want to right-censor this type of observation. In practice I achieve this in my simulation by randomly flipping the censor state of individual observations \sim Bernoulli(p)

In particular, my expectation is that censoring more values this way should result in higher uncertainty estimating parameters in the Weibull but that the collection of realizations should broadly follow the uncensored simulated data. What I am seeing is that the posterior realizations rapidly diverge from observed data as I increase the rate at which data is censored:

image
image1920×1520 110 KB

image
image1920×1520 135 KB

This data can reproduce the p=0.6 plot:
p_0_6_right_censored.csv (221 Bytes)

Should I be expecting this? Am I missing something about how right censored data should affect the posterior? That the rate of right censoring affects the result this significantly suggests to me there is something I don’t understand in my model as written. Any help is appreciated!

2

Your model assumes in the observed[n] == 0 case that you only know that the value falls in the interval, not that it’s lower than a min value or higher than a max value. Is that right? I’m used to censoring happening when values are below or above a threshold. Or in this case, out of the interval. Here the observations are in the interval.

Your model then assumes the observed[n] == 1 cases are not actually observed as the variable would seem to indicate, but that they’re known only to be above low_bound[n].

And I don’t see any actual observations—everything is “censored” in some way?

I’m afraid I didn’t understand the rest of the discussion about varying right censoring times as I didn’t see any right censoring in the model.

3

Bob,

You are correct in that I am assuming in the case of observed[n] == 0 that we know the value falls within the interval. A similar case in medical statistics would be “A patient previously in remission at their last check in is now testing positive for recurrence, so we know that the recurrence event happened sometime between their last check-in and now”. The event wasn’t “observed” in the traditional sense but we can put bounds on when it must have occurred.

For observed[n] == 1 being right censored, please forgive the notation. I took guidance and notation from this paper: https://ph02.tci-thaijo.org/index.php/thaistat/article/view/251054.

My understanding of right censoring matches “they’re only known to be above low_bound[n]”. I currently think everything with observed[n] == 1 is right-censored. The implementation of this is the same as in the user-guide on censored survival models. Is this mistaken?

My general thinking behind censoring implementation is that we are accumulating probability from where we know the event occurred. If we don’t have an observation, i.e., the patient survived the study, we can only say that the patient event occurred between t = end_of_study and t = inf as required by the definition of the cumulative hazard function. I’ve tried to model that by accumulating the lcdf from time t to inf via weibull_lccdf(t | alpha, sigma)

My understanding of interval censoring (and left censoring) follows similar arguments as above; we know where the event didn’t occur (outside the interval, or to the left of time t in the case of left censoring), so we accumulate probability from the region we do know the event must have occurred in, even if we don’t know exactly where it occurred.

The additional section on right censoring is related to having various times t where points could be censored. An example I found in literature similar to my specific domain problem is studying plant die-offs. The study had three fields where they were studying plant survival and ran out of funding to return to one of the fields. For two of the fields the right censoring time was t = t_1 due to the study ending while at the field they could not return to the censoring time was t=t_2 < t_1. This is similar to my problem as I have many different data points that have various points of right censorship, e.g., each “test” is a different “field” with one unit of observation. Sometimes things have their survival event observed, and sometimes we “run out of money” and can’t actually state where the unit of observation in that field saw the survival event, so that discrete field is right censored at a specific time t_i.

This is how I ended up with a data set that is entirely censored. Events interval-censored can be likened to “we returned to the field and found the plant dead, so it died sometime before the last visit”. Right censoring is either “we ran out of money and couldn’t return to the field, yielding some censoring time t_i” or "the test ended and is censored at time t_test_end.

Does this make sense? Apologies for the wall of text.

4

Your Stan code looks right now that I understand what’s going on. It’s just semantics in the first use of “censoring” that was confusing for me, so I wanted to make sure you really did intend inclusion in the interval in the first case and above the end point in the second.

That’s fine—you just have to record the bound for each observation, which you appear to do with low_bound and high_bound. That all looks fine.

Yes.

No. I would expect this model to be calibrated if you simulate from the prior correctly. This is true of any Bayesian model where the sampler can fit it. What I mean by this is that if you simulate parameters from the prior and data from the parameters, then fit the model with the data, you should get appropriate coverage of posterior intervals (e.g., 50% intervals should have 50% coverage, or more specifically binomial(0.5, N) should be covered with N parameters). Are you just simulating from the underlying distributions and then reporting the censoring? That’s the only way you’ll be able to get the appropriate proportion of censoring for the simulated parameters.