I'm trying to write a survival model of fossil species durations. In this case, the minimal possible duration for a species is 1. Also, the general idea in paleontology is that we are only observing a small amount of the actual species that existed. The fossil record then produces a natural form of truncation.

I have a large amount of durations. Some of these durations are right-censored because the species didn't go extinct. I'm fitting a Weibull distribution to the data.

For the uncensored observations, this is pretty easy. I increment the log probability by log Weibull at the observed duration (dur_unc) and then increment by the negative log complementary cummulative density function evaluated at the truncation point (1). This follows from page 248 in the Stan 2.5 manual. I'm using the simple T[,] method from Stan because it is much cleaner.

for(i in 1:N)    
  dur_unc[i] ~ weibull(alpha,sigma) T[1,];    

For the censored observations I'm integrating out the censored values by incrementing the log probability by the log complementary cummulative density function evaluated at the (censored) duration (dur_cen).


How do I handle censored data drawn from a truncated distribution like this? The uncensored example involves just subtracting the log ccdf which was the denominator from the truncated distribution. What is the correct denominator/thing to subtract after incrementing the log probability using the log ccdf already?

Hopefully this make sense.

This is a similar question, but for BUGS. It is mostly unanswered: Is it possible to model BOTH censoring and truncation in BUGS?


1 Answer 1


From the manual, 2.8.0, page 101, just swap the normal function for Weibull,

"One way to model censored data is to treat the censored data as missing data that is constrained to fall in the censored range of values. Since Stan does not allow unknown values in its arrays or matrices, the censored values must be represented explicitly, as in the following right-censored case."

    data {
      int<lower=0> N_obs;
      int<lower=0> N_cens;
      real y_obs[N_obs];
      real<lower=max(y_obs)> U;
    parameters {
      real<lower=U> y_cens[N_cens];
      real mu;
      real<lower=0> sigma;
} model {
      y_obs ~ normal(mu,sigma);
      y_cens ~ normal(mu,sigma);

So this code basically specifies the number (U) above which data are modelled as censored data.


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.