Stan Users Guide - 4.3 Censored Data

Modeling
1

I’m relatively new to STAN but have had success with most of the examples from the users guide so far. This week I’m playing around with censored data and ran into a bit of trouble. Section 4.3 of the user’s guide:

shows two approaches: estimating censored values versus integrating them out. I started by simulating some data in R as follows

set.seed(1983)
U <- 150
N_total <- 500
mu_true <- 130
sigma_true <- 25
y_true <- rnorm(N_total, mu_true, sigma_true)
y_obs <- ifelse(y_true < U, y_true, NA)
dat <- list(y_obs = y_obs[!is.na(y_obs)], 
            N_obs = sum(!is.na(y_obs)), 
            N_cens = sum(is.na(y_obs)), 
            U = U)

and then tried each of the two approaches, copying the STAN code from the users guide verbatim in each case. The first approach (estimating censored values) worked like a charm and the resulting inferences made sense. The second approach (integrating out), on the other hand, failed. Here’s the code:

data {
  int<lower=0> N_obs;
  int<lower=0> N_cens;
  array[N_obs] real y_obs;
  real<lower=max(y_obs)> U;
}
parameters {
  real mu;
  real<lower=0> sigma;
}
model {
  y_obs ~ normal(mu, sigma);
  target += N_cens * normal_lccdf(U | mu, sigma);
}

For every draw in every chain, I received the following error message:

Chain 1   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1   Stan can't start sampling from this initial value.
.
.
.
Chain 1 Initialization between (-2, 2) failed after 100 attempts.

If I understand correctly, the problem here comes from normal_lccdf() being evaluated at values of mu and sigma that result in a tail probability that equals zero to numerical precision.

To fix this, I tried specifying an initialization function for mu and sigma. In particular, I drew mu from a uniform on mu_lower and 2 * mu_lower where mu_lower is the mean of the uncensored data. I used the same approach for sigma. This worked perfectly. But I’m hoping to build up to much more complicated models with censoring, and it may not be easy to come up with a good initialization.

The “censored data as parameters approach” doesn’t have this problem, but the “integrating out” solution is more general. I ultimately want to work with censored, hierarchical Poisson regression models, where STAN’s inability to sample discrete parameters forces me to use the second approach.

I’d be grateful for any tips or suggestions for a less kludgy solution that might be helpful as I work my way up to more complicated models. Sorry if this is a bit vague!

2 Likes